Decay constant of Bs\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$B_s$$\end{document} and Bs∗\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$B^*_s$$\end{document} mesons from Nf=2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathrm{N_f}=2$$\end{document} lattice QCD

We report on a two-flavor lattice QCD estimate of the Bs\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$B_s$$\end{document} and Bs∗\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$B^*_s$$\end{document} leptonic decays parameterized by the decay constants fBs\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_{B_s}$$\end{document} and fBs∗\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_{B^*_s}$$\end{document}. In addition to their relevance for phenomenology, their extraction has allowed us to investigate whether the “step scaling in mass” strategy is suitable with Wilson–Clover fermions to smoothly extrapolate quantities of the heavy-strange sector up to the bottom scale. From the central value of fDs\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_{D_s}$$\end{document} quoted by FLAG at Nf=2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N_f=2$$\end{document} and our ratio fBsfDs\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{f_{B_s}}{f_{D_s}}$$\end{document}, we obtain fBs=215(10)(2)(-5+2)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_{B_s}=215(10)(2)(^{+2}_{-5})$$\end{document} MeV and fBs∗/fBs=1.02(2)(-0+2)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_{B^*_s}/f_{B_s}=1.02(2)(^{+2}_{-0})$$\end{document}.


Introduction
In the very active research of new effects in high-energy particle physics, flavour physics does play a key role at the so-called intensity frontier. Indeed, rare events are sensitive probes of New Physics (NP) scenarios with the exchange of extra particles in quantum loops with respect to what is known from the Standard Model (SM), However, theoretical uncertainties on hadronic quantities, for instance hadron decay constants, that encode the dynamics of QCD at large distance, severely weaken the constraints that are derived through the analysis of experimental data. Those hadronic constants cannot be reliably estimated in perturbation theory. b-quark physics is a particularly interesting place to search for NP effects and it has recently regained even stronger attention after experimental signs of several anomalies in B and B c decays. More precisely several ratios R D = (B→Dτ ν τ ) (B→K e + e − ) and R K * = (B→K * μ + μ − ) (B→K * e + e − ) show some discrepancy with SM expectations [1][2][3][4][5][6][7][8][9][10][11][12]. The three former might bring stringent constraints onbc currents, for instance mediated by the exchange of leptoquarks [13]. The fura e-mail: benoit.blossier@th.u-psud.fr (corresponding author) ther ratios R D ( * ) , under investigation at LHCb, will provide even more informations once, on the theory side, the hadronic matrix elements associated to B s → D ( * ) s are under comparable control by means of lattice QCD. Simulating the B s meson on the lattice is delicate as far as cut-off effects are concerned. Several strategies have been followed in the literature, including simulations of relativistic b-quarks using an action tuned so as to minimize discretization errors [14][15][16][17], the use of Non Relativistic QCD [18,19], performing computations in Heavy Quark Effective Theory (HQET) [20] and the extrapolation of simulation results obtained in the region between the charm quark mass m c and a mass ∼ 3m c to the physical b-quark mass [21,22]. As we plan to employ the latter approach to study B s decays with O(a) improved Wilson-Clover fermions, an intermediate step is to extract f B s and f B * s , in order to validate the method. The lattice QCD community has made a significant effort to compute f B s with N f = 2 [22,23], N f = 2 + 1 [14,15,18,24,25] and N f = 2 + 1 + 1 [19,[26][27][28]. Recently, the SU(3) symmetry breaking f B s / f B has been extracted at the physical point [29]. Concerning the spin-symmetry breaking ratio f B * s / f B s only 2 lattice results are available, both at N f = 2 + 1 + 1 [30,31]. Ratios f B * / f B and f B * s / f B s have been investigated with other methods than lattice simulations, i.e. constituent quark models [32,33] and QCD sum rules [34][35][36][37].
The paper is organized as follows: in Sect. 2 we recall what is the "step scaling in mass" strategy, in Sect. 3 we present the simulations details and our raw data, and in Sect. 4 we describe our analysis and comment the results. Finally we conclude in Sect. 5.

Step scaling in mass with Wilson-Clover fermions
The idea is to extract f B s ≡ In the following we focus on the two latter in the framework of N f = 2 Wilson-Clover fermions. With a given pion mass, through κ sea , the valence strange quark mass κ s and lattice spacing a, we consider the ratio where λ = m Bs m Ds heavy quark mass, and M P (a, κ sea , κ h 0 , κ s ) ≡ m D s , up to mistuning effects. For later usage it is convenient to redefine r P as , that will appear later in the paper, are the matching coefficients between the QCD cur- where μ Q is a scale related to the heavy quark mass m Q , for instance the heavy-light pseudoscalar meson mass M P (Q, q). C stat A and C stat V are known up to 3-loop of perturbation theory [38][39][40][41]. 1 r P is independent of μ f because it involves the renormalization group equation of C stat A integrated from M P (a, κ sea , κ h i , κ s )) to M P (a, κ sea , κ h i+1 , κ s ). Thanks to scaling laws in HQET, it is expected that r P has a simple expansion in the inverse heavy quark mass defined in a specific renormalization scheme, for instance the pole mass [42]. But, in the case of Wilson-Clover fermions and contrary to the case of twisted-mass fermions, there is so far no straightforward relation between the pole quark mass and the bare quark mass, through the renormalization group invariant (RGI) quark mass, if the quark is significantly heavier than the charm. Indeed, using the series of RGI masses m RGI already determined [43] and m RGI c known after the tuning of κ c , we define the improved RGI mass m RGI we can have 1+b m am h i < 0 because the improvement coefficient b m is negative [44]. Negative RGI masses are of course not physical. The issue would be solved by adding the O(a 2 ) term in the definition of the RGI mass, which is unfortunately unknown. That is why we have decided to consider the inverse of pseudoscalar heavy-strange pseudoscalar meson masses M P (κ h i , κ s ) ≡ M H s (i) as the parameter expansion of r P and we define the steps as sequential ratios that should be constant with the regulator a and the sea quark mass κ sea . They are the analogous of the ratios of RGI quark

