Perturbative corrections to $B \to D$ form factors in QCD

We compute perturbative QCD corrections to $B \to D$ form factors at leading power in $\Lambda/m_b$, at large hadronic recoil, from the light-cone sum rules (LCSR) with $B$-meson distribution amplitudes in HQET. QCD factorization for the vacuum-to-$B$-meson correlation function with an interpolating current for the $D$-meson is demonstrated explicitly at one loop with the power counting scheme $m_c \sim {\cal O} \left (\sqrt{\Lambda \, m_b} \right ) $. The jet functions encoding information of the hard-collinear dynamics in the above-mentioned correlation function are complicated by the appearance of an additional hard-collinear scale $m_c$, compared to the counterparts entering the factorization formula of the vacuum-to-$B$-meson correction function for the construction of $B \to \pi$ from factors. Inspecting the next-to-leading-logarithmic sum rules for the form factors of $B \to D \ell \nu$ indicates that perturbative corrections to the hard-collinear functions are more profound than that for the hard functions, with the default theory inputs, in the physical kinematic region. We further compute the subleading power correction induced by the three-particle quark-gluon distribution amplitudes of the $B$-meson at tree level employing the background gluon field approach. The LCSR predictions for the semileptonic $B \to D \ell \nu$ form factors are then extrapolated to the entire kinematic region with the $z$-series parametrization. Phenomenological implications of our determinations for the form factors $f_{BD}^{+, 0}(q^2)$ are explored by investigating the (differential) branching fractions and the $R(D)$ ratio of $B \to D \ell \nu$ and by determining the CKM matrix element $|V_{cb}|$ from the total decay rate of $B \to D \mu \nu_{\mu}$.


Introduction
Precision calculations of the semileptonic heavy-to-heavy B → D ν decays are indispensable for a stringent test of the CKM matrix element |V cb | exclusively and for a deep understanding of the strong interaction mechanism in the heavy-quark system from both QCD and heavyparticle effective field theories. The long-standing tension between the exclusive and inclusive determinations of |V cb | [1] as well as the topical R(D) ≡ BR(B → Dτ ν τ )/BR(B → D ν) anomaly [2] necessitates further QCD calculations of B → D form factors with an increasing accuracy. Employing the heavy-quark expansion (HQE) technique, systematic studies of B → D form factors near zero recoil were performed including both perturbative QCD corrections and subleading power contributions (see [3] for a review). Very recently, the unquenched lattice-QCD calculations of B → D ν form factors near zero recoil were reported by the FNAL/MILC Collaboration [4] using the asqtad-improved fermions for the light valence quarks and the improved Wilson ("clover") fermions for the heavy valence quarks and by the HPQCD Collaboration [5] independently applying the NRQCD action for bottom and the highly improved staggered quark action for charm quarks together with N f = 2 + 1 MILC gauge configurations. However, extrapolating the HQE and lattice calculations of B → D form factors to the entire physical kinematic range can be only achieved with certain parametrizations for the momentum dependence of form factors, which introduce an additional source of theoretical uncertainties for the determination of |V cb |. A direct QCD computation of the form factors of B → D ν at large recoil is therefore highly in demand to provide complementary information on the hadronic dynamics for a better understanding of the form-factor shapes. The major objective of this paper is to carry out perturbative QCD corrections to the two-particle contributions to B → D form factors at large recoil from the light-cone sum rules (LCSR) with B-meson distribution amplitudes (DA) following the techniques developed in [6][7][8][9], in an attempt to extend the leading-order calculation of the corresponding vacuum-to-B-meson correlation function performed in [10].
Applying the light-cone operator-product expansion (OPE) and parton-hadron duality, the LCSR approaches have been proven their usefulness in computing both the local and non-local hadronic matrix elements for the theory description of hadronic processes with large momentum transfer. Demonstrating QCD factorization for the vacuum-to-B-meson correlation function at leading power in Λ/m b is evidently the first step in the construction of sum rules for B → D ν form factors. To this end, we need to establish the power counting scheme for the distinct momentum scales involved in the correlation function under consideration. In contrast to the correlation function for the heavy-to-light form factors, the appearance of an additional energy scale (the charm-quark mass) further complicates perturbative QCD factorization for the vacuum-to-B-meson correlation function with an interpolating current for the D-meson. Motivated by the mass hierarchy between the bottom and charm quarks numerically, we will adopt a novel power counting scheme m c ∼ O √ Λ m b , as implemented in the analysis of the inclusive semileptonic B → X c ν decays [11,12] with shape functions of the B-meson, instead of the popular counting scheme m c ∼ O(m b ) widely used in the heavyquark physics (see, for instance [13,14]). It is then apparent that the on-shell bottom-quark field, the charm-quark field appeared in the interpolating current for the D-meson and the external light-quark field can be identified as hard, hard-collinear and soft modes in QCD are computed manifestly at one loop. The (partial) NLL resummation improved sum rules for the two-particle contributions to B → D form factors will be also derived here and they constitute the main new results of this paper. In Section 4 we further calculate the threeparticle contributions to the LCSR of B → D ν form factors at tree level, which will be shown to contribute only at subleading power in Λ/m b with the power counting scheme adopted here. Phenomenological implications of the newly derived sum rules for the decay form factors of B → D ν will be explored in Section 5 with different models for B-meson DA, including determinations of the shape parameters for the vector and scalar B → D form factors with both the z-series expansion [32] and the Caprini-Lellouch-Neubert (CLN) parametrization [33], the differential branching fractions of B → D ν and R(D) defined in the above, as well as the extraction of the CKM matrix element |V cb |. We will conclude in Section 6 with a summary of main observations and a perspective on the future refinements. We collect some useful results of loop integrals for evaluating the vacuum-to-B-meson correlation function at one loop, present the spectral representations for the convolution integrals essential to construct the (partial) NLL LCSR of B → D form factors and display the lengthy coefficient functions entering the "effective" B-meson DA, absorbing the hard-collinear corrections at one loop, in Appendices A, B and C, respectively.

