Leptonic D_s decays in two-flavour lattice QCD

We report on a two-flavour lattice QCD study of the D_s and D_s^* leptonic decays parameterized by the decay constants f_{D_s} and f_{D_s^*}. As the phenomenology in the D_s sector seems very promising in the next years with the experiments LHCb and Belle II, it is worth putting a big effort in lattice computations regarding its non-perturbative QCD contributions. Before examining more challenging processes such as hadron-hadron transitions, a natural first step is to address some basic aspects in the context of leptonic decays, where systematic uncertainties from excited state contaminations and cut-off effects in the computation of charmed meson decay matrix elements can be investigated in a more straightforward setting.


I. INTRODUCTION
After 60 years of joint effort by theoretical and experimental communities in particle physics, the Standard Model (SM) offers a complete picture of fundamental interactions up to the electroweak scale. In the quark sector, for instance, which is spread in 3 families, charged weak decays are mediated by a left-handed current with the exchange of a W boson, flavour changing neutral currents are forbidden at tree-level by virtue of the so-called Glashow-Iliopoulos-Maiani (GIM) mechanism [1], and CP violation shows up because the Jarlskog invariant J is different from zero [2]. On the electroweak side, the presence of a scalar field with a non-zero vacuum expectation value (VEV) induces a spontaneous breaking of the SU (2) W × U (1) Y electroweak symmetry into U (1) EM , which manifests in the Higgs mechanism [3][4][5]: Starting from a complex iso-doublet of scalar fields, 3 degrees of freedom are absorbed to give masses to the W and Z bosons such that a single field remains, the Higgs boson H. With a mass of 125.09 (24) GeV (as measured along the data taking of LHC Run 1 [6]) and a VEV of 246 GeV, its coupling to vector bosons is quadratic in their masses, while its coupling to quarks is linear in their masses.
A well-known issue with the SM Higgs is, however, that the quartic term in the Higgs Lagrangian generates for the Higgs mass m H a quadratic divergence with the hard scale of the theory, related to the so-called hierarchy problem [7,8]. There are several New Physics scenarios that are supposed to cure this caveat (amongst the others) of the SM. One class of proposed scenarios beyond the SM is characterized by a minimal extension of the Higgs sector. They contain 2 complex scalar iso-doublets Φ 1 and Φ 2 , which after the spontaneous breaking of the electroweak symmetry lead to 2 charged particles H ± , 2 CP-even particles h (an SM-like Higgs) and H plus 1 CP-odd particle A (see [9,10] and references therein for nice reviews). In these scenarios, quarks are coupled to a charged Higgs through a right-handed current. It is particularly this feature that has received a lot of attention recently, because several tests of lepton-flavour universality have shown some hints of an anomaly with respect to SM expectations, especially for the  [11][12][13]: Semi-leptonic decays with a τ lepton in the final state can have a non-SM contribution from the exchange of a right-handed current that is not helicity-suppressed by the mass of the charged lepton. In view of the highly promising perspective that Belle II uses a part of its integrated luminosity to run at the energy of Υ(5S) in order to accumulate B s pairs, it then might be very valuable to investigate equivalent ratios R D ( * ) s , where just the spectator quark in the aforementioned processes is changed. Therefore, on the theory side, the non-perturbative hadronic properties of the B s , D s and D * s mesons involved, which are accessible through lattice QCD, have to be under very good control. This paper represents the very first step in this program, and it reports on an estimate by N f = 2 lattice QCD of the ratio of charm-strange leptonic decay constants in the vector and pseudoscalar channels, f D * s /f Ds , which quantifies spin-breaking effects in heavy-strange mesons. While a large number of results for f Ds is available in the literature [14][15][16][17][18][19][20][21][22][23][24][25][26][27], there are so far only a few of them for f D * s [28][29][30]; moreover, this study serves as an opportunity to learn how to identify and alleviate systematic effects such as contaminations from excited states in correlation functions and potentially large lattice artifacts, which are typically encountered in simulations of D s meson systems.

