Neutral Meson Mixing in Unquenched Lattice QCD
Abstract
We study and mixing in unquenched lattice QCD employing the MILC collaboration gauge configurations that include , and sea quarks based on the improved staggered quark (AsqTad) action and a highly improved gluon action. We implement the valence light quarks also with the AsqTad action and use the nonrelativistic NRQCD action for the valence quark. We calculate hadronic matrix elements necessary for extracting CKM matrix elements from experimental measurements of mass differences and . We find , MeV and MeV. We also update previous results for decay constants and obtain MeV, MeV and . The new lattice results lead to updated values for the ratio of CKM matrix elements and for the Standard Model prediction for with reduced errors. We determine and .
pacs:
12.38.Gc, 13.20.Fc, 13.20.HeHPQCD Collaboration
I Introduction
The mass differences and between the “heavy” and “light” mass eigenstates in the neutral meson system have now been measured very accurately leading to the possibility of a precise determination of the ratio of two CabibboKobayashiMaskawa (CKM) matrix elements pdg (); cdf (); d0 (). This ratio is an important ingredient in fixing one of the sides of the “Unitarity Triangle”, and hence plays a crucial role in consistency checks of the Standard Model. Reaching the goal of determining from the experimental ’s, however, requires theory input on hadronic matrix elements of certain fourfermion operators sandwiched between the and states. Information on such hadronic matrix elements can only be obtained if one has control over the strong interactions, QCD, in the nonperturbative domain.
Set  size  

