Theory determination of B̄ →

We carry out an analysis of the full set of ten B̄ → D(∗) form factors within the framework of the Heavy-Quark Expansion (HQE) to order O ( αs, 1/mb, 1/m 2 c ) , both with and without the use of experimental data. This becomes possible due to a recent calculation of these form factors at and beyond the maximal physical recoil using QCD light-cone sum rules, in combination with constraints from lattice QCD, QCD three-point sum rules and unitarity. We find good agreement amongst the various theoretical results, as well as between the theoretical results and the kinematical distributions in B̄ → D(∗){e−, μ−}ν̄ measurements. The coefficients entering at the 1/mc level are found to be of O(1), indicating convergence of the HQE. The phenomenological implications of our study include an updated exclusive determination of |Vcb| in the HQE, which is compatible with both the exclusive determination using the BGL parametrization and with the inclusive determination. We also revisit predictions for the lepton-flavour universality ratios RD(∗) , the τ polarization observables P D(∗) τ , and the longitudinal polarization fraction FL. Posterior samples for the HQE parameters are provided as ancillary files, allowing for their use in subsequent studies.


I. INTRODUCTION
The decaysB → D −ν andB → D * −ν with = e, µ, τ are of great phenomenological interest for several reasons. First, the decays with light leptons in the final states are used to determine the CKM matrix element |V cb | in the Standard Model (SM). Second, New Physics (NP) scenarios -model-independently defined through the means of an Effective Field Theory (EFT) at low energies -are constrained by both the light-lepton modes and the ones involving a τ lepton. Third, the interplay between two heavy quarks provides a laboratory to study Heavy-Quark Effective Theory (HQET) and the Heavy-Quark Expansion (HQE) of the relevant hadronic matrix elements, to further our understanding of Quantum Chromodynamics (QCD). Interestingly, there are presently two tensions between theory predictions and the corresponding experimental measurements: the so-called V cb puzzle, i.e. the difference between the value for |V cb | as extracted from inclusive vs. exclusive modes, and a significant deviation from lepton-flavour universality in ratios of τ over µ and e modes [1].
The inference of phenomenological parameters such as |V cb | or the EFT Wilson coefficients from experimental measurements of branching ratios and kinematical distributions inB → D ( * ) −ν decays requires knowledge of the relevant hadronic matrix elements. The latter are commonly described by a set of ten independent hadronic form factors, which parametrize the strong-interaction dynamics in these modes as functions of the four-momentum transfer q 2 . The determination of these form factors requires nonperturbative methods, such as lattice QCD or QCD sum rules.
Until recently, the available theoretical calculations were insufficient to fully determine these form factors independently of experimental data; instead, the form factor shapes and |V cb | were fitted together to the light-lepton modes. However, this approach requires the assumption of absence of NP in these modes; this does not seem appropriate, given the anomalies not only in b → cτ ν data, but also in b → sµ + µ − modes, since models accommodating both anomalies commonly also modify the couplings to light leptons in charged-current transitions. Furthermore, these fits were based on a HQE up to O(1/m c,b ). Recently, the determination ofB → D ( * ) form factors has advanced, due to both experimental and theoretical improvements: on the experimental side, the Belle collaboration has released three measurements of the kinematical distributions of the modes in question, including correlations between bins [2][3][4]. BaBar has performed the first analysis of the four-fold differential rate [5], however, these data are not yet available in a form that could be used in our analysis. On the theory side, two lattice determinations of twoB → D form factors at finite recoil became available [6,7]. In addition, a second lattice calculation for oneB → D * form factor at zero recoil was published [8]. Moreover, a recent light-cone sum rule (LCSR) calculation [9] provides for the first time information on all form factors parametrizing matrix elements of the basis of dimension-six operators, including those appearing only in connection with NP effects. This calculation is complementary to the presently available lattice calculations in that it is applicable at q 2 0, while lattice calculations so far have been carried out at q 2 8 GeV 2 , and only for a subset of form factors. For those form factors where lattice data are available, the LCSR calculation therefore acts as an anchor for what would otherwise be an extrapolation of the lattice data based on heavy-quark symmetry relations. For all other form factors this is the only direct calculation available. The release of the unfolded Belle data has made it possible, for the first time, to analyze the spectra ofB → D ( * ) −ν with different approaches for the form factors, while most previous experimental analyses provided their results in terms of parameters of the CLN parametrization, which includes the aforementioned expansion up to O(1/m c ). The BGL parametrization, on the other hand, provides a model-independent parametrization of form factors based on unitarity and analyticity [10], neither expanding in 1/m b,c nor in α s . Assuming the convergence of the latter expansions, clearly the results obtained from either approach should coincide asymptotically. Analyses of the recent Belle measurements using the BGL and CLN approaches yielded the following observations: • BGL fits to the unfoldedB → D data employing also the recent lattice results work very well and yield a value for |V cb | in perfect agreement with the value from inclusive decays [11]. The CLN parametrization, while yielding a similar value for |V cb |, is not sufficiently flexible to accommodate the experimental and lattice data at the same time, indicating the importance of higher-order corrections [11].
• TheB → D experimental and lattice data can be combined in a HQE framework including consistently the correlations due to the parameters in the leading and subleading Isgur-Wise functions, at the cost of introducing partial 1/m 2 c corrections [12,13].
• Comparisons between the BGL and CLN parametrizations using the unfolded Belle 2017 data with hadronic tag [3] show a surprisingly large difference between the values for |V cb |, with the value extracted using the BGL parametrization again compatible with the one from inclusive decays [14][15][16]. The central values of such a fit violate expectations based on heavy-quark symmetry strongly [17], which is not the case, however, once information from the (at the time available) LCSR calculations [18] is included, at the price of a slightly lower increase of |V cb | [14,19].
• No such parametrization dependence is found when employing the recent untagged Belle results [4], but the value of |V cb | extracted from the combined 2017/2018 Belle data remains ∼ 2σ smaller than the one from inclusive decays [20].
Given these results, it is fair to say the V cb puzzle is reduced, but not fully resolved yet. The difficulties in fitting theB → D data and the large differences in the analysis of the tagged Belle data strongly motivate an analysis of higher-order corrections in the HQE framework.
The outline of this article is as follows: in section II we revisit the heavy-quark symmetry relations forB → D ( * ) form factors, with an emphasis on terms that have been generally omitted so far. In section III we combine all available theory information on these form factors, demonstrating the necessity of including additional terms compared to previous treatments. We analyze various scenarios with different classes of inputs in order to probe their mutual compatibility; we provide the fit results for form factors and quantities of interest like R(D ( * ) ) in the viable scenarios. We also apply our extracted form factors in fits to the available experimental data in the context of the SM, and show how the inclusion of the additional terms resolves the deviation in the extracted values of |V cb | in previous fits using the BGL or CLN parametrization. We summarize our results in section IV.

