Lattice QCD matrix elements for the Bs-Bsbar width difference beyond leading order

Christine T. H. Davies, Judd Harrison, 1 G. Peter Lepage, Christopher J. Monahan, 5, 6 Junko Shigemitsu, and Matthew Wingate a (HPQCD Collaboration) b SUPA, School of Physics and Astronomy, University of Glasgow, Glasgow, G12 8QQ, United Kingdom DAMTP, University of Cambridge, Cambridge CB3 0WA, United Kingdom Laboratory of Elementary Particle Physics, Cornell University, Ithaca, NY 14853, United States Institute for Nuclear Theory, University of Washington, Seattle, WA 98195-1550, United States Physics Department, College of William and Mary, Williamsburg, VA 23187, United States Thomas Jefferson National Accelerator Facility, Newport News, VA 23606, United States Department of Physics, Ohio State University, Columbus, OH 43210, United States Predicting the B s −B̄ s width difference ∆Γs relies on the heavy quark expansion and on hadronic matrix elements of ∆B = 2 operators. We present the first lattice QCD results for matrix elements of the dimension-7 operators R2,3 and linear combinations R̃2,3 using nonrelativistic QCD for the bottom quark and a highly improved staggered quark (HISQ) action for the strange quark. Computations use MILC ensembles of gauge field configuations with 2+1+1 flavors of sea quarks with the HISQ discretization, including lattices with physically light up/down quark masses. We discuss features unique to calculating matrix elements of these operators and analyze uncertainties from series truncation, discretization, and quark mass dependence. Finally we report the first Standard Model determination of ∆Γs using lattice QCD results for all hadronic matrix elements through O(1/mb). The main result of our calculations yields the 1/mb contribution ∆Γ1/mb = −0.022(10) ps −1. Adding this to the leading order contribution, the Standard Model prediction is ∆Γs = 0.092(14) ps −1.

