The $V_{cb}$ puzzle: an update

We discuss the impact of the recent untagged analysis of ${B}^0\rightarrow D^{*}l\bar{\nu}_l$ decays by the Belle Collaboration on the extraction of the CKM element $|V_{cb}|$ and provide updated SM predictions for the $b\to c\tau\nu$ observables $R(D^*)$, $P_\tau$, and $F_L^{D^*}$. The value of $|V_{cb}|$ that we find is about $2\sigma$ from the one from inclusive semileptonic $B$ decays, and is very sensitive to the slope of the form factor at zero recoil which should soon become available from lattice calculations.

Lattice calculations for the B → D * ν channel have not yet reached the same level of maturity, and there are published results only for the form factor at zero recoil, which implies an extrapolation of experimental data, until recently performed by the experimental collaborations using the Caprini-Lellouch-Neubert (CLN) parametrization [10]. Two years ago, the Belle collaboration published a preliminary tagged analysis with results for the q 2 and angular distributions of B → D * ν [11]. These unfolded distributions allowed, for the first time, independent fits to the form factors and |V cb | using different parametrizations, with the surprising result that |V cb | could vary by as much as 6% between the CLN and the Boyd-Grinstein-Lebed (BGL) [12] parametrizations, with the latter lifting |V cb | to agreement with (2) [13][14][15]. The CLN and BGL parametrizations are both built on the analytic properties of the form factors, which together with the operator product expansion applied to correlators of two hadroniccb currents allow us to constrain them significantly.
These constraints can be made stronger using HQET relations between B ( * ) → D ( * ) form factors, known at O(1/m, α s ), supplemented by QCD sum rules for the subleading Isgur-Wise functions, see [10,16] and references therein. Ref. [10] (CLN) provides a simplified parametrization which includes these strong constraints, but the uncertainties related to missing higher order contributions and to the simplified parametrization have never been fully addressed before 2017. In [17] two of us implemented the strong unitarity bounds taking into account recent lattice calculations as well as conservative estimates of the theoretical uncertainties, and showed that they reduce the BGL vs CLN discrepancy in |V cb |, but do not eliminate it.
As emphasized in [13], the large observed parametrization dependence could be specific to the only data set available at that time for the B → D * ν differential distributions. What was certainly clear from the 2017 analyses was the pivotal role of precise information on the form factors in the small recoil region, and in particular arXiv:1905.08209v2 [hep-ph] 13 Jun 2019 their zero-recoil slope.
Several new relevant lattice calculations have appeared since 2016, showing that the lattice community has recognized the crucial importance of the form factor determination for a resolution of the puzzle. First, the HPQCD collaboration has computed the B → D * ν form factor at zero recoil [18], confirming the result by FNAL-MILC, although with less precision. Their recent results for B s → D * s ν [19], obtained using the same action for the b quark as for the lighter ones and thereby circumventing a large source of systematic uncertainty, appear very promising. Preliminary results for the B → D * ν form factors at non-zero recoil have been presented by the JLQCD [20] and FNAL/MILC [21] collaborations. JLQCD has shown that large deviations from HQS in R 1 (a form factor ratio, see below) seem to be excluded, while FNAL/MILC has shown blinded results for all the four B → D * ν form factors, albeit without complete systematic uncertainties. In particular, the slope of the combination of form factors that enters the differential rate in w, F(w), appears to be determined with good accuracy at zero recoil. Still concerning form factor calculations, improved and updated results at q 2 ≤ 0 with Light Cone Sum Rules (LCSR) have been presented in [22].
On the experimental side, a significant effort has been devoted to reanalyzing Babar and Belle data. The Babar collaboration recently published a new tagged analysis of B → D * ν using the BGL parametrization [23], performing an unbinned four-dimensional fit. The reported value of |V cb | is (38.4 ± 0.2 ± 0.6 ± 0.6) × 10 −3 , but, unfortunately, their data are not yet available for independent analyses.
The purpose of this note is to discuss the latest development, namely a new untagged analysis of B → D * ν by the Belle collaboration [24], and to provide a first assessment of the global resulting situation.
The differential B → D ( * ) ν data also provide bounds on new physics. Indeed, the known differential distributions already place stringent bounds on all possible single-mediator models of new physics [25], and none of these models can contribute to solve the V cb puzzle in a significant way [25,26], see also [27]. Even though it is unlikely to signal new physics, the |V cb | puzzle is still very important because i) it may be a signal that something is not yet understood in either the inclusive or exclusive analysis, with possible implications on the interpretation of R(D ( * ) ) as well, and ii) it limits the accuracy of |V cb | which affects FCNC studies in an important way.