Lattice calculation details
We have performed our analysis from the CLS ensembles made of N f = 2 nonperturbatively O(a)-improved Wilson-Clover fermions [45,46] and the plaquette gauge action [47].
In Table 1  . The values we find for κ c are close to what is quoted in [50] where the tuning was realised thanks to a constraint on cut-off effect magnitude for the ratio of PCAC masses m PCAC c /m PCAC s . We have used the same procedure as in [49] to compute the statistical error at finite a and in the continuum limit, to compute stochastic all-to-all propagators and to reduce the contamination by excited states on 2-pt correlators by solving a 4 × 4 Generalized Eigenvalue Problem (GEVP) with one local and 3 Gaussian smeared interpolating fields. In our application of the step scaling in mass strategy we have chosen K = 6 steps. Similarly to [49] we have extracted the relevant matrix elements from projected correlators along the fixed ground state generalized eigenvectors v (1) take the values (4a, 3a) at β = 5.2, (4a, 3a) at β = 5.3 and (6a, 5a) at β = 5.5. We collect in Tables 10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25 of the Appendix the whole set of raw data we need in our analysis, i.e. ratios of pseudoscalar and vector heavy-strange mesons masses for 2 subsequent heavy bare quark masses, ratios of pseudoscalar and vector heavy-strange meson decay constants for 2 subsequent heavy bare quark masses and PCAC quark masses m hs defined by: Parameters of the simulations: bare coupling β = 6/g 2 0 , lattice resolution, hopping parameter κ, lattice spacing a in physical units, pion mass, number of gauge configurations, bare strange and charm quark masses , C A 0 P and C P P are axialpseudoscalar and pseudoscalar-pseudoscalar 2-pt correlation functions of the heavy-strange meson defined by the bare quark masses (κ h , κ s ) and c A is the improvement coefficient of the axial bilinear of Wilson-Clover fermions determined in [51]. We show in Fig. 1 the effective mass for two heavy-strange mesons, one at the first step scaling in mass and the other at the next to last step. For very heavy quarks the signal deteriorates quickly. It explains why we fixed shorter interval ranges to extract the hadronic properties, as indicated in Tables 10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25.