II. FORM FACTORS FORB → D ( * ) TRANSITIONS
The hadronic matrix elements forB(p) → D ( * ) (k) semileptonic transitions can be expressed in terms of ten independent form factors, which are scalar functions of the four-momentum transfer q 2 ≡ (p − k) 2 . A common basis of form factors arises from the following definitions: ForB → D, one commonly defines In the above, f + is the vector form factor, f T is the scale-dependent tensor form factor arising only in NP scenarios (its definition corresponds to the one in Ref. [21]), and f 0 doubles as the scalar form factor: The matrix elements of the remaining axial and pseudoscalar currents are zero by virtue of QCD conserving parity.
ForB → D * , one commonly defines where η denotes the D * polarization vector, V the vector form factor, and A 1,12 are the axial form factors. Note that the relative sign between our eq. (4) and the decomposition in ref. [22] arises from the different definition of the Levi-Civita tensor: we use ε 0123 = +1. Moreover, in the decomposition above A 12 correspond to longitudinal polarizations of the emitted virtual W , which is more convenient (e.g. when inferring form factors from lattice QCD) than parametrizations involving the form factor A 2 , see e.g. [22]. The function A 0 doubles as the pseudo-scalar form factor, whereas the matrix element of the scalar current vanishes by virtue of QCD conserving parity.
Exact relations at q 2 = 0 between some of the form factors ensure the absence of unphysical singularities in eq. (1) and eq. (5). These relations read: A further exact relation arises due to algebraic identities involving the Lorentz structures σ µν and σ µν γ 5 [22]: Further approximate relations arise from the HQE of the hadronic matrix elements. These relations, the parametric models involved, and theoretical inputs needed for the subsequent statistical analyses are the subject of the remainder of this section.