Oscillations of neutral mesons into their antiparticles have been important phenomena in the study of quark flavor. This flavor-changing, neutral mixing is absent in the Standard Model (SM) at the classical level; appearing at one-loop level it is suppressed by two powers of Fermi's constant G F relative to hadronic and quark mass scales. A few observables are, to a high level of precision, sensitive only to short-distance physics. Within the Standard Model, predictions for these are reliably calculable because the dominant contribution comes from top-quark loops; there is no significant contribution from intermediate-state hadronic physics. Prime examples are: in the mixing of strange mesons K 0 −K 0 , the indirect CP-violating ratio ϵ K and, in beauty mesons B 0 d,s −B 0 d,s , the mass differences (equivalently, the oscillation frequencies) ∆M d,s . Precise experimental measurements of these, together with accurate Standard Model predictions, constitute stringent tests of the SM description of quark flavor.
Beyond these observables are others, where contributions from hadronic intermediate states must be included. For example, mixing in neutral charm mesons D 0 −D 0 has significant long-distance contributions due to differing flavor structure from K 0 and B 0 mixing. Predictions for heavy meson and baryon lifetimes also require theoretical treatment of long-distance effects.
The B 0 s −B 0 s width difference ∆Γ s is another example and is the focus of this paper. Unlike the mass difference, which comes from the real part of the mixing a M.Wingate@damtp.cam.ac.uk b http://www.physics.gla.ac.uk/HPQCD amplitude, the width difference comes from the imaginary part which, by the optical theorem, describes the decays to real final states, primarily b → ccs decays. Consequently, ∆Γ s is dominantly due to processes with intrinsic charm and we expect it to be insensitive to new physics. Comparison between theory and experiment is a test of the theoretical methods involved. Agreement here is a necessary condition for trusting these methods to reliably yield a SM prediction for quantities where new physics could contribute more prominently, e.g. in D 0 −D 0 mixing.
Predicting the SM width difference ∆Γ s requires the determination of matrix elements of a nonlocal product of effective operators H ∆F =1 eff , with charm and up quarks in the virtual loops. Direct calculation of (0)}|B 0 s ⟩ using lattice QCD is not presently feasible. Therefore, an additional theoretical approximation is necessary in order to obtain a Standard Model prediction for ∆Γ s , namely the heavy quark expansion (HQE). This expansion makes use of the large b-quark mass compared to the c-quark and other scales in the problem such as Λ QCD , and approximates the imaginary part of the matrix elements above by a power series in 1/m b , composed of matrix elements of local, ∆F = 2 operators such as those appearing in H ∆F =2 eff [1]. (See Ref. [2] for a recent review of the HQE applied to ∆Γ s .) Matrix elements of the leading, dimension-6 operators in H ∆F =2 eff have been calculated using lattice QCD with increasing precision, motivated by their impact on predictions for ∆M d,s . Results are now available from several groups [3][4][5][6] (see Ref. [7] for a review and Refs. [8,9] for sum rule calculations). The precision of these determinations has become good enough that matrix elements of higher-dimension operators are needed in order to re-duce the SM uncertainty in ∆Γ s . Current estimates for the higher-dimension matrix elements come from the vacuum saturation approximation [2], although sum rules have also been applied [10]. The lack of any full QCD calculation of dimension-7 operators is the leading uncertainty in the Standard Model determination of ∆Γ s [2,11].
In this paper, we provide results of the first complete lattice QCD calculation needed for ∆Γ s through O(1/m b ). The goal here is to replace order-of-magnitude estimates based on the vacuum saturation approximation with first principles calculations including a quantitative analysis of errors. In this first step, we neglect O(α s ) corrections to the dimension-7 operators. Including these corrections involves a technically challenging perturbative calculation which is left as a future project.
The convention we use for the dimension-6 operators O 1−5 is as in Ref. [6]. At higher order in the HQE, one needs matrix elements of the following operators as well as the color-rearranged partnersR 1 ,R 2 , and R 3 . The perturbative coefficients α 1 and α 2 are given in Refs. [12][13][14]. Using Fierz identities and neglecting terms at higher order in 1/m b we havẽ In nonrelativistic lattice QCD, matrix elements of both sides of (2) hold exactly. We carry out our calculations using gauge field configurations generated by the MILC Collaboration [15][16][17]. These include the effects of 2 + 1 + 1 flavors of sea quarks using the HISQ fermion action [18,19]. We use five separate ensembles, those labeled sets 1, 3, 4, 6, and 7 in [6]. (Also see the Appendix.) Two of the ensembles have all of the quark masses tuned to be close to their physical values, e.g. m s /m l ≈ 27; these have lattice spacings of 0.12 and 0.15 fm. The other three ensembles span 3 lattice spacings from 0.09 to 0.15 fm, with unphysically large light-quark masses corresponding to pion masses of about 300 MeV, m s /m l ≈ 5. Our use of three lattice spacings allows us to estimate discretization errors, and the computations done with unphysical light quark masses gives us information with which to correct any slight quark mass mistunings.
The lattice actions used are the same as in our recent study of the dimension-6 operator matrix elements [6]. Correlation functions are computed using the HISQ action for the strange quark; the valence quark mass is tuned to be closer to the physical strange mass than the value which was used for the sea strange quark. The nonrelativistic QCD (NRQCD) action [20,21] is used for the bottom quark. The action coefficients and valence quark parameters are the same as in recent work [6,21,22]. Because the determination of ⟨R 2 ⟩ and ⟨R 3 ⟩ here will have an O(α s ) uncertainty due to tree-level matching between lattice and continuum regularization schemes, we only need a fraction of the statistics used in Ref. [6]. We will occasionally refer to Ref. [6] as the high-statistics companion to this work. Throughout this paper we will use the abbreviated notation ⟨·⟩ ≡ ⟨B s | · |B s ⟩.
Let us examine a unique feature of computing ⟨R 2 ⟩ and ⟨R 3 ⟩. In the rest frame of the heavy quark, only the temporal component of the following bilinear is important at where the sign is determined by whether the temporal derivative acts on an outgoing heavy quark or an incoming heavy antiquark. Γ represents either γ µ (1 − γ 5 ) (R 2 ) or 1 − γ 5 (R 3 ). Using the equation of motion for the strange quark, and neglecting contributions of O( ms m b ), , and similarly for R 2,3 . In order to implement the derivative operator, the calculation of the associated three-point correlation functions will require new strange quark propagators, beyond those used in [6], with point-split inversion sources at the operator location. We use a symmetric difference operator in each of the spatial dimensions as the source for the inversion of the HISQ Dirac matrix on Coulomb-gaugefixed configurations. Note we use a hat on an operator when we wish to call attention to the 4-quark operator computed directly on the lattice, in distinction to a linear combination such as a renormalized, matched, or subtracted operator.
As discussed below, we will also need matrix elements of the dimension-6 operatorsÔ 1 andÔ 2 . At no additional cost, we recompute these here and check that they agree with the high-statistics study [6].
Matrix elements are obtained from fits to two-and three-point correlation functions using methods developed in Refs. [6,[23][24][25]. Some details have appeared in [26]. (Further details are provided in the Appendix.) The prediction of a matrix element of a higher dimension operator using lattice NRQCD is complicated by mixing with lower dimension operators [27]. The presence of the lattice cutoff a means that the matrix elements ⟨R 2,3 ⟩ will contain contributions from ⟨Ô 1,2 ⟩ of the order O(α s /(am b )). We have used lattice perturbation theory to calculate the coefficients ξ ij which parametrize this mixing at one-loop level. Matrix elements of the subtracted operator, will have power-law mixing cancelled through O(α s ).
The coefficients ξ ij have not been calculated before. The procedure is a straightforward extension of Ref. [28], in particular Sec. IV.B. In this instance the derivative acts on the light quark propagator instead of the heavy quark propagator. Numerical values are tabulated in Table I. For the numerical value of the strong coupling constant, we use α V (2/a) (see Table I of [29], inferred from the work of [30,31]). Fig. 1 illustrates the effect of the subtraction (4). Matrix elements ofR i and R sub i are shown normalized to ⟨O 1 ⟩ [6]. The fact that the matrix element of the subtracted operator is 50% (R 2 ) to 70% (R 3 ) of the unsubtracted operator shows that this subtraction is significant. This is of comparable size to the effects seen in the 1/m b contributions to the B and B s decay constants [27,32] and matrix elements ⟨O j ⟩ [33]. Fig. 1 shows that, as expected, the subtraction leads to a reduction in a-dependence over the range of lattice spacings used.
We must estimate the uncertainty due to not knowing the O(α s ) MS-to-lattice matching of the R operators nor the O(α 2 s ) contributions from the dimension-6 operators. Both of these are suppressed by a power of α s compared to matrix elements of the two terms in (4). Therefore, we include these truncation errors in our results by multiplying our results by a noisy esti- where δ am b = 0 ± 1 is a Gaussian-distributed random variable, one for each of the 3 lattice spacings. This is our largest source of uncertainty.
Our calculations include numerical data obtained with all quark masses tuned close to their physical values. In order to include some data with unphysically large values for the up/down sea quark masses, we assume an analytic dependence on these masses. We also parametrize discretization errors, e.g. due to the gluon and staggered fermion actions, in powers of (aΛ QCD ) 2 taking Λ QCD = 500 MeV. Results presented here come from fits to where sea quark mass dependence is parametrized by x ′ ℓ = M 2 π /(2Λ χ ) 2 and x ′ s = (2M 2 K − M 2 π )/(2Λ χ ) 2 , with Λ χ = 1 GeV. We use the lattice masses aM π and aM K tabulated in [17]. For the valence strange quark x s = M 2 ηs /(2Λ χ ) 2 , the η s being a fictitious flavor-nonsinglet ss pseudoscalar meson, which is nevertheless well-defined in chiral perturbation theory. Values for aM ηs on these ensembles [34] are used to parametrize the difference between the input valence strange quark mass and the physical one, using for the "physical" value M ηs = 688.5(2.2) MeV [34,35]. We assume Gaussian priors of 0 ± 1 for the fit parameters, except for d 2 , which we take to be 0 ± 5. In the fits, we find |d 2 | ≈ 2.5 ± 2.0.
As one would expect, any sensitivity of the B s −B s matrix elements to the light sea quark mass is much smaller than our uncertainties. The slight mistunings in the sea or valence strange quark masses are not large enough to give a nonzero result for c s 1 parameters. Including terms quadratic in the x variables has no effect on the fit. Similarly, we obtain d 4 consistent with zero. Within errors the data is consistent with a mild a 2 dependence. Our main results are for the pair of matrix elements or for the linearly dependent, color-rearranged operator matrix elements (2) At the accuracy with which we work, these results can be interpreted as the MS-scheme results at the scale µ 2 = m b , with an uncertainty included to account for the fact that the lattice-continuum matching is tree-level. It is sometimes convenient to include a numerical factor which arises from the vacuum saturation approximation (VSA). We define 7 6 , and 5 6 for R 2 ,R 2 , R 3 , andR 3 , respectively. The "unprimed" bag factors B R k which are equal to 1 in the VSA include a mass factor such that . These B R k are the ones taken in recent phenomenological estimates [2,11,13] to be 1.0 ± 0.5 in the absence of a QCD calculation. We tabulate numerical results of our work in Table II. It turns out that the VSA expectation is a reasonable back-of-the-envelope estimate. (Note that while the bag factors depend on the definition of m pow b , the B ′factors do not.) Our results replace the rough estimates with a lattice QCD computation with all uncertainties quantified.
The matrix elements determined in Ref. [6] allow determination of the remaining 3 matrix elements The 1/m b contribution to ∆Γ s can be expressed as a linear combination of the matrix elements of the R operators, times perturbative coefficients γ k [1,13]. Writing ∆Γ 1/m b = −2Γ 12,1/m b cos ϕ 12 , we havẽ Here k is an index that runs over the 4 operators in (1) plus the 3 color-rearranged operators. The γ k (z) are related to the g k (z) of [13] by the numerical coefficients relating the matrix elements to the B ′ factors; additionally γ 1 andγ 1 include a factor ofm [38] and the leading order H ∆F =1 Wilson coefficients C Our result for (11) isΓ 12,1/m b = 0.0110(52) ps −1 which, given cos ϕ 12 = 1 to the precision relevant here, contributes to the width difference as The uncertainty in (12) is dominated by that of ⟨R 2 ⟩. From studies of the leading ∆Γ s term, we expect scale and scheme uncertainties here to similarly be at the 10% level, i.e. not significant compared to the present hadronic uncertainty.
The error in (13) is mostly due to the uncertainty in ∆Γ 1/m b ; its variance contributes approximately 60% to the total variance in (13). The next largest uncertainty, contributing 30%, comes from the coefficients in first term of (10). The variance of B 1 contributes 8% of the total. The HFLAV average of experimental measurements is ∆Γ s = 0.085(6) ps −1 [39]. This is in good agreement with our result (13) with half the uncertainty.
There remains more to do in order for the theoretical prediction to match the experimental precision. The next generation lattice calculation will require one-loop matching of lattice to MS regularization schemes in order to reduce the uncertainty in ∆Γ 1/m b . At the same time the work to determine the perturbative coefficients appearing in (10) through NNLO must be completed. First steps have already begun [14].   Table I of [22]; errors are statistical, NRQCD systematic, experiment respectively.

APPENDIX
In this Appendix we collect a few pertinent details from the literature cited in the main body, and we give a few numerical values for the short-distance coefficients entering the Standard Model calculation of ∆Γ 1/m b .
The parameters of the MILC gauge field ensembles we used are given in Table III [15][16][17], and the valence quark parameters we used in Table IV [6,21,22].
In order to determine the energies and amplitudes associated with the meson creation and annihilation operators, we perform multi-exponential Bayesian fits [23]. The fit functions for the two-point and three-point func-tions are, respectively, and, abbreviating terms containing any Y parameters, In practice, we fit to energy differences for all but the ground state energy E 0 . The X and Y parameters are the amplitudes for meson creation/annihilation. The labels a and b run over the 3 types of smearings used with the B s meson interpolating operators, a local operator and two Gaussian-smeared operators with different widths [29]. The V nn,ij are parameters related to the matrix elements of 4-quark operators. Not all fit parameters are wellconstrained by the data, so we introduce Bayesian priors with means and widths chosen as described in Ref. [26]. We take the two-point functions from earlier studies of B (s) decay constants [22,29]; these were obtained with 16 time sources on each gauge-field configuration, allowing precise determination of the ground state energies and decay amplitudes. Given that our determinations of the matrix element of R 2,3 will have O(α s ) truncation errors, we do not need such high statistics. The three-point correlators in this work come from using just 2 time sources on each configuration.
We find the best approach in this case is to use chained [25], marginalized fits [24]. The high-statistics two-point functions are fit using N 2pt = 5. The resulting E I and X a,i , central values and errors, are used as priors for the fits to the three-point functions. Excited state contamination is accounted for in the three-point functions by incorporating noisy estimates using the results of the two-point fits, then fitting using N 3pt = 1 [6,24].
While we perform separate fits to the correlation functions associated with eachR 2,3 andÔ i in order to study the relative contributions to the subtracted matrix elements, for the main results we perform fits to linear combinations of three-point functions so that the subtracted matrix element is directly extracted from the fit. This allows correlations to be propagated straightforwardly. Table V gives numerical values for the short-distance quantities used in combination with our matrix elements to determine ∆Γ 1/m b . Numerical values for the Wilson coefficients have been taken from Table 2 of Ref. [14]. For the charm quark mass we use the world average [41] of lattice results with 2 + 1 + 1 flavors of sea quarks [31,[41][42][43], and for the bottom mass and the ratio m b /m c we use the result from Ref. [43].