C1  2.645  0.005/0.050  0.005  677  4/2  
0.040  
C2  2.635  0.007/0.050  0.007  834  4/2  
0.040  
C3  2.619  0.010/0.050  0.010  672  4/3  
0.040  
C4  2.651  0.020/0.050  0.020  459  4/3  
0.040  
F1  3.701  0.0062/0.031  0.0062  547  4/2  
0.031  
F2  3.721  0.0124/0.031  0.0124  534  4/2  
0.031 
In this article we use lattice QCD methods to calculate the hadronic matrix elements that appear in the Standard Model to describe neutral meson mixing. Simulations are carried out on unquenched configurations created by the MILC collaboration milc1 (). These configurations include effects from vacuum polarization due to three light quark flavors, , and ( configurations, where equals the number of sea quark flavors). The and quark masses are set equal to each other. Table I lists the six different ensembles used, together with their characteristics such as the number of configurations, sea quark masses in lattice units, the valence quark masses employed for each ensemble, number of time sources and the number of smearings for the quark propagators. Information on the lattice spacing is presented in terms of the ratio , where is obtained from the static potential and has been calculated by the MILC collaboration for each of their ensembles r1 (). The bare and quark masses have been fixed already in previous simulations of the upsilon () and Kaon milc2 () systems. The MILC collaboration unquenched configurations have been created using the “fourth root” procedure to remove the four fold degeneracy of staggered fermions and some theoretical issues remain concerning the validity of this procedure. Considerable progress has been made, however, in addressing this important issue fourthrt0 () and several recent reviews fourthrt () summarize our current understanding of the situation. In this work we assume that physical QCD is obtained in the continuum limit, as implied by existing evidence.
In a previous article the HPQCD collaboration presented the first unquenched results for meson mixing parameters, based on simulations on two out of the above 6 MILC ensembles (sets C3 and C4) bsmixing (). In the present work we broaden considerably the scope of our studies of mixing phenomena. We generalize to include both and mixing and we use two sets of lattice spacings (the first 4 ensembles in Table I have fm and are called “coarse” whereas the last two have fm and are refered to as “fine” lattices). Furthermore we now employ smeared operators for the meson interpolating operators and even on those ensembles used previously we have doubled the statistics, by going from two to four time sources. Unquenched lattice calculations by other groups exist in the literature. Several years ago the JLQCD collaboration published studies of and mixing jlqcd () and the Fermilab/MILC collaboration has recently presented preliminary results based on the same AsqTad light quarks as in the present article, however using a different action for the quarks tod (); elvira ().
In the next section we summarize the formulas needed for analysis of meson mixing phenomena. We introduce the relevant fourfermion operators and describe how their matrix elements are parameterized and how they can be related to the CKM matrix elements and . We then discuss the lattice fourfermion operators used in the simulations and how they can be matched onto the operators in continuum QCD. In section III we describe our simulation data and the fitting procedures one must go through in order to extract the matrix elements of interest. Section IV focuses on chiral and continuum extrapolations and section V presents results for , and together with discussions of systematic errors. This section, section V, summarizes the main results of the present work for quantities most directly associated with mixing analysis. As part of our simulations, however, we have also accumulated more data on meson decay constants, and . Hence in section VI we update the results for these decay constants published previously in fbprl (); fbsprl (). Section VII presents a summary of the current work and a discussion of future directions in our program.
We conclude this introductory section with a comment on notation. The decay constants with are defined in eq.(2) below and are used together with appropriate bag parameters to parameterize fourfermion operator matrix elements. can of course be identified with , the decay constant of the charged B mesons, and can be measured through the latter meson’s leptonic decays. The meson, on the other hand, cannot decay leptonically via a single W boson and hence by itself is not a directly measurable quantity in the Standard Model. Although is the more physical quantity we use the notation throughout this article in order to facilitate uniform treatment of and mixing.
Ii Relevant Four Fermion Operators and Matching
Neutral meson mixing occurs at lowest order in the Standard Model through box diagrams involving the exchange of two bosons. These box diagrams can be well approximated by an effective Hamiltonian expressed in terms of fourfermion operators. More specifically, for calculations of in QCD one is interested in the operator with [VA] x [VA] structure,
(1) 
where and are color indices and are summed over. The symbol stands for either the or the quark. Working in the scheme, it is customary to parameterize the matrix element of between a and a state as,
(2) 
Here is the meson decay constant and its “bag parameter”. Factors such as ensure that in the “vacuum saturation” approximation. Given the definitions in (2) the Standard Model prediction for the mass difference is buras (),
(3) 
where depends on the quark and the boson masses and , is a perturbative QCD correction factor and the InamiLim function. is the renormalization group invariant bag parameter and at twoloop accuracy one has for the present case. From eq.(3) one sees that an experimental measurement of would yield directly the CKM matrix element combination provided the quantity is available. One also sees that the ratio can be obtained from,
(4) 
The goal is to evaluate the hadronic matrix element in eq.(2) using lattice QCD methods. Several steps are required in going from what is actually simulated on the lattice to the scheme quantities appearing in the continuum phenomenology formulas. One important step is to relate fourfermion operators in continuum QCD to operators written in terms of the heavy and light quark fields appearing in the lattice actions that we employ. Another crucial step will be to correct for the fact that simulations are carried out at nonzero lattice spacings and with light quark masses larger than the or quark masses in the real world. In the remainder of this section we address the first step, namely matching between the continuum QCD operator and its counterpart in the effective lattice theory that we simulate. The other step of chiral and continuum extrapolations will be discussed in section IV.
Our simulations are carried out using the improved staggered (AsqTad) quark action for the light quarks stagg () and the nonrelativistic (NRQCD) action for the heavy quarks nrqcd (). Matching through for the lattice action of this article was completed in reference fourfmatch (), where is the heavy quark mass. We refer the reader to that paper for details and just summarize the most important formulas here. In effective theories such as NRQCD one works separately with heavy quark fields that create heavy quarks () and with those that annihilate heavy antiquarks (). The operator that contributes to mixing at treelevel and that matches onto (1) at lowest order in has the form,
(5)  
As is well known, even at lowest order in 1/M there is a oneloop order mixing with another fourfermion operator,
(6)  
This is true both in NRQCD and in HQET. If one introduces an effective theory field,
(7) 
then and the QCD field are related by a FoldyWouthuysenTani (FWT) transformation. In particular,
(8) 
where the acts to the left. The FWT transformation determines the treelevel 1/M corrections to the fourfermion operators in the effective theory. For they come in as,
(9)  
Taking these corrections into account one can work through and finds the following matching relation,
(10) 
The matching coefficients , , and are listed (for ) in fourfmatch (). As explained there, the terms proportional to are needed to remove power law contributions in the matrix elements .
Iii Simulation Data and Fitting
The starting point for a lattice simulation determination of , with , or , is the calculation of the threepoint correlator,
One works with dimensionless operators which are kept fixed at the origin of the lattice. is an interpolating operator for the meson of smearing type “”, and spatial sums over and ensure one is dealing with zero momentum and incoming and outgoing states. The meson is created at time and propagates to time slice where it mixes into a meson. The meson then propagates further in time until it is annihilated at time . We have accumulated data for with on the coarse lattices and on the fine lattices. Given the well known properties of staggered light quarks, for fixed the threepoint correlator must be fit to
This ansatz allows for nonoscillatory and oscillatory contributions to the correlator (in practice we have worked with ). Not all the amplitudes etc. are independent due to symmetries. Similarly twopoint correlators are fit to,
The relation between the amplitudes or the and the matrix elements of the previous section can be identified as follows.
(14) 
The energy eigenstates in the numerator are taken to have conventional relativistic normalization and the factors in the denominator are needed to make up the difference between this continuum normalization and the one in the effective lattice theory. For the ground state contribution , and recalling that , one has,
(15) 
which includes the matrix element that we are interested in. Similarly for the 2ptfunctions one has,
(16) 
Using one then has,
(17) 
In order to assemble all the terms on the RHS of (II) we have tried two approaches. In the first approach we did separate fits for each of the operators , and and inserted their ground state matrix elements into (II). In the second approach we went through the analysis in the opposite order. Namely we first obtained the renormalized fourfermion operator at the threepoint function level by forming the appropriate linear combinations of the ’s, and then carried out fits to extract for the full renormalized threepoint function. Consistent results were obtained from the two methods. For our final analysis we adopted the second approach which we found to be more convenient in practice.
Our smearings consist of Gaussian smearings of the heavy quark propagator at both source and sink. In addition to point sources and sinks we use Gaussians with widths, in units of the lattice spacing, of 2.0 and 6.0 for sets C1 and C2 and one Gaussian each of width 6.0 for C3 and C4 and of width 8.0 for sets F1 and F2. To extract for our renormalized threepoint function we carry out simultaneous fits to an matrix of twopoint correlators (eq.(III) with ) ( corresponds to local, to first Gaussian etc.) and to the renormalized threepoint functions with . Bayesian fitting bayes () methods are employed to enable these complicated fits with large numbers of exponentials, i.e. of fit parameters. We fit to all data points within and for , and . We have used ranging between 4 to 16 and looked for consistency in fit results as the number of exponentials was increased. An example of fit results on one of the coarse lattices is shown in Fig.1. One sees that good and consistent results are obtained for . When becomes very large (in the case of Fig.1 ), errors tend to increase again indicating that our fit ansatz has become too complicated for the minimization routines to handle, given the amount of statistics that we have. Fig.2 shows an example for one of the fine lattices. Here we find good results for . In general we have relied on our and fits for all our ensembles.
We summarize fit results in Table II. The dimensionful quantities are given in units of . Errors include both statistical plus fitting errors and errors coming from uncertaintiy in which we take to be %. Note that we have also gone to the renormalization group invariant bag parameter .
Set  