A. Lattice set-up and analysis techniques
This work is based on a subset of the CLS ensembles [31], made of N f = 2 non-perturbatively O(a) improved Wilson-Clover fermions [32,33], together with the plaquette gauge action [34] for gluon fields, and generated using either the DD-HMC algorithm [35][36][37][38] or the MP-HMC algorithm [39]. We collect our simulation parameters in Table I. Two values of the lattice spacing, a β=5.5 = 0.04831 (38) fm and a β=5.3 = 0.06531(60) fm as determined from a fit in the chiral sector [40], are considered, with pion masses in the range [190 , 440] MeV. The bare (valence) strange quark mass had been tuned by imposing the ratios m 2 K /f 2 K and m 2 π /f 2 K to coincide with their physical values, where the scale is set by f K [41]. Irrespective of the sea quark mass, the bare mass parameter of the valence charm quark (represented by the hopping parameter κ c ) was fixed in [43] through a linear interpolation of m 2 Ds in 1/κ c to the physical D s meson mass (= 1968 MeV [42]), since this functional form turned out to be best supported by the actual data in the vicinity of the targeted m phys Ds . Eventually, statistical errors on each ensemble are estimated from the jackknife procedure 1 , whereas statistical errors on quantities calculated as results of a final joint chiral and continuum limit extrapolation to the physical point are obtained by the following bootstrap-inspired prescription: Create a large set of N event "event" vectors, the dimension of which is given by the number of CLS ensembles considered in our analysis (= 6 in the present case), and fill them component-wise with entries randomly chosen from the available sample of N bin[#id] jackknife-binned data per ensemble (#id=1,...,6). The statistical error of the result of any extrapolating fit is then estimated as the variance over N event such fits, taking these random vectors (which via their very composition can be understood as being drawn from the actual statistical distribution of the contributing raw ensemble data) as inputs.
Two-point correlation functions are evaluated in the standard manner by first expressing them as expectation values of a product of two quark propagators (indicated as [...] in the equation below), where the latter are computed using stochastic sources defined in a randomly chosen timeslice with spin dilution. In addition, the noise has been reduced by applying the one-end trick [44,45]. We study the two-point correlation functions where V is the spatial volume of the lattice, ... denotes the expectation value over gauge configurations, and the interpolating fieldscΓs are not necessarily local. Four Gaussian smearing levels for the quark fields s and c, including the case of no smearing, are considered to build a 4 × 4 matrix of correlators, from which we extract all O(a) improved hadronic quantities relevant here, after analyzing the associated generalized eigenvalue problem (GEVP) [46][47][48]. Solving the GEVP for the pseudoscalar-pseudoscalar and vector-vector matrix of correlators, we have constructed the corresponding projected correlators (as well as their symmetric counterparts by exchanging operators at the source and at the sink); as local quark bilinears, we here employ the composite fields These projections, together with their asymptotic behaviour as t/a 1, read where the various Z stand for the meson-to-vacuum matrix elements of the operators in the respective channels, which arise in the spectral decompositions of the correlation functions. The label "L" refers to a local interpolating field, while the sums over i and j run over the 4 smearing levels.
Upon including the appropriate O(a) terms in the definitions of the axial and vector currents, viz.
the forward (∂ ν ) and backward (∂ * ν ) difference operators act on the (foregoing large-t/a asymptotics of the) projected correlators as Then, together with the asymptotic behaviour of the correlation functions in eqs. (5) and (6), the lattice expressions for the pseudoscalar-and vector-to-vacuum matrix element of the renormalized axial and vector current, which are proportional to leptonic pseudoscalar and vector meson decay constants of interest, can be split into leading and O(a) improvement contributions according to: where p and λ i stand for the spatial meson three-momentum and the vector meson polarization, respectively.
The renormalization constants Z A and Z V , which multiply the (time component of the) axial current and the (spatial components of the) vector current, were determined non-perturbatively for two-flavour QCD with O(a) improved Wilson fermions in [49,50]. For the corresponding improvement coefficients c A , c V , b A and b V we have used non-perturbative estimates, if available, and perturbative formulas elsewhere [51,52]. Moreover, the mass dependent factors in eqs. (11) and (12) involve the matching coefficient Z between the average of bare quark masses defined via the axial Ward identity (also called PCAC quark mass) and its counterpart defined through the vector Ward identity (also called subtracted quark mass). On the lattice, and specifically written for the charm-strange sector considered in this work, these two definitions of m cs ≡ (m c + m s )/2 translate into and am VWI in practice, the local mass m AWI cs (t) exhibits an extended plateau over the timeslices sufficiently far from the boundaries of the lattice and therefore can be accurately determined with confidence as the plateau average over central timeslices. Coming back to the factor Z in the present N f = 2 case, non-perturbative numbers for it were taken from [41,53]. The uncertainties of the non-perturbative values of those various c-and b-coefficients and Z-factors, as far as they are available from the literature, were incorporated in the subsequent analysis, adding them in quadrature with the independent statistical errors of our quantities in question. Thereby, these effects are properly propagated into the final results, yet their contribution is subdominant compared to the combined statistical and lattice spacing errors (from the scale setting) discussed below.

B. Results
As a starting point of our analysis to extract the charm-strange pseudoscalar and vector meson masses and decay constants along the lines outlined in the last subsection, we have fitted the projected correlators C within a time range [t min , t max ] such that the statistical error δm stat (t min ) on the ground-state effective mass E 1 is significantly larger than the systematic error caused by contaminations from higher, excited states. As an estimate of the latter, we advocate a mass gap of To be on the conservative side, we always imposed δm stat (t min ) > 4∆m sys (t min ). This is also visualized semi-logarithmically in the left panel of Fig. 1 for a typical case, which demonstrates that, by virtue of our prescription to select the fitting interval for the ground-state extraction, the systematic effect from residual excited state contaminations within this interval is negligible by contrast with the statistical errors. Despite having solved a 4 × 4 GEVP, we notice that the third excited state is not well under control, because huge statistical errors on the corresponding effective mass are encountered. Nevertheless, based on the overall landscape of resulting energies from the low-lying states, we find our above guesstimate ∆ E ∼ 2 GeV to be safe. The upper bound t max of the particular fit interval per ensemble was fixed by an individual visual inspection of the quality of the effective mass plateaux; its actual influence on the final numbers is insignificant though, as long as one stays within the plateau region. Since it can be assumed (and was numerically confirmed on a few representative datasets) that ∆ P E ∼ ∆ V E holds to a good approximation for the mass gaps in the pseudoscalar and vector channels at fixed lattice spacing and pion mass, identical fit intervals were chosen for pseudoscalar and vector meson states. Lastly, as for the timeslice parameter t 0 in the GEVP analysis (see eqs. (5) and (6)), we have stuck to t 0 = 3a at β = 5.3 and t 0 = 5a at β = 5.5, after we had observed that the resulting energies and matrix elements do not change appreciably and only the statistical errors grow.
The right panel of Fig. 1 illustrates the time dependence of the effective masses of the D s and D * s mesons obtained from the gauge configuration ensemble F7. They were calculated by the formulas In Table II of the appendix we specify all masses and decay constants as they result from the application of the GEVP analysis, in conjunction with eqs. (5) - (12), to the CLS ensembles for the purpose of this study. After all, an approach to the physical point still amounts to perform the chiral and continuum limits. In order to extrapolate the masses and (the ratio of) decay constants from the GEVP analysis to this point comprising both limits simultaneously, we have employed a simple ansatz with a linear term in m 2 π and leading cutoff effects, in the non-perturbatively O(a) improved theory, proportional to a 2 : Note that our calculation, with two lattice spacings below 0.1 fm and (a max /a min ) 2 ≈ 1.8 > 1.4 fulfills the criteria of the second level (out of three) within the quality rating regarding the continuum extrapolation by the FLAG Working Group [54]. Since the datasets actually involve lattice spacings only, which both are even below 0.07 fm, and the underlying action and composite fields are improved and renormalized non-perturbatively, we consider the linear description of the lattice spacing dependence in eq. (16) to be adequate to estimate the systematic uncertainty related to residual discretization errors. This is also in accordance with earlier findings in the B meson sector building upon almost the same collection of CLS datasets [55,56]. The strange and charm (valence) quark mass hopping parameters κ s and κ c were tuned at every sea quark mass (given by κ sea ) such that m Ds (κ sea ; κ s , κ c ) = m phys Ds . As a consequence, the formulas inspired by partially quenched heavy-light meson chiral perturbation theory (HLMχPT) [57], which commonly serve as a guide to extrapolate lattice results of heavy-light meson masses and decay constants to the physical point, can be simplified to a leading-order (LO) chiral expression as done with our model ansatz, eq. (16). As a cross-check, we have also added a next-to-leading (NLO) HLMχPT term, which has the logarithmic form ∝ m 2 π log m 2 π with fit parameter X 4 , to this ansatz. The outcome of these LO and NLO fits to the joint chiral and continuum limits, together with the fit parameters, their errors and associated χ 2 /d.o.f.'s, are detailed in Table III of the appendix. As can be read off from it, the inclusion of a NLO contribution leads to a (in most cases noticeable) deterioration of the quality of the fits, and the corresponding fitted coefficient X 4 is always compatible with zero, while the physical results stay consistent. We thus do not account for the discrepancy between the LO and NLO chiral extrapolations as a separate piece in the systematic error, because there are not enough data points to claim sensitivity of our data to such an NLO term, which would be a prerequisite to justify a reliable NLO analysis.
Before presenting the final results, we still have to address the effect of a possible mistuning of κ c and κ s . As for κ c , it was fixed to match m phys  Table II then reveals that all data points for am Ds are fully compatible with am phys Ds -even more so upon accounting for the uncertainty on the lattice spacings from the scale setting paper [41] -, except for ensemble E5. To quantify the impact of mistuning on the results from this ensemble, we have performed additional calculations at a reasonably varied hopping parameter for the valence charm quark, κ c = 0.12735 instead of 0.12724. This yields a relative decrease by O(0.6%) from am Ds = 0.659(1) to 0.654(1), accompanied by changes in the other quantities (af Ds , am D * s , af D * s , f D * s /f Ds ) that are at most about 1/2 of their statistical errors. As for the effect from the mistuning of the hopping parameter of the valence strange quark, κ s (which is largely independent of the scale setting, because it was obtained via the ratio m 2 K /f 2 K in [41]), we have studied it through additional measurements of the relevant correlators for E5 at a shifted value 2 of κ s = 0.135827, which lies apart by three times the statistical error quoted for κ s = 0.135777 (17) in [41]. Again, also in this case, the corresponding changes in our mesonic observables are at the same level of about 1/2 of the respective statistical errors (or even smaller), with the exception of f D * s /f Ds that slightly drops to 1.23(2) but within its error remains consistent with the value f D * s /f Ds = 1.25 (2) in Table II. Therefore, since a marginal mistuning of κ c and κ s only affects the (coarser and least chiral) ensemble E5, but all in all turns out to be insignificant there, and is absent for all the other ensembles, it appears safe to neglect this contribution to the systematic error beyond the uncertainty on the lattice spacing, which of course is propagated into the overall errors quoted below.
We now come to the results of our analysis for the physical quantities under disposal. As can be inferred from vector channel, m D * s , nicely reproduces its experimental value of 2.112 GeV [42], with cutoff effects being limited to about 0.5% at β = 5.3. More precisely, we arrive at where the first error is of statistical nature, and the second one reflects the uncertainty in the lattice spacing induced by the scale setting [41]. Analogous physical point extrapolations of the pseudoscalar and vector meson decay constants, f Ds and f D * s , as well as of their ratio f D * s /f Ds are displayed in Fig. 3 and overall show an only quite mild dependence on the pion mass and the lattice spacing. In particular, cutoff effects on f Ds are limited to ∼ 1% at β = 5.3, while stronger scaling violations of the order of 7% are observed for f D * s . Hence, they also propagate with a contribution of about 6% into the total error on the ratio f D * s /f Ds . We quote as our main result where in case of this ratio the systematic error that stems from the uncertainty in lattice spacings is negligible on the level of precision here. Finally, as our estimate for f Ds at the physical point we obtain f Ds = 238(5)(2) MeV. Again, the first error is statistical, while the second one incorporates the uncertainty from the scale setting. This value is lower by about 1.8σ than the N f = 2 lattice QCD average quoted by the FLAG Working Group f FLAG,N f =2 Ds = 250 (7) MeV [22,54] 3 . Let us emphasize, however, that the agreement in f Ds is satisfactory, when comparing to the outcome of the independent two-flavour computation in [43], which uses almost the same CLS ensembles as in this work. There, the extraction of the D s meson decay constant follows from an expression, which combines the axial Ward identity (resp. PCAC) quark mass m AWI cs , cf. eq. (13), with the pseudoscalar-to-vacuum matrix element of the pseudoscalar density operator (the latter being obtained through a fit of the local pseudoscalar correlator C P L P L (t)) and leads to results consistent with the ones reported here.

C. Discussion
So far, there are only two lattice estimates of f D * s /f Ds for N f = 2, namely by ETMC [28] and us, and two other ones have been performed for N f = 2 + 1 by HPQCD [29] and N f = 2 + 1 + 1 by ETMC [30],  respectively. We summarize the various results in Fig. 4. In the past, owing to the apparent discrepancy between the two-flavour ETMC result of about 1.25 and the N f = 2+1(+1) determinations, it was thought that f D * s /f Ds could be a quantity, where a quite large quenching effect of the strange quark shows up, with an amount of ∼ 10% or even more. Our finding, employing a lattice discretization of the two-flavour theory different from the ETMC calculation in [28], tends however to point to the conclusion that this effect is significantly less pronounced; nonetheless, the trend that less spin-breaking effects are present, if more flavours are active, still remains to be visible when placing the result of our study to the circle of the other lattice estimates in Fig. 4.

III. CONCLUSION
In this paper we have reported on a two-flavour lattice QCD computation of the ratio of vector to pseudoscalar decay constants f D * s /f Ds , based on simulations with non-perturbatively improved Wilson fermions, satisfying non-perturbative renormalization and accounting for correlations between the two decay constants in the (statistical) error analysis. As an interesting lesson we can state that the quenching of the strange and the charm quark has an only moderate impact (of the order of ∼ 5%) on a ratio, which is quite of phenomenological pertinence in quantifying the rate of spin-symmetry breaking in heavy-strange  s transitions by means of lattice QCD, a strategy inspired by the "step scaling in mass" method of [59,60] will be adopted to extrapolate results to the B meson region. A significant difference, however, is that in the present case of Wilson-Clover regularization we cannot use renormalization group invariant (RGI) quark masses to impose lines of constant physics for continuum limit extrapolations. Indeed, owing to the lack of knowledge of the O(a 2 ) improvement coefficient b m , which would enter in the definition of the RGI quark mass in terms of the bare (vector Ward identity) quark mass, viz.
where H is denotes a heavy-strange meson made out of quarks with masses κ hi and κ s . By this we are able to access the B physics region in a controlled way, which constitutes a firm basis to proceed with our program of studying there (semi-)leptonic decays and lepton-flavour universality violations through lattice QCD.
In Table II we list our numerical results for m Ds , f Ds , m D * s and f D * s in lattice units, as well as for f D * s /f Ds , which were obtained via the GEVP analysis of the correlation functions evaluated on the CLS gauge field ensembles contributing to this work. Table III summarizes the outcome of the fits of the final results on m D * s , f Ds , f D * s and f D * s /f Ds in physical units to the ansatz (16) for the global fit modeling their approach to the chiral and continuum limits. The next-to-leading-order (NLO) fit type also incorporates a logarithmic term from HLMχPT (proportional to m 2 π log m 2 π with fit parameter X 4 ), on top of the leading-order (LO) expression in eq. (16).  Results of the leading-order (LO) and next-to-leading-order (NLO, i.e., also including a logarithmic HLMχPT-inspired term with fit parameter X4) joint chiral and continuum extrapolations to the ansatz in eq. (16). Note that the χ 2 /d.o.f.'s of the underlying fits are given in the second column, too.