Improving the theoretical prediction for the $B_s-\bar{B}_s$ width difference: matrix elements of next-to-leading order $\Delta B=2$ operators

We present lattice QCD results for the matrix elements of $R_2$ and other dimension-7, $\Delta B = 2$ operators relevant for calculations of $\Delta \Gamma_s$, the $B_s-\bar{B}_s$ width difference. We have computed correlation functions using 5 ensembles of the MILC Collaboration's 2+1+1-flavour gauge field configurations, spanning 3 lattice spacings and light sea quarks masses down to the physical point. The HISQ action is used for the valence strange quarks, and the NRQCD action is used for the bottom quarks. Once our analysis is complete, the theoretical uncertainty in the Standard Model prediction for $\Delta \Gamma_s$ will be substantially reduced.


Introduction
Mixing between particle and antiparticle states of neutral mesons has now been observed in K 0 , D 0 , B 0 , and B 0 s mesons. These mixings are due to couplings between generations of quark SU(2) L doublets after electroweak symmetry breaking. Since the leading-order weak process, represented by the "box" diagrams, is at 1-loop level, there is the chance that new heavy particles, beyond those in the Standard Model, could cause differences between Standard Model predictions and experimental results.
To a good approximation, three parameters are sufficient to describe neutral meson mixing: the moduli of the off-diagonal matrix elements of the 2 × 2 mass and width matrices, M and Γ, and their relative phase φ = arg(−M 12 /Γ 12 ). For the B 0 s system these parameters are constrained by experimental measurements of the B 0 s -B 0 s mass difference, width difference, and a flavour-specific CP asymmetry: ∆M s = 2|M s 12 | , ∆Γ s = 2|Γ s 12 | cos φ s , and a s fs = ∆Γ s ∆M s tan φ s .
In (1) and the remainder of this paper we use notation specific to B 0 s mixing. See Ref. [1] for a recent review of B 0 s mixing and references to a rich literature.
Speaker, e-mail: M.Wingate@damtp.cam.ac.uk The calculation of ∆Γ s within the Standard Model is summarized in Ref. [2]. Contributions to Γ s 12 come from matrix elements of the non-local product of 2 ∆B = 1 effective Hamiltonians The contributions from charm and up quarks in the intermediate states depend on the corresponding CKM matrix elements; i.e. Γ s At the present level of accuracy, only the CKM-leading contribution from Γ cc 12 is important. Direct lattice QCD calculation of matrix elements such as B s |T |B s is not generally feasible due to the difficulty in correctly treating all intermediate states. 1 However, one can employ an operator product expansion known as the heavy quark expansion (HQE). Order-by-order in Λ QCD /m B , one relates the matrix elements of nonlocal operators to a series of matrix elements of local ∆B = 2 operators. Using the most advantageous basis [2] the charm-charm loop contribution to Γ s 12 is given by with the next order in the HQE given bỹ A full basis of dimension-6 ∆B = 2 operators can be written as At higher order in the HQE, one needs matrix elements of the following operators Matrix elements of the Q i operators (5) have long been calculated using lattice QCD; unquenched results for B s mixing appear in [4][5][6][7]. Until this work there have been no calculations of R 2 and R 3 matrix elements. In phenomenological analyses [2,8] the vacuum saturation approximation was used, allowing a 50% uncertainty. Sum rule estimates suggest these matrix elements should be within a few percent of the vacuum saturation approximation values [9], although the VSA predictions depend sensitively on the value of the b-quark pole mass.

Method
We use MILC's HISQ ensembles, which include sea quark effects of degenerate up and down quarks and physical-mass strange and charm quarks [10,11]. We use the HISQ action for the valence s quark and the NRQCD action for the b quark. The 5 ensembles include 3 distinct lattice spacings which we respectively refer to as fine (F), coarse (C), and very coarse (VC). For each of these spacings we use configurations with dynamical pion mass of approximately 300 MeV, and for the coarse and very coarse spacings, we used the physical ensembles which have pion mass approximately 130 MeV. Table 1 lists specific input parameters and the lattice spacings as determined from the Υ(2S − 1S ) splitting [12,13].
In carrying out the calculation of B s |R i |B s , with i = 2, 3, we need not compute all 4 terms in the Lorentz dot product (6). The temporal derivative acting on the b field is O(m b ):b ← D 0 = ±m bb γ 0 , the sign depending on whether we have an outgoing b quark or incomingb antiquark. Thus we can write Applying the equations of motion, iγ 0 D 0 s = ( γ · D)s, we arrive at The lattice calculation of B s |R 2,3 |B s proceeds just as for the Q j matrix elements, except for the need to have a derivative operate on the strange quark at the operator. In addition to needing a staggered propagator g(y, z) computed from local source [14] K we need propagators from a point-split source (k = 1, 2, 3) Naive quark propagators are constructed from staggered propagators via where Ω(x) = 3 µ=0 (γ µ ) x µ /a . Since we will need to sum over spatial directions, we require 4 strange quark inversions on each configuration and for each source location.

