Extraction of $B_s \to D^{(*)}_s$ form factors from ${ N_f}=2$ lattice QCD

We report on a two-flavour lattice QCD study of the $B_s \to D_s$ and $B_s \to D^*_s$ transitions parameterized, in the heavy quark limit, by the form factors ${\cal G}$, and $h_{A_1}$, $h_{A_2}$ and $h_{A_3}$, respectively. In the search for New Physics through tests of lepton flavour universality, $B_s$ decay channels are complementary to the $B$ decays widely studied at $B$ factories and LHCb, while on the theory side they can be better controlled than the $B_c$ and $\Lambda_b$ decays. The purpose of this exploratory two-flavour study is, in preparation for future analyses of lattice QCD simulations with ${N_f>2}$ and physical quark-masses, to gain experience on a suitable method for a lattice extraction of form factors associated with $b \to c$ currents that may yield tighter control over systematic effects like contamination from excited states and cut-off effects. We obtain the zero-recoil values ${\cal G}^{B_s \to D_s}(1)=1.03(14)$ and $h^{B_s \to D^*_s}(1)=0.85(16)$.


Introduction
Though the Standard Model of particle physics (SM) has shown to describe the fundamental interactions up to the electroweak scale pretty well, experimental measurements sometimes give results that were totally unexpected. In particular, a recent and spectacular example in flavour physics phenomenology is the so-called B anomalies in the test of lepton flavour universality in b → c semileptonic decays. The ratio of decay widths R Hc ≡ Γ(H b →Hcτ ντ ) Γ(H b →Hc ν ) =e,µ has been considered in cases (H b , H c ) = (B, D), (B, D * ) and (B c , J/ψ) and compared with theoretical expectations. A discrepancy of 2σ has been observed for R D [1][2][3][4][5][6][7][8][9][10], 3σ for R D * [8][9][10][11][12][13][14][15] and ∼ σ for R J/ψ [16,17]. A bunch of New Physics models have been advocated to explain those discrepancies, together with other anomalies in b → s transitions, like models with extra doublets of Higgs fields [18] or with leptoquarks [19]. Further measurements are led to explore more b → c processes, such as Λ b → Λ c lν l or B s → D ( * ) s lν, and are expected to come soon at LHCb and Belle 2. A theoretical effort has been undertaken to extract the form factors associated with B s → D ( * ) s lν, but the number of results is still limited [20][21][22][23][24][25]. Our aim here is to study whether Wilson-Clover fermions, in combination with the stepscaling in mass method [20,26], provide a suitable lattice regularisation to get reliable results for these processes, as far as cut-off effects and contamination by excited states are concerned. While our current exploratory work focuses on N f = 2, any phenomenologically relevant application of the methods investigated here has to be done in lattice QCD with N f = 2 + 1(+1) flavours.
The paper is organised as follows: in Section 2 we highlight the strategy we have come up with, in Section 3 we give the simulation details and collect our raw data. In Section 4 we present our analysis and comment on the results. Section 5 contains the conclusion.

Strategy
The starting point is to consider the hadronic matrix element D s (p Ds )|V µ |H s (p Hs ) , with V µ =cγ µ h. Its Lorentz structure decomposition reads D s (p Ds )|V µ |H s (p Hs ) = f + (q 2 )(p Ds + p Hs ) µ + f − (q 2 )q µ , (2.1) with the quantities f + and f − familiar as form factors depending on q = p Hs − p Ds . We will be interested in the physics case H s ≡ B s . In phenomenology it is convenient to introduce the scalar form factor f 0 defined by: When the H s meson is put at rest, the recoil w = p Ds ·p Hs m Hs m Ds ≡ v Hs · v Ds has the simple expression w = E Ds /m Ds . Our goal is to extract the form factor f + (B s → D s ) at zero recoil, w = 1. However, except in the elastic case D s → D s where one finds trivially that f + (w = 1) = 1, this kinematical point is not directly accessible to a lattice measurement. It is only by extrapolating in w that one can obtain results at zero recoil. With q = −(θ, θ, θ), w ∼ m 2 Ds + 3θ 2 /m Ds , we have  h ± (1), f ± (1) and f 0 (1) can then be extracted by doing a polynomial extrapolation in w − 1. In the following we will concentrate on estimating of G Bs→Ds (w = 1).
To obtain results at the b quark mass (while trying to keep cut-off effects under control) we have decided to apply the strategy of step-scaling in masses [20,26] that was already adapted to the Wilson-Clover regularisation in [27]. Steps in RGI heavy quark mass cannot be used at this stage of knowledge as far as the Wilon-Clover regularisation is concerned. Indeed without a still unknown O(a 2 ) term in an improvement factor, the RGI quark mass becomes negative, hence unphysical, for quark masses greater than ∼ 2m c and the lattice spacings at our disposal.
The idea is to consider ratios of G(1) at 2 consecutive heavy-strange meson masses m Hs(i+1) and m Hs(i) in the step-scaling mass chain: where K is the number of steps. Thus we have G lat (m Hs(i−1) ) is the ratio of successive G's on the given ensemble. By construction, G(1, m Ds ) is equal to 1. It will be a useful check that our numerical data obey that constraint. Taking this into account, we have G(1, m Bs ) = K−1 i=0 σ i G . Concerning the decay H s → D * s , a convenient framework is again HQET. Taking into account parities under C, P and T symmetries, the Lorentz structure decomposition of the hadronic matrix element D * and the normalisation equation of the polarisation vectors of a vector meson of mass m V and momentum k, According to our step-scaling procedure, these quantities can finally be expressed as: with

