Standard Model prediction of the $B_c$ lifetime

Applying an operator product expansion approach we update the Standard Model prediction of the $B_c$ lifetime from over 20 years ago. The non-perturbative velocity expansion is carried out up to third order in the relative velocity of the heavy quarks. The scheme dependence is studied using three different mass schemes for the $\bar b$ and $c$ quarks, resulting in three different values consistent with each other and with experiment. Special focus has been laid on renormalon cancellation in the computation. Uncertainties resulting from scale dependence, neglecting the strange quark mass, non-perturbative matrix elements and parametric uncertainties are discussed in detail. The resulting uncertainties are still rather large compared to the experimental ones, and therefore do not allow for clear-cut conclusions concerning New Physics effects in the $B_c$ decay.


Introduction
The B c meson is the lightest state with both "naked" beauty and charm. As such it is stable against both strong and electromagnetic decay. Its weak decay can proceed through three distinct mechanisms: eitherb-quark decay, c-quark decay orb-c annihilation. Experimentally, the B c lifetime has been measured by LHCb [1,2] and CMS [3] with a world average of [4]: τ exp Bc = 0.510 (9) ps , which corresponds to a decay width Since both valence quarks in the B c meson are heavy, the state is similar in structure to it's quarkonium cousins, the η c and η b pseudoscalar mesons, the lightest members of the J/ψ and Υ towers of states. This circumstance allows for an effective treatment in terms of Non-Relativistic QCD (NRQCD), where (much as in the case of heavy quarkonium) the anti-quark corresponding to the valence quark, or the quark corresponding to the valence anti-quark, is integrated out at the respective scale. This method has been used to estimate the lifetime of the B c meson in the Standard Model (SM) [5][6][7].
The most precise value has been obtained in [6] to be τ Bc = 0.52 ps, with an uncertainty from varying the input charm quark 1 mass that results in 0.4 ps < τ Bc < 0.7 ps, corresponding to Γ Bc = 1.92 +0.58 −0.49 ps −1 , exclusive of other sources of uncertainty that include, among others, ±22% from scale uncertainty in the perturbative calculation. Similar results were found in other OPE calculations [5,7] as well as using QCD sum rules [8] or potential models [9]. A comparison of the different predictions can be found in [10].
The experimental measurement in Eq. (1) has a much smaller uncertainty than the theory prediction in Eq. (3). This motivates a reinvestigation of the SM prediction from over 20 years ago with the goal of improving the theoretical precision, and the eventual hope of a precision comparable to the experimental one.
Renewed interest in the B c lifetime has arisen because it is susceptible to New Physics (NP) effects. Consequently a more precise SM prediction allows to place stronger constrains on NP models. In particular, experimentally measured deviations from SM expectations in the semileptonic decays B → Dτ ν, B → D * τ ν and B c → ψτ ν suggest NP contributions to the quark level process b → cτ ν [11][12][13][14][15][16][17]. These so-called R(D), R(D * ) and R(B c ) anomalies can be accounted for by several extensions of the SM. If the NP is realized as an effective pseudoscalar interaction, that is, a four-fermion interaction involving a pseudoscalar hadronic bilinear times a leptonic S −P bilinear, then the B c lifetime is especially effective in placing constraints on its strength [18,19]. This type of interaction is often found in models of NP proposed in the interpretation of the R(D ( * ) ) anomalies, such as the two-Higgs-doublet model (2HDM) and leptoquarks. The B c lifetime constraint rules out any of the 2HDM interpretations of R(D * ), including the type III versions of the model that contain general Yukawa couplings to the different fermions [20], which can explain simultaneously R(D) and R(D * ), [21][22][23][24][25][26], going beyond the more restricted set of models [27][28][29][30] on which BaBar [11,12] and Belle [14] place constraints. Leptoquark models have also been proposed to explain these anomalies and they are susceptible to these bounds when they generate sizable pseudo-scalar operators [18,31].
In this work we follow the OPE methodology of Ref. [6] (henceforth "BB"). There are several ways in which we can improve the result of that work. First and foremost, we use and compare three mass schemes to eliminate pole (or "on-shell") masses in favor of well defined masses and therefore eliminate renormalon ambiguities that arise in the on-shell scheme. The largest source of uncertainty in the calculation of the width is from the pole masses, since the width scales as the fifth power of these. BB use ad hoc values for b and c pole masses. We perform the calculation in the MS mass scheme, the "Upsilon" scheme [32,33], and a "meson" scheme that incorporates some aspects of the Upsilon scheme. We will expand on these below but for now we point out that the Upsilon scheme organizes the perturbative expansion in a manner that not only eliminates renormalon ambiguities from the expression for the decay width but in addition is empirically seen to have better convergence than other schemes. Moreover, the masses of the Υ and J/ψ particles, used as inputs, are determined accurately. Second, for the one-loop calculation of the subprocess b → ccs BB uses the result of Bagan et al [34,35]. Later, Krinner, Lenz and Rauh (henceforth KLR) inferred the presence of typos in the analytic expressions of Bagan et al [36], because those expressions do not produce the numerical results and graphs they present. KLR compute the 1-loop corrections anew and agree with the numerical results of Bagan et al, inferring thus the presence of typos. We have used, and verified partially, the results of KLR. Third, BB neglect, and we include, the contribution of penguin operators, that are formally of the same order as other contributions retained in the calculation albeit numerically small. Fourth, we use better determined input parameters, such as the strong coupling constant α s (M Z ), the CKM matrix elements, and the non-perturbative B c decay constant, f Bc , for which there is now a lattice calculation [37], among others. Fifth, we use spin symmetry to relate some of the non-perturbative matrix elements appearing in the calculation.
In addition to updating and improving on the BB result, we have tried to clarify several issues in their presentation. Among these we clarify below the need to go up to fourth order in velocity in the NREFT, the role of spin symmetry in the calculation of matrix elements, the relative (un)importance of various corrections, the distinction between quarks and anti-quarks in the NREFT as well as the interpretation and precise value of the momenta p b ± p c that enter the operators in the OPE for weak annihilation (WA) and Pauli interference (PI) diagrams.
The rest of the paper is organized as follows: In Sec. 2 the different mass schemes used for the calculation are discussed. In Sec. 3 we outline the Effective Field Theory approach used at the electroweak (EW) scale involving the effective Hamiltonian. The matching onto NRQCD is performed in Sec. 4. The Operator Product Expansion (OPE) is discussed in Sec. 5. Section 6 describes the computation of relevant matrix elements and in Sec. 7 the numerical computation and analysis of the uncertainties for the SM prediction is performed before we summarize our results and offer some observations in 4 the conclusions, in Sec. 8.

