$B \to D^*\ell\nu_\ell$ semileptonic form factors from lattice QCD with M\"obius domain-wall quarks

We calculate the form factors for the $B \to D^*\ell\nu_\ell$ decay in 2+1 flavor lattice QCD. For all quark flavors, we employ the M\"obius domain-wall action, which preserves chiral symmetry to a good precision. Our gauge ensembles are generated at three lattice cutoffs $a^{-1} \sim 2.5$, 3.6 and 4.5 GeV with pion masses as low as $M_\pi \sim 230$ MeV. The physical lattice size $L$ satisfies the condition $M_\pi L \geq 4$ to control finite volume effects (FVEs), while we simulate a smaller size at the smallest $M_\pi$ to directly examine FVEs. The bottom quark masses are chosen in a range from the physical charm quark mass to $0.7 a^{-1}$ to control discretization effects. We extrapolate the form factors to the continuum limit and physical quark masses based on heavy meson chiral perturbation theory at next-to-leading order. Then the recoil parameter dependence is parametrized using a model independent form leading to our estimate of the decay rate ratio between the tau ($\ell = \tau$) and light lepton ($\ell = e,\mu$) channels $R(D^*) = 0.252(22)$ in the Standard Model. A simultaneous fit with recent data from the Belle experiment yields $|V_{cb}| = 39.19(91)\times 10^{-3}$, which is consistent with previous exclusive determinations, and shows good consistency in the kinematical distribution of the differential decay rate between the lattice and experimental data.

R(D * ) = 0.252 (22) in the Standard Model.A simultaneous fit with recent data from the Belle experiment yields |V cb | = 39.19(91)×10−3 , which is consistent with previous exclusive determinations, and shows good consistency in the kinematical distribution of the differential decay rate between the lattice and experimental data.