Computational setup
We have used in our analysis gauge field configuration ensembles generated by the CLS effort with N f = 2 non-perturbatively O(a) improved Wilson-Clover fermions [28,29] and the plaquette gauge action [30]. The parameter κ s was taken from [31] and κ c from [32]. They have been used in previous works eg. [27,33]. We have injected momenta to quarks by imposing isotropic twisted-boundary conditions in space, i.e. ψ(x+L e i ) = e iθ i ψ(x) with twisting angles θ i , i = 1, 2, 3. The kinematics for 2-pt. and 3-pt. correlation function is depicted in Figure 1, where strange, charm and heavy flavours are denoted by s, c and h, respectively. Only the charm quark carries the non-zero momentum. To reduce some of the O(a 2 ) effects we have averaged our results over positive and negative momenta. Six pairs ± θ, in addition to the zero-momentum point, help us to extrapolate the various form factors associated with H s → D s and H s → D * s from data available in the recoil region 1 ≤ w 1.1.
In Table 1 we collect the parameter specifications of the ensembles. Three lattice spacings a β=5.5 = 0.04831 (38)     masses are the same as in [27]. Statistical errors have been computed employing the "Γ-method" [35], based on estimating autocorrelation functions 1 . As in [27] all-to-all propagators were estimated stochastically with U (1) spin-diluted walltime noise sources. We have also reduced the contamination by excited states on 2-pt. and 3-pt. correlators by solving a 4 × 4 Generalized Eigenvalue Problem (GEVP) with one local and 3 Gaussian smeared interpolating fields. In the mass step-scaling strategy outlined above, we set K = 6. In Table 2 we collect the twisting angles used to inject momenta to the charm quark. The phase shift was applied isotropically in all spatial directions. Our analysis involves the following 2-pt. correlation functions evaluated on the gauge field ensembles at hand: where quark labels are as above, i and j specify smearing labels and spatial coordinates (summed over) are suppressed. The entering quark bilinears are defined as: P hs =ψ h γ 5 ψ s , P θ cs =ψ θ c γ 5 ψ s and V θ cs, k =ψ θ c γ k ψ s . Now we continue with the 3-pt. correlation functions: where the improved currents read ). Accordingly, we have to solve the three GEVPs Then we project the 2-pt. and 3-pt. correlation functions onto the generalised eigenvector chosen at a given time t fix , b ≡ v(t fix , t 0 ): (3.4) The asymptotic behaviour of the 2-pt. functions is known to be given bỹ Finally, the desired matrix elements may be computed from the large-time asymptotics of suitable ratios of the foregoing 2pt.-and 3pt.-correlation functions asC where the label (b) refers to bare hadronic matrix elements. Later, it will be convenient to note Upon mass dependent O(a) improvement and multiplicative renormalisation the physical matrix elements are

Extraction of hadronic matrix elements
As described in the previous section, we have solved 4 × 4 GEVP systems, except for the most chiral ensemble G8 for which we had to restrict the analysis to the 2 × 2 most smeared part of the correlators matrices, because the data are too noisy. We have set t 0 /a = 3 and t fix /a = 6. We vary the parameters t fix , t 0 and the time ranges used to fit the 2-pt. correlation functions and ratios of 3-pt. and 2-pt. correlation functions, in order to estimate the systematic errors on the hadronic matrix elements. In practice, we found that the alteration of the parameters t 0 and t fix does not influence the results significantly. Moreover we have observed that the GEVP solutions are in large part compatible with choosing the most smeared source for our interpolating fields. These findings are illustrated in Figure 2. The source-sink time separation t i is equal to T /2.
We have obtained hadron masses am Hs by fitting the plateau of the effective mass data am eff (t) coming from the generalised eigenvalues. 2 The fit range  [t min , t max ] is chosen after a visual inspection of the plateau. The data are strongly correlated, therefore the error of the plateau is not much smaller than the error of the individual data points. For this reason it was more important to exclude points outside the plateau than to maximise the available statistics. The same is done to evaluate the hadronic matrix elements kl along the formulae above.
In more detail, the extraction of the matrix elements involves the following elements 1 We apply symmetries, solve the GEVP and project the correlators following eq. (3.4).