A. Heavy-Quark Expansion and models
The combination of heavy-quark spin symmetry and heavy-quark flavour symmetry permits to relateB ( * ) (v) → D ( * ) (v ) matrix elements with each other in a simultaneous expansion in the strong coupling α s and the inverse pole masses 1/m Q , where Q = b, c is the quark flavour. The coefficients of this HQE -up to kinematical and combinatorial factors -are the Isgur-Wise functions, which depend exclusively on the recoil parameter w ≡ v · v . For convenience, the expansion is commonly expressed in terms of dimensionless quantities ε Q ≡ Λ/2m Q , where Λ arises in the HQE of the heavy meson masses.
We begin by adopting the power counting ε b ∼ ε 2 c ∼ α s /π ∼ ε 2 for the HQE. Consequently, when expanding up to O(ε 2 ), we need to account for all leading-order radiative and subleading-power corrections, as well as partial subsubleading-power corrections. Higher powers in our expansion or mixed terms are assumed to be negligible. The HQE is well known, and we follow ref. [12] closely in our analysis. By virtue of our power counting, any form factor F discussed in section II can be expressed in terms of ten independent functions: ξ, the leading Isgur-Wise (IW) function; χ 2,3 and η, the subleading IW functions; and 1,...,6 , the subsubleading IW functions at order ε 2 c as introduced in ref. [23]; see also appendix B 3 for more details. Each of these functions depends on the recoil parameter w. In the complex half plane Re w ≥ 1, the form factors, the HQET Wilson coefficients, and the IW functions are free of singularities due to QCD dynamics. Singularities of kinematical origin can always be removed by redefining the form factors. Consequently, for f being any of the ten IW functions considered here, we expand it around w = 1: 1 Following ref. [10,24], we can further trade the variable w for z(w), which correctly captures the analytic properties of the matrix elements, i.e. it develops a branch cut corresponding to the B ( * ) D ( * ) pair production at w ≤ −1. While z(w) is a small expansion parameter in the semileptonic phase space, absent further modifications to the form factors as discussed in ref. [10] we cannot generally expect small coefficients in an expansion in z. We proceed to expand each monomial (w − 1) k in eq. (10) around z (1), where the maximum order in z − z(1) depends on our concrete parameter models discussed later. In this way, we keep the benefits inherent to parametrizing in z, while conserving at the same time the physical meaning of the fit parameters f (n) (1) as derivatives of the IW functions at the zero-recoil point. In this setup, we follow ref. [13] closely.
Both HQET and the HQE of the heavy meson masses provide us with some information on the parameters arising in the HQE of the hadronic matrix elements at hand. The remaining ones need inferring from theoretical or experimental inputs. For our statistical analyses we define fit models, which vary only in our choice of the order to which the different Isgur-Wise functions are expanded in z. All our models include all ten Isgur-Wise functions above; the expansion up to 1/m 2 c is not only preferable from the point of view of precision, but, as mentioned in the introduction, necessary, given the available lattice data at w = 1. Employing the recent LCSR results [9] allows to include all subsubleading IW functions, which is an improvement compared to ref. [13], where only the functions 1,2 could be included. The models used in this work are denoted as k/l/m, where the numbers k, l and m have the following meaning: k : the order to which the leading IW function is expanded in z around z = 0; l : the order to which the subleading IW functions are expanded; and m : the order to which the subsubleading IW functions are expanded.
We keep all purely kinematical powers of w, i.e. terms that arise when relating the form factors to the IW functions in the HQE. Within the scope of this work we will discuss the models 2/1/0 and 3/2/1.
We emphasize that increasing the maximum order in the z expansion from one fit model to another can always be expressed in terms of a non-zero shift to the new parameter appearing in the higher-order term. As an example, consider increasing the order of the z expansion from a 2/l/m to a 3/l/m model. The 2/l/m models can be recovered in the 3/l/m model by assigning ξ (1) = −ξ (1)/2 − 3ξ (1)/64.