C1  1.430(21)  1.193(27)  1.199(29) 
C2  1.442(16)  1.248(35)  1.155(33) 
C3  1.382(21)  1.179(21)  1.172(23) 
C4  1.413(18)  1.263(22)  1.119(22) 
F1  1.353(17)  1.138(28)  1.189(26) 
F2  1.334(20)  1.193(27)  1.118(24) 
Iv Chiral and Continuum Extrapolations
The lattice data presented in Table II are for simulations with and quark masses larger than in the real world and need to be extrapolated to the physical point. Reaching this physical point also involves taking the lattice spacing limit. We use staggered chiral perturbation theory (SChPT) schpt1 (); schpt2 (); schpt3 () augmented by further general discretization correction terms to carry out the simultaneous chiral and continuum extrapolations. Continuum heavy meson chiral perturbation for and mixing was developed in detlin (); kamenik1 () including for the partially quenched case. These formulas were generalized recently to nexttoleading order SChPT by Bernard, Laiho and Van de Water schpt4 () and generously made available to us prior to publication. We use the following fit ansatz,
(18) 
stands for the chiral log contributions and includes the staggered light quark action specific taste breaking terms. The factor of comes about since was calculated for the square, namely for . We use the notation and for the sea and quark masses respectively, and (or ) for the valence quark masses. The second bracket parametrizes further discretization corrections that are expected to come in at and . We have also tried adding more analytic terms with higher powers of quark masses.
includes the coupling which has not been measured experimentally. However, based on Heavy Quark Effective Theory (HQET) arguments, is believed to be close to an analogous coupling in the meson system for which some experimental information is available. The latter coupling is estimated to be between gddpi (). As we discuss below, we have carried out two types of fits, one where we did a whole sequence of fits with varying between but where this coupling was kept fixed during each individual fit. In the second type of fit we let the coupling float and be one of the fit parameters. Both types of fits favored smaller values with , however as long as fit results were quite insensitive to its exact value.
For the ratio we use,
(19) 
Here we have split up,
(20) 
and then let become one of the fit parameters. In Fig.3 we show a simultaneous fit to the six entries in the last column of Table II. The green and blue curves are the curves from this fit appropriate to the coarse and fine lattice data points respectively and the red curve is the “continuum” curve obtained by retaining the fitted values for , and and turning off the and correction terms plus the taste breaking contributions inside and . One sees that within our statistical and fitting errors of %, there is consistency between the three curves. In other words, we see almost no statistically significant lattice spacing dependence in this ratio. At the physical point the difference between the green and blue curves is 1.8%, which reduces to 1.3% if the green curve is adjusted and corrected for having a sea quark mass on the coarse lattices that is about 20% too large. One might be surprised that the magenta curve lies below the blue curve. This comes about because the various discretization effects inside and in the & terms can have different signs and come in with different relative weights between the coarse and fine lattices. All these effects come in at the % or less level, and are hence too small to allow us to disentangle one from the other in a meaningful way. The fit shown in Fig.3 has bayes2 () and gives .
Fig.4 shows chiral & continuum extrapolation curves for using the fit ansatz of eq.(IV). Again the green and blue full curves are the fit curves for the coarse and fine lattice data respectively, and the dotted lines show the error bands around these central curves. Turning off the and contributions and the taste breaking terms inside leads to the red curve which can be followed down to the physical point. In contrast to the situation for the ratio , here, with , one finds a noticeable shift between the coarse and fine lattice points. The difference between the green and blue curves is a 5.5% effect. Going from the fine (blue) curve to the red continuum extrapolated curve is a 4% shift, which is also the size of the chiral & continuum extrapolation error at the physical point. The fit in Fig.4 has .
Finally, in Fig.5 we show results for , where for the simultaneous fit to all the data points. Here the difference between the green and blue curves is a 6% effect and between the blue and red curve a 5.7% effect. These shifts are slightly larger than but similar to those for in Fig.4. In both cases the discretization effects we are seeing in are larger than the naive expectation of a leading correction of which would be % or % on the coarse or fine lattices respectively. It was hence very important to have simulations results at more than one lattice spacing and carry out an explicit continuum extrapolation. Fortunately, for the important ratio these discretization corrections cancel out to a large extent, as expected and as we have already verified in Fig.3.
V Main Results and Error Budget
Table III gives our error budget for the main uncertainties in the three quantities, , and .
source of error  