Mass schemes
The largest source of uncertainty in the calculation of the inclusive rate of the B c decay is in the value of the pole mass. As we will review below, the calculation of the rate, Γ Bc , is based on an OPE of the two point function of the effective Hamiltonian, cf Eqs. (33) and (34). The leading term in the calculation corresponds to Γ b + Γ c , the sum of decay rates of b and c quarks computed perturbatively as if the quarks were not bound to the B c meson. The rate Γ q is given in terms of the pole mass, m q , as an expansion in powers of α s where κ is a constant independent of m q , x and α s chosen so that the tree-level rate has f 0 (0) = 1. For reasons that will be explained momentarily, we have introduced the power counting parameter = 1 that is set to unity at the end of the calculation. The coefficients of the expansion, f i (x) are functions of the ratios of final state pole masses, m q , to the mass of the decaying quark, x = m q /m q . Pole masses are convenient for the perturbative calculations, but are beset by both computational and conceptual difficulties: their perturbative expansion is poorly convergent and suffers from a renormalon ambiguity. Moreover, the roots of the implicit equation that relates them to short distance (Lagrangian) masses are complex for the light quarks. When expressed in terms of the pole mass, the perturbative rate also suffers from a renormalon ambiguity. Remarkably, eliminating the pole mass in favour of well defined (e.g., short distance) masses, gives a perturbative expansion of the rate that is free of renormalon ambiguities [38][39][40][41]. For each choice of well defined mass the manner in which renormalon ambiguities cancel in the rate suggests how to organize the perturbative expansion. We refer to any one choice of mass and of reorganization of the perturbative series as a "mass-scheme". When the well defined quark masses are chosen as the modified minimally subtracted masses, m q , the expansion of the pole mass takes the form so that the renormalon free expansion for the rate for the flavour transition q → q takes the form wherex = m q /m q . As we see, in this "MS scheme" the expansion in is equivalent to the perturbative expansion in powers of α s . By contrast, in the "Upsilon scheme" the power counting in does not correspond to powers of α s . Up to small non-perturbative corrections the q pole mass is given in terms of that of the mass Mq q of some quarkoniumqq state by For example, if in Γ b→cud we neglect the light quark masses and use this scheme for both b and c quarks, then where now X = Mc c /Mb b . While any good choice of well defined masses yields a well behaved perturbative expansion in the sense that it is free of renormalon ambiguities, different mass-schemes may differ in how rapidly the expansion converges (in the asymptotic expansion sense). Were we able to compute to high order, all mass-schemes would give the same numerical value for the rate (up to a small higher order term). But computations of rates are available only to low orders in perturbation theory, often only including 1-loop corrections to the leading, tree level term. In practical term it is best to choose a scheme that converges most rapidly. In the sub-sections below we present and compare the three schemes in which we perform calculations.
As stated above, the leading term in the OPE for the decay width of B c mesons is the perturbative Γ b + Γ c . The corrections to those are expressed as products of nonperturbative matrix elements and Wilson coefficients. The latter are perturbatively computed as functions of pole masses. In the calculations below we use the same scheme choice for these sub-leading terms as for the leading ones.

The MS mass-scheme
This scheme uses the running renormalized Lagrangian masses m b (µ) and m c (µ) evaluated at a sufficiently high renormalization scale µ. In this scheme the expansion parameter simply counts powers of α s , making its use particularly simple. The expansion of m b,c in terms of the MS masses is known to third order in α s [42,43], but we only need it to first order: where we have retained explicitly the power of used in organizing the perturbative expansion, corresponding to Eq. (5). It is convenient to cast the c-quark decay rate in terms of the running mass evaluated at itself, m c (m c ), and the b-quark decay rate in terms of both quark masses evaluated at the shorter distance scale, m b (m b ) and m c (m b ). The value of the masses, m b (m b ) and m c (m c ), is reported by the PDG with 1%-2% accuracy. Because the rates scale as the fifth power of the mass, this results in an uncertainty of 5%-10%. In addition, as discussed below, the convergence of the perturbative series is slow in the case of semileptonic b and c decays, and there is no reason to suspect it is any better for non-leptonic decays. 6

The Upsilon scheme
The Upsilon expansion was introduced by Hoang, Ligeti and Manohar in Refs. [32,33] (henceforth "HLM") to address the largest source of uncertainty in the calculation of the inclusive rate of semileptonic B and D meson decays. The calculation of the rate for semileptonic B decays is based on an operator product expansion of the two-point function of hadronic charged currents in terms of operators in the Heavy Quark Effective Theory (HQET) [44]. The expansion is naturally given in inverse powers of the pole mass of the decaying heavy quark, and the resulting rate is proportional to the fifth power of the pole mass. In the Upsilon scheme the masses are chosen to be the well measured 1S masses, 2 Mb b = m Υ and Mc c = m J/ψ , resulting in a negligible uncertainty in the decay rate from the value of these masses. The expansion of the pole masses in terms of m Υ and m J/ψ has been determined to fourth order in α s , with corrections that start at order α 2 s [45,46]: where β 0 = 11 − 2 3 n f is the coefficient of the first term in the QCD β-function, C F = (N 2 c − 1)/2N c = 4/3, and we have omitted the known term of order α 4 s . For quarkonium systems the quark potential suffers from a renormalon ambiguity, as does the pole mass, but the resulting quarkonium state energy is free from ambiguities and hence the onium mass is a good candidate for an unambiguous, well defined mass [47,48].
The organization of the expansion is unusual. The parameter = 1 is inserted to indicate the order in the expansion. For the expansion of the rate in Eq. (4) the power of matches that of α s . But, as seen above, the parameter in Eq. (7) carries one additional power of α s , i.e., n α n+1 s . As seen in Eq. (8) the leading correction to the rate includes order α s terms from the perturbative corrections to the rate, and order α 2 s from the perturbative expansion of the masses. This is dictated by the requirement that the renormalons cancel in the expression for the rate.
In the implementation of this scheme for non-leptonic decays HLM expand the Wilson coefficients of the electroweak effective Hamiltonian in powers of α s ln(m W /m b ) and truncate in the power of α s to which they perform the calculation of the rate. We do not adopt this prescription and instead retain the full resummed values of the Wilson coefficients. The reason is that this series does not contribute to the cancellation of renormalons, as can be seen from the fact that M W is an arbitrary parameter. This is the case for the other mass schemes as well.

Light quarks
The available 1-loop calculations give decay rates in terms of pole masses. At L-loops the pole mass m q is determined implicitly, as a root of the equation where the MS-mass, m q (µ), and the running coupling, α s (µ), are computed at L-loops [42]. 3 Beyond 1-loop this relation has only complex roots for small input MS-mass; that is, the relation involves necessarily ln(ln(m q /Λ QCD )) which is complex for the light quarks, q = u, d, s. The effect of including the light quark masses is expected to be small, so one may get around this difficulty by working in the massless limit; chiral symmetry guarantees, perturbatively, the vanishing of the pole mass. Yet, the strange quark is sufficiently heavy as to have a significant effect on charm decays. Therefore, in principle Eq. (14) has to be solved iteratively in order to give m q in terms of m q (µ) at a sufficiently high scale µ, where perturbation theory is still valid. In this work however we only use the 1-loop result given in Eq. (9) together with the results from Monte Carlo simulations of QCD on the lattice that reliably determine m q (µ) at µ = 2 GeV. For the expansion in we use the power counting indicated in Eq. (9), regardless of the scheme used for the heavy quarks. The numerical impact of the non-zero strange quark mass will be discussed below. We use a non-zero strange mass for charm decays and neglect m s inb-decays as well as in WA and PI. Clearly the non-vanishing mass tends to decrease rates as it restricts, if only marginally, the phase space available to the decay products. But this effect is very small in decays of the heavier b quark and even smaller in WA. Because we do not consistently include a non-zero strange quark mass in all B c decay channels, we give our results using a vanishing strange quark mass in all B c decay channels, and then, in addition, we list the results of decay channels from charm decay anew with non-zero strange quark mass.

Nomenclature
For clarity we summarize our definitions of "schemes" as applied to the computation of the B c lifetime in the rest of this work: MS For all partial rates we use m b and m c given in terms of m b (µ b ) and m c (µ c ) respectively.
Meson For all partial rates we use Eq. (11) to give m c in terms of m b and m B − m D , and we use the Upsilon expansion to give m b in terms of m Υ .
Upsilon The contributions to the total width that arise fromb decays, WA and PI are computed as in the meson scheme, but for those from c decays the c-quark pole mass is given by the Upsilon expansion in terms of m J/ψ .
The contributions to the total width arising fromb decays, WA and PI are computed with , while for c decays we use µ c = m c (m c ). These choices are motivated by the typical total energy released in each of these processes. In all three schemes the light quarks, q = u, d, s are assumed to be massless, except for the strange quark for which we use m s in terms of m s (µ c ) when computing partial B c decay widths from c-quark decays.