The framework
We briefly review the LCSR approach to compute B → D form factors with B-meson DA following the strategy presented in [15]. The vacuum-to-B-meson correlation function adopted to construct the sum rules for the form factors f + BD (q 2 ) and f 0 BD (q 2 ) is defined as Π µ (n · p,n · p) = i d 4 x e ip·x 0|T {q(x) n γ 5 c(x),c(0) γ µ b(0)} |B(p B ) , where p B ≡ m B v indicates the four-momentum of the B-meson and p refers to the fourmomentum carried by the interpolating current of the D-meson. We will work in the rest frame of the B-meson and introduce two light-cone vectors n µ andn µ satisfying n·v =n·v = 1 and n ·n = 2. We will further employ the following power counting scheme at large hadronic recoil, where q = p B − p is the transfer momentum and Λ is a hadronic scale of order Λ QCD . For the space-like interpolating momentum p, the leading power contribution to the correlation function (1) at tree level can be achieved by evaluating the diagram in figure  1 with the light-cone OPE and we obtain b q q p c Figure 1: Diagrammatical representation of the two-particle contributions to the vacuum-to-B-meson correlation function Π µ (n · p,n · p) defined in (1) at tree level.
where ω c = m 2 c /n·p and the convolution integral in the second step corresponds to the spectral representation of the considered correlation function. It is apparent that (3) reproduces precisely the result for the corresponding correlation function defined with a pion interpolating current in the m c → 0 limit.
The light-cone DA of the B-meson in the coordinate space are defined as follows [34,35]: where b v indicates the effective bottom-quark field in heavy-quark effective theory (HQET), the soft gauge link is given by and φ ± B (ω, µ) are obtained from the Fourier transformation ofφ ± B (t, µ). The static B-meson decay constantf B (µ) can be further expressed in terms of the QCD decay constant Employing the standard definitions for the decay constant of the D-meson and for B → D transition form factors it is straightforward to write down the hadronic dispersion relation for the correlation function at leading power in Λ/m b , where ω s is the effective threshold parameter for the D-meson channel. Applying the parton-hadron duality approximation for the hadronic dispersion integral and performing the Borel transformation with respect to the variablen · p lead to the LCSR for the two B → D form factors at tree level where the upper limit of the ω -integration is consistent with ω 0 (q 2 , s D 0 ) derived in [10] at leading power in Λ/m b keeping in mind the replacement rule s D → n · p ω s . In contrast to the sum rules for B → π form factors, the power counting rules of ω c , ω s and ω M now read Based upon the canonical behaviour for the B-meson DA φ − B (ω ), one can readily deduce the power counting of f + BD (q 2 ) at large hadronic recoil as O (Λ/m b ) 3/2 from the tree-level sum rules (9), different from the scaling law f + BD (q 2 ) ∼ O(1) [14] obtained with the counting scheme m c ∼ m b . In addition, the two form factors f +,0 BD (q 2 ) obtained in the above respect the large-energy symmetry relation as discussed in [35] and the symmetry breaking effects can arise from both perturbative QCD corrections to the short-distance coefficient functions and the subleading power contributions from the higher-twist DA of the B-meson.
3 Two-particle contributions to the LCSR at O(α s ) The purpose of this section is to demonstrate QCD factorization for the two-particle contribution to the vacuum-to-B-meson correlation function (1) Figure 2: Diagrammatical representation of the two-particle contribution to the vacuum-to-B-meson correlation function Π µ (n · p,n · p) defined in (1) at O(α s ).
at leading power in HQE. We will compute the hard matching coefficients C i,k (i = ±, k = n orn) and the jet functions J i,k at next-to-leading order (NLO) in α s manifestly and perform QCD resummation for the large logarithms in the hard functions at NLL accuracy with the RG evolution in momentum space. Since the perturbative matching coefficients must be independent of the external partonic states in the OPE calculation, we will choose the initial state to be of the minimal quark and gluon degrees of freedom |b(p b )q(k) for the convenience.

Perturbative matching functions at NLO
The NLO hard and jet functions entering the factorization formula (11) can be extracted by computing the leading power contributions of the one-loop QCD diagrams displayed in figure  2 from the hard and hard-collinear regions. It is evident that the soft contributions to the above-mentioned QCD diagrams at leading power in Λ/m b will be cancelled out precisely by the infrared subtraction terms following the presentation [15] and we will mainly focus on the hard and hard-collinear contributions of the one-loop diagrams in figure 2 in the following.

Weak vertex diagram
We are now ready to compute the one-loop correction to the weak vertex diagram displayed in figure 2 (a) where the external bottom and light quarks are taken to be on the mass-sell with the momenta m b v and k (with k 2 = 0). Employing the power counting scheme for the external momenta it is straightforward to identify that the leading power contributions to the scalar integral come from the hard, hard-collinear and soft regions of the loop momentum. Evaluating the hard contribution from the weak vertex diagram yields (15) where the two hard functions are given by with r = n · p/m b . The hard-collinear contribution to the weak vertex diagram can be further computed by expanding (12), in the hard-collinear region, at leading power in Λ/m b Applying the results of loop integrals presented in Appendix A leads to − ln (1 + r 1 ) + 1 + ln 2 µ 2 n · p (ω −n · p) + 2 ln µ 2 n · p (ω −n · p) − 2 ln(1 + r 1 ) ln µ 2 n · p (ω −n · p) + 1 + ln 2 (1 + r 1 ) with and ω ≡n · k. It is apparent that our result of J (weak) −,n reproduces the hard-collinear contribution to the weak vertex diagram, displayed in (29) of [15], for constructing the sum rules of B → π form factors in the m c → 0 (i.e., r 1 → 0) limit. Proceeding in a similar manner, we can extract the soft contribution to the weak vertex diagram by expanding (12) in the soft region which is precisely the same as the soft subtraction term defined by the convolution integral of the partonic DA of the B-meson at NLO in α s , calculable from the Wilson-line Feynman rules, and the tree-level short-distance function (see [15] for more discussion). We are then led to conclude that the soft QCD dynamics of the weak vertex diagram is indeed characterized by the B-meson DA at leading power in Λ/m b , independent of the renormalization scheme.

