Theory determination of $\varvec{\bar{B}\rightarrow D^{(*)}\ell ^-\bar{\nu }}$ form factors at $\varvec{\mathcal {O}(1/m_c^2)}$

Abstract

We carry out an analysis of the full set of ten $\bar{B}\rightarrow D^{(*)}$ form factors within the framework of the Heavy-Quark Expansion (HQE) to order $\mathcal {O}\left( \alpha _s,\,1/m_b,\,1/m_c^2\right)$, 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 $\bar{B}\rightarrow D^{(*)}\lbrace e^-,\mu ^-\rbrace \bar{\nu }$ measurements. The coefficients entering at the $1/m_c^2$ level are found to be of $\mathcal {O}(1)$, indicating convergence of the HQE. The phenomenological implications of our study include an updated exclusive determination of $|V_{cb}|$ 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 $R_{D^{(*)}}$, the $\tau$ polarization observables $P_\tau ^{D^{(*)}}$, and the longitudinal polarization fraction $F_L$. Posterior samples for the HQE parameters are provided as ancillary files, allowing for their use in subsequent studies.

1 Introduction

The decays ${\bar{B}}\rightarrow D\ell ^-{\bar{\nu }}$ and ${\bar{B}}\rightarrow D^{*}\ell ^-{\bar{\nu }}$ with $\ell =e,\mu ,\tau$ 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 $\tau$ 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 $\tau$ over $\mu$ 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 in ${\bar{B}} \rightarrow D^{(*)}\ell ^-{\bar{\nu }}$ 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.

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\rightarrow c\tau \nu$ data, but also in $b\rightarrow s\mu ^+\mu ^-$ 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 ${\mathcal {O}}(1/m_{c,b})$. Recently, the determination of ${\bar{B}}\rightarrow 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 two $\bar{B}\rightarrow D$ form factors at finite recoil became available [6, 7]. In addition, a second lattice calculation for one $\bar{B}\rightarrow 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\lesssim 0$, while lattice calculations so far have been carried out at $q^2\gtrsim 8~\mathrm{GeV}^2$, and only for a subset of form factors. For the $B\rightarrow D^*$ form factor $A_1$, 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 $B\rightarrow D^*$ form factors this is the only direct calculation available. In view of the irreducible systematic uncertainties of this method our analysis should be revisited once precise lattice QCD results become available for $B\rightarrow D^*$ from factors.