Comparison of schemes
We have already indicated that both the Upsilon and the meson schemes have the advantage that the input masses are extremely well known. This gives them a clear advantage over the MS scheme. However, one should keep in mind that the relation giving the pole mass in terms of the 1S mass is subject to poorly known non-perturbative corrections. The authors of HLM estimate the non-perturbative mass shift δm as δm ∼ a 3 Λ 4 QCD , where a is the Bohr radius of the 1S. For the Upsilon state they estimate this correction as δm ∼ 15 MeV for Λ QCD = 350 MeV (δm ∼ 60 MeV for Λ QCD = 500 MeV). The leading contribution to the mass shift is from the gluon condensate α s G 2 µν , and has been calculated in Refs. [51,52]: This is estimated as δm ≈ 60 MeV in Ref. [45]. HLM uses δm ∼ 100 MeV as a conservative estimate of the error in the pole mass relation. This seems to result from an abundance of caution. We do not include the condensate correction in our calculations. Were we to include it in the calculations in the Upsilon scheme the uncertainty in the charm quark mass would render it useless. Instead we assume it can be neglected. This can be checked a posteriori, and much like the implicit assumption that quark-hadron duality can be used for the implementation of the OPE method we view this as an assumption that can be validated either by further theoretical progress or by experimental results.
Turning to the nature of the expansions, we present calculations of semileptonic decay rates of B and D mesons to compare some of the schemes. Using = 1 to indicate orders of the expansion, the rate for B → X u eν in the Upsilon scheme behaves as 4 1 − 0.115 − 0.030 2 + · · · while for the MS scheme 1 + 0.30 + 0.20 2 + · · · For the numerical estimate we have used α s = 0.223, corresponding to the running coupling evaluated at µ = m b (m b ). The ellipses stand for higher orders of as well as non-perturbative corrections. These estimates seem to indicate more rapid convergence of the Upsilon scheme than for the MS scheme. For comparison we also give the expansion for the rate in terms of the pole mass: The expansion in the Upsilon scheme seems to be approaching a limit quickly, as each term in the expansion is an order of magnitude smaller than the previous. By contrast the expansion in the MS scheme displays much slower convergence. One expects of course that after a sufficient number of terms are included, both Upsilon and MS scheme expansions will give equivalent rates. But in practice only a small number of terms in the series can be computed and it is most practical to use the scheme that converges fastest. We note that the convergence in terms of the pole mass is similar to that of the MS scheme. This should not suggest to use the pole mass as a scheme: as mentioned above, both the mass and rate expansion are ambiguous. For B → X c eν we find the following expansions 1 − 0.10 − 0.03 2 + · · · (meson/Upsilon scheme) In addition, the reader may find in HLM similar estimates for B → X u,c τ ν, and B → X c(s+d) (that only retain the BLM part of the 2 term). In all cases it is apparent that the expansion in the Upsilon scheme converges faster than in the MS scheme and faster than for the (a priori ill defined) rate in terms of pole masses.

Effective Hamiltonian
We employ the standard Effective Field Theory approach where the heavy SM particles (top quark, Higgs and EW gauge bosons) are integrated out at the EW scale and matched onto the following effective Hamiltonian [55,56]: 5 with the current-current operators and the QCD-penguin operators Here P R,L = 1 2 (1 ± γ 5 ) and α, β denote colour indices. We have used V tb V * ts ≈ −V cb V * cs , which holds to high precision. The sums run over the active quarks at the given scale.
To compute the relevant partial decays rates at one-loop, one has to combine one-loop matrix elements together with Wilson coefficients resulting from a one-loop matching and two-loop running calculation. In our approach we only work up to this order for the operators Q 1,2 but not for the QCD penguin operators. This simplification is justified by the fact that the operators Q 3,4,5,6 have much smaller Wilson coefficients. For convenience we report in Tab. 1 the LO Wilson coefficients of the operators in eq. (15) at the scales of the b-and c-quark.
The semileptonic channels also contribute to the decay rate of the B c meson. Integrating out the W at the EW scale leads to the following effective charged-current operators: where the sum runs over all lepton flavours. At tree-level in the SM with the chosen normalization the Wilson coefficients are equal to one for all lepton flavours, C = 1. These operators do not run under QCD due to current conservation and have a small running under QED [57]. Therefore, we only consider the tree-level matching at the electroweak scale and neglect RG effects for these operators.
With a slight abuse of notation, we use the same symbols to denote operators and Wilson coefficients of the effective Hamiltonian for charm hadronic decay, with the current-current operators Because these operators carry four separate flavours, there is no QCD mixing from Penguin operators. For semileptonic decays we have the analogue of (18), As described in Sec. 5 below, the effective Hamiltonians in Eqs. (15), (18), (19) and (21) are then used together with an OPE to obtain the B c lifetime.