D-meson vertex diagram
We turn to compute the one-loop QCD correction to the D-meson vertex diagram displayed in figure 2 Inspecting the scalar integral and the Dirac algebra of (20) with the power counting scheme (13), one can verify that the leading power contributions to the D-meson vertex diagram only arise from the hard-collinear and soft regions of the loop momentum. As already discussed in [15], it turns out to be more apparent to compute the loop integrals directly and then to expand the obtained expression up to the leading power in Λ/m b , instead of applying the strategy of regions. Employing the results of loop integrals collected in Appendix A and the light-cone projector of the B-meson in momentum space derived in [35,36] yields where we have introduced the following jet functions with Here, we have taken advantage of the fact that the soft contribution to the D-meson vertex diagram vanishes in dimensional regularization (a similar observation already made in [15]). Several remarks on the resulting expressions for the QCD correction to the D-meson vertex diagram presented in (21) are in order.
• It is evident that the linear term in the charm-quark mass corresponds to the power enhanced effect compared to the remaining contributions due to the power counting . This observation can be readily understood from the fact that the strange-quark mass effect in B → K form factors is suppressed by m s /Λ, but not power of m s /m b , compared to the leading power contribution in HQE [17], and the charm-quark mass effect would naturally generate corrections to B → D form factors with a scaling factor of m c /Λ ∼ O( m b /Λ) in light of our power counting scheme. (This clearly does not imply an expansion of m c /Λ in the QCD calculation of the correlation function (1).) To develop a better understanding of such power-enhancement mechanism, we inspect the linear charm-quark mass term in (20) before performing the loop integration As the loop integral in the hard-collinear region induces a power-enhanced factor O (m b /Λ), one can then identify the scaling behaviour of the charm-quark mass term as O( m b /Λ) with respect to the tree-level contribution. Furthermore, one can investigate the charmquark mass effect directly in the context of B → D form factors applying the QCD factorization approach. A straightforward calculation of the spectator interaction diagram with a hard-collinear gluon exchange between the energetic charm quark and the light spectator quark yields where φ D (u) is the twist-2 DA of the D-meson on the light cone. It is evident that the second convolution integral in (28) suffers from the end-point divergence and the above-mentioned charm-quark mass effect should be identified as the non-factorizable contribution to B → D form factors.
• Since there is no power enhanced contribution proportional to the charm-quark mass at tree level, the NLO jet function J (D) +,n (n · p) must be infrared finite to validate the factorization formula (11) for the considered correlation function.
• Setting m c → 0 (namely r 2 → 0), our results of J (D) −,n (n · p) and J (D) +,n (n · p) will be reduced to the corresponding hard-collinear contributions from the pion vertex diagram as displayed in (38) of [15].
• Applying the Wilson-line Feynman rules, one can easily verify that the soft contribution to the D-meson vertex diagram, at leading power in Λ/m b , is cancelled out precisely by the corresponding soft subtraction term.

Wave function renormalization
We proceed to compute the self-energy correction to the intermediate charm-quark propagator shown in figure 2(c). With the expressions of loop integrals collected in Appendix A, we obtain which is apparently free of soft and collinear divergences. The resulting jet function reads where the ultraviolet divergence of J (wf c, 2) −,n will be subtracted after the charm-quark mass renormalization dependent on the subtraction scheme. We will employ the MS renormalization scheme for the charm-loop mass, which is appropriate for a Lagrange parameter entering the short-distance matching function in the OPE calculation, and more discussion on the renormalization schemes of the charm-quark mass can be found in [12].
It is straightforward to compute the matching coefficients from the wave function renormalization of the external quarks where Φ bq indicates the partonic DA of the B-meson defined in (12) of [15], the tree-level hard kernel T is given by

Box diagram
The one-loop QCD correction to the box diagram displayed in figure 2(d) can be further computed as Taking advantage of the power counting scheme (13), one can identify that the leading power contributions to the box diagram come from the hard-collinear and soft regions of the loop momentum. Expanding the QCD expression of the box diagram (36) in the hard-collinear region to the leading power in Λ/m b yields Employing the results of loop integrals collected in Appendix A and the momentum-space projector of the B-meson leads to where (26). One can readily verify that the resulting jet functions J (box) −,n and J (box) +,n reproduce the hard-collinear contribution to the one-loop box diagram for the vacuum-to-B-meson correlation function with a pion interpolating current as presented in (52) of [15], in the m c → 0 limit.
The soft contribution to the one-loop box diagram at leading power in Λ/m b can be further extracted from (36) as follows which is precisely the same as the soft subtraction term Φ (1) bq,box ⊗ T (0) µ computed with the Wilson-line Feynman rules. We then conclude that the soft dynamics of the vacuum-to-Bmeson correlation function under discussion is indeed parameterized by the light-cone DA of the B-meson in HQET.

