Neutral B Meson Mixing in Unquenched Lattice QCD

Neutral Meson Mixing in Unquenched Lattice QCD

Elvira Gámiz Department of Physics, University of Illinois, Urbana, IL 61801, USA    Christine T. H. Davies Department of Physics & Astronomy, University of Glasgow, Glasgow, G12 8QQ, UK    G. Peter Lepage Laboratory of Elementary Particle Physics, Cornell University, Ithaca, NY 14853, USA    Junko Shigemitsu Department of Physics, The Ohio State University, Columbus, OH 43210, USA    Matthew Wingate Department of Applied Mathematics and Theoretical Physics, University of Cambridge, Cambridge CB3 0WA, UK.

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 .

12.38.Gc, 13.20.Fc, 13.20.He

HPQCD 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 Cabibbo-Kobayashi-Maskawa (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 four-fermion 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
C2 2.635 0.007/0.050 0.007 834 4/2
C3 2.619 0.010/0.050 0.010 672 4/3
C4 2.651 0.020/0.050 0.020 459 4/3
F1 3.701 0.0062/0.031 0.0062 547 4/2
F2 3.721 0.0124/0.031 0.0124 534 4/2
Table 1: Details of configurations employed. is the number of time sources used per configuration and the number of smearings on the heavy propagator including the unsmeared local case. All quark masses are given in the MILC collaboration normalization convention with . Errors in are estimated to be at the 0.5% level.

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 four-fermion 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 four-fermion 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 four-fermion 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 four-fermion operators. More specifically, for calculations of in QCD one is interested in the operator with [V-A] x [V-A] structure,


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,


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 (),


where depends on the quark and the boson masses and , is a perturbative QCD correction factor and the Inami-Lim function. is the renormalization group invariant bag parameter and at two-loop 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,


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 four-fermion 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 tree-level and that matches onto (1) at lowest order in has the form,


As is well known, even at lowest order in 1/M there is a one-loop order mixing with another four-fermion operator,


This is true both in NRQCD and in HQET. If one introduces an effective theory field,


then and the QCD field are related by a Foldy-Wouthuysen-Tani (FWT) transformation. In particular,


where the acts to the left. The FWT transformation determines the tree-level 1/M corrections to the four-fermion operators in the effective theory. For they come in as,


Taking these corrections into account one can work through and finds the following matching relation,


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 three-point 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 three-point correlator must be fit to

This ansatz allows for non-oscillatory and oscillatory contributions to the correlator (in practice we have worked with ). Not all the amplitudes etc. are independent due to symmetries. Similarly two-point correlators are fit to,

The relation between the amplitudes or the and the matrix elements of the previous section can be identified as follows.


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,


which includes the matrix element that we are interested in. Similarly for the 2pt-functions one has,


Using one then has,


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 four-fermion operator at the three-point function level by forming the appropriate linear combinations of the ’s, and then carried out fits to extract for the full renormalized three-point 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.

Figure 1: Fit results for versus the number of exponentials for one of the coarse ensembles, Set C2 with .
Figure 2: Same as Fig.1 for a fine ensemble, Set F1 with .

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 three-point function we carry out simultaneous fits to an matrix of two-point correlators (eq.(III) with ) ( corresponds to local, to first Gaussian etc.) and to the renormalized three-point 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 .

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)
Table 2: Fit results for in units of and for the dimensionless ratio . Errors in the last column are statistical + fitting errors. Those in the second and third columns include additional errors coming from the 0.5% uncertainty in .

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 next-to-leading 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,


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,


Here we have split up,


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 .

Figure 3: Chiral and continuum extrapolation of the ratio . Errors shown are statistical plus fitting errors. The physical point is at , where is the valence strange quark mass.

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 .

Figure 4: Chiral and continuum extrapolation of . Errors on the data points are statistical plus fitting errors combined with uncertainty in . The dashed curves correspond to the error bands about the central green and blue full lines. The physical point is at , where is the valence strange quark mass.

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.

Figure 5: Same as Fig.4 for versus .

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 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
Table 3: Errors in % for , and .
Figure 6: at the physical point versus , the number of fit parameters. The numbers below the data points give . For we give results for two different types of term added, a term proportional to the quark mass squared (right data point) or a further discretization correction (left data point).

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.


and using upsilon (),


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 two-point and three-point functions, such as in eqns.(III) and (III), that enabled us to extract the B-mixing 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 heavy-light 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,


Just as with the four-fermion operators of section II, matching is required between the heavy-light current in continuum QCD and currents made out of quark fields of the effective lattice theory. This matching has been carried out at the one-loop order for NRQCD/AsqTad currents in pert1 () based on formalism developed in pert2 ().


The heavy-light currents in the effective theory are defined as,


with and


The matching coefficients and are given in pert1 (). Note that the matching for the heavy-light current includes contributions at and hence is more accurate than the matching in (II) for the four-fermion operator.

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)
Table 4: Fit results for in units of and for the ratio . Errors are as described in Table II.

We have evaluated the two-point functions,