The release of the unfolded Belle data has made it possible, for the first time, to analyze the spectra of ${\bar{B}}\rightarrow D^{(*)}\ell ^-{\bar{\nu }}$ 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 ${\mathcal {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 $\alpha _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 unfolded $\bar{B}\rightarrow 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].

• The $\bar{B}\rightarrow 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_c^2$ 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 $\sim 2\sigma$ 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 the $\bar{B}\rightarrow 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 Sect. 2 we revisit the heavy-quark symmetry relations for ${\bar{B}}\rightarrow D^{(*)}$ form factors, with an emphasis on terms that have been generally omitted so far. In Sect. 3 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 Sect. 4.

2 Form factors for $\varvec{\bar{B}\rightarrow D^{(*)}}$ transitions

The hadronic matrix elements for $\bar{B}(p)\rightarrow 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 \equiv (p - k)^2$. A common basis of form factors arises from the following definitions: For $\bar{B}\rightarrow D$, one commonly defines

$$\begin{aligned}&\langle D(k)| \bar{c} \gamma ^\mu b |\bar{B}(p)\rangle \nonumber \\&\quad = \left[ (p + k)^\mu - \frac{M_B^2 - M_D^2}{q^2} q^\mu \right] f^{B\rightarrow D}_+(q^2)\nonumber \\&\qquad + \frac{M_B^2 - M_D^2}{q^2} q^\mu f^{B\rightarrow D}_0(q^2), \end{aligned}$$

(1)

$$\begin{aligned}&\langle D(k)| \bar{c} \sigma ^{\mu \nu } b |\bar{B}(p)\rangle \nonumber \\&\quad = \frac{2i}{M_B+M_D}(k^\mu p^\nu -p^\mu k^\nu )f_T(q^2,\mu ), \end{aligned}$$

(2)

with $\sigma ^{\mu \nu }=\frac{i}{2}[\gamma ^\mu ,\gamma ^\nu ]$. 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:

$$\begin{aligned} \langle D(k)| \bar{c} b |\bar{B}(p)\rangle&= \frac{M_B^2 - M_D^2}{m_b - m_c}\, f^{B\rightarrow D}_0(q^2). \end{aligned}$$

(3)

The matrix elements of the remaining axial and pseudoscalar currents are zero by virtue of QCD conserving parity.

For $\bar{B}\rightarrow D^*$, one commonly defines

$$\begin{aligned}&\langle D^*(k, \eta )| \bar{c} \gamma ^\mu b |\bar{B}(p)\rangle \nonumber \\&\quad = -\epsilon ^{\mu \nu \rho \sigma } \eta ^*_{\nu }(k)\, p_\rho \, k_\sigma \frac{2\, V(q^2)}{M_B + M_{D^*}}, \end{aligned}$$

(4)

$$\begin{aligned}&\langle D^*(k, \eta )| \bar{c} \gamma ^\mu \gamma _5 b |\bar{B}(p)\rangle \nonumber \\&\quad = i \eta ^*_{\nu } \left\{ 2 M_{D^*} A_0(q^2) \frac{q^\mu q^\nu }{q^2}+ 16\frac{M_B M_{D^*}^2}{\lambda } A_{12}(q^2) \right. \nonumber \\&\qquad \left[ 2 p^\mu q^\nu - \frac{M_B^2 - M_{D^*}^2 + q^2}{q^2} q^\mu q^\nu \right] \nonumber \\&\qquad + (M_B + M_{D^*})\, A_1(q^2) \left[ g^{\mu \nu } + \frac{2(M_B^2 + M_{D^*}^2 - q^2)}{\lambda } q^\mu q^\nu \right. \nonumber \\&\qquad \left. \left. - \frac{2(M_B^2 - M_{D^*}^2 - q^2)}{\lambda } p^\mu q^\nu \right] \right\} ,\end{aligned}$$

(5)

$$\begin{aligned}&\langle D^*(k, \eta )| \bar{c} \sigma ^{\mu \nu } b |\bar{B}(p)\rangle \nonumber \\&\quad = i \eta ^*_{\alpha }\epsilon ^{\mu \nu }\!{}_{\rho \sigma }\left\{ -\left[ \left( (p+k)^\rho -\frac{M_B^2-M_{D^*}^2}{q^2}q^\rho \right) g^{\alpha \sigma }\right. \right. \nonumber \\&\qquad \left. + \frac{2}{q^2} p^\alpha p^\rho k^\sigma \right] T_1(q^2) \nonumber \\&\qquad -\left. \left( \frac{2}{q^2}p^\alpha p^\rho k^\sigma -\frac{M_B^2-M_{D^*}^2}{q^2}q^\rho g^{\alpha \sigma }\right) T_2(q^2)\right. \nonumber \\&\qquad \left. +\frac{2}{M_B^2-M_{D^*}^2}p^\alpha p^\rho k^\sigma T_3(q^2)\right\} . \end{aligned}$$

(6)

where $\eta$ 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 $\varepsilon ^{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,

$$\begin{aligned} \langle D^*(k, \eta )| \bar{c} \gamma _5 b |{\bar{B}}(p)\rangle&= -2i M_{D^*}\,\frac{\eta ^* \cdot q}{m_b + m_c} A_0, \end{aligned}$$

(7)

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 Eqs. (1) and (5). These relations read:

$$\begin{aligned} \begin{aligned} f_+(q^2 = 0)&= f_0(q^2 = 0), \\ A_0(q^2 = 0)&= \frac{M_B + M_{D^*}}{2 M_{D^*}}\,A_1(q^2 = 0) \\&\quad - \frac{M_B - M_{D^*}}{2 M_{D^*}}\,A_2(q^2 = 0). \end{aligned} \end{aligned}$$

(8)

A further exact relation arises due to algebraic identities involving the Lorentz structures $\sigma ^{\mu \nu }$ and $\sigma ^{\mu \nu }\gamma _5$ [22]:

$$\begin{aligned} T_1(0) = T_2(0). \end{aligned}$$

(9)

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.

2.1 Heavy-quark expansion and models

The combination of heavy-quark spin symmetry and heavy-quark flavour symmetry permits to relate $\bar{B}^{(*)}(v) \rightarrow D^{(*)}(v')$ matrix elements with each other in a simultaneous expansion in the strong coupling $\alpha _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 \equiv v\cdot v'$. For convenience, the expansion is commonly expressed in terms of dimensionless quantities $\varepsilon _Q \equiv \overline{\Lambda }/ 2 m_Q$, where $\overline{\Lambda }$ arises in the HQE of the heavy meson masses.

We begin by adopting the power counting $\varepsilon _b\sim \varepsilon _c^2\sim \alpha _s/\pi \sim \varepsilon ^2$ for the HQE. Consequently, when expanding up to ${\mathcal {O}}(\varepsilon ^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 Sect. 2 can be expressed in terms of ten independent functions: $\xi$, the leading Isgur–Wise (IW) function; $\chi _{2,3}$ and $\eta$, the subleading IW functions; and $\ell _{1,\dots ,6}$, the subsubleading IW functions at order $\varepsilon _c^2$ as introduced in Ref. [23]; see also appendix B3 for more details. Each of these functions depends on the recoil parameter w. In the complex half plane ${\text {Re}}{w} \ge 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$:¹

$$\begin{aligned} \begin{aligned} f(w)&= \sum _{k=0}^K \frac{f^{(k)}}{k!}\, \left( w - 1\right) ^k. \end{aligned} \end{aligned}$$

(10)

Following Ref. [10, 24], we can further trade the variable w for

$$\begin{aligned} z(w)\equiv z(w,1)\quad \mathrm{with}\quad z(w,a) = \frac{\sqrt{1+w}-\sqrt{2}a}{\sqrt{1+w}+\sqrt{2}a}, \end{aligned}$$

(11)

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 \le -1$. While z(w) is a small expansion parameter in the semileptonic phase space, in absence of 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_c^2$ 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 $\ell _{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 $\xi '''(1) = -\xi ''(1) / 2 - 3 \xi '(1) / 64$.

2.2 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::

For $\bar{B}\rightarrow 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 \ge 1$.² We use correlated pseudo data points from both studies: seven from [6] and five from [7]. Note that at $w = w_{\text {max},D}$ the form factors fulfill an equation of motion that reduces the number of data points per study. For $\bar{B}\rightarrow D^*$ the HPQCD and FNAL/MILC collaborations have independently determined the form factor $h_{A_1}$ at $w = 1$ [8, 25], averaged by FLAG to $h_{A_1}(w=1)=0.904\pm 0.012$ [26].

QCDSR::

The subleading IW functions $\chi _{2,3}$ and $\eta$ 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 data points in total.

LCSR::

At $w \gtrsim 1.5$ the $\bar{B}\rightarrow 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 of $\bar{B}\rightarrow D^*$ form factors at four phase space points, albeit with two caveats: the form factors are available only at $w \ge w_\text {max}\simeq 1.5$ and the $\bar{B}\rightarrow 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 data points.

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 $J_{0^+} \equiv \bar{c}b$, $J_{0^-} \equiv \bar{c}\gamma _5 b$, $J_{1^-}^\mu \equiv \bar{c}\gamma ^\mu b$, and $J_{1^+}^\mu \equiv \bar{c}\gamma ^\mu \gamma _5 b$. 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 $\chi _{L^P}$ and ${\tilde{\chi }}_{L^P}$, where the tilde indicates subtraction of known one-body contributions [24]. Updated values for these four quantities have been recently provided in Refs. [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 1.

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$, $\alpha _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\ge l\ge 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_Q^2$ contributions are only fully included for $n\le m$, $1/m_Q$ contributions for $n\le l$ and leading-order contributions for $n\le 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\ge 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\rightarrow D^{(*)}\tau \nu$ observables?

3 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 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].

3.1 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 $\chi ^2\sim 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 $\chi ^2=22.87$ for 38 dof. This model therefore represents the minimal viable fit model.³ 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 $\chi ^2$ in the best-fit point by $\sim 13$, at the expense of 10 additional parameters. Details for these fits are given in Table 1.

Table 1 Summary of the goodness of fit in terms of the $\chi ^2$ values at the best-fit point for all combinations of fit models and datasets. The largest single pull arises from the QCDSR constraint on $\chi _3'(1)$ in the 2/1/0 fit model with $\chi ^2 \simeq 4$

Using samples of the posteriors of the fits to both models, we produce posterior predictive distributions for all $\bar{B}\rightarrow D^{(*)}$ form factors, including $f_T$. The median values and $68\%$ probability envelopes for each form factor are shown in Fig. 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~\mathrm{GeV}^2\le q^2\le 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$. Note, however, that this behaviour is not statistically significant.

• Neither of the two models is able to simultaneously fit all the nominal theory constraints plus the LCSR constraints on the ${\bar{B}}\rightarrow D$ form factor $f_T$, which have been excluded from the nominal fit. This is not too surprising, given the inability to determine the effective threshold parameter for $f_T$ and subsequent approximation with the threshold parameter of $f_+$ as discussed in Ref. [9]. This approximation and the additional systematic uncertainty assigned due to this approximation do not yield $f_T$ results that are compatible with the remaining form factors within the HQE; we therefore do not consider the LCSR results for $f_T$ reliable. We emphasize that $f_T$ is the only form factor affected by this issue and that the constraints on this form factor do not enter any of our results.

• The unitarity bounds for the scalar and pseudoscalar channels are essentially saturated for the central values of our fits, implying that they are effective form factor constraints. Due to the heavy-quark relations, they constrain all form factors, not just the (pseudo-)scalar ones. On the other hand, since the posterior distributions for these bounds are broad (the $68\%$ intervals include values down to 0.28 and 0.31 for the scalar and pseudoscalar bound, respectively), this observation does not imply any significant tension with the physical bound, even when allowing for additional intermediate states. The bounds for the vector and axial-vector channels are not effective in constraining the model parameters, since their posterior intervals are safely below unity.

Fig. 1

The full set of $\bar{B}\rightarrow D^{(*)}$ form factors as a function of $q^2$ are used to showcase our results for the nominal fit model 3/2/1 (orange lines and areas), in comparison to the minimal viable fit model 2/1/0 (light blue lines and areas). For both models we show the central values and $68\%$ probability envelopes from posterior-predictive distributions of the respective fits. The lattice constraints used in the fits are shown as green data points. The LCSR constraints used in the fits are shown as purple data points. The superseded LCSR results not used in the fits are shown as red data points for comparison, only

Fig. 2

The four 1D kinematical probability distributions arising from the full 4D differential decay rate $\bar{B}\rightarrow D^*\{e^-,\mu ^-\}\bar{\nu }$ next to one of the 1D kinematical distribution in $\bar{B}\rightarrow D\{e^-,\mu ^-\}{\bar{\nu }}$ 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

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 $\tau$ polarizations $P_\tau ^{D^{(*)}}$ in $\bar{B} \rightarrow D^{(*)}\tau ^-{\bar{\nu }}$ decays, and the longitudinal polarisation fraction $F_L$ in $\bar{B}\rightarrow D^*\tau ^-{\bar{\nu }}$ decays. We obtain

(12)

The longitudinal $D^*$ polarisation for electrons is obtained as $F^e_L = 0.540 \pm 0.013$, which is well compatible with the Belle measurement of $0.56 \pm 0.02$ [34]. We also produce posterior-predictive distributions for the $\bar{B}\rightarrow D^{(*)}\{e^-,\mu ^-\}{\bar{\nu }}$ branching ratios for both of our fit models. Their summaries in form of mean value, standard deviations and correlations are collected in Table 3. We emphasize again that these results are independent of experimental data and improve on previous works using the HQE by including higher orders in $1/m_c$ and z explicitly, thereby providing a reliable error estimate.

Table 2 Best-fit point for the parameters of the 3/2/1 model in a simultaneous fit to theory constraints and all available experimental measurements. Uncertainty ranges presented here are meant for illustrative purpose only, and should not be interpreted a standard deviations due to non-Gaussianity of the joint posterior

Table 3 Posterior predictions for the branching ratios of $\bar{B}^0\rightarrow \lbrace D^+,D^{*+}\rbrace \{e^-,\mu ^-\}{\bar{\nu }}$ decays in units of $|V_{cb}|^{2}$, as well as the values of $|V_{cb}|$ extracted from the isospin-averaged branching ratios. Throughout we use the HFLAV world averages [1] for the determination of $|V_{cb}|$. For the columns marked 2017 and 2018, the values in parentheses are obtained by using only the respective Belle measurements of the branching ratios [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

3.2 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 are $\bar{B}\rightarrow D^{(*)}\ell ^-{\bar{\nu }}$ kinematical distributions published by the Belle collaboration [2,3,4] and the world averages for the branching fractions [1]. The $\bar{B}\rightarrow D\ell ^-{\bar{\nu }}$ distribution $P_D(w)$ from Ref. [2], and the four $\bar{B}\rightarrow D^*\ell ^-{\bar{\nu }}$ distributions $P_{D^*}(w)$, $P_{D^*}(\chi )$, $P_{D^*}(\cos \theta _{D^*})$, and $P(\cos \theta _\ell )$ from Ref. [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 Fig. 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(\chi )$, $P(\cos \theta _\ell )$ and $P(\cos \theta _{D^*})$ in $\bar{B}\rightarrow D^*\{e^-,\mu ^-\}{\bar{\nu }}$. The agreement between our predictions and the experimental measurements for $P(\cos \theta _\ell )$ 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 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 1. Our nominal best-fit point, obtained in the 3/2/1 model, is presented in Table 2. Due to the non-Gaussianity of the posterior we refrain from providing linear correlations. We note that the slopes of the subsubleading IW functions $\ell _i'(1)$ are all compatible with zero at $\simeq 68\%$ probability.

Our predictions for the ${\bar{B}}\rightarrow D^{(*)}\tau ^-{\bar{\nu }}$ observables including the experimental information from the light-lepton modes read:

(13)

The predictions for ${\bar{B}}\rightarrow D\tau ^-{\bar{\nu }}$ remain unchanged, which is not surprising, given the high precision of the available ${\bar{B}}\rightarrow D$ form factor constraints. For the three ${\bar{B}}\rightarrow D^*\tau ^-{\bar{\nu }}$ observables we find a shift of $\sim 0.5\sigma$ and a reduction in the uncertainties by 0.003, which is particularly significant for $R(D^*)$. Again these results improve on previous works using the HQE by including higher orders in $1/m_c$ and z explicitly. Furthermore, the latest experimental measurements are included. As a result, the uncertainties are on the same level as some of the previous estimates using the HQE, while including more conservative error estimates.

Finally, we produce posterior predictions for the integrated branching ratios of $\bar{B}\rightarrow D^{(*)}\ell ^-{\bar{\nu }}$ decays in units of $|V_{cb}|^2$. We choose to present our results for the $\bar{B}^0$ mode only. Our results are listed in the top half of Table 3. 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 3. The isospin average of the necessary branching ratios, expressed as branching ratios of the $\bar{B}^0$ mode, are:

$$\begin{aligned} \begin{aligned}&{\mathcal {B}}(\bar{B}^0\rightarrow D^+\{e^-,\mu ^-\}{\bar{\nu }}) = (2.24 \pm 0.07)\%, \\&{\mathcal {B}}(\bar{B}^0\rightarrow D^{*+}\{e^-,\mu ^-\}{\bar{\nu }}) = (5.11 \pm 0.10)\%. \end{aligned} \end{aligned}$$

(14)

Our nominal result for the exclusive determination of $|V_{cb}|$, obtained by combining all available theoretical and experimental information, is:

$$\begin{aligned} \left| V_{cb}^\text {excl}\right| = \left( 40.3 \pm 0.8\right) \times 10^{-3}. \end{aligned}$$

(15)

Its agreement with the individual values from ${\bar{B}}\rightarrow D^{(*)}\ell ^-{\bar{\nu }}$ is excellent. Averaging the two exclusive determinations with the inclusive one [35], we find $|V_{cb}|=(41.3\pm 0.5)\times 10^{-3}$, where the three values are compatible at the $1.2\sigma$ level.

4 Summary and outlook

In this work we carry out a comprehensive analysis of the full set of ${\bar{B}}\rightarrow D^{(*)}\ell ^-{\bar{\nu }}$ form factors. The basis of our analysis is the Heavy-Quark Expansion (HQE) up to ${\mathcal {O}}(\varepsilon ^2)$, where our power-counting is defined as $\varepsilon _b\sim \varepsilon _c^2\sim \alpha _s/\pi \sim \varepsilon ^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\rightarrow c$ transitions with light leptons. Our work focuses on two applications: precision predictions for ${\bar{B}}\rightarrow D^{(*)}\tau ^-{\bar{\nu }}$ observables and accurate determinations of $|V_{cb}|$ from ${\bar{B}}\rightarrow D^{(*)}\lbrace e^-,\mu ^-\rbrace {\bar{\nu }}$ 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 ${\mathcal {O}}(1/m_c^2)$. We find that the expansion in $1/m_Q$ is well-behaved, similar to what has been found in a recent analysis of $\Lambda _b\rightarrow \Lambda _c$ form factors [36]. Based on these findings we expect the terms at ${\mathcal {O}}(\varepsilon ^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 for ${\bar{B}}\rightarrow D^{(*)}\tau ^-{\bar{\nu }}$ 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\sigma$ level. Predictions with and without the use of experimental inputs are given in Eqs. (13) and (12), 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\sigma$ 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

$$\begin{aligned} \left| V_{cb}^\text {excl}\right| = \left( 40.3 \pm 0.8\right) \times 10^{-3}. \end{aligned}$$

The upcoming lattice analyses of four of the ${\bar{B}}\rightarrow D^*$ form factors at nonzero recoil [37,38,39] will benefit our approach and help to determine the HQE parameters to even higher precision.

Appendices

Appendix A: Statistical treatment of the unitarity bounds

We consider a positive semi-definite function $B(\vec {\theta }\,)$, where $\vec {\theta }$ denotes the parameters of interest. In the context of our work, B corresponds to any of the previously discussed strong unitarity bounds. Moreover, we have $B(\vec {\theta }\,) \le \chi _B$. The probability density of the parameters $\vec {\theta }$ can then be expressed as

$$\begin{aligned} P(\vec {\theta }\,) = \frac{N(\vec {\theta }\,)}{\int d\vec {\theta }\, N(\vec {\theta }\,)}, \end{aligned}$$

(A1)

where we use

$$\begin{aligned} N(\vec {\theta }\,) = \int _{-\infty }^{+\infty } d\chi _B\,\mathbb {1}\left[ B(\vec {\theta }\,) \le \chi _B\right] \, P(\chi _B). \end{aligned}$$

(A2)

In the above $\mathbb {1}$ denotes the indicator function, and $P(\chi _B)$ is the probability density for the upper bound on B. For the case that $\chi _B$ is precisely known, $P(\chi _B)$ can be approximated as a Dirac delta and consequently $P(\vec {\theta }\,)$ is given by the indicator function. For our application, the upper bounds carry a significant theoretical uncertainty, which we model as independent Gaussian distribution for each bound, centered around a value $\mu _B$ and with a standard deviation of $\sigma _B$. In that case we obtain:

$$\begin{aligned} -2 \ln N(\vec {\theta }\,) \simeq {\left\{ \begin{array}{ll} 0 &{} B(\vec {\theta }\,) \le \mu _B\\ \left( \frac{B(\vec {\theta }\,) - \mu _B}{\sigma _B}\right) ^2 &{} \text {otherwise} \end{array}\right. }. \end{aligned}$$

(A3)

Appendix B: Form factor definitions

1.1 1. $\bar{B}\rightarrow D$

The $\bar{B}\rightarrow D$ form factors have been defined in Eqs. (1)–(2). Their translation to the heavy-quark matrix elements defined in Ref. [12] reads

$$\begin{aligned} f_+&= \frac{1}{2\sqrt{r}}\,\left[ (1 + r) h_+ - (1 - r) h_-\right] , \end{aligned}$$

(B1)

$$\begin{aligned} f_-&= \frac{1}{2\sqrt{r}}\,\left[ (1 + r) h_- - (1 - r) h_+\right] ,\end{aligned}$$

(B2)

$$\begin{aligned} f_0&= f_+ + \frac{1 + r^2 - 2 r w}{1 - r^2}\, f_-,\end{aligned}$$

(B3)

$$\begin{aligned} f_T&= \frac{1 + r}{2 \sqrt{r}}\, h_T, \end{aligned}$$

(B4)

where $r \equiv M_D / M_B$.

1.2 2. $\bar{B}\rightarrow D^*$

The form factors for $\bar{B}\rightarrow D^*$ are defined in Eqs. (4)–(6). The factors in Eq. (6) have been chosen such that the conventions of Ref. [22] for the following contractions are recovered:⁵

$$\begin{aligned}&\langle D^*(k, \eta )| \bar{c} q^\nu \sigma _{\mu \nu } b |{\bar{B}}(p)\rangle \nonumber \\&\quad = -2 i \epsilon _{\mu \nu \rho \sigma } \eta ^{*\nu }p^\rho k^\sigma T_1(q^2), \end{aligned}$$

(B5)

$$\begin{aligned}&\langle D^*(k, \eta )| \bar{c} q^\nu \sigma _{\mu \nu }\gamma _5 b |{\bar{B}}(p)\rangle \nonumber \\&\quad = \eta ^*_\alpha \left\{ \left[ (M_B^2-M_{D^*}^2)g^{\alpha }_{\mu }-p^\alpha (p+k)_\mu \right] T_2(q^2)\right. \nonumber \\&\qquad \left. +p^\alpha \left[ q_\mu -\frac{q^2}{M_B^2-M_{D^*}^2}(p+k)_\mu \right] T_3(q^2)\right\} . \end{aligned}$$

(B6)

Translating again to the heavy-quark form factors in [12], we find⁶

$$\begin{aligned} V&= \frac{1 + r}{2 \sqrt{r}} h_V,\end{aligned}$$

(B7)

$$\begin{aligned} A_0&= \frac{1}{2 \sqrt{r}} \left[ (w + 1) h_{A_1} + (r w - 1) h_{A_2} + (r - w) h_{A_3}\right] ,\end{aligned}$$

(B8)

$$\begin{aligned} A_1&= \frac{\sqrt{r}\,(1 + w)}{1 + r} h_{A_1},\end{aligned}$$

(B9)

$$\begin{aligned} A_{12}&= \frac{1 + w}{8 \sqrt{r}}\left[ (w - r) h_{A_1} - (w - 1) r h_{A_2} - (w - 1) h_{A_3}\right] ,\end{aligned}$$

(B10)

$$\begin{aligned} T_1&= -\frac{1}{2\sqrt{r}}\left[ (1 - r) h_{T_2} - (1 + r) h_{T_1}\right] ,\end{aligned}$$

(B11)

$$\begin{aligned} T_2&= \frac{1}{2\sqrt{r}}\left[ \frac{2 r (w + 1)}{1 + r} h_{T_1}- \frac{2r(w - 1)}{1 - r} h_{T_2}\right] ,\end{aligned}$$

(B12)

$$\begin{aligned} T_3&= \frac{1}{2\sqrt{r}}\left[ (1 - r) h_{T_1} - (1 + r) h_{T_2} + \left( 1 - r^2\right) h_{T_3} \right] , \end{aligned}$$

(B13)

where $r \equiv M_{D^*} / M_B$.

1.3 3. Details on the HQE for the form factors

A generic heavy-quark form factor h is expanded in $\alpha _s/\pi$, $\varepsilon _Q \equiv \bar{\Lambda }/2m_Q$, and z as follows:

$$\begin{aligned} h(z)= & {} \xi (w){\hat{h}}(w) \nonumber \\= & {} \xi (w)\left( a + \frac{\alpha _s}{\pi }b\left[ {\hat{C}}_i(w)\right] + \varepsilon _b \, c_b\left[ {\hat{L}}_i(w)\right] \right. \nonumber \\&\left. + \varepsilon _c \, c_c\left[ {\hat{L}}_i(w)\right] + \varepsilon _c^2 \, d\left[ {\hat{\ell }}_i(w)\right] \right) _{w = w(z)}\, \end{aligned}$$

(B14)

where $a=1,0$, depending on whether the form factor vanishes in heavy-quark limit or is proportional to the leading IW function $\xi (w)$, and b, $c_b$, $c_c$, and d represent linear functionals of their arguments. The $\alpha _s$ corrections are written in terms of the leading order results for the HQET Wilson coefficients ${\hat{C}}_i(w)$; see [40] for a review. The $\mathcal {O}\left( \varepsilon _Q\right)$ corrections fulfill equations of motion that reduce the number of independent functions $\hat{L}_i$ from six to four. These four are commonly denoted as ${\hat{\chi }}_{1,2,3}(w)$ and ${\hat{\eta }}(w)$, where heavy-quark spin symmetry allows ${\hat{\chi }}_1(w)$ to be absorbed into $\xi (w)$. Hence, only three independent functions remain [40]. Similarly, equations of motion relate the functions at $\mathcal {O}\left( \varepsilon _Q^2\right)$ with each other. However, the set of independent functions is comprised by more than six functions [23]. Since we are only interested in the $\varepsilon _c^2$ corrections, it suffices to define the first six IW functions at $\mathcal {O}\left( \varepsilon _Q^2\right)$ as the six functions arising at $\mathcal {O}\left( \varepsilon _c^2\right)$. Finally, all functions are expanded in z according to Eq. (10), to the order given in the corresponding model.