The NLL hard and jet functions
Now we are in a position to present the one-loop hard and jet functions entering the factorization formula (11) with resummation of the large logarithms by solving the corresponding RG equations in momentum space. Putting different pieces together, we first derive the renormalized hard coefficients including the O(α s ) corrections and the renormalized jet functions at one-loop accuracy are given by It is apparent that the hard functions C −,n and C −,n are identical to the perturbative matching coefficients of the weak currentq γ µ (1 − γ 5 ) b from QCD onto SCET, as discussed in [15], due to the power counting rule of the charm-quark mass displayed in (2). However, the hardcollinear functions entering the factorization formula (11) turn out to be significantly more involved, due to the massive charm-quark effect, than the counterpart jet functions for the massless hard-collinear quark.
To verify the factorization-scale independence of the correlation function Π µ at O(α s ), we make use of the RG evolution equation of the charm-quark mass [37,38] it is then straightforward to write down where the three terms in the bracket arise from the RG running of the charm-quark mass, the jet function J −,n and the hard function C −,n , respectively. Applying the one-loop evolution where the renormalization kernel H (1) − (ω, ω , µ) can be found in [39,40], we can readily compute the last term of (48) as Substituting (50) into (48) immediately leads to indicating that the correlation function computed from the factorization formula (11) is indeed independent of the renormalization scale to the one-loop accuracy.
We will proceed to perform the summation of parametrically large logarithms in the perturbative expansion of the hard matching coefficients at NLL by applying the corresponding RG equations in momentum space. Following the arguments of [41], we will not aim at summing over logarithms of µ/µ 0 , with µ 0 being a hadonic scale of O(Λ), from the RG evolution of the B-meson DA φ − B (ω, µ), since we will take the factorization scale µ as a hard-collinear scale µ hc ∼ O √ Λ m b which is quite close to the hadronic scale µ 0 numerically. Employing the RG equations for the hard coefficient C −,n (n · p, µ) and the static B-meson decay constant where the various anomalous dimensions can be inferred from [41]. To achieve the resummation improved factorization formula for the considered correlation function at NLL accuracy, the cusp anomalous dimension Γ cusp (α s ) needs to be expanded at the three-loop accuracy, while the soft anomalous dimensions γ(α s ) andγ(α s ) need to be expanded up to two loops. The NLL resummation improved expressions for C −,n andf B can be further computed as where the explicit expressions of the evolution functions U 1 and U 2 can be found in [20]. It is then a straightforward task to deduce the (partial) NLL resummation improved factorization formula for the correlation function Π µ × C +,n (n · p, µ) J +,n (n · p, ω, µ) n µ + C +,n (n · p, µ) J +,n (n · p, ω, µ)n µ φ + B (ω, µ) where µ h1 and µ h2 are the hard scales of O(m b ) and the factorization scale µ needs to be taken at a hard-collinear scale O( √ m b Λ).

The NLL LCSR for B → D form factors
We are now ready to derive the (partial) NLL resummation improved sum rules for the vector and scalar B → D form factors. To this end, it is mandatory to work out the spectral representation of the factorization formula for Π µ obtained in the above, which turns out to be a nontrivial task compared to the construction of the B-meson LCSR for B → π form factors [15]. Applying the dispersion representations of convolution integrals collected in Appendix B yields where we have introduced the "effective" DA Φ eff i,k (ω , µ) (i = ±, k = n orn) absorbing the hard-collinear QCD corrections to the correlation function −,n (ω, ω ) with P indicating the principle-value prescription. The tedious expressions of the coefficient functions ρ (i) −,n in (59) are collected in Appendix C. Applying the standard strategy to construct the sum rules for hadronic transition form factors by matching (55) and (8) with the aid of the parton-hadron duality approximation and the Borel transformation leads to Figure 3: Diagrammatical representation of the three-particle contributions to the vacuumto-B-meson correlation function Π µ (n · p,n · p) defined in (1) at tree level.
which serves as the master formula for the two-particle contributions to B → D form factors obtained in this paper. It is evident that the symmetry-breaking corrections to the formfactor relation (9) arise from both hard and hard-collinear fluctuations as displayed in the last line of (60), and the power enhanced contribution due to the charm-quark mass effect preserves the large-recoil symmetry relation. The symmetry relations of B → D form factors at large hadronic recoil could be also obtained in SCET with massive collinear quark fields by generalizing the discussion for heavy-to-light form factors [42][43][44][45] and we leave a systematic investigation of the SCET formulation of B → D form factors for future work.
4 Three-particle contributions to the LCSR at tree level The objective of this section is to compute the tree-level contribution to the sum rules of B → D form factors from the three-particle B-meson DA. To this end, we first need to demonstrate QCD factorization for the three-particle contributions to the correlation function (1) with space-like interpolating momentum, which can be achieved by evaluating the diagram displayed in figure 3 with the aid of the light-cone OPE.
Employing the massive quark propagator near the light cone from the background gluon field method [46][47][48] with G µν = g s T a G a µν , it is then straightforward to derive the factorization formula where the coefficient functionsρ i,k (u, ω, ξ) are given bỹ The appeared three-particle DA of the B-meson can be defined by the position-space matrix element on the light cone [49,50] 0|q where we have neglected the soft Wilson lines to maintain the gauge invariance of the string operator on the left-hand side. The following conventions are further introduced in (63) for convenience. Our current understanding of the model independent properties of the three-particle quark-gluon B-meson DA is limited to the twist-3 DA Φ 3 (ω, ξ) ≡ Ψ A (ω, ξ) − Ψ V (ω, ξ) which was shown to be completely integrable in the large N c limit and to be solvable exactly at one loop [51]. Investigating perturbative QCD constraints on the higher-twist B-meson DA from the OPE analysis in the partonic picture and from the renormalization evolution equations in momentum (or "dual") space, along the lines of [52,53], would be interesting for future work. Expressing the factorization formula for the three-particle contributions to the correlation function Π µ, 3P in a dispersion form with respect to the variablen · p and implementing the continuum subtraction with the aid of the parton-hadron duality relation as well as the Borel transformation lead to the following sum rules where the explicit expressions of the coefficient functions F 3P,n and F 3P,n are It is evident from (66) that the large-recoil symmetry breaking effect between the vector and scalar B → D form factors can be induced by the higher-twist contributions from the three-particle B-meson DA, already at tree level.
We are now ready to derive the power counting behaviour of the three-particle correction to B → D form factors at tree level. Applying the canonical behaviour of the three-particle B-meson DA in the end-point region [7] we can verify that the three-particle contribution to the sum rules of B → D form factors (66) is counted as O (Λ/m b ) 5/2 in the heavy quark limit and is indeed suppressed by a factor of Λ/m b compared to the two-particle contribution at tree level as presented in (9). However, it needs to point out that higher twist effects from the three-particle gluon-quark DA can give rise to the leading power contributions to B → D form factors at O(α s ) on account of the reasonings of [42,54] and a transparent demonstration of this observation by computing the three-particle contribution to the sum rules at NLO in α s directly is in high demand on both conceptual and phenomenological aspects. The final expressions for the form factors of B → D ν can be obtained by adding up the two-particle and the three-particle contributions together where the detailed expressions of f +,0 BD, 2P (q 2 ) and f +,0 BD, 3P (q 2 ) are presented in (60) and (66), respectively.