Extrapolation to the physical point
Once we have obtained the hadronic matrix elements, we are in principle able to determine the heavy quark symmetry form factors h + , h − and h A 1,3 . Unfortunately the statistical quality of our data is not sufficient to reliably calculate h A 3 . Therefore we will restrict ourselves to the two former quantities at non zero recoil and on h A 1 (w = 1). We have tried to extrapolate h + and h − to w = 1 separately, and compute G(1) from it, as well as to extrapolate G(w) directly. The two extrapolations are mostly compatible, because the range in w is very small. We have performed linear and quadratic extrapolations in (w − 1). As quadratic fit coefficients are consistent with zero for almost all ensembles, we work with the result from the linear  extrapolations for the remainder of the analysis. In Figure 7 we show extrapolations of G to zero recoil for a selection of ensembles. The next-to-last step is to extrapolate Σ i G = G(1,a,m 2 π ,m Hs(i) ) G(1,a,m 2 π ,m Hs(i−1) ) to the physical point, where i denotes the mass step-scaling step corresponding to m Hs(i) = λ i m Ds . We have used the following fit ansatz: (4.1) The fit parameters and χ 2 /d.o.f. are collected in Table 3.
Parametrisations with additional terms were also studied. In particular, one might include a mistuning term 1 − m Hs λ i m Ds to account for the fact that the heavy meson masses are not tuned exactly equally on each ensemble.
However we decided to only include the terms in eq. (4.1)), because the data are so noisy that the fits can not reliably resolve these terms. In fact, as seen in Table  3, only very few non-constant terms are statistically significant. Cut off effects are limited to 5%, while a pion mass dependence is almost absent. Then, we can write Unfortunately, no clear trend in m Hs dependence is seen, as shown in the left   panel of Figure 8. That is why we have fitted the ratios σ i G by a constant σ G . We get σ G = 1.005 (23).
Moreover, we have checked, whether our data obey the constraint G(1, m Ds ) = 1. To do so, we have performed the following extrapolation: The left panel of Figure 9 shows that our continuum extrapolation G(1, m Ds ) is fully compatible with 1 within the limited statistics. As a final result for G we quote We have also performed a combined fit for the ratio using all data points simultaneously to explore the mass dependence and the effects of mistuning between the different ensembles:      where r mis is defined as am Hs λ i am Ds and m Hs denotes the physical meson mass. While the mistuning is not significant, Σ G 3 = 0.69(71), an inclusion of this term makes the fit less stable. The last term does contribute Σ G 4 = −0.0248(91), but the ratios (analogous to Figure 8) are still compatible with a constant. We obtain a final result for this method by computing Coincidently we get the same result as we get with the individual fits. The extrapolation of h A 1 to the physical point is completely analogous to that of G. We have decided to use individual fits as we did for the extrapolation of G. First we use the ansatz (1) as shown in the right panel of Figure 9. Then the ratios of successive mass-step-scaling steps are fitted with (4.8) The fit results are listed in Table 4. Again, the pion mass dependence is very weak, below 1%. It is less straightforward to conclude about cut off effects. They seem to be more significant than for Σ i G , as large as 12%. As for Σ i G , adding a mistuning term 1 − m Hs λ i m Ds has made the fits unstable. In the same way as before, from the ratios at the physical point and in the continuum limit, a fit to a constant σ h A 1 of the ratios σ i h A 1 is performed, because again no clear trend in m Hs is observed as shown in the right panel of Figure 8. The result for h A 1 can now be calculated according to eq. (2.24) (23)) 6 = 0.85 (16). (4.9)