Perturbative matching
One-loop matching between the lattice theory and the continuum MS renormalization schemes has been carried out for the Q i operators, including tree-level 1/m b corrections [15]. The 1-loop mixing between operators is parametrized by the ρ i j matrix and 1/m b corrections are given byQ sub i,1 , which are of the form 1 Because derivatives are implemented as finite difference operators the 1 aQ i mix withQ i,1 ; this can similarly be computed in perturbation theory. We define a subtracted operator which gives a more accurate determination of the next-to-leading contribution: The coefficients ρ i j and ζ i j are tabulated in [15]. A similar subtraction is done here for the R 2,3 operators: Values for ξ i j are given in Table 2. We use the α V values as tabulated in [16]. Note that we have not carried out the 1-loop matching between lattice and MS schemes forQ sub i,1 orR sub i . Therefore our results for their MS matrix elements will have an O(α s ) systematic uncertainty.

Fits to correlation functions
On each of the 1000 or so configurations in the 5 ensembles listed in Table 1, we created strange quark propagators with inversion sources on 2 timeslices per configuration -except for the VC5 ensemble where we weighed the benefits of doubling the number of sources per configuration. 2 We calculated 3-point functions with local B s andB s sinks as well as Gaussian smeared sinks with 2 radii. The smearing was done with the links fixed to Coulomb gauge. The 2-point correlators are taken from earlier work where 16 sources per configuration were used [13].
Correlator data are fit to functions of the form 2 We concluded that increased statistics were not sufficiently beneficial to warrant the cost of doubling the data set on other ensembles.
using the corrfitter package [17]. The oscillating states in (16) and (17) appear due to oppositeparity temporal doublers present in staggered quark formulations. In (17), t is the temporal distance between the initial state interpolating operator and the 4-quark operator and T is the distance between the initial and final state interpolating operators. Values used in the fits presented here are given in Table 3. The Gaussian priors for the fit amplitudes and energies were set as follows. We first performed 2exponential fits (N = 2 exponentials in each parity channel) to the 2-point data using wide priors and t min 1.2 fm. From the output of this fit we took the ground state energy and amplitude, multiplied their uncertainties by 10, and used this as the prior means and half-widths for subsequent fits. For the excited states, we took the energy splittings to be O(aΛ QCD ) ± 50% and the amplitudes to be 0 ± 1.
After fixing the priors for the energies and 2-point amplitudes, we performed N = 3 fits to 3-point correlator data with t min ≈ 1.0 fm and 2 large values of T using 0 ± 1 for the priors on the V fit parameters. This gave an order-of-magnitude estimate for the ground state amplitude. In subsequent fits we set the prior on V nn,00 to be the fit result ±50 − 100%; for the amplitudes in the oscillating terms, we used standard deviations of 100 − 400% of the results from the preliminary fits.
In the fits presented below, we found that convergence was improved by first fitting the 2-point correlator data and using the results as priors for the fits to the 3-point correlators. In most cases the difference between these "chained" fits and fully simultaneous fits is not significant [18]; however, there were some cases where the simultaneous fits failed to converge.
In Fig. 1 we show preliminary results of fits to the 3-point amplitudes V nn,00 determining the R 2 and R 3 matrix elements. We observe results which give consistent results once enough exponentials are included to account for excited state contributions to the correlators. In order to obtain this, it was necessary to impose an SVD cut of 0.001 on the correlation matrix whose smaller singular values are not well-determined by the data. Figure 2 shows preliminary results vs. a 2 for matrix elements of R 2 and R 3 after subtraction (14), on the VC5, C5, and F5 ensembles. The fits on the physical point ensembles VCp and Cp are not as far along in the process of being checked. We are still assessing fitting uncertainties and ensuring results are robust against different fitting choices. What we present here are the results of separate fits to correlators for the operators R 2 , R 3 , Q 1 , and Q 2 . Once we have finished investigating these fits, we will form the linear combinations of correlators, configuration-by-configuration, which will allow us to determine matrix elements of R sub 2 and R sub 3 directly.

Outlook
We presented preliminary results for B s |R 2 |B s and B s |R 3 |B s on 5 ensembles spanning a range of 3 lattice spacing and including 2 physical mass ensembles. We are presently verifying stability of fit   results. The statistical precision may be improved by performing fits to the linear combinations of correlators directly yielding the 1-loop subtracted matrix elements. The results from different ensembles will then enable an assessment of discretization and quark-mass tuning effects. We expect the dominant uncertainty to be due to the O(α s ) difference between lattice and continuum regularization schemes. This will be the first time these matrix elements have been computed using lattice QCD.