II. FITTING THE NEW DATA
Like in Ref. [11], the authors of Ref. [24] provide onedimensional distributions in the variables w, cos θ l , cos θ v , and χ, with 10 bins for each variable. Unlike Ref. [11], however, they do not provide unfolded distributions, but give the efficiencies and response functions necessary to fold theoretical predictions in order to get predictions for the yield in each bin. They also perform binned fits to their data with both CLN and BGL parametrizations and find very similar results for |V cb |, (38.4 ± 0.2 ± 0.6 ± 0.6)×10 −3 and (38.3±0.3±0.7±0.6)×10 −3 , respectively (here the errors are statistical, systematic, and due to the lattice form factor at zero recoil).
These results are very different from those based on [11]: there is no sign of parametrization dependence in |V cb |. A first possibility is then that the two datasets are incompatible. To investigate this point, we take the unfolded distributions given in [11] and fold them in the experimental environment of Ref. [24], propagating the errors 2 . The bin by bin comparison of the yields for the sum of electrons and muons is shown in Fig. 1. Despite visible deviations in a couple of bins, the visual impression is of general compatibility. We will discuss the issue in a more quantitative way when we will present fits including both datasets. Fig. 1 shows clearly that the 2018 data are considerably more precise than the 2017 ones.
It is useful to recall at this point that in the BGL framework the four form factors relevant in B → D * ν, f, g, F 1 and F 2 (the latter entering only for = τ ), are written as power series in the variable z = ( The outer functions φ i (z) and the Blaschke factors P i (z) are given, for instance, in [13,17]. 3 The z-expansion is truncated at order n i , which may differ among different form factors. The weak unitarity constraints on the parameters a i k are ng k=0 (a g k ) 2 < 1, they ensure a rapid convergence of the z-expansion over the whole physical region, 0 < z < 0.056. The Belle analysis [24] employs n f = 1, n g = 0, n F1 = 2, which we denote by BGL (102) . 4 We have performed fits to the 2018 data of [24] using the information provided in that paper. A BGL (102) 2 The two analyses differ slightly in the endpoint for w, which leads to negligible differences in the angular bins but has a visible impact on bin 10. 3 The P i (z) depend on the masses of the Bc resonances below the lowest threshold. We use Table III of [17] updating the masses of the second 0 − state to 6.871GeV and of the second 1 − state to 6.910GeV due to their experimental discovery [28][29][30]. This has a negligible impact on all of our results. The normalization of the outer function has been updated in [13,17]. 4 Notice that a F 1 0 is fixed by the value of a f 0 , cfr. Eq. (9) of [13]. Note further that our notation BGL (n 1 ,n 2 ,n 3 ) corresponds to BGL n 2 +1,n 1 +1,n 3 in Ref. [31]. fit without systematic errors and using the same inputs roughly reproduces the results reported in the paper, with |V cb | about 1% higher, see Table I. The same applies to fits performed on the electron and muon data separately. Notice that, unless specified, we always perform constrained fits subject to Eq. (5). 5 As expected, the inclusion of systematic errors and correlations has a considerable impact on the result of the fits, and it generally increases the central value of |V cb |, see Table I. Of course, a fit based on statistical uncertainties only, is likely to be biased, because it gives too much weight to bins with larger systematic uncertainties, like those at small w which depend on the soft pion reconstruction.
On the other hand, the high degree of correlation of the systematic errors for the angular bins suggests some caution. We will perform specific tests on the stability of the fit later on. For the moment, we observe that the complete covariance matrix does not show correlations exceeding 0.94 and that the correlations are generally slightly higher than in the 2017 analysis. An important point concerns the implementation of the systematic uncertainties, which Ref. [24] provides as relative uncertainties. It is well-known that computing the systematic uncertainty as a fraction of the yield in each bin can lead to a bias (the D'Agostini effect [32]), which however can be avoided by expressing the systematic uncertainty as a fraction of the predicted yield. Indeed we observe a significant shift in |V cb | due to this effect, see Table I. There is a residual ambiguity depending on the form factors employed to predict the yields, but it is numerically very small, as the main effect is related The |V cb | error always includes the lattice uncertainty. In the second column stat stands for "only statistical errors", naive stands for "systematic errors as fractions of the yield", in all other cases we take the D'Agostini effect into account (see text). All fits are with weak unitarity constraints.
to employing a prediction that is not subject to fluctuations. In Table I and in the following we will always compute the systematic errors from predictions based on the form factors obtained in a fit where the systematic errors are a fraction of the yields, unless explicitly stated.
In BGL fits the power of z at which the series in (4) is truncated is potentially important for the extraction of |V cb |. For instance, the optimal choice of (n f , n g , n F1 ) in BGL fits has been recently discussed in [31]. We believe that, generally speaking, the problem has a simple solution: the optimal truncation of the z-expansion oc-curs when adding more terms does not change the result of the fit in any relevant way. Eqs. (5) together with 0 ≤ z ≤ 0.056 guarantee the convergence of such a procedure. Although this may imply adding (almost) redundant parameters subject to Eq. (5), it is crucial for determining the uncertainty of |V cb | in a reliable way. We illustrate the point by comparing the BGL (102) fit with a BGL (222) fit, having checked that nothing changes by adding even more parameters 6 . The total uncertainty increases from 0.9 to 1.4 10 −3 , which we think is the correct uncertainty of |V cb | in a BGL fit. The argument that a certain parameter can be dropped because the fit is unable to constrain it effectively is, in this particular case, ill-conceived. The BGL parametrization is not modelindependent if one arbitrarily drops parameters. From now on we will limit ourselves to BGL (222) fits only.
Ref. [31] also mentions the risk of overfitting. Imposing at least weak unitarity, which is avoided in [31], minimizes this risk and is completely safe, because the unitarity constraints (5) are very far from being saturated by the B → D * channel alone, see [17].
Let us now consider a fit to the combined 2017 [11] and 2018 [24] Belle datasets. Unlike the previous fits, where we were comparing directly with [24], we now employ the FLAG average for the form factor at zero-recoil, h A1 (1) = 0.904(12) [9]. The complete results of this fit are given in Table II: they show a marked increase in the minimal χ 2 /dof , implying some tension between the 2017 and 2018 data. Nevertheless, the combined fit still has an excellent p-value of ∼ 24%. In Fig. 2 we compare our fit result for η 2 EW |V cb | 2 |F(w)| 2 with the two Belle data sets. In order to show the data points of Ref. [24] in the same plot with those of Ref. [11], we employ an effective binby-bin rescaling factor obtained by comparing yields and binned differential branching fraction in the case of our best fit. We have performed a few checks on the stability of this fit: first, we have removed a few bins from the 2018 analysis, aiming at eliminating the strongest systematic correlations, and we did not observe any relevant change in |V cb |. If we remove all angular bins we get almost the same central value with larger uncertainty: |V cb | = (39.8 ± 1.5)10 −3 . The w bins are more important for the determination of |V cb |, and the first two in particular; the fit prefers lower |V cb | only if we remove the first two w bins of both 2017 and 2018 analyses, otherwise it is almost unchanged.
As mentioned in the introduction, the weak unitarity constraints of Eq. (5) can be made stronger using additional information related to Heavy Quark Symmetry. In Table II we report the results of a fit that adopts the strong unitarity bounds derived in [17]. 7 Interestingly,  the results do not differ significantly from the fit with weak unitarity bounds, in contrast to analogous fits to the 2017 data only [17]. It seems that the new and more precise data bring the fit naturally closer to the physical region, even in the absence of strong unitarity bounds. Another feature of the fits to the 2017 data presented in [13,17] was that the vector form factor g(z) grew with z (or decreased with q 2 ). This behaviour is unphysical and led to strong deviations from the HQET expectation [17]. The fits in Table II do not show this pathological behaviour. Like in Refs. [13,17] we also study the inclusion of LCSR results at q 2 = 0 in the fits, employing a recent updated analysis [22]: there seems to be excellent compatibility and |V cb | is basically unchanged, both with weak and strong unitarity bounds, see Table I.
As discussed in the Introduction, at the moment we are unable to include the recent Babar results [23] in our fit, and all previous Babar analyses report results only in the CLN parametrization. However, the total B → D * ν branching fraction is essentially independent  Table II including weak and strong unitarity bounds, and to a fit with LCSR inputs and weak bounds, in order of decreasing uncertainty. The respective HQET estimates are between the dashed lines.
of the parametrization employed. We have checked that including the previous Babar results for the total branching fraction leaves the reference fits of Table II almost  unaffected. As already done in [13,17] we compare our results with expectations based on NLO HQET, supplemented by QCD sum rules [16], for which we use conservative error estimates [17]. In particular, we show results for the two ratios of form factors Previous fits to the 2017 Belle data showed a marked discrepancy of R 1 with HQET, likely due to the unphysical behaviour of g discussed above. The plot in Fig. 3 shows the predictions based on the fits of Table II (left). The uncertainties of the fit with weak unitarity constraints are as large as those in the fit to 2017 data only, but now there is everywhere good agreement with HQET. The large uncertainty in R 1 at low and high recoil is due to low sensitivity to g(z) in these two regions. Using the strong constraints, the uncertainty in those two regions decreases, and it becomes even smaller when using LCSR results in the fit, see Fig. 3. As mentioned above, better knowledge of the form factors in the small recoil region would improve significantly the determination of |V cb |. In [13] this was explicitly illustrated with a fit where we assumed that a future lattice calculation would provide the slope of F(w) at zero recoil. Here, we repeat the exercise taking inspiration from the preliminary plots shown in [21] (although the results are still blinded, the slope of F(w) depends only marginally on the blinding factor). We adopt dF/dw| w=1 = −1.40 (7). The 5% uncertainty appears a realistic goal for the calculations currently in progress. This value shows some tension with the 2018 data at small recoil, but while its inclusion in the fit increases the total χ 2 by about 4.4, see Table I, it does not compromise the overall quality of the fit. At the same time |V cb | increases by over 1σ. This simple exercise does not anticipate in any way the final results of the FNAL/MILC collaboration; its only purpose is to illustrate the potential impact of lattice calculations on the fit. Inclusion of strong unitarity bounds and of LCSR does not change this picture, see Table I.
Finally, let us comment on the binning chosen in [11,24] for the angular variables. It is known that the single angular differential rates have a very simple form that can be parametrized in terms of only 3 parameters in the cases of θ l and χ, and only 2 parameters in the case of θ v , even beyond the SM. Rather than using 10 highly correlated bins, completely integrated over q 2 , taking the first few moments of cos θ l,v , sin θ l,v or their analogue in χ in q 2 bins would enhance the sensitivity of the analysis, a point emphasized also in Ref. [23].

III. SEMITAUONIC DECAYS
The results presented in the previous Section allow us to provide predictions for three quantities related to semitauonic decays: we update our predictions for R(D * ) (the ratio of semitauonic to light lepton widths) and for the τ polarization asymmetry P τ [17], and we compute the longitudinal polarization fraction of the D * , F D * L . There is a new form factor that enters semitauonic decays, the pseudoscalar form factor, which is unconstrained by the present experimental data and whose calculation on the lattice has not yet been completed. Here, to constrain its values we follow the third method employed in [17]: it is based on a kinematic relation linking it to F 1 at maximum recoil and on the use of an HQET relation with conservative uncertainties at zero recoil. We obtain where we use weak unitarity only and no LCSR input. In comparison with [17] the error for R(D * ) is reduced by 20% (but remains larger than in [15,16]) and the central value is about 1σ lower. The discrepancy of our SM prediction for R(D * ) with the experimental world average 0.295(11)(8) [1] is therefore now 2.8σ. On the other hand, our prediction for P τ is almost unchanged, and of course agrees with the experimental measurement P τ = −0.38(51)(21) by Belle [33]. Our new F D * L prediction is in good agreement with previous estimates [34][35][36] and 1.4σ from the recent experimental measurement F D * L = 0.60(8)(4) by the Belle collaboration [37].  [3,6,13] and this work, based on the quoted lattice QCD and experimental results. We show here the results obtained with weak unitarity constraints and no LCSR input only.

IV. CONCLUSIONS
In this paper we have studied the impact of a new Belle untagged analysis of the B → D * ν decay. When we analyse the data of Ref. [24] we obtain a value of |V cb | 2.1% higher than reported there, with an uncertainty about 50% larger. Including in the fit the previous tagged analysis by Belle [11] we get |V cb | = (39.6 +1.1 −1.0 ) × 10 −3 (11) which still differs from the inclusive determination by about 1.9σ and is in excellent agreement with the determination from B → D ν, see the overview that we provide in Fig. 4. We find that the inclusion of strong unitarity bounds and of LCSR results at maximum recoil in the fit does not change the central value of |V cb |, although it helps constraining the individual form factors. As a byproduct of our analysis, we provide in Eqs. (8)(9)(10) updated predictions for R(D * ), P τ , and F D * L .
We also show that higher values of |V cb | may still be compatible with the available data. Indeed, preliminary results of lattice calculations suggest a slope of the relevant form factor F(w) at zero recoil steeper than expected from the experimental data. We have shown that if such a high value for the slope were confirmed, |V cb | extracted from a global fit to B → D * ν data would agree with the inclusive determination. In other words, it is lattice QCD that will decide the eventual fate of the |V cb | puzzle.