Numerical analysis
Having at our disposal the resummation improved sum rules for B → D form factors, we are in a position to explore the phenomenological consequences of perturbative corrections to the two-particle contributions and higher-twist effects from the three-particle quark-gluon Bmeson DA at tree level. Employing the obtained theory predictions for B → D form factors, we will further present our results for the semileptonic B → D ν ( = µ, τ ) decay observables, including the invariant-mass distributions of the lepton pair, the topical R(D) ratio, and the CKM matrix element |V cb |.

Theory input parameters
We will first specify the theory inputs entering the LCSR for B → D form factors shown in (60) and (66). The fundamental nonperturbative quantities describing the soft strong interaction dynamics encoded in the correlation function (1) are the B-meson light-cone DA defined in (4) and (64). A detailed account of the current understanding towards perturbative and nonperturbative aspects of the two-particle B-meson DA was already presented in [15], following which we will consider two models of the leading twist DA φ + B (ω, µ 0 ) proposed in [34,55] where the shape parameter ω 0 can be converted to the inverse moment of the leading-twist B-meson DA λ B (µ 0 ). The renormalization scale evolution of λ B (µ) at one loop can be derived from the Lange-Neubert equation The inverse-logarithmic moment at a hadronic scale µ 0 = 1 GeV will be taken as σ 1 (µ 0 ) = 1.4 ± 0.4, determined from a QCD sum rule analysis [55]. The higher-twist two-particle DA of the B-meson φ − B (ω, µ 0 ) will be determined from the QCD equations of motion in the heavy quark limit [49] which was also demonstrated to be valid from a non-relativistic toy model manifestly at NLO in the strong coupling constant [39]. With regarding to the three-particle quark-gluon DA of the B-meson, we will employ the exponential model inspired from the canonical behaviour predicted by the HQET sum rules at tree level [7] Ψ where the normalization constant defined by the matrix element of the chromoelectric operator was estimated as λ 2 E (µ 0 ) = (0.03 ± 0.02) GeV 2 [57], from the two-point QCD sum rules including higher-order perturbative and nonperturbative corrections. We already implemented the approximation λ E = λ H in the above expressions, following [7], which is supported by the nonperturbative QCD calculations numerically [34,57].
With the matching relation (6) the static B-meson decay constantf B (µ) can be traded into the QCD decay constant f B , whose values will be taken from lattice QCD simulations f B = (192.0 ± 4.3) MeV [58] with N f = 2 + 1. Likewise, we will adopt the intervals for the Dmeson decay constant in QCD from lattice simulations f D = (209.2 ± 3.3) MeV [58] with N f = 2 + 1. In addition, we will employ the MS bottom quark mass m b (m b ) = (4.193 +0.022 −0.035 ) GeV [59] from non-relativistic sum rules and the MS charm quark mass m c (m c ) = (1.288 ± 0.020) GeV [60] from relativistic sum rules, employing the quark vector correlation function computed at O(α 3 s ). Following [15,41], the hard scales µ h1 and µ h2 are taken to be equal and will be varied in the interval [m b /2 , 2 m b ] around the default value m b , and the factorization scale is chosen as 1.0 GeV ≤ µ ≤ 2.0 GeV with the default value µ = 1.5 GeV.
The determination of the sum rule parameters ω M and ω s can be achieved by implementing the standard procedure described in [15] and applying the same strategies leads to which are in agreement with the intervals adopted in [10,61].