I. INTRODUCTION
The B → D * ℓν ℓ semileptonic decay, where ℓ = e, µ, τ and ν ℓ represents the corresponding neutrino, plays a key role in stringent tests of the Standard Model (SM) and searches for new physics.The channels associated with light leptons ℓ = e, µ provide a determination of the Cabibbo-Kobayashi-Maskawa (CKM) matrix element |V cb |, which is a fundamental parameter of the SM.On the other hand, the τ channel is expected to be a good probe of new physics, since, for instance, it is expected to strongly couple to the charged Higgs boson predicted by supersymmetric models.Indeed, there is an intriguing ≳ 3 σ tension between the SM and experiments on the decay rate ratio describing the lepton-flavor universality violation [1].
However, there has been a long-standing tension between the |V cb | determinations from exclusive and inclusive decays (B → Dℓν ℓ + D * ℓν ℓ + D * * ℓν ℓ + • • • ).The Heavy Flavor Averaging Group reported that analysis of the inclusive decay yields [1] |V cb | × 10 3 = 41.98 (45), (2) which shows a ≲ 4 σ (9 %) tension with the determination from the exclusive decays ( It has been argued [2] that such a tension can be explained by introducing a tensor interaction beyond the SM, which, however, largely distorts the Z → bb decay rate measured precisely by experiments.Therefore, it is likely that the |V cb | tension is due to our incomplete understanding of the theoretical and/or experimental uncertainty.
The largest theoretical uncertainty in the exclusive determination comes from form factors, which describe non-perturbative QCD effects to the decay amplitude through the matrix elements where w = vv ′ is the recoil parameter defined by four velocities v = p/M B and v ′ = p ′ /M D ( * ) , and ϵ ′ is the polarization vector of D * , which satisfies p ′ ϵ ′ = 0. Note that, here and in the following, kinematical variables (momentum, velocity, polarization vector, but not space-time coordinates) with and without the symbol "′" represent those for D * and B, respectively.
Lattice QCD is a powerful method to provide a first-principles calculation of the form factors with systematically improvable accuracy [3,4].However, until recently, only h A 1 (1) at the zero-recoil limit w = 1 had been calculated in unquenched lattice QCD [5][6][7].In the previous determination of |V cb |, therefore, other information of the form factors in addition to |V cb | had to be determined by fitting experimental data of the differential decay rate to a theoretical expression [1,[8][9][10].While a model independent parametrization of the form factors by Boyd, Grinstein and Lebed (the BGL parametrization) is available [11], its large number of parameters due to the poor knowledge on the form factors makes the fit unstable [12][13][14][15] even with experimental data of differential distribution with respect to all kinematical variables, namely the recoil parameter and three decay angles, from the KEKB/Belle experiment [16, 17] [80].
One of the largest problems in the lattice study of B meson physics is that the simulation cost rapidly increases as the lattice spacing decreases.In spite of recent progress in computer performance and development of simulation algorithms, it is still difficult to simulate lattice cutoffs a −1 sufficiently larger than the physical bottom quark mass m b,phys to suppress O((am b ) n ) discretization errors to, say, the 1 % level.For the time being, a practical approach to B physics on the lattice is to employ a heavy quark action based on an effective field theory, such as the heavy quark effective theory (HQET), to directly simulate m b,phys at currently available cutoffs a −1 ≲ m b,phys .Another straightforward approach is the so-called relativistic approach, where a lattice QCD action is used also for heavy quarks but with their masses sufficiently smaller than the lattice cutoff to suppress discretization effects.The Fermilab/MILC reported the first lattice calculation of all form factors at zero and non-zero recoils [18] using the Fermilab approach for heavy quarks, namely a HQET interpretation of the anisotropic Wilson-clover quark action [19], and staggered light quarks.
It was recently followed by the HPQCD Collaboration with a fully relativistic approach [20], where the highly improved staggered quark action was used both for light and heavy flavors.
In this paper, we present an independent study with a fully relativistic approach using the Möbius domain-wall quark action [21] for all relevant quark flavors.This preserves chiral symmetry to good precision at finite lattice spacing, reducing the leading discretization errors to second order, i.e.O((am b ) 2 ).We can also use chiral perturbation theory in the continuum limit (a = 0) to guide the chiral extrapolation to the physical pion mass.Note also that we do not need explicit renormalization of the weak axial and vector currents on the lattice to calculate relevant form factors, because the renormalization constants of these currents coincide with each other and are canceled by taking appropriate ratios of relevant correlation functions [22].Our preliminary analyses and discussions on these features can be found in our earlier reports [23][24][25].
This paper is organized as follows.After a brief introduction to our choice of simulation parameters, the form factors are extracted from correlation functions at the simulation points in Sec.II.These results are extrapolated to the continuum limit and physical quark masses in Sec.III.In Sec.IV, we generate synthetic data of the form factors and fit them to a model-independent parametrization in terms of the recoil parameter to make a comparison with the previous lattice and phenomenological studies as well as to estimate the SM value of R(D * ).We also perform a simultaneous fit with experimental data to estimate |V cb |.Our conclusions are given in Sec.V.

II. FORM FACTORS AT SIMULATION POINTS
A. Gauge ensembles Our gauge ensembles are generated for 2+1-flavor lattice QCD, where dynamical up and down quarks are degenerate and dynamical strange quark has a distinct mass.We take three lattice spacings of a ≃ 0.044, 0.055 and 0.080 fm [26], which correspond to the lattice cutoffs a −1 ∼ 2.453(4), 3.610(9) and 4.496(9) GeV, respectively.They are determined using the Yang-Mills gradient flow [27].The tree-level improved Symanzik gauge action [28][29][30] and the Möbius domain-wall quark action [21] are used to control discretization errors and to preserve chiral symmetry to a good precision.We refer the interested reader to Refs.[21,26,31] for the five dimensional formulation of the quark action and our implementation.Its four-dimensional effective Dirac operator is given by where m q represents the quark mass.We employ the kernel operator , where D W (−M ) is the Wilson-Dirac operator with a negative mass with M = 1 and with stout smearing [32] applied three successive times to the gauge links.With our choice of the staple weight (ρ = 0.1), the smearing radius ∼ 1.5 a is at the scale of the lattice spacing and vanishes in the limit of a → 0. We therefore expect that the link smearing does not change the continuum limit of the B → D * form factors.It may induce additional discretization effects, which, however, do not turn out to be large in our continuum and chiral extrapolation in Sec.III.The polar decomposition approximation in Eq. ( 6) is realized for the sign function ϵ [H M ] in the five-dimensional domain-wall implementation.With this choice, the residual quark mass, which is a measure of the chiral symmetry violation, is suppressed to the level of 1 MeV at a ≃ 0.080 fm, and 0.2 MeV or less on finer lattices, which are significantly smaller than the physical up and down quark masses [26].
With reasonably small values of the fifth-dimensional size N 5 = 8 -12, the computational cost is largely reduced from that of our previous large-scale simulation [33,34] using the overlap formulation that also preserves chiral symmetry [35].
We employ the Möbius domain-wall action for all relevant quark flavors.Our choices of   [36,37] for our chiral extrapolation of simulation results to the physical pion mass without introducing terms to describe discretization effects.We take a strange quark mass m s close to its physical value fixed from The charm quark mass m c is set to its physical value fixed from the spin-averaged charmonium mass (M ηc +3M J/Ψ )/4 [38].
Depending on the lattice spacing a, we take three to six bottom quark masses m to avoid large discretization errors.We note that, with chiral symmetry, the leading discretization error is reduced to O((am Q ) 2 ).
As we will see in the next section, a and m Q dependences of the form factors are mild, and the extrapolation to the continuum limit and physical quark masses is under good control.
Our simulation parameters are summarized in Tables I and II.
The spatial lattice size L = N s a is chosen to satisfy the condition M π L ≳ 4 to suppress finite volume effects (FVEs).At our smallest M π on the coarsest lattice, we also simulate a smaller physical lattice size M π L ∼ 3 to directly examine the FVEs.
The statistics at the lattice cutoffs a −1 ≃ 2.5, 3.6 and 4.5 GeV, are 5,000, 2,500 and 2,500 Hybrid Monte Carlo trajectories with a unit trajectory length of 1, 2 and 4 times the molecular dynamics time unit, respectively, which is increased to take account of the longer auto correlation toward the continuum limit.On each ensemble, we take N conf configurations in an equal trajectory interval to calculate correlation functions relevant to the B → D * ℓν decay.Our choice of N conf is listed in Table I.More details on our simulation setup can be found in Ref. [26].

B. Ratio method
The B → D * matrix elements can be extracted from the three-and two-point functions through their ground state contribution where O Γ is the weak vector (V µ ) or axial (A µ ) current, and we omit the symbol "′" on the which are used in a correlator ratio (13) below.
The three-point functions are calculated by using the sequential source method, where the total temporal separation ∆t+∆t ′ between the source and sink operators is fixed and the temporal location of the weak current is varied.In order to control the contamination from the excited states, we repeat our measurement for four values of the source-sink separation ∆t + ∆t ′ listed in Table I.
The statistical accuracy of the two-and three-point functions are improved by averaging over the spatial location of the source operator x src as indicated in Eqs.(7) and (8).To this end, we employ a volume source with Z 2 noise.At M π ≲ 300 MeV, we repeat our measurement over two values of the source time-slices t src = 0 and N t /2 to take the average of the correlators.Namely N tsrc in Eqs. ( 7) and ( 8) is 2 at M π ≲ 300 MeV, and 1 otherwise.
The correlators with non-zero momentum p are also averaged over all possible p's based on parity and rotational symmetries on the lattice.
The number of measurements on each ensemble is given as N conf N tsrc in Table I.As mentioned above, correlation functions for each configuration are averaged over N tsrc values of the source time-slice t src .Then, simulation data of N conf configurations are divided into 50 bins: namely, the bin size is two (one) configurations for β = 4.17 (4.35 and 4.47).The statistical error is estimated by the bootstrap method with 5,000 replicas.Figure 1 shows the bin size dependence of the relative statistical error of a three-point function at our largest cutoff, where the topological charge Q changes much less frequently leading to larger autocorrelation compared to the coarser lattices [26].Our choice of the bin size on this finest lattice is one configuration, and we do not observe significant increase of the statistical error toward larger bin sizes.This suggests that our bin size is sufficiently large partly due to our choice of larger unit trajectory length towards the continuum limit.We also note that the error estimate is stable when we vary the number of bootstrap replicas as well as when we employ the jackknife method instead.
As observed in our previous simulation in the trivial topological sector [33], the local topological fluctuation are active even if we fix the "global" topological charge Q, namely the whole-volume integral of the topological charge density.In Ref. [39], we demonstrated that the topological susceptibility χ t is calculable within the fixed topology setup.Our result for χ t shows quark mass dependence consistent with ChPT, and its value extrapolated to the chiral limit is in good agreement with that from our simulation in the ϵ-regime [40,41].
References [42,43] argued that the bias due to the freezing of the global topology can be considered as FVEs suppressed by the inverse space-time volume.In our study of the B → πℓν decay [26], we confirmed that such FVEs on the pion effective mass are well suppressed, so that its |Q| dependence is not significant.Figure 2 shows the statistical average of a B → D * three-point function calculated for each |Q|.We do not observe any significant |Q| dependence of the average at |Q| ≤ 4. At larger |Q|, we do not have enough data for a reliable error estimate, but individual data plotted by crosses are consistent with the averages at |Q| ≤ 4. We therefore conclude that the topology freezing effect is insignificant within our statistical accuracy.
To precisely extract the form factors, we construct ratios of the correlation functions, in which unnecessary overlap factors and exponential damping factors cancel for the groundstate contribution [22].The statistical fluctuation is also expected to partly cancel in such a ratio.With B and D * mesons at rest, the vector current matrix element (4) vanishes, and the axial vector one ( 5) is sensitive only to the axial vector form factor h A 1 at zero recoil, which is the fundamental input to determine |V cb |.In order to precisely determine h A 1 (1), we employ a double ratio where the D * polarization is chosen as ϵ ′ = (1, 0, 0, 0) along the polarization of the current inserted, i.e.A 1 .The left panel of Fig. 3 shows our results for this ratio at our smallest We find that all data are well described by the following fitting form including effects from the first excited states with χ 2 /d.o.f.≲ 1 for a wide range of the values of ∆t and ∆t ′ .Here ∆M B(D * ) represents the energy difference between the B(D * ) meson ground state and the first excited state of the same quantum numbers, and is set to the value estimated from a two-exponential fit to the two-point function C B(D * ) (∆t; 0).The fit range is chosen such that all the data of R 1 (∆t, ∆t ′ ) with the source-to-current (current-to-sink) separation ∆t (′) equal to or larger than a lower cut ∆t min are included irrespective of the source-sink separation ∆t + ∆t ′ .As shown in the left panel of Fig. 4, our result for h A 1 (1) is stable against the choice of the lower cuts ∆t min and ∆t ′ min .
The right panels of Figs. 3 and 4 for our largest a −1 also show ground state dominance with ∆t + ∆t ′ ≳ 23a ∼ 1 fm and stability of h A 1 (1) against the choice of the fit range, respectively.The situation is similar for other simulation parameters.In order to extract form factors at non-zero recoils, the D * momentum and polarization vector need to be chosen appropriately.The axial vector matrix element ( 5) is sensitive only to h A 1 with a polarization vector ϵ ′ = (1, 0, 0) and a momentum p ′ ⊥ that satisfies ϵ ′ * v = 0. We note that, here and in the following, three dimensional polarization vector ϵ is accompanied by its temporal component to be fixed from the convention ϵ ′ * p ′ ⊥ = 0, which is assumed for Eqs. ( 4) and (5).Then, the recoil parameter dependence of h A 1 may be studied from the following ratio Here we use the two-point function with the local sink C D * sl (10) as in our study of the K → πℓν ℓ decay [44].The vector form factor h V can be extracted from the following ratio through the vector current matrix element (4) Here we use two different polarization vectors, ϵ ′ = (1, 0, 0) and ϵ ′′ = (0, 1, 0), and momenta to share the same recoil parameter w, and v ′′ ⊥j = p ′′ ⊥j /M D * .We emphasize that renormalization factors of the weak vector and axial currents cancel thanks to chiral symmetry preserved in our simulations.
in these ratios among all simulated values of the source-sink separation ∆t + ∆t ′ ≳ 1 fm.A similar ground state dominance is also observed at other simulation points.We carry out a simultaneous fit using a fitting form that takes account of the first excited state contribution as to extract the form factor ratio h A 1 (w)/h A 1 (1) from the ground state contribution R ′ 1 (p ′ ⊥ ), and a similar form for R V (p ′′ ⊥ , p ′ ⊥ ; ∆t, ∆t ′ ) to extract h V (w)/h A 1 (w).Our data are well described by this fitting form with χ 2 /d.o.f.≲ 1.From these form factor ratios and h A 1 (1) from R 1 , we calculate h A 1 (w) and h V (w) at simulated values of w.
There is a 2 σ bump of R ′ 1 for ∆t+∆t ′ = 40 (black triangles) around ∆t = 28 (∆t ′ = 12) in Fig. 5. Since the global topological charge does not fluctuate very much at the corresponding cutoff, one may expect bump(s) in the exponential decay of relevant correlators due to instantons frozen at the temporal separation ∆t (′) for a certain period of our simulation time.However, Fig. 7 shows that there is no statistically significant bump of the relevant correlators at ∆t ∼ 28 and ∆t ′ ∼ 12.As mentioned above, the local topological fluctuation is active even when the global topological charge is fixed.We attribute 1 -2 σ bumps of correlator ratios in Figs. 3,   The axial vector form factors h A 2 and h A 3 can be extracted from the axial vector matrix element (5) with the D * spatial momentum p ′ ̸ ⊥ not perpendicular to the spatial polarization vector.The matrix element of A 1 , for instance, has non-zero sensitivity to h A 1 and h A 3 .We use the following ratio to extract h where ) appears, since the overlap factor of D * depends on whether the polarization vector is perpendicular to the momentum, even if While r Z can be estimated from individual fit to the relevant two-point functions, we employ a ratio  the numerator is also sensitive to We plot our results for R 3 and R 2 on our finest lattice in Figs. 8 and 9, respectively.
The double ratio R 1 enables us to calculate h A 1 (1) with 1 % statistical error or better, as it only involves correlators with zero momentum and the statistical fluctuation partially cancels between the numerator and denominator.At non-zero recoil, the numerator of , and the statistical error of h A 1 and h V is typically a few %.This is not the case for h A 2 and h A 3 with our setup where the B meson is at rest.
The numerator of R 3 involves h A 1 and h A 3 , and the typical precision for h A 3 is 10 -20 %.
The statistical error increases to ≳ 40 % for h A 2 , since the central value is suppressed by heavy quark symmetry and R 2 is a linear combination of all axial form factors.
While discretization effects to the wave function renormalization are not small [45], these also cancel in the ratios.

III. EXTRAPOLATION TO CONTINUUM LIMIT AND PHYSICAL QUARK MASSES A. Choice of fitting form
Our results for the form factors obtained at various simulation parameters are plotted as a function of the recoil parameter w in Fig. 10.We find that the results with different heavy and light quark masses are all close to each other.It also turns out that the discretization effects are not large.
We extrapolate these results to the continuum limit and physical quark masses based on next-to-leading order heavy meson chiral perturbation theory (NLO HMChPT) [36,37].The extrapolation form is generically written as which comprises the constant term, chiral logarithm at w = 1 and leading corrections in small expansion parameters with  ), we use its explicit form in Ref. [37].For the D * Dπ coupling, we take a value of g D * Dπ = 0.53 (8), which was employed in the Fermilab/MILC paper [6] and covers previous lattice estimates [46][47][48][49][50].This choice does not lead to a large systematic uncertainty, because the chiral logarithm is suppressed by the small mass splitting squared ∆ 2 D .
In Fig. 11, we plot h A 1 and h V at our smallest value of w as a function of M 2 π .The M 2 π dependence is rather mild, possibly because there is no valence pion in this decay.We note that there is a non-trivial concave structure at small M 2 π due to the opening of the D * → Dπ decay.The effect estimated in NLO HMChPT is, however, not large compared to our statistical accuracy.
For our chiral extrapolation, we employ the so-called ξ expansion in ξ π and ξ ηs : namely, m s , respectively, as well as quark-mass dependent f π for Λ χ instead of the decay constant f in the chiral limit.This is motivated by our observation [51] that the ξ-expansion shows better convergence of the chiral expansion for light meson observables compared to the x-expansion with x = 2Bm q /(4πf ) 2 , where B gives the lowest order relation M 2 π = 2Bm ud .
The extrapolation form explicitly includes the one-loop radiative correction η X for the matching of the heavy-heavy currents between QCD and HQET [52][53][54].Remaining heavy quark mass dependence of the form factors is then parametrized by polynomials of The dependence is mild as shown in Fig. 12.Our results are well described by the linear terms c Q ϵ Q and c am Q ξ am Q in Eq. ( 19).We note that, for h A 1 , the term c Q ϵ Q is modified as c Q (w − 1)ϵ Q to be consistent with Luke's theorem [55], and we add the quadratic The discretization effects in the form factors also turn out to be mild with our choice of parameters, a −1 ≳ 2.5 GeV and m Q ≤ 0.7a −1 .These are described by an expansion, i.e.
linear in the m Q independent and dependent parameters ξ a and ξ am Q , respectively.Figure 12 shows that the m Q dependent discretization error is significant in our results.However, the net m Q dependence is described reasonably well with our fitting form (19).We note that the m Q dependence for h A 2 and h A 3 is insignificant due to their large relative uncertainty.As we will mention later, fit results do not change significantly if we exclude data at 0.5 ≤ am Q ≤ 0.7 in our continuum and chiral extrapolation.In this study, we first extrapolate the form factors to the continuum limit and physical quark masses, and then parametrize the w dependence by using the model independent BGL parametrization, because it may not be simply combined with the chiral and heavy quark expansions.Concerning the w dependence, therefore, the fit form (19) with an expansion in w − 1 is used to interpolate our results.Figure 10 shows that our data do not show strong curvature in the simulation region of 1 ≤ w ≲ 1.1, and can be well described with only up to the quadratic term.
Numerical results of our continuum and chiral extrapolation are summarized in Table III.
The polynomial coefficients turn out to be O(1) or less with our choice of the dimensionless expansion parameters (20), which are normalized by appropriate scales.Many of them are consistent with zero reflecting the mild dependence on the lattice spacing and quark masses mentioned above.As a result, our choice of the fitting forms (19) results in small values of χ 2 /d.o.f.≲ 0.2, which, however, are only a rough measure of the goodness of the fits, because we employ approximated estimators of the covariance matrix as discussed in the next subsection.

Correlated fit
With limited statistics, the covariance matrix C of our form factor data has exceptionally small eigenvalues with large statistical error.In this study, we test two methods to take the correlation into account in the combined continuum and chiral extrapolation.First, as in the study of B → πℓν ℓ [26], we introduce a lower cut of the eigenvalue (singular value) λ cut .
Effects of the exceptionally small eigenvalues are suppressed by replacing eigenvalues smaller than λ cut by λ cut in the singular value decomposition of C as where u k represents the eigenvector corresponding to the eigenvalue λ k .Here eigenvalues are sorted in ascending order (λ k < λ k+1 ) and n cut satisfies λ ncut < λ cut < λ ncut+1 .We choose λ cut such that all eigenvalues below λ cut have statistical error larger than 100 %.The second method is the so-called shrinkage estimate of the covariance matrix [56,57], which amounts to replacing C ij by With 0 ≤ ρ ≤ 1, the shrinkage is an interpolation between the uncorrelated (ρ = 1) and correlated (ρ = 0) fits.An optimal value of ρ can be estimated by minimizing a loss function [57].
We observed, however, that results of the continuum and chiral extrapolation do not strongly depend on the choice of λ cut and ρ.In this study, we therefore employ the singular value cut for our main results, since this is a conservative approach, which only tends to increase the final error.The shift of the fit results obtained with the shrinkage estimator is treated as a systematic error due to the estimate of C.
With the approximated estimators of C, χ 2 /d.o.f. is only a rough measure of the goodness of the fit.The impact is shown in Fig. 13.We observe that, in the case of our analysis, χ 2 /d.o.f.tends to be slightly smaller than that for the unquenched fit, and rather stable against the choice of the estimator (singular value cut, or shrinkage) and parameter ρ.While C is not invertible for the correlated fit (with ρ = 0 for the shrinkage estimator), we may expect that χ 2 /d.o.f.does not change drastically and serves as a "rough" measure of the goodness of the fit even with the approximated estimators of C. We also note that, in this paper, the approximated estimators are used only for the continuum and chiral extrapolation of form factors due to the large size of C with different choices of the quark masses (including m Q for the relativistic approach), lattice cutoff and recoil parameter.shows that with the shrinkage estimator as a function of the parameter ρ in Eq. ( 22).The error is estimated by the bootstrap method.