Non-Relativistic QCD
The b-quark and c-antiquark in theB c meson can be well described in NRQCD. There are several advantages in utilizing this effective theory. First, it organizes the computation in an expansion in powers of the relative velocity v = | v| of the heavy quarks bound in the meson. Second, the expansion makes explicit additional approximate "spin" symmetries for the b and c quarks separately. These then give relations among matrix elements that hold even non-perturbatively (in the QCD coupling expansion). Third, separate conservation of b and c numbers fixes, non-perturbatively, the values of the matrix elements of the leading operators in the OPE that determines the semi-inclusive partial lifetimes B c → X cc and B c → X bb . The corrections arise from the order in perturbation theory to which the Wilson coefficients are computed, and from the higher order terms in the OPE. The matrix elements of the latter can be estimated with, say, potential models; to the extent that the expansion in v is well behaved, the potentially large uncertainty in the sub-leading matrix elements may still translate in a manageable uncertainty in the total rate. We briefly review elements of NRQCD.
In our treatment the quark NRQCD fields are still given in terms of Dirac spinors, much like is commonplace in HQET. The 4-velocity of the quarkonium state is u µ , with u 2 = 1. We take away the fast oscillation in the near on-shell evolution of the field of the QCD heavy quark, Q(x), of mass m, and furthermore project out the positive energy components Ψ + (that correspond to the particle annihilation operator), where The equation of motion, (i / D − m)Q = 0 can then be used to solve for the "small component" where for any vector a µ the spatial component in the quarkonium restframe is a µ ⊥ = a µ − (u · a)u µ . Using this in the QCD Lagrangian for the quark gives The NRQCD Lagrangian is obtained by expanding the Lagrangian in Eq. 22 in powers of iu · D/2m and truncating the expansion. The characteristic scale of the derivatives on the field Ψ + is 1 2 mv 2 for the energy, iu · D, and mv for the spatial momentum, D ⊥ . Hence, counting powers of v one has iu · D ∼ D 2 ⊥ /2m ∼ v 2 . In addition, the integral d 3 xΨ + Ψ + ∼ 1 will be concentrated over the volume ∼ (M v) −3 of a quarkonium state of mass M , and it follows that Ψ + ∼ v 3/2 . Counting rules for all fields and derivatives follow from these and the additional field equations [58]. 6 In particular, for the gauge field in Coulomb gauge, g s u · A ∼ v 2 , g s A ⊥ ∼ v 3 , which determines the velocity-scaling of the chromo-electric and magnetic fields (in the quarkonium restframe) g s E ∼ v 3 and g s B ∼ v 4 . Finally, since the typical momentum is the inverse Bohr radius, one has α s (M v) ∼ v. Expanding the Lagrangian and retaining the lowest order (∼ v 5 ) one has The NRQCD treatment of antiquark fields is analogous. Now one has where X − is an antiquark creation operator (containing only negative energy components). It follows that In order to compare these results to calculations that use the non-relativistic 2component spinor notation, and to use estimates of matrix elements that use quark potential models, we recast them in the rest frame, u = (1, 0). In the Dirac basis of γ-matrices, γ 0 = σ 3 ⊗ 1, γ = iσ 2 ⊗ σ, the matrices (1 + / u)/2 and (1 − / u)/2 project out the upper and lower two components of the 4-component spinor, which we denote as ψ q and χ q , respectively. Then, for the quark we have and for the antiquark While ψ q is only an annihilation operator, χ q is only a creation operator. It is convenient to rewrite the Lagrangian for the antiquark in terms of the annihilation spinor In terms of this the anti-quark Lagrangian is where it should be noted that the covariant derivatives involve the generators −T aT , corresponding to those of anti-triplets. The lowest order quark and antiquark Lagrangian in NRQCD is symmetric under unitary transformations of the components of the Dirac spinors Ψ + and X − that preserve the conditions These transformations form an internal symmetry group isomorphic to SU (2) × SU (2), and hence the heavy quark spin component of total angular momentum is separately conserved, precisely as expected from non-relativistic quantum mechanics. This is the additional spin-symmetry alluded to above.
The spin symmetry has considerable implications relating matrix elements of composite operators. An important example is for the matrix elements of the four-fermion operators entering our computation. For gamma-matrices Γ and Γ spin symmetry implies where Ψ (c) − stand for the Ψ + and X − fields of the c and b quarks. This relates all of the matrix elements of these 4-fermion operators to a single invariant "reduced matrix element" Q . It also gives for arbitrary Dirac matrices Γ c and Γ b . Higher order terms in the Lagrangian are readily incorporated. The first corrections to L 0 are of order v 7 (relative order v 2 ). Operators with powers of iu · D are difficult to simulate on the lattice, and therefore it is conventional (but unnecessary) to eliminate them from the higher order corrections by means of field transformations. For example, the terms are eliminated in favour of the ones without temporal derivatives by the transformation which is chosen such that temporal derivatives in L appear only in commutators with spatial derivatives, thus: In order to use the result of this expansion in inverse powers of m together with an expansion in powers of the relative velocity v, it is necessary to display explicitly the (restframe) electric and magnetic fields.
in the quarkonium restframe), this gives, indicates the derivative acting only on the field strength tensor. The last line in (28) is from the square of the commutator, the last term in (27). Therefore, at order v 7 The dimensionless coefficients c n , n = 1 . . . , 4 are functions of α s and are determined by standard EFT matching procedures. 7 The tree-level calculation above gives c n = 1 + O(α s ), for all four coefficients. Beyond tree level other operators may appear. Using dimensional analysis, and imposing symmetries, the additional possible operators can be listed. Reparametrization invariance implies that the coefficient of the second term in the Lagrangian in (23) remains unrenormalized (relative to the first term). It also imposes relations on the coefficients of higher dimensional operators. In particular c 1 = 1 and At the next order in the expansion in inverse powers of m one has from the expansion of the Lagrangian in Eq. (22), and various terms containing a single power of iu · D and four powers of i / D ⊥ that arise from the product of the lower order terms in Eq. (22) and those in the change of field variables in (26). Terms containing a single power of iu·D can be combined into commutators so as to eliminate single temporal derivatives, iu · D, at order 1/m 4 , provided one further changes field variables via On the other hand, commutators of D ⊥ with u · D alone cannot remove all temporal derivatives from the Lagrangian at order 1/m 4 . 8 The conserved Noether current associated with quark number that follows from L 0 is Then evaluating at p = p = M Bc u and using spin symmetry 9 and an analogous expression for the antiquark. This expression holds non-perturbatively in α s and receives corrections of order v 2 from the symmetry breaking terms in L 1 . This is the basis for the statement above that quark number conservation gives, nonperturbatively, the values of the matrix elements of the leading operators in the OPE that determines Γ(B c → X cc ) and Γ(B c → X bb ). Furthermore we have This has a simple physical interpretation: the right hand side is a vector, but can only depend on the 4-velocity u, which however has no perp-component (no spatial component in the B c rest-frame).

Operator product expansion
To obtain the lifetime of the B c meson, the optical theorem is used to relate the decay width to the forward scattering of the B c meson: where the transition operator is given by the absorptive part of the time-ordered product: 8 Although higher temporal derivatives are acceptable in effective field theories, as mentioned above, they are better avoided in lattice simulations of the quantum field theory. That these order (iu · D) 3 terms cannot be removed may have escaped attention because they come in at order v 11 , or relative to the lowest order, they constitute a relativistic correction of order v 6 . 9 In the QCD case, defining the quark form factor f (q 2 ) by B c (p )|j µ |B c (p) = f (q 2 )(p + p ) µ , where j µ =cγ µ c and q = p − p , charge conservation gives f (0) = 1.
Invoking quark-hadron duality one expects that for a large energy release an OPE can be performed to express the transition operator as a series of local operators of increasingly higher dimensions with coefficients suppressed by the large energy released, corresponding in the case at hand to the heavy quark masses, m q . The calculation of the transition operator is organized by separating the contributions for b and c decays, and those of WA and PI terms, Each of these is an expansion of operators, with T b,c starting at dimension 3 while T WA,PI starting at dimension 6. In NRQCD, these correspond to expansions starting at order v 3 and v 6 , respectively. Retaining the contributions T WA,PI is physically important. It would appear then that for consistency we need a full calculation to order v 6 which, since α s ∼ v, should include corrections of order α 3 s to the coefficients of the dimension-3 operators. However, for WA and PI the 2-body final state is enhanced relative to the 3-body phase space of single quark decay by 16π 2 which is numerically ∼ α s (m b ) −3 , and suppressed by the probability that the quarks in the B c meet at a point, controlled by the wave-function at the origin, |Ψ(0)| 2 ∝ f 2 Bc , where f Bc is the B c -meson decay constant. Therefore we consider the expansions of T b,c and T WA,PI independently and carry out each to some fixed order.
It must be noted that, since we use NRQCD to organize the calculation in powers of the relative velocity u, the OPEs are reorganized. For example, the 4-fermion operator X − and the magnetic moment operator Ψ (c) + are of mass dimension 6 and 5 respectively but of order v 6 and v 7 , respectively. We follow Ref. [6] in including some terms of order v 7 in the expansion, merely to explore their significance. We do, in contrast, retain all operators of order v 6 in the expansion of T WA,PI . In particular, we retain the operator Ψ (c) + which we count as of order v 6 since it is equivalent to a combination of 4-fermion penguin operators Q 3 -Q 6 which we keep in our calculations. The expansion of T b,c is done consistently to order v 4 . The first non-perturbative effect in the calculation of T b,c comes in at order v 5 and we retain this together with the next non-perturbative effect, of order v 7 , to explore the relevance of non-perturbative corrections. These corrections are particularly important for charm decays, namely roughly of order 20%. By comparison the omitted perturbative corrections are expected to be of order α 2 s (m c ) ∼ (0.37) 2 . As inferred from the discussion above, we choose to expand in terms of operators in NRQCD. This has several advantages, if only conceptually. First, it organizes the expansion systematically. To see this consider that in the QCD expansion, an operator with (∂/m q ) n q, where q is a heavy quark, is not suppressed even though it carries arbitrary powers of the large mass in the denominator. In NRQCD i∂/m q is replaced by v + i∂/m q the rest-frame velocity plus a derivative that corresponds to the residual momentum, a small quantity by the NRQCD power counting. Second, as we will see, the QCD expansion contains non-local operators, but not so the NRQCD expansion. This is often ignored by writing Wilson coefficients as functions of the momenta, p b,c , of the b, c quarks in the corresponding operators, but one should keep in mind that these just stand for derivatives acting on the fields in the operators. The Wilson coefficients often contain negative powers of these momenta, that is, they are non-local operators. In NRQCD, p −2 = (m q u + k) −2 has a local expansion in terms of the derivative k → i∂. The third advantage is that spin symmetry allows for vast simplifications which are absent in QCD. And finally, if using NRQCD there is no non-trivial perturbative matching correction of the leading OPE operator, as we now explain.
The perturbative matching calculation involves on-shell, or near on-shell external quarks. The leading terms in the OPE consist of dimension 3 operators of the form Ψ + ΓΨ + (or the analogous anti-quark operators) with Γ a Dirac matrix. The matching calculation consists of evaluating the matrix element in single quark states | p.s of the transition operator in (34) on the one hand, and the same matrix element of the OPE on the other, and then fixing the Wilson coefficients by imposing equality between the two calculations: In order to compute C (3) i at n-th order in perturbation theory one must compute the matrix elements on the right-and left-hand sides of this equation. However in NRQCD (but not in QCD) the dimension 3 operators on the right-hand side are all protected from radiative corrections at zero momentum transfer because they are related by spin symmetry to the conserved current (31) and just as in the case of mesons, the quark form factor has F (0) = 1. 10 Incidentally, spin-symmetry allows us to treat the dimension 3 operators as a single one, and this is standard practice: by choosing an operator in Eq. (32) that has a unit matrix element the imaginary part of the Wilson coefficient takes the value of the perturbative quark decay width.
With this understanding we can now write The computation of the C (i) q as well as that of T WA,PI will be discussed in the following subsections.