Predictions for B → D form factors
We will turn to discuss the choice of the inverse moment λ B (1 GeV) which serves as an important source for theory uncertainties, prior to presenting the sum rule predictions for the form factors of B → D ν. Albeit with the intensive investigations of determining λ B theoretically, the current constraints on λ B are still far from satisfactory, due to the emerged tension between the NLO sum rule predictions and the implications of hadronic B-meson decay data from QCD factorization. In particular, the subleading power contributions in HQE can give rise to sizeable impact on the determination of λ B numerically. This has been demonstrated explicitly by extracting λ B from the partial branching fractions of B → γ ν with the power suppressed effects estimated from the dispersion approach [20,62]. Following the above argument, it is very plausible that the unaccounted subleading power contributions to the LCSR for B → π form factors can yield significant corrections to the fitted values of λ B (µ 0 ) [15] in addition to the systematic uncertainty induced by the parton-hadron duality approximation. In order to be insensitive to the unconsidered effects in the sum rule determinations [15], we will perform an independent determination of λ B (µ 0 ) by matching the B-meson LCSR of the vector B → D form factor at q 2 = 0 to the lattice-QCD calculation with an extrapolation to the large recoil region using the z-series parametrization 1 . Taking f + BD (0) = 0.672 ± 0.027 [4] as an input and implementing the above-mentioned matching procedure lead to which differ from the intervals of λ B (µ 0 ) obtained from matching two different versions of sum rules for the vector B → π form factor [15], however, are comparable to the values determined with distinct QCD approaches [55,63]. In the following we will take φ + B,I (ω, µ 0 ) as our default model for the illustration purpose and the systematic uncertainty induced by the model dependence of φ + B (ω, µ 0 ) will be taken into account in the final predictions for the form factors of B → D ν. Now we are ready to explore the numerical features of the LCSR predictions for B → D form factors including the (partial) NLL resummation for the two-particle contributions and the subleading power corrections from the three-particle quark-gluon DA. To demonstrate the reliability of the sum rule calculations, we first display the dependence of f + BD (0) on the sum rule parameters M 2 and s 0 in figure 4. It is apparent that the variations of the Borel parameter and the effective threshold within the intervals (76) only bring about negligible influence on the sum rule predictions for the vector form factor of B → D ν. We further present the separate contributions to f + BD (0) from the two-particle and the three-particle contributions in figure 4 in an attempt to understand the subleading power corrections from the higher Fock states. It can be readily observed that the tree-level contributions from the three-particle Bmeson DA turn out to be of minor importance numerically, approximately O(1%), compared to the two-particle contributions. However, it is worthwhile to mention that the smallness of the three-particle contribution at leading order (LO) in perturbative expansion does not imply the insignificant impact of the subleading power contributions in the B → D ν decay amplitude in general, due to the yet unaccounted power suppressed effects induced by the off light-cone corrections to the nonlocal matrix element defining the two-particle B-meson DA, by the subleading power corrections to the perturbative coefficient functions entering the factorization formula (11) and by the additional contributions generated by new momentum regions (or equivalently, new field modes in the language of SCET) when applying the method of regions to the evaluation of loop integrals involved in perturbative corrections to the QCD amplitude.
We proceed to investigate the impact of perturbative corrections to the short-distance functions and resummation effects for the parametrically large logarithms on predicting the form factors of B → D ν. It is evident from figure 5 that NLO QCD corrections to the perturbative matching coefficients can reduce the tree-level prediction for f + BD (0) by approximately 10% at µ = 1.5 GeV and the (partial) NLL resummmation effect can enhance the NLO QCD calculation by an amount of 3% numerically at the same value of µ. Both the NLO and NLL predictions exhibit weaker dependencies on the factorization scale when compared to the LO result. Inspecting the different origins of perturbative QCD corrections displayed in figure  5 shows that the one-loop hard-collinear correction turns out to be more pronounced than the corresponding hard correction, with the inverse moment λ B (µ 0 ) = 570 MeV, at q 2 ≥ 0, highlighting the significance of computing the NLO jet functions accomplished in this paper. We further plot the theory predictions for the ratio [f + BD,2P (q 2 )] NLL /[f + BD,2P (q 2 )] LL in figure . Bottom: Breakdown of the form factor f + BD (0) from the two-particle and from the three-particle contributions. The two-particle contribution to f + BD (0) indicated by the dashed curve is almost indistinguishable from the total result represented by the solid curve, due to the negligible effect from the three-particle B-meson DA as shown by the dot-dashed curve. Figure 5: Left: Factorization scale dependence of the two-particle contributions to the form factor f + BD (0) computed from the sum rules (60). The curves labelled with "LL", "NLO" and "NLL" are obtained from the resulting predictions at leading-logarithmic (LL), NLO and (partial) NLL accuracy. Right: Breakdown of the two-particle contributions to f + BD (0) from the LL effect, from the NLL hard correction ("Hard") and from the NLO hard-collinear correction ("Jet"). 6 with uncertainties from the variations of both the hard and hard-collinear scales within the intervals given above. To develop a better understanding of the sensitivity of B → D form factors on the inverse moment λ B , we also present the LO, NLO and (partial) NLL sum rule predictions for the leading power contributions to f + BD (q 2 ), in figure 6, in a wide range 300 MeV ≤ λ B (µ 0 ) ≤ 700 MeV. We can readily observe that the sum rule predictions for f + BD (q 2 ) increase steadily with the reduction of λ B , in analogy to the observation for the radiative leptonic B → γ ν decay in [20].
As already explained in [10], the light-cone OPE for the vacuum-to-B-meson correlation function (1) can be only justified near the maximal recoil to fulfill the power counting rules n · p m c Λ and m c ∼ O( √ n · p Λ). In addition, QCD factorization for the correlation function (1) is fully applicable at space-like momentum transfer on the basis of the power counting analysis. It is then distinctly that the B-meson LCSR for the form factors of B → D ν derived in (70) can be trusted at q 2 min ≤ q 2 ≤ q 2 max = 2 GeV 2 , where a moderate value q 2 min = −3 GeV 2 will be employed (see also [61]) in the following analysis for the sake of adopting the same intervals of the sum rules parameters shown in (76). In order to access the information of B → D form factors in the whole kinematic region, we need to extrapolate the LCSR predictions obtained above toward large momentum transfer with a certain parametrization for the form factors. Following the arguments of [32,61], we will take advantage of the z-series parametrization, in line with the analytical properties and perturbative QCD scaling behaviours of B → D form factors, which can be readily introduced by mapping the cut complex q 2 -plane onto the unit disk |z(q 2 , t 0 ) ≤ 1| according to the conformal transformation [64] Here, the parameter t + = (m B + m D ) 2 is determined by the onset of the branching cut in the |B D channel, and t 0 < t + is an auxiliary parameter determining the value of q 2 mapped to the origin in the z-plane. An optimal choice of t 0 = (m B − m D ) 2 will be implemented in the numerical computation to achieve a narrow interval of |z|.
Taking into account the near-threshold behaviour from angular momentum conservation leads to the suggested parametrization [32,61] for the vector form factor, where m B * c = (6.330 ± 0.009) GeV [65], and the z-series expansion will be truncated at N = 2 in the practical matching procedure. Along the same vein, we will employ the following parametrization Figure 7: The transfer momentum dependence of the form factor f + BD (q 2 ) (left panel) and of the form factor f 0 BD (q 2 ) (right panel) predicted from the LCSR calculations, including the twoparticle contributions at (partial) NLL accuracy and the tree-level three-particle contributions, with an extrapolation toward large q 2 applying the z-series parametrization. The pink, blue and green curves correspond to theory predictions from this work, from the lattice QCD calculations by the HPQCD Collaboration [5], and from a combined fit of the BarBar and Belle data as well as the HPQCD and FNAL/MILC calculations [68]. Theory uncertainties for all the calculations are indicated by the shaded regions.
for the scalar form factor, where m B (0) c = (6.420 ± 0.009) GeV [5], and we will only keep the terms up to O(z) in the z-expansion. An alternative parametrization of B → D form factors proposed in [66] (see also [67] for a recent discussion) including the outer function and the Blaschke factor will not be considered here following the reasonings of [32].
It is now a straightforward task to implement the matching procedure described above with the aid of the LCSR predictions at q 2 min ≤ q 2 ≤ q 2 max and the z-series parametrization for the determination of the momentum-transfer dependence in the entire kinematic region. To achieve a better accuracy for the resulting shape parameters b 1 andb 1 , we will further employ the synthetic data points for the form factors f +,0 BD (q 2 ) at q 2 = 8.47 GeV 2 , 10.05 GeV 2 and 11.63 GeV 2 from the FNAL/MILC Collaboration [4] in the numerical fitting. We present the yielding predictions for the q 2 dependence of f +,0 BD (q 2 ) in the physical kinematic range 0 ≤ q 2 ≤ (m B − m D ) 2 with theory uncertainties in figure 7 where recent determinations from the lattice QCD simulation combined with a similar z-parametrization by the HPQCD Collaboration [5] and from a joint fit of the BaBar and Belle data combined with the lattice calculations [68] are also shown for a comparison. It is apparent that our predictions for the B → D ν form factors are in good agreement with those displayed in [5,68] within the theory uncertainties. We further collect the obtained results for the shape parameters, with numerically important uncertainties, entering the z-series expansion in Table 1, where the predictions for form factors at q 2 = 0 are also presented for completeness. It turns out that the dominant theory uncertainties for the shape parameters arise from the variations of the hard scales µ h1 (2) and from the errors in the determination of ω 0 . The relatively more precise predictions for the form factors of B → D ν at high q 2 are mainly due to the high precision lattice date points from the FNAL/MILC Collaboration [4], whose accuracy is even significantly better than that computed by the HPQCD Collaboration [5]. Moreover, the LCSR predictions for f +,0 BD (q 2 ) appear to grow faster with the increase of the momentum transfer squared compared to those obtained from the lattice calculations [5].  Table 1: Summary of the predicted shape parameters and the normalization constant entering the z-series parametrizations (79) and (80) for the B → D ν form factors with dominant uncertainties from the variations of theory inputs.
Another popular parametrization for the form factors of B → D ν taking into account the constraints from the heavy quark symmetry, proposed in [33], was also extensively employed in the phenomenological applications. Including the subleading power corrections to the heavyquark relations and implementing the dispersive analysis yield the following one-parameter representations for the vector and scalar form factors [33,67] where we have introduced the following conventions Apparently, this z parameter is equivalent to z(q 2 , t 0 ), defined by (78), It needs to point out that the ratio f 0 BD (z)/f + BD (z) is fully determined in HQET including both perturbative corrections to the leading Wilson coefficients and subleading power contributions computed from QCD sum rules. Matching the LCSR calculations for the form factors of B → D ν at q 2 min ≤ q 2 ≤ q 2 max onto the CLN parametrization (82) leads to f + BD (z = 0) = 1.22 ± 0.02 , ρ = 1.07 +0.08 −0.11 , Figure 8: The differential decay rates for B → Dµν µ (pink band) and B → Dτ ν τ (blue band) as a function of the momentum transfer squared q 2 , where the experimental data points for B → Dµν µ (purple squares) from the Belle Collaboration [71] are also presented for a comparison.
which are well consistent with the fitted values obtained in [5], albeit with the comparably large theory uncertainties for the slop parameter ρ.