Discussion
To our knowledge, only three lattice results for G Bs→Ds (1) are quoted in the literature. The ETM Collaboration, in their analysis of ensembles with N f = 2 twistedmass fermions, gets G Bs→Ds (1) = 1.052 (46), using the step-scaling in mass method defined through RGI quark masses [20]. The HPQCD Collaboration has analysed ensembles with N f = 2 + 1 staggered fermions, using the non-relativistic framework to regularise the b quark. The result reads G Bs→Ds (1) = 1.068(40) [21]. They have also analysed ensembles with N f = 2 + 1 + 1 staggered fermions, regularising the b quark with the heavy HISQ action [24]: they get G Bs→Ds (1) = 1.070(42) 3 . Our result is in the same ballpark as the three previous lattice computations, with a significantly larger error. 4 A possible explanation is that we have set a large source-sink time separation > ∼ 2 fm. This conservative choice helps to reduce the contamination from excited states with the caveat that the signal deteriorates faster for correlators build from heavy quarks.
Let us emphasize that a gain in statistics at each individual point of the stepscaling in mass strategy would have a significant impact to reduce the error on the final result, because the individual errors multiply with each other along the steps.
Three lattice results have been quoted in the literature for h (1) by HPQCD Collaboration at N f = 2 + 1 and N f = 2 + 1 + 1. Using non-relativistic b quark, again, they get h (1) = 0.906 (20), using the Fermilab action to regularize the b quark [39]. Two more recent preliminaries studies by this collaboration state h B→D * (1) leads to the conclusion that heavy quark mass dependence on this form factor is smaller than 10%, as it is for G Hs→Ds (1).

Conclusion
In the present paper we have reported an N f = 2 lattice QCD computation of the form factors G(1) and h A 1 (1) associated to the B s → D  from suitable correlators, we have clearly demonstrated that solving the GEVP is beneficial to reduce their contamination by excited states. Yet, only smearing D ( * ) s and B s interpolating fields is already a good starting point.
Step-scaling in masses helps to control cut-off effects originating from the heavy quark regularised on the lattice with the Wilson-Clover action. The twisted-boundary condition technique proves to be very helpful to extrapolate G(w) at zero recoil. It allows us to consider kinematics near-zero recoil so that a linear extrapolation in (w − 1) is fully satisfactory. As we were conservative in setting the source-sink time separation larger than 2 fm in 3pt. correlation functions, we end up with statistical errors ∼ 2.5% in the heavy-strange meson masses sequence of form factor ratios. This precision is not yet sufficient to be sensitive to the heavy quark mass dependence and thus also entails the still considerably large overall errors of the form factors from the present pilot study.
From the findings of our study, it can be deduced that a significant improvement in precision can be expected from a reduced, less conservative choice for the sourcesink time separation and focusing on the most smeared sources, on top of generally higher statistics of the individual ensembles of gauge field configurations that enter the analysis. However, rather than pursuing this for the two-flavour theory which is known to have only limited impact on precision phenomenology, we intend to take these aspects into account in a future calculation of B s → D ( * ) s form factors using lattice QCD ensembles with N f = 2 + 1 dynamical quarks in the medium term. In fact, thanks again to large-scale simulations carried out by the CLS (Coordinated Lattice Simulations) cooperation of researchers 5 , a rich landscape of such (2 + 1)flavour gauge field configuration ensembles is already available [42][43][44], which are to be regarded as natural successors to the ones investigated in the present work. These ensembles employ non-perturbatively O(a) improved Wilson-Clover fermions, together with the tree-level Symanzik-improved gauge action, and cover lattice spacings close to the continuum limit as well as the physical pion mass point. Since they also offer statistics substantially higher than underlying our N f = 2 analysis, a computationally feasible application of the methods studied (and its extensions suggested) here to the N f = 2 + 1 case, including additional form factors of the B s system, now appears realistic. Therefore, even though the cost of large-volume calculations via the step-scaling in mass strategy remains challenging, an outcome of phenomenologically interesting target precision of a few percent may be accomplished from this approach on a reasonable time scale.
Finally, let us also mention the recently proposed, so-called stabilised Wilson fermions as a new avenue for QCD calculations with (improved) Wilson-type fermions [45]. First simulations of lattice QCD with 2 + 1 dynamical quark flavours within this formulation [46], which amounts to replace the original clover term with an exponentiated version of it, already hint at improved properties such as, apart from stabilising simulations in large volumes towards small pion masses, an overall good continuum scaling and milder relative quark mass dependent cutoff effects. Hence, this framework brings into sight another promising direction for estimating strong interaction effects in precision weak interaction phenomenology through lattice QCD that we are also envisaging for future applications of the techniques explored in the present study.      [8,27] 1.23825(51) 0.100172 [7,23] 1.46300(53) 0.088934 [8,20] 1.72617(57)        [7,27] 1.08818(74) 0.10593 [8,22] 1.28429(94) 0.096211 [5,18] 1.51684(82)