C
c : Free c-quark decay We turn first to the leading order term in the expansion in (35). C Writing Q(x) in terms of Ψ + (x) and performing the change of field variables in (26) and (30), and retaining only terms up to order u 7 , one obtains This can be used for T c and an analogous expansion in terms of X − for T b . When we need to evaluate matrix elements of the right hand side of Eq. (36) we use the equations of motion to eliminate iu · D in favour of −(iD ⊥ ) 2 /2m, we go to the restframe of the meson to decompose the field strength into chromo-electric and magnetic components and use 2-component spinors: Below we neglect the spin-orbit coupling ψ † q g s σ·( E × D− D× E)ψ q , since it has a vanishing matrix element.
The coefficient C c is given by the decay rate of the charm quark as if it were unbound, Throughout we approximate the light quarks, u, d as massless; as explained in Sec. 2.4 we compute both for m s = 0 and with a non-vanishing s-quark pole mass given in terms of the running MS mass. The partial decay rate for c → sud at m s = 0 and in the absence of running of Wilson coefficients was first calculated by Guberina, Peccei and Rückl [62]. Altarelli, Curci, Martinelli and Petrarca first computed 2-loop running and computed the 1-loop corrections to the decay rate anew using a subtraction scheme common to both calculations [63], and a missing ln(µ/m) term was corrected by Altarelli and Petrarca [64]. Their result was confirmed by Buchalla [65]. Hokim and Pham computed the α s corrections to the decay rate for arbitrary final quark masses, including only the effect of Q 2 [66]; their results are trivially adapted to compute the α s corrections to the decay rate if the only operator were Q 1 . The contribution to the rate from Q 1 -Q 2 interference, for one massive and two massless quarks in the final state, as is the case in, e.g., c → sud and in b → cūd, was given by Bagan, Ball, Braun and Gosdzinsky [34]. Cabibbo and Maiani computed the α s corrections to the semileptonic decay neglecting the charged lepton mass; the dependence on the final state quark mass was computed numerically [67]. Finally, Nir gives Γb →c ν , or equivalently, Γc →s ν , for massless leptons [68]. The analytic expression for the 1-loop correction to the semileptonic rate can be inferred from the work of Hokim and Pham cited above. We note in passing that, although not stated explicitly, it may be inferred that these works provide the rates in terms of quark pole (on-shell) masses [68]. We assume this is also the case for the non-leptonic decay rates. As explained above, given a mass scheme one has a perturbative expansion of the pole mass in terms of a well defined mass. Consistency requires that when expressed in terms of the well defined mass, the rate can be expanded (and truncated) to the appropriate order, i.e., α n s with n = 2 in the Upsilon and meson schemes, and n = 1 in the MS scheme. We compute the quark decay rates at the one-loop order for the semi-leptonic decays and also for the contributions of the operators Q 1,2 in H ef f in the case of hadronic decays. For the penguin operators, Q 3−6 , however we only take compute the rate to leading order.
Using the decay rate Γ c to compute the coefficient C ofQQ in the OPE is only appropriate for theΨ + Ψ + term in the expansion ofQQ in terms of NRQCD fields, that is, the first term on the right side of Eq. (36). The coefficients of each of the remaining terms should be computed from individually matching them: in the OPE one should really write where the ellipsis stand both for other terms in the matching ofQQ and for the higher order terms in Eq. (35). At lowest order C q , with n = 1, 2, . . . for operators that are present in the expansion (36) and vanish for operators that arise in theQQ expansion only at or above 1-loop order. Reparametrization invariance requires that some of these operators come in fixed combinations, and this is reflected in exact relations between some of these coefficients, e.g., C In our calculations we have used the leading order expression for Γ q (q = c,b) for the coefficients C (3,n) q , and the next to leading expression for Γ q for C (3) q . Since the subleading operators in (36) are of order v 2 and higher, this truncation introduces uncertainties of order v 2 α s ∼ v 3 .
Throughout we approximate the light quarks, u, d and s as massless; the effect of a non-vanishing strange quark mass will also be estimated. The perturbative calculation to order α s of individual (anti)-b quark decay rates can be found in the literature as detailed above, including the full m c -dependence. In addition, Ref. [35] gives Γb →cτ ν including both m c and m τ -dependence. For the transitionb →cud the expressions for the decay rate are given in [34], who also estimated the case of a single massless and a massive pair of quarks in the final state, which is the case for b → ccs. KLR corrected a misprint in the latter 11 and included additional effects, most importantly that of penguin operators and of writing the rate in terms of a well defined mass. As in KLR, we treat the penguin Wilson coefficients as of sub-leading order; while this is not formally correct, it is practical, since the full calculation of the decay rate at 1-loop is unavailable. We also follow KLR in including also Q 1,2 contributions to the rate through "penguin diagrams", and the contribution of the chromomagnetic operator,bσ µν G A µν T A b. 11 The full decay rate has first been reported in [70]. As pointed out by KLR, there are several misprints in [70], resulting in a much lower decay rate for b →ccs. Also, KLR reports complex values for the contributions of vertex corrections to the decay rate. We have determined that the contribution to the rate is from the real part of those results (i.e., not twice the real part, nor the imaginary part). We have verified this by an explicit recalculation of some of the vertex graphs and the corresponding real emission graphs. In addition we have verified that with this interpretation the resulting width has a finite limit as the gluon (IR regulator) mass approaches zero.

C
The coefficient C governs the contribution of the chromomagnetic operator to the total decay rate Γ Bc . Writing Q(x) in terms of Ψ + (x) and performing the change of field variables in (26) and (30), and retaining only terms up to order v 7 , one obtains This can be used for T c and an analogous expansion in terms of X − for T b . When we need to evaluate matrix elements of the right-hand side of Eq. (41) we go to the restframe of the meson to decompose the field strength into chromo-electric and magnetic components and use 2-component spinors: The chromomagnetic moment operator, with coefficient C above is of order v 7 in the NR expansion, so to this order there are no additional operators resulting from the field redefinition. An electromagnetic moment operator has been ignored since it is further suppressed by the smallness of the fine structure constant.
The coefficient C consists of several contributions: where the normalization constant sets the scale for tree-level decay rate of theb quark. The subscript of the "phase space" factors 12 P i denote the particles in the loop and N i denote the following combinations of Wilson coefficients: Neglecting all fermion masses except for the charm and tau mass one finds the following phase space factors [6]. The light contributions are given by The semi-leptonic mode is given by In the limit of massless tau leptons (x τ → 0) Eq. (48) reduces to Eq. (46). The contributions from Q 1 and Q 2 are given by and for the second insertion by Again, these expressions are given as a function of pole masses, which have to be written in terms of well defined masses. At this level we only have expressions at zeroth order in α s , and for consistency we truncate the expansion of pole masses at zeroth order as well, e.g., m b = m Υ /2 in the Upsilon scheme. The c-quark analogue of the chromomagnetic operator involving b-quarks is given by