Phenomenological implications
Having in our hands the theory predictions of the two form factors f +,0 BD (q 2 ), we are ready to explore their phenomenological implications on the semileptonic B → D ν decays. The differential decay rate of B → D ν in the rest frame of the B-meson can be computed as 0.02 ± 0.01 Table 2: Theory predictions for the partial decay rates of B → Dµν µ compared with the Belle measurements from [71].
where | p D | = λ (m 2 B , m 2 D , q 2 )/(2 m B ) with λ(a, b, c) = a 2 + b 2 + c 2 − 2ab − 2ac − 2bc is the magnitude of the three-momentum of the D-meson and originates from the short-distance QED corrections to the four-fermion operator responsible for the B → D ν decays [69,70]. The differential q 2 distributions for B → D ν obtained with the form factors f +,0 BD (q 2 ) displayed in Table 1 are plotted in figure 8, including also the recent experimental measurements for the combination of B + →D 0 e + ν e , B 0 → D − e + ν e , B + →D 0 µ + ν µ and B 0 → D − µ + ν µ from the Belle Collaboration [71]. It is evident that our predictions for the q 2 shape of B → Dµν µ are in excellent agreement with the experimental data bins. Furthermore, we collect the numerical results for the (normalized) partial decay rates of B → D ν in Tables 2 and 3 with selections of the q 2 bins identical to that from the Belle and BaBar measurements [71,72].  Table 3: Theory predictions for the partial decay rates of B → Dτ ν τ and for the binned distributions of ∆R(t 1 , t 2 ) defined in (90). The semileptonic B → D form factors obtained by fitting the experimental data with the Boyd-Grinstein-Lebed parametrization [66] are employed for the recent calculations presented in [74].
In particular, our predictions for the total decay width of B → Dµν µ in units of 1/|V cb | 2 can be obtained straightforwardly ∆Γ µ (0, 11.63 GeV 2 ) = 6.06 +1. 18 −0.92 × 10 −12 GeV , with all separate uncertainties from variations of the theory inputs added in quadrature, from which the exclusive determinations of |V cb | are achieved  [4], HPQCD [5] and from a joint fit [67] of the available experimental data and the lattice calculations including the updated unitarity bounds. Finally, we turn to compute the differential distributions of the celebrated ratio where most of the hadronic uncertainties from the B → D ν form factors are cancelled out. Inspecting the obtained results in Table 3 indeed implies an incredibly precise one-percent accuracy of ∆R(t 1 , t 2 ), albeit with the implementation of much less accurate predictions for the hadronic form factors shown in Table 1. Our predictions for the binned distributions of ∆R(t 1 , t 2 ) are also compatible with the recent determinations reported in [74], employing the B → D ν form factors extracted from a joint fit of the experimental data and two recent lattice calculations [67]. We further present our predictions for the ratio of the total branching fractions of two semileptonic decay channels which coincides with the previous determinations [4,5,67,[74][75][76][77] in the Standard Model (SM) at the 1 σ level and needs to be compared with the HFAG average value R(D) HFAG = 0.403 ± 0.040 ± 0.024 [2]. We mention in passing that the relatively high theory uncertainty of R(D) in (91), approximately 8%, can be traced back to the uncancelled hadronic uncertainties for determining the partial branching fraction of B → Dµν µ in the phase-space region 0 ≤ q 2 ≤ m 2 τ .

Concluding discussion
In this paper we have presented perturbative QCD corrections to the semileptonic B → D ν form factors with the power counting scheme m c ∼ O √ Λ m b , at leading power in Λ/m b , employing the LCSR with the two-particle B-meson DA. QCD factorization for the vacuumto-B-meson correlation function (1) was demonstrated explicitly at one loop applying the diagrammatic factorization approach. Due to the appearance of a new hard-collinear scale m c , the resulting jet functions turn out to be more complex than the counterparts in the evaluation of the correlation function for constructing the sum rules of B → π for factors. Taking advantage of the evolution equation of the B-meson DA φ − B (ω, µ), factorization-scale independence of the correlation function (1) was verified at O(α s ) with the obtained hard and hard-collinear functions. The (partial) NLL resummation improved sum rules for B → D form factors (60) derived with the dispersion representations in Appendix B constitute the main new ingredients of this paper. The subleading power contributions from the three-particle quark-gluon B-meson DA were also computed from the same LCSR method at tree level. In the light of the canonical behavious of the three-particle DA of the B-meson from the QCD sum rule analysis [7], the power suppressed three-particle corrections were demonstrated to invalidate the large-recoil symmetry relation between the vector and scalar B → D form factors.
We proceeded to explore the phenomenological implications of the resulting sum rules for the B → D ν form factors applying two nonperturbative models for the B-meson DA φ + B (ω, µ 0 ) inspired from the QCD sum rule calculations [34,55]. The perturbative QCD corrections from the two-particle DA were found to generate an approximately O(10%) shift to the LL predictions for f + BD (q 2 ) with the default theory inputs, and the one-loop hardcollinear corrections appear to have a more profound influence at q 2 ≥ 0 numerically when compared with the corresponding hard corrections. Moreover, the subleading power effects from the three-particle B-meson DA were shown to be insignificant numerically. The z-series expansion fulfilling the analytical properties of B → D form factors was further employed to extrapolate the (partial) NLL LCSR predictions toward the low recoil region. In addition, we presented theory predictions for the binned distributions of the B → D ν decay rates and of the ratio ∆R(t 1 , t 2 ) by applying our determinations of the form factors f +,0 BD (q 2 ). Matching the predicted results for the normalized decay width of B → Dµν µ and the recent experimental measurements from Belle [71] and BaBar [73] led to the extracted values of |V cb | at O(10%) accuracy as displayed in (89). Our predictions for the B → D form factors also gave rise to the determinations of R(D) presented in (91), confronted with the HFAG average value in [2].
Further developments of the B-meson LCSR approach for computing the form factors of B → D ν can be pushed forward in distinct directions. First, perturbative QCD corrections to the LCSR (66) from the three-particle B-meson DA can be carried out for a complete understanding of the leading power contributions in the heavy quark limit. In doing so, the one-loop evolution equations for the remaining high-twist three-particle DA of the B-meson in (64) are in demand and they have been recently worked out in [78]. Second, computing yet higher-order QCD corrections to the two-particle contributions (60) are of both conceptual and phenomenological interest for exploring the renormalization properties of the B-meson DA φ ± B (ω, µ) at two loops (e.g., the eigenfunctions and the analytical structures of renor-malization kernels) and for bringing down the still sizeable perturbative uncertainties of the theory predictions displayed in Table 1. Third, the present strategies can be readily applied to calculate QCD corrections to the semileptonic B → D * ν form factors based upon the LCSR with the B-meson DA (with additional attention to the renormalization prescription of γ 5 in dimensional regularization), allowing for a comprehensive analysis of the full angular distributions B → D * (→ D X) τ (→ Y ν τ )ν τ with X = (π, γ) and Y = ( ν, π) as discussed in [79]. To conclude, precision QCD calculations of semileptonic B-meson form factors with analytical QCD approaches will continually provide us with a deeper insight into the strong interaction dynamics of heavy quark decays and into the general properties of effective field theories.
Here, D = 4 − 2 , the integration measure is defined as follows and we also introduce the conventions . (109)

B Spectral representations
We collect the dispersion representations of various convolution integrals entering the (partial) NLL resummation improved factorization formula (54), for the sake of constructing the sum rules of B → D form factors given by (60). It needs to point out that we have validated each spectral function by verifying the corresponding dispersion integral numerically.