Extraction of f B s
We have performed extrapolations to the physical point by doing a global fit analysis. However, in a preparatory stage, we restrict our analysis to a given step in heavy mass i and study the pion mass and the cut-off dependence of r P (i).
Here we ignore the mistuning effects because our goal is to determine how large are the discretisation effects. We show in Fig. 2 the extrapolation in a 2 and m 2 π of r P (i). We observe very big cut-off effects for the 4 th ratio (10%) and the 5 th ratio (17%). Hence we are not quite confident in using the ensembles at a = 0.075 fm and, most probably, those at a = 0.065 fm as well, at the final stage of our analysis. That is why we prefer, in the combined fit analysis, to exclude the data at the fourth and the fifth heavy quark mass at the lattice spacings 0.075 fm and 0.065 fm, Then, we have used the following fit ansatz: Table 2 the fit parameters and χ 2 /d.o. f.. 2 We show in Fig. 3 the dependence of r P on 1/M H s found with the fit formula (4). We retrieve the parametrically large cut-off effects ∝ (a M H s ) 2 on r P , justifying our decision to exclude some data of our 2 coarsest lattices in the analysis. It is reassuring that the pion mass dependence is found to be of the order of a few % and that it is numerically a sub-leading effect: by construction of r P , pion mass effects are expected to vanish. Our way to define r P is such that, according to the HQET scaling law telling that lim . With our value of λ = 1.18, the limit is expected to be 0.92, in excellent agreement with our fit parameter 1 + r P 0 = 0.94 (2). Then, in the continuum, we interpolate r P at the 6 points m D s λ i+1 to get a set of 6 ratios r P (i): We recall that r P (i) is independent of the renormalization scale μ f and that we have f Hs (i) . We collect in Table 3 values of r P at the reference points The last step is straightforward: f B s / f D s is obtained by a series of products: The effect on the statistical error of the correlation among the different terms of the product of r P (i) is taken into account by the mean of computing the errors described in [49].
To address the systematic error, we have performed two other fits: -fit(A), adding to (4) a "next to leading order" chiral contribution in m 2 π ln(m 2 π ) -fit(B), adding to (4) a contribution in 1/m 2 H s (i) to count for a higher order in the heavy quark expansion -fit(C): fit (4) but using the matching coefficient C stat A at NLO Other fits with extra terms in (a M H s ) 2 or in a 3 give non reliable results. We collect the corresponding fit parameters and χ 2 /d.o. f. in Table 4 and we get A fourth source of systematics can be included by propagating the uncertainty on raw data if we change t min → t min + 2a to extract plateaus. In this case we get the following result: We collect the corresponding fit parameters in Table 5.
Adding together the different sources of systematics, we obtain where the first error is statistical and the second error counts for the systematic error.
Concerning the D s meson decay constant, an update of the analysis reported in [49], now that we have the additional coarsest ensembles A5 and B6, meaning a third lattice spacing at our disposal, gives where the first error is statistical and the second error comes from the uncertainty on the lattice spacings. As in the previous analysis, a next to leading order contribution to the chiral fir destabilises the fit with chiral fit parameters compatible with zero.
Then we get where the first error is the statistical error, the second one counts for the systermatic error on f D s while the third error corresponds to the systematic error on f B s / f D s . FLAG has recently made an update collection of lattice estimates of f B s [52]. Our estimate of f B s using the step scaling in mass strategy is compatible with the value obtained by the ALPHA Collaboration f B s = 224(14) MeV [23] by a     computation, performed over almost the same CLS ensembles as in this paper, of hadronic matrix elements in the framework of HQET with a non-perturbative matching of the HQET parameters with QCD. It is 2σ lower than the result reported by the ETM Collaboration [22] with N f = 2 twisted-mass fermions defined at maximal twist. The fact that we cannot constrain the static limit of the ratio r P to be equal to 1, due to mistuning effects of the heavy quark mass, explains a part of that discrepancy. The second source of discrepancy is the presence of large a 2 (a M H s ) 2 cut-off effects in our data while they are numerically absent in ETMC data. Having to take them into account necessarily increases the uncertainty in extrapolation to the continuum limit because more parameters are required to described the data.

Extraction of f B * s / f B s
To extract f B * s / f B s we have performed an alternative analysis to the one discussed in the previous subsection. We have examined the ratios As the HQET anomalous dimension of the axial and vector static-light operator are the same, applying the renormalization group equation makes R f * independent of the renormal-ization scale μ f . To extrapolate to the physical point we have used the following fit ansatz: We can impose the static limit constraint lim M Hs →∞ R m * = lim M Hs →∞ R f * = 1 because those ratios are free of heavy quark mistuning effect. We collect in Table 6 the corresponding fit parameters and we obtain We show in Fig. 4 the extrapolation to the physical point of R m * and R f * .
To estimate the systermatic error, we have performed fits (A') and (B'), that read -fit(A'): add to (14) and (15) a contribution in m 2 π ln(m 2 π ) -fit(B'): add to (14) a contribution in 1/m 3 H s -fit(C'): fit (15) but using matching coefficients C stat A and C stat V at NLO  (15) and (C') is using the expression (15) and the formulae at NLO of the matching coefficient C stat  Other terms in the fit lead to unstable and reliable results. We collect the fit parameters in Tables 7 and 8 and we obtain: As for f B s / f D s , we have counted for the systematic error coming from the change t min → t min + 2a in the plateaus extraction. We get the following results: Fit parameters are collected in Table 9.
Eventually, we quote where the first error is statistical and the second error counts for the systematic error estimated by the fits (A') and (B') and contamination from excited states. There is no reason why the ratio

Comment
We collect in Fig. 5 the lattice QCD estimates of f B s at N f = 2 [22,23], with the corresponding FLAG average [52] and those of f B * s / f B s [30], [31]. Of course the fact that we get f B * s / f B s > 1 while the 2 other lattice QCD results read f B * s / f B s < 1 is puzzling. However, a computation performed by the ETM Collaboration with N f = 2 dynamical quarks indicated the hierarchy f B * / f B > 1, f B * / f B = 1.050(16) [54]. So, a plausible explanation for the observed tension is the effect of the quenching of the strange quark in the spinbreaking contribution of the heavy quark symmetry to the ratio f B * s / f B s . It might be of the same order of magnitude as in f D * s f D s but with a more important qualitative impact because we examine a region of paramaters closer to the symmetric point f H * / f H = 1. In that respect, studies of this ratio with N f = 2 + 1 ensembles are welcome.