C
and it contributes in the following way to the transition operator T c : 13 with the tree-level decay rate and the phase space factors [6] with x s = (m s /m b ) 2 . The Wilson coefficient combinations are defined in Eq. (45) but we have indicated that they are evaluated at a scale µ c that may be different than for C (5) b , since the choice of charm scale µ c ∼ O(m c ) is more appropriate.

Pauli interference
For the dimension-six contributions to the transition operator in Eq. (35) we write, as is customary, the full contribution coming from the operator insertions rather than the individual coefficients (C (6) q ) i . This allows for a more compact presentation of the results [6,7,71]: where p − = p b − p c and z − = m 2 c /p 2 − . We have summed over s and d quarks, neglected their masses and used |V cd | 2 + |V cs | 2 ≈ 1.
The notation p − = p b − p c should be understood in the operator sense, that is, as derivatives acting on the b and c fields. Recall that we have chosen to match the OPE directly to NRQCD; by a slight abuse of notation we have denoted the fields in Eq. (56) as in the full theory, but it should be understood that they really stand for EFT fields, i.e., c → Ψ (c) − . One advantage of matching directly to the EFT is that p − = (m b −m c )u+(k b −k c ) where the derivatives k b −k c correspond to residual momenta, and are clearly sub-leading in the NR expansion. While some authors use

Bc
[7], in our approach the leading term in the NR expansion is unambiguous and the higher order terms correspond to matrix elements of well defined operators. Another advantage of matching to the effective theory is that one may then avoid non-local operators in the expansion. Note that there are inverse powers of p 2 − implicit in the functions of z − in Eq. (56). By matching to the EFT the inverse powers give well defined expansions in terms of local operators, 1/p 2

Weak annihilation
Also for the WA contributions we will report the full contribution to the transition operator instead of just the coefficients (C (6) q ) i . The semi-leptonic contribution to the transition operator reads [6]: with p + = p b + p c and z τ = m 2 τ /p 2 + . As was the case for p − in the PI computation, the notation p + = p b + p c here should be understood in the operator sense, that is, as derivatives acting on the b and c fields. In addition, since we have chosen to match the OPE directly to NRQCD the fields indicated should be those for NRQCD; by a slight abuse of notation, and in order to make the long expression more legible, we have denoted the fields in Eq. (57) and below in Eq. (58) as those of the full theory, but it should be understood that they really stand for fields in the effective theory, i.e., c → Ψ (c) − . One advantage of matching directly to the EFT is that p + = (m b + m c )u + (k b + k c ) where the derivatives in k b + k c correspond to residual momenta, and are clearly sub-leading in the NR expansion, and the leading term is well defined in terms of pole masses and the meson 4-velocity.
The hadronic WA decay gives with z + = m 2 c /p 2 + . We have summed over s and d quarks, neglected their masses and used |V cd | 2 + |V cs | 2 ≈ 1.

Matrix elements
As mentioned above, we denote 2-component spinor fields by the lowercase counterparts of their 4-component Dirac forebearers, Ψ + → ψ q and X − → χ and ψq = iσ 2 (χ † ) T . The first correction in Eq. (35) can be estimated using potential models [72] where T is the expectation value of the kinetic energy computed in potential models. Note that the first equality above is interpreted as (m c v c ) 2 = (m b v b ) 2 which is useful in estimating the NR-quark velocities, v b and v c . For our calculations we estimate the matrix element of D 4 as the square of that of D 2 , from Eq. (59), thus: The leading matrix elements for the chromomagnetic operator are given by and the corresponding matrix elements for the charm quark are obtained by the replacement m b ↔ m c . The wave function at the origin, Ψ(0), relates these to other physical quantities that can be determined from Monte Carlo simulations of NRQCD on the lattice: The matrix elements of four quark operators are all related by spin symmetry per Eqs. (24) and (25). We have Using spin symmetry, we find for the new matrix elements for the penguin operators that enter the calculation of PI: The colour-crossed matrix elements are obtained from these by replacing B Bc → 3B Bc . For WA the additional matrix elements are given by: The "bag parameters" B Bc and B Bc have been chosen so that B Bc = B Bc = 1 in the vacuum insertion approximation. 14 The vacuum insertion approximation can be justified in the large N c limit, so the errors incurred in using these expressions with B Bc = B Bc = 1 are of the order O(1/N c ) and O(v). 15 For the numerical estimates below we adopt this approximation and then analyze the error incurred in the computation of the lifetime by varying B Bc and B Bc away form unity. Eventually, calculations of these matrix elements in Monte Carlo simulations of NRQCD on the lattice will remove this uncertainty. In terms of these parameters WA and PI contributions to the width are obtained from the matrix elements of the transition operators in Eqs. (58) and (56), respectively, yielding: and where z ± = (m b /m c ± 1) −2 .

Numerical Analysis
In this section we present the results for the B c decay width in the MS, the meson and the Upsilon scheme. For our analysis we use the input values summarized in Tab The QCD coupling constant in the QCD corrections are calculated using the 1-loop beta function. As explained in previous sections, the QCD corrections are carried out to 1st order in α s , with the Wilson coefficients computed to NLL order. The running of the Wilson coefficients C 1,2 incorporates analytically the effect of the 2-loop beta function, while that of C 3−6 is computed only at LL; since C 3−6 are very small, the numerical effect of this approximation is negligible in the total rate. As is well known, the consistent counting for resummation of logs involves 2-loop beta functions and anomalous dimensions in the running and 1-loop matching and matrix elements. Hence it is appropriate to include only a 1-loop running α s (µ) in the matrix elements. In particular, including the effects of 2-loop running in the matrix elements only increases the µ dependence of the final results for partial and total decay widths. The renormalization scale µ has been chosen differently for different partial widths. In calculations ofb decays and for WA and PI we use µ = m b (m b ), while for c-decays we use µ = m c (m c ). 16 For the semileptonicb-decays and for WA the m τ -dependence is taken into account, whereas the light leptons are assumed to be massless. Furthermore we neglect the light quark masses. In particular, the strange mass, m s , is neglected in theb-decays.
It is worth repeating that for the computation we take into account QCD corrections truncated to order α s and carry out the non-relativistic expansion up to v 7 (relative order v 4 since the leading order is v 3 ) as presented in the previous sections. Since the power counting in NRQCD has v ∼ α s this is not fully consistent. However, the numerically important effects of WA and PI come in first at order v 6 . Roughly, these effects are amplified by a factor of 16π 2 from the 2-body vs 3-body decay phase space, and suppressed by a relative order v 3 factor of (f Bc /M Bc ) 2 , and 16π 2 (f Bc /M Bc ) 2 = 0.73. Corrections to 3-body decays of relative order v 2 and v 4 are included, so their numerical effect can be estimated and analyzed. We do not include any QCD corrections to WA and PI, since these would correspond to higher-order velocity terms which are neglected in our counting.