Fitting form
The highest order is quadratic in the recoil parameter w − 1 expansion and heavy quark expansion for h A 1 , and linear for other expansions including chiral and finite a expansions.
The systematic error due to this choice is estimated by repeating the continuum and chiral extrapolation without one of the highest order terms, if the corresponding coefficient is poorly determined, otherwise by adding a yet higher order term.
Possible higher order corrections in the chiral expansion are also examined by testing the x-expansion.We also test a multiplicative form to examine the potential effect of cross terms.We find that these alternative fits yield a change of the fit results well below the statistical error.

Inputs
The lattice cutoff a −1 , physical meson masses M π , M K , M η b , and the D * Dπ coupling g D * Dπ are inputs to the continuum and chiral extrapolation.Since the masses have been measured very precisely, we take account of the uncertainty due to a −1 and g D * Dπ by repeating the extrapolation with each input shifted by ± 1 σ.
We set M π and M K to those in the isospin limit as in our determination of a [26].To examine the isospin breaking effects to the form factors, we test the continuum and chiral extrapolation to M π = M π 0 and M π ± as well as M K = M K 0 and M K ± .We also test the heavy quark expansion parameter replaced with ϵ Q = Λ/2M B and compare the extrapolations using The difference in the form factor can be considered as an estimate of the strong isospin breaking effect.It turns out to be small as suggested by the mild quark mass dependence of the form factors.We include this in the systematic uncertainty so that our synthetic data for the form factors can be used for both the neutral and charged B meson decays.