stat. + chiral extrap.  2.3  4.1  2.0 
residual extrap.  3.0  2.0  0.3 
uncertainty  
uncertainty  2.3  2.3  — 
uncertainty  1.0  1.0  1.0 
and tuning  1.5  1.0  1.0 
operator matching  4.0  4.0  0.7 
relativistic corr.  2.5  2.5  0.4 
Total  6.7  7.1  2.6 
We explain each entry in Table III in turn.

statistics and chiral extrapolations: These are the errors shown on the “physical points” in Figs.3, 4 and 5 and are outputs from our chiral & continuum extrapolation fits.

residual extrapolation error: It is necessary to list this error separately since the degree to which the red curves in the above figures actually correspond to the true continuum limit depends on how well one has modelled discretization errors in our simulations. In other words one needs to assess the error in the fit ansatz for the continuum extrapolation (we assume the chiral extrapolation is handled sufficiently accurately by Staggered ChPT) and this turns out to be a nontrivial task.
On the one hand the data appears to be consistent with the fit ansätze of eqs.(IV) and (IV). We have tried adding further terms and found that fit results shifted by an amount less than, and in most cases much less than, the “statistical + chiral extrapolation” errors. Fig.6 shows results for the chirally and continuum extrapolated versus the number of parameters in the fit ansatz. For one has discretization corrections only through the term and one sees that a good fit to the data points at the two lattice spacings cannot be obtained. However, once , good fits are achieved and results and their errors are stable with respect to changes in .
On the other hand we know that due to our use of the NRQCD action to describe the heavy quark, coefficients such as are in general complicated functions of (although does include a constant piece coming from the light quarks).
We have approached this complicated situation in the following way. We interpret the red “continuum” curves in Figs.3, 4 and 5 as the curves one would get after taking care of all discretization errors coming from the light quark and the glue sectors. Then under “residual extrapolation uncertainty” one would include errors coming from the heavy quark action. The leading such error in our calculations is of multiplied by some function of which we initially take to be of leading to an additional uncertainty of % in . This would correspond to a standard power counting estimate of discretization errors where one takes coefficients of order one in higher order corrections. We have opted to be slightly more conservative in our power counting assessment and apply a factor of 1.5 rather than 1.0 for those cases where the “physical” (the magenta) points deviate by more than one from the fine lattice (blue) curve. By “” we mean here the “statistical + chiral extrapolation” errors. For we have multiplied the error for by a factor of .
Since the second row in Table III gives an assessment of the uncertainty in our extrapolation rather than an estimate of the full discretization error, we believe the procedure outlined here to fix it is a reasonably conservative one. 
uncertainty: follows from the 1.5% error in current determinations of the physical value for .