B. Theory constraints
With the definition of the models at hand we proceed to the available theoretical calculations of the hadronic matrix elements as well as theoretical bounds on the parameter space derived from dispersion relations [10,24]. The individual pieces of theory information entering the likelihood are: Lattice: ForB → D the HPQCD and FNAL/MILC collaborations have, independently from each other, determined the vector form factor f + and the scalar form factor f 0 at several values of the recoil parameter w ≥ 1. We use correlated pseudo data points from both studies: seven from [6] and five from [7]. Note that at w = w max,D the form factors fulfill an equation of motion that reduces the number of observations per study. ForB → D * the HPQCD and FNAL/MILC collaborations have independently determined the form factor h A1 at w = 1 [8,25], averaged by FLAG to h A1 (w = 1) = 0.904 ± 0.012 [26].
QCDSR: The subleading IW functions χ 2,3 and η have been studied within three-point QCD Sum Rules [27][28][29]. These sum rules have been used to infer the normalization and slope of the subleading IW functions at w = 1, yielding five observations in total.
LCSR: At w 1.5 theB → D ( * ) form factor can be accessed using LCSRs with B-meson Light-Cone Distribution Amplitudes (LCDAs) [18]. These results have been superseded by an updated analysis [9], which includes for the first time all two-particle and three-particle LCDAs in a consistent twist-expansion up to twist 4 [30]. Moreover, the recent analysis provides for the first time information about the shape of the complete set ofB → D * form factors at four phase space points, albeit with two caveats: the form factors are available only at w ≥ w max 1.5 and theB → D form factor f T could not be extracted as part of the same approach as the other form factors.
The first point requires attention in the context of the z expansion, since it increases the maximal value of |z| to values larger than encountered in other studies. The second point dissuades us from including f T in the analysis; instead, we choose to predict f T within our approach and compare it with the prediction of ref. [9]. Following the introductory discussion concerning exact relations between some of the form factors, we arrive at a total number of 33 observations.
Beyond the likelihood, we include further information on the hadronic matrix elements. This additional information is expressed as so-called unitarity bounds [10,24]. In the context of the HQE, it is convenient to adopt the approach of ref. [24]. We consider the bounds for the currents For all currents J L P we derive the bounds in terms of the full set of ten independent IW functions present in our models. The results of the perturbative OPE calculations for the bounds are denoted as χ L P andχ L P , where the tilde indicates subtraction of known one-body contributions [24]. Updated values for these four quantities have been recently provided in ref. [11,14,19] based on ref. [31]. The general problem of how to include positivity bounds in a Bayesian fit and our approach to solve it is discussed in appendix A.
In the construction of the unitarity bounds, a choice must be made to which order n in z the bound is formulated. Using BGL-like form factors, this coincides with the order to which the form factors are expanded. However, the treatment of the unitarity bounds in the context of the HQE is non-trivial. The reason for their complexity arises from the simultaneous expansion in 1/m Q , α s , and z. As indicated above, we expand the IW functions on different levels in the 1/m Q expansion to different orders in z, according to their combined power-counting, i.e., a z 3 contribution might be relevant for the form factors when entering the leading IW function, while such a term in a subsubleading IW function is expected to be negligible. Hence we choose generally k ≥ l ≥ m. However, in a BGL setup the unitarity bounds are written as quadratic forms of the BGL coefficients without explicit factors of z. Therefore the relative importance of higher-order contributions in z is larger in these bounds than in the form factors themselves. Consequently, the treatment of these higher-order contributions is important. Specifically, 1/m 2 Q contributions are only fully included for n ≤ m, 1/m Q contributions for n ≤ l and leading-order contributions for n ≤ k. Since particularly 1/m c contributions can be large, and the terms in the 1/m Q expansion are not necessarily positive, the order n for the unitarity bounds should be chosen to be at most n = l, with l ≥ 1.
The combination of the described constraints allows to include higher-order contributions in the HQE for the full set of form factors. Within our determination of these contributions we pose the following questions: • Are the various theoretical constraints mutually compatible in the context of the HQE? If yes, what is the minimal k/l/m model that achieves a good description?
• In case of a successful combined fit, what are the phenomenological consequences with respect to |V cb | and predictions for B → D ( * ) τ ν observables?

III. STATISTICAL ANALYSES
The numerical and statistical results presented in the following have been obtained by means of two completely independent implementations. One of these is publicly available as part of the EOS software [32], which has also been 2/1/0 3/2/1 3/2/1 3/2/1 3  used to prepare all of the following numerical results and plots. The posterior samples used to produce these results and plots are available as ancillary files [33].