Finite volume effects (FVEs)
Our spatial lattice size satisfies M π L ≥ 4. At our smallest pion mass on the coarsest lattice, we also simulate a smaller size L ≃ 3.0M −1 π to directly examine FVEs.Comparison between the two volumes at that simulation point does not show any significant deviation: for instance, h A 1 at w = 1 is consistent within statistical error of ≲ 1 %.We also repeat the continuum and chiral extrapolation by replacing the form factors at the smallest M π with those on the smaller L. The change in the fit results is included in the systematic uncertainty, but turns out to be well below the statistical uncertainty.We conclude that finite volume effects are not significant in our simulation region of M π ≳ 230 MeV.

Accuracy of continuum and chiral extrapolation
Figure 14 shows the estimate of various uncertainties for the form factors in the continuum limit and at the physical quark masses as a function of w.As discussed in Sec.II, the axial (vector) matrix element is exclusively sensitive to h A 1 (w) (h V (w)) with an appropriate choice of the D * meson polarization vector and momentum.Combined with a precise determination at zero-recoil, the total uncertainty of h A 1 is about 2 % in our simulation region of w, and the largest uncertainty comes from statistics and discretization effects.On the other hand, either D * or B needs to have non-zero momentum to determine h V , which is, therefore, less accurate.The total accuracy is ≲ 7 % with 3 -4 % statistical and ∼ 5 % discretization errors.
Since the simulated bottom quark masses are not far below the lattice cut-off (m Q ≤ 0.7a −1 ), we repeat the continuum and chiral extrapolation only using the data with m Q ≤ 0.5 a −1 .The change in the fit results is well below the statistical uncertainty.We therefore conclude that our continuum extrapolation in the relativistic approach is controlled.as in Fig. 10.The green bands show these ratios in the continuum limit and at physical quark masses reproduced by the continuum and chiral extrapolation of form factors (19).
On the other hand, the axial form factors, h A 2 and h A 3 , have larger uncertainty, i.e.
40 % and 15 %, respectively, dominated by the statistical error.Towards a more precise determination, it would be advantageous to simulate moving B mesons so as to have matrix elements more sensitive to these form factors.

Previous determinations of |V cb | often employed a HQET-based parametrization by
Caprini, Lellouch and Neubert of h A 1 and form factor ratios where higher order HQET corrections are expected to be partially canceled [58].Figure 15 shows that these ratios also have mild parameter dependences and can be safely extrapolated to the continuum limit and physical quark masses.We, however, do not employ the modeldependent CLN parametrization in the following analysis.

IV. FORM FACTORS AS A FUNCTION OF THE RECOIL-PARAMETER
A. Fit to lattice data From the continuum and chiral extrapolation in the previous section, we generate synthetic data of the form factors in the relativistic convention, with which the BGL parametrization was originally proposed.A kinematic invariant k = (t + − q 2 ) (t − − q 2 ) /4t is defined with the maximum value of the momentum transfer q 2 max = t − = (M B − M D * ) 2 , which corresponds to the zero recoil limit, and the threshold t + = (M B + M D * ) 2 for the W → BD * production channel.Instead of a + and a − , it is common to use linear combinations These are related to the form factors in the HQET convention as We note that there are two kinematical constraints to be satisfied at the minimum (w = 1) and maximum values of w (w max = (1 + r 2 )/2r with m ℓ = 0) With these form factors, we can write the helicity amplitudes for a given helicity of the intermediate W boson, and the differential decay rate with respect to w in simple forms H 0 (w) = 1 and From the continuum and chiral extrapolation presented in Sec.III, synthetic data of f , g and F 1,2 are calculated at reference values of the recoil parameter w ref through the relations ( 29) - (32).We take w ref = 1.025, 1.060 and 1.100, which span the whole simulated region of w, since our continuum and chiral extrapolation interpolates h A {1,2,3} and h V in w.Three values are chosen, because the statistical covariance matrix of the synthetic data develops ill-determined eigenvalues with four or more values of w ref .This might be partly because, in our continuum and chiral extrapolation, we parametrize the w dependence of each form factor using a quadratic form with three independent parameters.The numerical values together with their covariance matrix are presented in Table IV.
The synthetic data are fitted to the model-independent BGL parametrization where the z parameter maps the whole semileptonic kinematical region 0 a small parameter region of z, 0 ≤ z ≲ 0.06.
The Blaschke factors TABLE IV: Synthetic data of form factors f , F 1,2 and g in the relativistic convention.The central value and total error are in the second column, and the correlation matrix is in the third to fourteenth columns.
. GeV units.The input-A is a compilation of slightly old experimental measurement [59] and lattice QCD [60,61] and phenomenological [62,63] studies.A different phenomenological estimate [64] for the (axial) vector resonances is used for input-B, where the pseudo-scalar masses are set from a recent edition of the Particle Data Book [65] (n = 0, 1) and another phenomenological study [66] (n = 2).We employ input-A in our main analysis, whereas input-B is used to test the stability of the BGL parametrization against the choice of the input.See the text for more details.
. factor out the pole singularities outside the semileptonic region but below the threshold t + due to the resonances for a given channel J P of bc mesons.The explicit form is given as where is the z-parameter corresponding to the k-th resonance (in term of q 2 rather than w).In our main analysis, we employ resonance masses denoted as "input-A" in Table V.This choice consists of slightly old experimental measurement [59] (M pole,{0,1} for J P = 0 − ) and estimates based on lattice QCD [60,61] (M pole,0 for 1 + and M pole,{0,1} for 1 − ) and quark potential models [62,63] (other M pole,n 's).This input has been used in some of previous phenomenological analysis using the BGL parametrization [9,14] and lattice calculations of form factors [18,20].Another choice ("input-B" in Table V) in the literature is used to estimate the systematics uncertainty due to the use of input-A (see below).
The outer function ϕ F (z) is chosen such that a constraint on the expansion coefficients  45) - (48).
Perturbative estimate (upper line) is employed in our main analysis.
. [12,68,69]  derived from unitarity of the theory has a simple form The explicit form is given as where n I = 2.6 is the number of the spectator quarks with a correction due to SU(3) breaking.

By χ
{T,L} J P (0), we denote the derivative of the vacuum polarization function of the weak current at zero momentum q 2 = 0 for a given channel J P and polarization {T, L}.As in the literature, we replace χ T 1 − by χT 1 − , from which the one-particle state contribution due to B * c is subtracted for a better saturation of the unitarity bound (44).While a lattice QCD estimate is available [67], we employ a perturbative estimate [12,68,69] summarized in Table VI.
We note that our choice of the resonance masses (input-A) and derivatives are the same as previous lattice studies to allow a direct comparison of the BGL parametrization with them.
It is known that the unitarity bound (44) is not well saturated by the B → D * ℓν ℓ channel, and hence is rather weak.In this study, we do not impose this bound, but confirm that fit results satisfy it within the error.
Figure 16 shows the regular part of the axial form factor P f (z)ϕ f (z)f (z), which is to be expanded in z.We observe rather mild z dependence of all form factors.While there is no clear sign of curvature in our simulation region, we employ a quadratic form N f = N g = N F 1 = N F 2 = 2 to safely suppress the truncation error.Among twelve coefficients, constants a F 1 ,0 and a F 2 ,0 are fixed to fulfill the kinematical constraints (33) at w = 1, where z = 0 and the linear and quadratic terms in the z-parameter expansion (39) vanish, and (34) at w max , where z ∼ 0.05 is well below unity, and hence the constant term gives rise to a large contribution to F 2 .To be consistent with the unitarity bound (44), we impose the bound for all form factors F = g, f, F 1 and F 2 .Table VII shows fit results, which describe the synthetic data with an acceptable value of χ 2 /d.o.f.= 0.51.
We test the stability of the fit with different choices of the input and reference values of w.
Input-B in Table V of the resonance masses for the Blaschke factors consists of the pseudoscalar masses M pole,{0,1} from recent experimental measurements [65], and other resonance masses from phenomenological studies different from input-A: namely, Ref. [66] for J P = 0 − and n = 2, and Ref. [64] for J P = 1 {+,−} .We note that these vector and axial masses have been used in a phenomenological analysis [10] and Belle analyses [16,17] of the B → D * ℓν decay.For the derivatives of the vacuum polarization functions χ T,L J P , we test the lattice QCD estimate in Ref. [67].Shift of the expansion coefficients a {g,f,F 1 ,F 2 },{0,1} in units of their uncertainty is summarized in Table VIII, where we test four choices of the input: namely, i) vector and axial resonance masses from Ref. [64], ii) ground and first excited pseudo-scalar masses from Ref. [65], iii) second excited pseudo-scalar mass from Ref. [66], and iv) lattice estimate of χ T,L J P .The shift of a {g,f,F 1 ,F 2 },{0,1} is at the 2-σ level or typically well below it.A large shift (∼ 4σ) in a f,0 with the lattice QCD estimate of χ T,L J P can be attributed to about 20% shift of χ T 1 + for the axial channel.As suggested in Eq. ( 46), this simply leads to rescaling    VII.We list their shift in units of their uncertainty due to different choices of input and reference values w ref .We omit a F 1 ,0 , which is proportional to a f,0 due to the kinematical constraint (33), as well as a {g,f,F 1 ,F 2 },2 for quadratic terms, which show tiny shift due to their large relative uncertainty.Note that a g,n is independent of the pseudo-scalar resonance mass input.The maximum (minimum) value of w ref is denoted as w ref,max (w ref,min ).The quadratic fit (N {f,g,F 1 ,F 2 } = 2) is impossible with two values of w ref , with which we compare coefficients of the linear fit (N {f,g,F 1 ,F 2 } = 1).
.  Table VIII shows that the shift of the coefficients a {g,f,F 1 ,F 2 },{0,1} is well below 1 σ.Their relative accuracy is slightly worse (typically by 5 -10 %) with the above choices of w ref suggesting that our main choice is reasonably good.We confirm the unitarity bound for the axial vector channel as  (10), suggest that the unitarity bound is rather weak as mentioned above.
The BGL parametrization and synthetic data of the form factors are plotted in Fig. 17.
We observe reasonable consistency with the Fermilab/MILC results [18] except the vector form factor g(w) near zero-recoil and the axial form factor F 2 extrapolated to large recoils.This can be seen in a more detailed comparison of the expansion coefficients in Fig. 18, where Fermilab/MILC (triangles down) [18].We also plot results from a recent phenomenological analysis of the Belle data [16] as open black squares [14].Since all these studies impose the kinematical constraint (33) leading to a F 1 ,0 ∝ a f,0 , we omit a panel for a F 1 ,0 .
there is slight tension in a g,0 , a g,1 and a F 1 ,1 .Recent HPQCD results show better consistency with ours, except for a F 1 ,1 and a F 2 ,0 .
In the same figure, we also plot a phenomenological estimate [14] obtained from the Belle data [16] combined with a lattice input of the FLAG average of h A 1 (1) ∝ f (1) [70].Note that this analysis employs the same BGL parametrization with the input (M pole,n and χ {T,L} J P ) and orders (N f , N g , N F 1 ), which are the same as our main analysis.The phenomenological result of a f,0 is well determined and shows good consistency with lattice results, simply because it is mainly fixed from the lattice input.Other poorly determined coefficients suggest the importance of lattice calculations to obtain controlled BGL parametrization.
On the other hand, the slope a F 1 ,1 is well determined from the phenomenological analysis of the experimental data, which are precise at large recoils.This study observes a reasonable agreement in a F 1 ,1 and, hence, in the differential decay rate as discussed in the next subsection (see Fig. 19).We note, however, that the Fermilab/MILC and HPQCD studies reported tension in the w dependence of the differential decay rate with experiment.
Since F 2 describes the contribution suppressed by m 2 ℓ to the decay rate (38), its expansion coefficients are poorly constrained by the experimental data for the light lepton channels B → D * ℓν ℓ (ℓ = e, µ).The lattice determination of F 2 is, therefore, helpful towards a precision new physics search using the τ channel.Through numerical integration of Eq. (38) for ℓ = e, µ and τ with fit results in Table VII (13), respectively, and to be compared with the experimental average 0.295( 14) [1].For a more precise and reliable estimate of R(D * ), it is helpful to simulate larger recoils as well as resolving the tensions in a F 2 ,{0,1} shown in Fig. 18.

B. Fit including Belle data
In this subsection, we carry out a simultaneous fit to our lattice and Belle data to estimate The differential decay rate with respect to w and three decay angles θ l , θ v and χ is given as with the helicity amplitudes given in (35) and (36).For a simultaneous fit with the Belle data for B 0 → D * − ℓ + ν ℓ (ℓ = e, µ) in Ref. [16], we set m ℓ = 0, and include the Coulomb factor  and dΓ/dχ by analytically integrating over two other decay angles and by numerically integrating over w via Simpson's rule with 800 steps within each bin as in the Belle paper [16].
Our results for the expansion coefficients are summarized in Table IX where the first error comes from our main fit with input-A in Table V [14,15] due to D'Agostini effects [71] and the strong systematic correlation of experimental data.That analysis is left for future work.with previous lattice studies.However, the expansion coefficients for g, F 1 and F 2 show significant tension among lattice studies.In particular, our slope for F 1 is slightly larger and consistent with a phenomenological analysis of the Belle data.As a result, our data show better consistency with the Belle data in the differential decay rates.
High statistics simulations on finer lattices would be helpful in resolving this tension.
Simulating larger recoils can lead to a better estimate of R(D * ) as well as more detailed comparison with experiment in a wide range of the recoil parameter to search for new physics.In order to resolve the |V cb | tension, it would also be very helpful to study the inclusive decay on the lattice for more direct comparison of the exclusive and inclusive analyses in the same simulations [72].While it is an ill-posed problem to calculate the relevant hadronic tensor from discrete lattice data on a finite volume, a joint project of the JLQCD and RBC/UKQCD Collaborations along a strategy proposed in Ref. [73,74] is in progress [75][76][77] towards a realistic study of B → X c ℓν ℓ .
momentum variable and the argument ϵ ′ for the D * two-point function C D * (∆t; p ′ , ϵ ′ ) for simplicity.The meson interpolating field is denoted by O P (P = B, D * ), where Gaussian smearing is applied to enhance their overlap with the ground state Z P (p) = ⟨P |O † P ⟩.The B meson is at rest (p = 0) throughout this paper.The w dependence of the form factors is studied by varying the three-momentum of D * as |p ′ | 2 = 0, 1, 2, 3, 4 (in this paper, we denote the momentum on the lattice in units of 2π/L).To this end, we also calculate D * two-point functions with the local sink operator O † D * ,lcl

FIG. 3 :
FIG. 3: Double ratio (11) to estimate h A 1 (1) as a function of temporal location of weak currents ∆t.The open symbols show all data at different choices of the source-sink separations ∆t + ∆t ′ , whereas the filled symbols are data used in the fit (12) to extract the ground state matrix element shown by the horizontal black band.The thick (thin) curves are fit curves inside (outside) the fit range.The left panel shows the data at our smallest pion mass M π ≃ 230 MeV (a −1 ≃ 2.5 GeV, m Q = 1.25m c ), whereas the right panel is for our largest cutoff a −1 ≃ 4.5 GeV (M π ≃ 300 GeV, m Q = 1.25 4 m c ).All data are shifted along the horizontal axis so that the mid-point ∆t = ∆t ′ is located at the center of the panel (∆t = 10 and 20 in the left and right panels, respectively).The symbols are further but slightly shifted along the horizontal axis for clarity.

Figures 5 2 | 10 β 2 FIG. 6 :
Figures 5 and 6 show an example of our results for R ′ 1 and R V at our largest cutoff a −1 ≃ 4.5 GeV and m Q = 1.25 4 m c .Around the mid-point ∆t ∼ ∆t ′ , we observe good consistency 5 -6, 8 -9  to large statistical fluctuations of correlators at the largest temporal separation ∆t+∆t ′ = 40, and accidental anti correlation of the statistical fluctuations.Since these data have larger statistical error than those at smaller ∆t + ∆t ′ , these bumps do not change results of the fits (12) and (15) to extract the ground state contribution.

FIG. 10 :
FIG. 10: Form factors as a function of recoil parameter w.The top-left, top-right, bottom-left and bottom-right panels show results for h A 1 , h A 2 , h A 3 and h V , respectively.The black, blue and red symbols show data at a −1 ∼ 2.5, 3.6 and 4.5 GeV, whereas the open, pale-shaded, filled and dark-shaded ones are at M π ∼ 500, 400, 300 and 230 MeV, respectively.Symbols with different shapes are obtained with different bottom quark masses.(See the legend in the plots.)The green bands show the form factors extrapolated to the continuum limit and physical quark masses.The dark and pale shaded bands represent their statistical and total uncertainties, respectively.

FIG. 11 :
FIG. 11: Axial vector form factor h A 1 (left panel) and vector form factor h v (right panel) as a function of M 2 π .Both panels show data at m Q = 1.25 2 m c and the smallest value of w: namely, w = 1.0 for h A 1 and ∼ 1.205 for h V , respectively.For the latter, w is not exactly equal to 1.025 due to the slight difference in L among simulated lattice spacings as well as discretization effects to w = E * D /M * D itself.The black, blue and red symbols and lines show data and fit curve at a −1 ∼ 2.5, 3.6 and 4.5 GeV, respectively.The green dashed line shows the form factors extrapolated to the continuum limit and physical strange and bottom quark masses.The dark and pale shaded green bands represent the statistical and total error, respectively.

FIG. 12 :
FIG.12: Same as Fig.11but as a function of HQET expansion parameter ϵ Q .To cancel nontrivial m Q dependence due to the matching factor η {A 1 ,V } , we plot a ratio to this factor, which is fitted to the HMChPT-based polynomial form(19) in ϵ Q and ξ am Q .

FIG. 13 :
FIG. 13: Values of χ 2 /d.o.f. for our continuum and chiral extrapolation of form factors.The left panel shows χ 2 /d.o.f. with the singular value cut for the covariance matrix, whereas the right panel

2 FIG. 14 :
FIG.14: Uncertainties of form factors as a function of w.The top-left, top-right, bottom-left and bottom-right panels show uncertainties of h A 1 , h A 2 , h A 3 and h V , respectively.The thick solid and dashed lines show the statistical and total systematic uncertainties, while the thin lines show a breakdown of the latter.Just for clarity, we omit small systematic errors below 0.2 % (h A 1 ), 4 % (h A 2 ), 1 % (h A 3 ) and 0.5 % (h V ) in the whole simulated region of w.

FIG. 16 :
FIG.16: Regular part of axial form factor P f (z)ϕ f (z)f (z) as a function of z.The two error bars in this figure and Fig.17show the statistical and total errors, respectively.

FIG. 17 :
FIG. 17: BGL parametrization (green band) of the synthetic data (green circles) in the whole semileptonic region 1 ≤ w ≲ 1.5.The top-left, top-right, bottom-left and bottom-right panels show f (w), F 1 (w), F 2 (w) and g(w), respectively.We also plot the parametrization by Fermilab/MILC[18] by the open black bands.The stars in the panels for F 1 and F 2 represent the values fixed from the kinematical constraints(33) and(34), respectively.

FIG. 19 :
FIG. 19: Differential decay rate with respect to recoil parameter, dΓ/dw (top-left panel), and three decay angles, dΓ/d cos[θ ℓ ] (top-right panel), dΓ/d cos[θ v ] (bottom-left panel) and dΓ/dχ (bottomright panel).We divide dΓ/dw by a phase space factor √ w 2 − 1 to make its w dependence milder and monotonous.The symbols show Belle results, whereas the green band is reproduced from our simultaneous fit to our lattice and Belle data.

TABLE II :
Bare heavy quark masses in lattice units used to calculate relevant two-and three-point functions.The smallest value corresponds to the physical charm mass fixed from the spin averaged which are stable against the choice of the fit range: we do not observe ∆t , because less data are available for the fit.From R 2 , R 3 and h A 1 extracted above, we can calculate h A 2 and h A 3 . min

TABLE V :
Two sets of input resonance masses with quantum number J P .All values are in

TABLE VI :
Derivatives of vacuum polarization functions entering other functions (

TABLE VII :
Expansion coefficients and their correlation for BGL parametrization (39) of form factors g, f , F

TABLE VIII :
Stability of BGL parametrization coefficients a {g,f,F 1 ,F 2 },{0,1} given in Table This also means that our results for a {g,f,F 1 ,F 2 },{0,1} in Table VII should be used with our choice of the input, namely input-A for resonance masses in Table V and perturbative estimate of χ T,L J P in Table VI, otherwise the uncertainty suggested in Table VIII should be taken into account.In Table VIII, we also test three different choices of the reference values w ref : namely, of all a {f,F 1 },{0,1,2} for the corresponding channel.Also for other coefficients, the changes in the Blaschke factors and outer functions due to the different inputs are largely absorbed into the shift in a {g,f,F 1 ,F 2 },{0,1,2} to reproduce the synthetic data of the form factors.As a result, the shift in physical observables, R(D * ) and |V cb | in this paper, is well below 1 σ as we will see later.
resort to phenomenological models nor experimental data.The second error comes from the test of the stability of the BGL parametrization summarized in Table VIII.It is mainly comes from the shift in w ref min and w ref max , but smaller than statistical and other systematic errors of the BGL parametrization.This estimate can be improved as

Table VI ,
(3)erturbative χ and w ref = 1.025, 1.060, 1.100.The second error is obtained by testing the stability of the fit against the different choices of the input and w ref as in TableVIII.This is in good agreement with the previous determination from the exclusive decay(3).In contrast to previous lattice studies, as shown in Fig.19, we observe consistency between our lattice and Belle data leading to an acceptable value of χ 2 /d.o.f.∼ 0.90.Towards the resolution of the |V cb | tension, we need a more detailed analysis taking care of, for instance, bias in |V cb | discussed in Refs.