Results
The results obtained for each individual channel as well as for the total decay width are collected in Tab. 3 and can be compared to the results obtain previously by BB shown in the first column. For theb-decays the values of all three schemes are significantly smaller than those obtained by BB. Since BB use an on-shell (OS) scheme, with masses declared as having some particular values, a direct comparison is difficult. For one thing, the OS corrections enhance the partial decay widths by up to 21% and 11% for b decays, and where the first number on the right-hand side is the OS value and the second number is the OS correction. The third number corresponds to non-perturbative corrections which are about 4-5% of the partial decay width, which is consistent with the NRQCD counting v 2 ∼ α 2 s ≈ 0.04. We only estimate the effect of a non-vanishing strange quark mass on charm decays, that comprise 60-70% of the total width. Table 4 gives the partial decay rates for c decays including the effect of a non-zero strange quark mass. The partial inclusive width for c decays includes also a contribution of the doubly Cabibbo suppressed c → dus channel, which is simply estimated as |V cd V us /V cs V ud | 2 Γ c→sud . The strange quark mass effect is expected to be suppressed inb-decays relative to c-decays by a factor of ∼ (m c /m b ) 2 ∼ 0.1, and in any case its computation requires a calculation of QCD corrections to the decay b → ccs with m s = 0 which is not available. 18 The total decay width of the B c listed in Tab. 4 includes the contributions ofb-decays and of WA and PI listed in Tab. 3. As a consistency check we confirm that the semi-leptonicb-decay agrees with determination of V cb . The PDG's B 0 partial semileptonic width is 0.068 ± 0.002 ps −1 .
The central values of the perturbative contributions to b → ceν in Tab. 3 are 0.055 ps −1 and 0.071 ps −1 for the MS scheme, and for the meson and Upsilon schemes, respectively. To compare with data these must be corrected for non-perturbative effects.  Table 3: Results for the partial decay rates in ps −1 . The results from BB in [6] in the second column are compared to the results in this paper given in the MS, meson and Upsilon schemes. The hadronic WA and PI contributions of our results include the QCD penguin contributions, which were neglected in BB.
non-perturbative correction is of order 1/m 2 b in the HQET expansion, and is given, as a fraction of the total rate, in terms of the non-perturbative matrix elements in HQET, λ 1,2 , by −(0.12λ 1 + 0.28λ 2 ) GeV −2 , in the meson and Upsilon schemes [32,33]. Using λ 1 = −0.27 ± 0.14 GeV 2 and λ 2 = 0.12 GeV 2 , this gives 0.071 ± 0.001 ps −1 in the meson and Upsilon schemes, where the range indicates only the uncertainty in λ 1 , exclusive of uncertainties in the perturbative calculation (discussed below).

Uncertainties
There are various sources of uncertainty in the calculated width. Use of quark-hadron duality in the formulation of the OPE method introduces an irreducible uncertainty that is, moreover, difficult to quantify. We have little to say about this other than to remind the reader that for semileptonic decays the situation is much more favourable since one may compute the width in terms of an OPE for Euclidean momentum region (that is, for imaginary time) [44]. In effect, in absence of new physics effects, this calculation can  provide a test of the validity of the assumption of quark-hadron duality. Similarly, in the meson and Upsilon schemes we have neglected the nonperturbative correction to the pole mass, and again the calculation may be used to provide a test of this assumption.

Perturbative expansion and QCD-scale uncertainty
The leading contributions to the width of the B c are from the perturbatively calculatedb and c quark decays. For example, forb decays, the perturbative calculation of the partial widths gives The uncertainty from omitting higher order terms in the perturbative expansion is readily estimated as order (α s (µ b )) 2 ∼ 4% and (α s (µ c )) 2 ∼ 10% forb and c decays,   to which the term with the explicit factor of α s ln(µ) in the NLO decay rate is added, displaying cancellation of scale dependence to order α s . The dotted-red line shows the difference between the blue and green lines, that is, the NLO decay rate sans the term with an explicit factor of ln(µ).
respectively. This is reflected in the uncertainty introduced by the arbitrary choice of renormalization scale µ in the calculations. For example, the fractional change in the width for a fractional change in scale, (∆Γ/Γ)/(∆µ/µ), centered at µ = µ b and µ = µ c for b and c decays, respectively, is −0.36, −0.20 and −0.11 for the MS, meson and Upsilon schemes. It is dominated in all cases by the µ-dependence of hadronic decays, and among these, of c decays. The scale dependence formally cancels out to order α s , so the residual dependence is of order α 2 s . The variation with respect to ln(µ) of the leading order is of order α s , and this is canceled exactly by the terms involving α s ln(µ) at NLO. This is shown in Fig. 1 that displays the scale dependence of the NLO calculations of Γ(b → cud) and Γ(c → sud) in the MS scheme. The dashed (orange) line shows the result of the leading order calculation, displaying stronger scale dependence than the result of the NLO calculation, displayed in the solid (blue) line. The latter is the sum of the dasheddot (green) line and the dotted (red) line. This split is intended to demonstrate the approximate cancellation of the µ dependence that comes in at LO from the Wilson coefficients and through the running mass and NLO terms with an explicit α s ln(µ): this combination is shown with the dashed-dot (green) line. Both green and red line should be flat up to corrections of order α 2 s . We can estimate the uncertainty due the scale dependence from the range of variation in the figures, about ±17% and ±35% forb and c decays, respectively. Since these hadronic decays dominate the total width of the B c , one can estimate the total uncertainty by weighing these by their relative importance to the total width, [(0.202)17% + (0.743)35%]/(0.202 + 0.743) = 31%. This can be compared with the fractional change in scale, (∆Γ/Γ)/(∆µ/µ), using ∆µ/µ ≈ ∆ ln µ which gives ∆Γ/Γ = −0.37 ln 2 = −0.26.
Repeating this calculation for the other schemes we infer an uncertainty from order α 2 s /scale dependence of 26%, 14% and 8% in the MS, meson and Upsilon schemes.

Non-relativistic expansion and Non-perturbative uncertainties
Additional uncertainties are introduced by the truncation of the non-relativistic expansion. As already explained we take into account QCD corrections truncated to order α s while carrying out the non-relativistic expansion up to relative order v 4 . Since the power counting in NRQCD has v ∼ α s this is not fully consistent. However, WA and PI come in first at order v 6 and need to be retained because they are numerically important because the 2-body phase space gives an amplification by a factor of 16π 2 . Corrections to 3-body decays of relative order v 2 and v 4 are included, so their numerical effect can be estimated and analyzed. We do not include any QCD corrections to WA and PI, since they would correspond to higher-order velocity terms which we are neglecting. Furthermore, to have a fully consistent calculation at O(v 3 ), three-loop corrections to the rate are required, consisting of three-loop corrections to the matrix elements, a three-loop matching calculation as well as four-loop running. Such corrections are however not available at present and are expected to be smaller than the WA and PI contributions, since they are not enhanced by 16π 2 from the 2-body phase space. Table 5 shows the fractional change in the width per fractional change in individual parameters. This table can be used to adjust the estimated partial widths in Tab. 3 by a change in the input parameters in Tab. 2. The second column in Tab. 5, "∆p/p", gives our best guess of the fractional uncertainty in the parameter p, so the fractional uncertainties ∆Γ/Γ are estimated by taking the product of the second column with the entries in column 3-5. It is seen that the uncertainty introduced by the poor knowledge of the matrix elements that go into the non-relativistic expansion is ≤ 2% for all entries. The reason is that, in fact, the non-relativistic corrections arising from non-perturbative matrix elements (that is, other than from the QCD-perturbative expansion of the free quark decays) are small.
We have already seen in (69) that the non-relativistic corrections are small, of order 4% for Γb →cu(s+d) . The corrections are slightly larger for the otherb decay modes, up to ∼ 10% for the semileptonic modes, but since the width is dominated by the hadronic modes, the non-relativistic correction to Γb →c is 4% and 5% in the MS and meson/Upsilon schemes, respectively. Similarly, we find the non-relativistic correction to Γ c→(s+d) is 6%, 7% and 10% in the MS, meson and Upsilon scheme. The small magnitude of the corrections both suggest the error from truncating the expansion is small, and the uncertainty introduced from the limited knowledge of the matrix elements is small, as expected from the previous paragraphs.
The rate of convergence of the non-relativistic expansion can also be seen by estimating the matrix element of the bilinearQQ, that enters at leading order in the OPE, expressed in terms of NRQCD fields, as given in (37). We find, in the Upsilon scheme, the expansion: where the first line is forbb and the second forcc, and the first, second and third corrections correspond to the terms D 2 , g s σ · B and g s D · E, and D 4 , respectively. As expected the non-relativistic expansion converges faster for b quark decay than for c quarks. We remind the reader that for the coefficient of the correction terms we have not included radiative corrections to the coefficient C (3) q in the OPE; see the discussion at the end of Sec. 5.1. This truncation represents a correction ∼ 6% in thecc case.