A. Fits to only theory constraints and SM predictions
The minimal model fulfilling all criteria laid out in the previous section is the 1/1/0 model. We find this model to provide a bad fit to the available theory constraints, with χ 2 ∼ 560 in the best-fit point for 39 degrees of freedom (dof). We therefore proceed to fit the theory constraints with the 2/1/0 model, which yields an excellent fit with χ 2 = 22.87 for 38 dof. This model therefore represents the minimal viable fit model. 2 Following the discussion in the previous section, it is important to account for systematic uncertainties inherent to the HQE by increasing the order of the z expansions by one. The corresponding model, 3/2/1, reduces the χ 2 in the best-fit point by ∼ 13, at the expense of 10 additional parameters. Details for these fits are given in table I.
Using samples of the posteriors of the fits to both models, we produce posterior predictive distributions for all B → D ( * ) form factors, including f T . The median values and 68% probability envelopes for each form factor are shown in figure 1, together with data points illustrating the theory constraints where applicable. We make the following observations: • As expected, the uncertainty bands are systematically broader in the 3/2/1 model.
• For the form factors f 0 , A 1 and T 2 we observe that model 2/1/0 produces a local minimum for −15 GeV 2 ≤ q 2 ≤ 0, where the LCSR constraints are available. This does not conform to the usual expectation in a dispersive picture: Far below the production threshold and sub-threshold poles it should be possible to approximate the dispersive integrals for the form factors with a single effective pole, leading to a monotonically falling form factor with decreasing q 2 .
• Neither of the two models is able to simultaneously fit all the nominal theory constraints plus the LCSR constraints on theB → D form factor f T . This is not surprising, given the different framework used for its prediction compared to the other form factors, as discussed in ref. [9].
In addition, we use the posterior samples for the 3/2/1 fit model to produce posterior-predictive distributions in the SM for the LFU ratios R D and R D * , the τ polarizations P D ( * ) τ inB → D ( * ) τ −ν decays, and the longitudinal polarisation fraction F L inB → D * τ −ν decays. We obtain We also produce posterior-predictive distributions for theB → D ( * ) {e − , µ − }ν branching ratios for both of our fit models. Their summaries in form of mean value, standard deviations and correlations are collected in table III.

B. Challenging measurements and extraction of |V cb |
We apply the form factors obtained in the previous subsection to the available experimental information to perform phenomenological studies with high accuracy. Specifically, we confront our predictions with the measured spectral information and extract |V cb |, assuming the SM. Our extraction of |V cb | to subsubleading power in the HQE is the first of its kind.
The publicly available experimental results areB → D ( * ) −ν kinematical distributions published by the Belle collaboration [2][3][4] and the world averages for the branching fractions [1]. TheB → D −ν distribution P D (w) from ref. [2], and the fourB → D * −ν distributions P D * (w), P D * (χ), P D * (cos θ D * ), and P (cos θ ) from ref. [3] 3 are unfolded of detector effects by the Belle collaboration. The data presented in ref. [4] are still folded, and the necessary information for the unfolding process is provided in the publication.
In a first step, we compare in figure 2 our posterior predictions for the kinematical PDFs with the experimental results. Both of our fit models yield visually indistinguishable posterior predictions for the three angular distributions P (χ), P (cos θ ) and P (cos θ D * ) inB → D * {e − , µ − }ν. The agreement between our predictions and the experimental measurements for P (cos θ ) is visibly worse than the excellent agreement for the remaining two angular distributions. However, we find that our predictions for these three distributions are considerably more precise than the experimental results. We therefore conclude that the latter do not further constrain the form factor parameters within our two models; we hence abstain from using them in the following. However, we find that the results for the distributions P D (w) and P D * (w) do have the potential to further constrain the form factor parameters.
In a second step, we fit the HQE expressions for the form factors simultaneously to the previously discussed theory constraints and different sets of publicly available experimental results for P D (w) and P D * (w). These sets are: only P D * (w) from the 2017 data, only P D * (w) from the 2018 data, and the combination of all experimental results for P D ( * ) (w). For all these sets we find that the simultaneous fits show excellent agreement between the theoretical . .   [3,4]. The row for the combined |V cb | takes into account correlations in the posterior predictions, but not the small and unpublished correlations in the HFLAV world averages.
constraints and the experimental PDFs. A summary of the goodness of fit in the best-fit points of all considered sets is presented in table I. Our nominal best-fit point, obtained in the 3/2/1 model, is presented in table II. Due to the non-Gaussianity of the posterior we refrain from providing linear correlations. We note that the slopes of the subsubleading IW functions i (1) are all compatible with zero at 68% probability.
Our predictions for theB → D ( * ) τ −ν observables including the experimental information read: While the predictions forB → Dτ −ν remain unchanged, we find a shift of ∼ 0.5σ for the threeB → D * τ −ν observables. This is not surprising, given the high precision of the availableB → D form factor constraints.
Finally, we produce posterior predictions for the integrated branching ratios ofB → D ( * ) −ν decays in units of |V cb | 2 . We choose to present our results for theB 0 mode only. Our results are listed in the top half of table III. We then proceed to extract the value of |V cb | from the isospin averages of the respective branching ratios. Our results for |V cb | are listed in the bottom half of table III. The isospin average of the necessary branching ratios, expressed as branching ratios of theB 0 mode, are: Our nominal result for the exclusive determination of |V cb |, obtained by combining all available theoretical and experimental information, is: Its agreement with the individual values fromB → D ( * ) −ν is excellent. Averaging the two exclusive determinations with the inclusive one [34], we find |V cb | = (41.3 ± 0.5) × 10 −3 , where the three values are compatible at the 1.2σ level.