uncertainty in : we carried out fixed coupling chiral fits for the range and looked at the spread in the results at the physical point. For couplings larger than 0.6, starts to deteriorate.

tuning of strange and bottom quark masses: The largest mistuning, which occurs in the sea strange quark mass on the coarse lattices, has been corrected for when calculating fit results at the physical point and residual effects have been estimated by varying this adjusted value for . Errors due to uncertainty in the valence quark mass have been assessed by comparing as one goes from valence quark mass down to and errors coming from mistuning of have been estimated from the dependence of decay constants studied in fbsprl ().

operator matching and relativistic corrections: These two sources of error are intimately intertwined and again how to separate the two is not clear cut. As indicated in eq.(II), our matching for has been carried out up to correction of and . In Table III we have listed the first correction under “operator matching” and the latter correction under “relativistic corrections”. And again the errors for are reduced by a factor of 1/6 relative to those for the two non ratio quantities.
Using central values coming from the physical (red) points in the figures and the errors summarized in Table III, we can now present our main results.
(21) 
and using upsilon (),
(22) 
(23) 
where the first error comes from statistics + chiral extrapolation and the second is the sum of all other systematic errors added in quadrature. From the individual , q=s or d, one obtains a ratio of 1.231(58)(21) which is consistent with (21) however with larger errors. The result for in eq.(22) is consistent with but more accurate than our previously published value of MeV bsmixing ().
Vi Updates on , and and Estimates of Bag Parameters
The numerical simulations of twopoint and threepoint functions, such as in eqns.(III) and (III), that enabled us to extract the Bmixing parameters of the previous section also provide information necessary to determine and meson decay constants and . Decay constants are defined through the matrix element of the heavylight axial vector current between the meson state and the hadronic vacuum. Using the temporal component and working in the heavy meson rest frame one has,
(24) 
Just as with the fourfermion operators of section II, matching is required between the heavylight current in continuum QCD and currents made out of quark fields of the effective lattice theory. This matching has been carried out at the oneloop order for NRQCD/AsqTad currents in pert1 () based on formalism developed in pert2 ().
(25)  
The heavylight currents in the effective theory are defined as,
(26)  
(27)  
(28) 
with and
(29) 
The matching coefficients and are given in pert1 (). Note that the matching for the heavylight current includes contributions at and hence is more accurate than the matching in (II) for the fourfermion operator.
Set  

C1  1.261(12)  1.085(14)  1.162(14) 
C2  1.246(11)  1.073(14)  1.162(12) 
C3  1.236(12)  1.071(14)  1.155(14) 
C4  1.248(16)  1.128(17)  1.107(20) 
F1  1.175(13)  0.990(22)  1.188(20) 
F2  1.180(13)  1.047(16)  1.120(11) 
We have evaluated the twopoint functions,