Parametric and numerical uncertainties
Additional uncertainties are introduced parametrically, from the uncertainty in the input parameters. The parametric uncertainty from α s is subsumed into the scale uncertainty discussed in Sec. 7.2.1 above. Except for |V cb | and the non-perturbative parameters that enter the non-relativistic expansion, the input parameters in Tab. 2 have negligible uncertainties compared to those that have already been discussed. The uncertainty from |V cb | can be read-off Tab. 5 that gives a fractional uncertainty in Γ Bc slightly smaller than the fractional error in |V cb |: ∆Γ/Γ ∆|V cb |/|V cb |.
For the MS scheme there are additional parametric uncertainties from the input value of the quark masses. The resulting uncertainty is largely due to the fifth power dependence of the partial widths on the mass of the decaying quark. The conservative estimate of the uncertainty assumes that both quark masses deviate in the same direction from the central value quoted in Tab. 2. In Tab. 5 we show the result of varying m b (m b ) and m c (m c ) independently, and in the final result we (over)estimate the uncertainty by adding these linearly, rather than by varying them together.
Additional uncertainties are introduced in the numerical integration of various QCD corrections and the need to compute at non-zero but sufficiently small gluon mass. We have verified that the error introduce by the numerical integration and zero gluon mass extrapolation is much below a percent.

Strange quark mass
The non-vanishing of the strange quark mass reduces the c quark decay rate, relative to the rate for m s = 0, by 7%, 7% and 6% in the MS, meson and Upsilon schemes. A naive estimate of the effect of the strange mass on theb-quark decay is obtained as this fraction times (m c /m b ) 2 ∼ 0.1 of theb-quark width, or ∆Γ b 0.003 ps −1 . More dramatic is the uncertainty derived from the lesser well determined input, m s (2 GeV), listed in Tab. 2. Estimating this as ±10% of the 6-7% reduction in the rate for c-quark decay, this results in an uncertainty ∆Γ c ∼ 0.01 ps −1 .

Conclusions
The result for the decay rate of the B c meson, displaying the contributions by decay channel and for each of the three mass schemes introduced in Sec. 2, is presented in Tab. 3 for m s = 0, and in Tab. 4 for m s = 0 in decay channels involving c-quark decay. In Sec. 7.2 we gave a detailed accounting of uncertainties in those calculated partial widths. We summarize our results: neglecting light quark masses, including that of the strange quark, our final result for the total width of the B c meson in the MS-, mesonand Upsilon-schemes is: where we indicate the uncertainties due to scale dependence (µ), non-perturbative effects (n.p.), the value of |V cb |, and for the MS scheme the input values of the masses (m). As discussed above, for m s = 0, the central values of the c-decay widths get reduced by about 7%, leading to the following total decay rates: For the scale uncertainty we have used Tab. 5 with ∆µ/µ = ln (2). The uncertainty due to non-perturbative effects captures the result of the non-relativistic truncation as well as the uncertainty in the matrix elements. We add linearly the absolute values of the errors from rows 4 -9 (parameters T through f Bc ) using for ∆p/p the values indicated on the second column, and to this we add linearly the absolute value of the estimate of a non-relativistic truncation error. The latter is estimated as a fraction, equal to the quark velocity v b ≈ 0.19 or v c ≈ 0.60, of the non-relativistic correction to theb and c decays, respectively, that has been included in the partial widths recorded on Tabs. 3 and 4. Linearly adding these uncertainties is the most conservative way to proceed. Adding them instead in quadratures leads to the smaller uncertainties ±0.05 ps −1 , ±0.11 ps −1 and ±0.12 ps −1 for both massless and massive strange quarks in the MS, meson and Upsilon scheme, respectively. For either massless or massive strange quarks the three different schemes are consistent with each other within their respective uncertainties. The difference between the values in the meson and Upsilon scheme can mainly be traced back to the different charm mass used in the charm decays. The wide range spanned by the widths calculated in the three schemes, i.e., the strong scheme dependence, however calls for an improvement of the SM prediction. The main uncertainty, from scale dependence, can be reduced by going one order higher in the perturbative expansion, i.e., using NNLL Wilson coefficients, as well as computing the matrix elements up to two loops. While one should see convergence of the various schemes on a single value for the perturbative partial widths of the underlyinḡ b and c decays as the calculation is carried out to higher order in α s , we expect the Upsilon scheme to converge fastest. Evidence for this was presented in Sec. 2.6 using the known 2-loop results for the semileptonic decays of b and c quarks. This is also substantiated by the milder scale dependence of the result in the Upsilon scheme compared to the other schemes. The non-perturbative uncertainty is dominated by our estimate of the order v 3 corrections. These arise solely from order α s corrections to the Wilson coefficients of order v 2 operators. Additional perturbative corrections to Wilson coefficients will have a marginal effect on the overall precision. More precise determinations of the relevant nonperturbative parameters (for example from lattice QCD) might however become important in the future.
It is interesting to note that branching fractions predicted by the three schemes are in good agreement with each other, as shown in Tab. 6. This may be interpreted as evidence that the dominant factor in the differences between scheme predictions is from the sensitive dependence on masses of theb and c quarks, which differ vastly among the schemes. We have not attempted to estimate uncertainties in the branching fractions, which however are expected to be significantly smaller in the ratios than in the individual partial widths because of cancellation of correlated uncertainties.
Previous calculations of the lifetime of B c mesons yield results closer to the experimental measured value than we present here. However, the total width depends sensitively on the choice of masses for the b and c quarks. An arbitrary choice can be made to yield a result close to the experimental value. We have, instead, computed systematically, eliminating the pole mass in favour of well defined masses and truncating the perturbative expansion consistently. This guarantees cancellation of renormalon ambiguities, that remain present if calculating in terms of ad hoc values for the on-shell masses. We find it is interesting to note that the Upsilon scheme results have at once the smallest uncertainty and give partial widths that are seemingly too large. In particular the semileptonic partial width Γ c→(s+d)eν is much larger than in both the two other schemes and the result of BB. However, as we have remarked earlier, it is this higher value that gives excellent agreement with the decay width for inclusive semileptonic decays of D mesons. Similarly, the Upsilon and meson schemes give excellent agreement with the inclusive width for semileptonic B meson decay.
The calculations presented depend on the assumption that quark-hadron duality validates the OPE method. Similarly, for the Upsilon scheme, and to lesser extend for the meson scheme, the nonperturbative contribution to the pole masses have been neglected. These assumptions need to be validated through other means. Were the predictions to fail the conclusion could well be that one or both of these assumptions are not valid.
Finally, the initial motivation for this computation, namely to distinguish NP contributions in the B c life-time remains challenging, and further progress on the theory side has to be made before clear-cut conclusions concerning new-physics effects in Γ Bc can be drawn.