IV. SUMMARY AND OUTLOOK
In this work we carry out a comprehensive analysis of the full set ofB → D ( * ) −ν form factors. The basis of our analysis is the Heavy-Quark Expansion (HQE) up to O(ε 2 ), where our power-counting is defined as ε b ∼ ε 2 c ∼ α s /π ∼ ε 2 . By determining the coefficients of this expansion from all available theoretical constraints we are able to predict the full set of form factors with high precision. This allows for their consistent and accurate use in a variety of phenomenological applications without the assumption of absence of NP effects in b → c transitions with light leptons. Our work focuses on two applications: precision predictions forB → D ( * ) τ −ν observables and accurate determinations of |V cb | . .    The four 1D kinematical probability distibutions arising from the full 4D differential decay rateB → D * {e − , µ − }ν next to one of the 1D kinematical distribution inB → D{e − , µ − }ν decays. We show all available constraints published by the Belle collaboration along side the 68% probability envelopes based on the two fit models fitted to only theory inputs. Note that for the 2018 Belle results we have unfolded the results ourselves; see the text for details. fromB → D ( * ) {e − , µ − }ν decays.
We find excellent agreement between the various theoretical constraints on the relevant hadronic matrix elements. The minimal viable fit model is found to be the 2/1/0 model, where the numbers refer to the order in the z expansion of the leading, subleading, and subsubleading Isgur-Wise (IW) functions, respectively. To account for systematic uncertainties inherent to the HQE, we increase the order of the z expansion for all three sets of IW functions, which defines our nominal fit model. Within our analysis we pay particular attention to the subsubleading terms in the HQE. In previous analyses it turned out to be necessary to include at least two such terms. Our analysis is the first to include the full set of IW functions at the order O(1/m 2 c ). We find that the expansion in 1/m Q is well-behaved, similar to what has been found in a recent analysis of Λ b → Λ c form factors [35]. Based on these findings we expect the terms at O(ε 3 ) to be negligible at the present level of precision. This assumption should be revisited once more precise theoretical and experimental information becomes available.
Our predictions forB → D ( * ) τ −ν observables benefit from the improved treatment of the HQE. This is reflected by significantly smaller uncertainties compared to previous analyses, while staying compatible at the 1σ level. Predictions with and without the use of experimental inputs are given in eq. (12) and eq. (11), respectively.
Our determinations of |V cb | from D and D * final states are mutually compatible and also compatible with the inclusive determination at the 1.2σ level. Unlike for CLN analyses, we find no tension with the BGL determinations. Our nominal result for |V cb | using all exclusive experimental inputs reads The upcoming lattice analyses of four of theB → D * form factors at nonzero recoil [36][37][38] will benefit our approach and help to determine the HQE parameters to even higher precision.
as independent Gaussian distribution for each bound, centered around a value µ B and with a standard deviation of σ B . In that case we obtain: