The B\to\pi\ell\nu semileptonic form factor from three-flavor lattice QCD: A model-independent determination of |V_{ub}|

The semileptonic form factor from three-flavor lattice QCD: A model-independent determination of

Jon A. Bailey Fermi National Accelerator Laboratory, Batavia, Illinois, USA    C. Bernard Department of Physics, Washington University, St. Louis, Missouri, USA    C. DeTar Physics Department, University of Utah, Salt Lake City, Utah, USA    M. Di Pierro School of Computer Sci., Telecom. and Info. Systems, DePaul University, Chicago, Illinois, USA    A.X. El-Khadra Physics Department, University of Illinois, Urbana, Illinois, USA    R.T. Evans Physics Department, University of Illinois, Urbana, Illinois, USA    E. D. Freeland Liberal Arts Department, The School of the Art Institute of Chicago, Chicago, Illinois, USA    E. Gamiz Physics Department, University of Illinois, Urbana, Illinois, USA    Steven Gottlieb Department of Physics, Indiana University, Bloomington, Indiana, USA    U.M. Heller American Physical Society, One Research Road, Ridge, New York, USA    J.E. Hetrick Physics Department, University of the Pacific, Stockton, California, USA    A.S. Kronfeld Fermi National Accelerator Laboratory, Batavia, Illinois, USA    J. Laiho Department of Physics, Washington University, St. Louis, Missouri, USA    L. Levkova Physics Department, University of Utah, Salt Lake City, Utah, USA    P.B. Mackenzie Fermi National Accelerator Laboratory, Batavia, Illinois, USA    M. Okamoto Fermi National Accelerator Laboratory, Batavia, Illinois, USA    J. N. Simone Fermi National Accelerator Laboratory, Batavia, Illinois, USA    R. Sugar Department of Physics, University of California, Santa Barbara, California, USA    D. Toussaint Department of Physics, University of Arizona, Tucson, Arizona, USA    R.S. Van de Water Fermi National Accelerator Laboratory, Batavia, Illinois, USA
July 16, 2019

We calculate the form factor for -meson semileptonic decay in unquenched lattice QCD with 2+1 flavors of light sea quarks. We use Asqtad-improved staggered light quarks and a Fermilab bottom quark on gauge configurations generated by the MILC Collaboration. We simulate with several light quark masses and at two lattice spacings, and extrapolate to the physical quark mass and continuum limit using heavy-light meson staggered chiral perturbation theory. We then fit the lattice result for simultaneously with that measured by the BABAR experiment using a parameterization of the form factor shape in which relies only on analyticity and unitarity in order to determine the CKM matrix element . This approach reduces the total uncertainty in by combining the lattice and experimental information in an optimal, model-independent manner. We find a value of .

12.15.Hh, 12.38.Gc, 13.20.He

Fermilab Lattice and MILC Collaborations

I Introduction

The semileptonic decay is a sensitive probe of the heavy-to-light quark-flavor changing transition. When combined with an experimental measurement of the differential decay rate, a precise QCD determination of the form factor allows a clean determination of the Cabibbo-Kobayashi-Maskawa (CKM) matrix element with all sources of systematic uncertainty under control. In the Standard Model, the differential decay rate for this process is


where is the momentum transferred from the -meson to the outgoing leptons. The form factor, , parameterizes the hadronic contribution to the weak decay, and must be calculated nonperturbatively from first principles using lattice QCD.

A precise knowledge of CKM matrix elements such as is important not only because they are fundamental parameters of the Standard Model, but because inconsistencies between independent determinations of the CKM matrix elements and -violating phase would provide evidence for new physics. Although the Standard Model has been amazingly successful in describing the outcome of most particle physics experiments to date, it cannot account for gravity, dark matter and dark energy, or the large matter-antimatter asymmetry of the universe. Thus we know that it is incomplete, and expect new physics to affect the quark-flavor sector to some degree, although we do not know a priori what experimental and theoretical precision will be needed to observe it.

The determination of from semileptonic decay relies upon the assumption that, because the leading Standard Model decay amplitude is mediated by tree-level -boson exchange, it will not be significantly affected by new physics at the current level of achievable precision. Recently, however, hints of new physics have appeared in various regions of the heavy-quark flavor sector such as -asymmetries in  Lin et al. (2008), constraints on from neutral meson mixing and 1-loop penguin-induced decays Lunghi and Soni (2008), and the phase of the -mixing amplitude Aaltonen et al. (2008); Abazov et al. (2008); Bona et al. (2008). The unexpected inconsistency most relevant to our new lattice QCD calculation of the form factor and is the current “ puzzle” Dobrescu and Kronfeld (2008). The HPQCD Collaboration’s lattice QCD calculation of the -meson leptonic decay constant  Follana et al. (2008) disagrees with the latest results from the Belle, BABAR, and CLEO experiments Abe et al. (2007); Aubert et al. (2007a); Pedlar et al. (2007); Ecklund et al. (2008); Stone (2008) at the 3- level, although HPQCD’s determinations of the masses and and the decay constants , , and all agree quite well with experimental measurements Amsler et al. (2008); Eisenstein et al. (2008). Furthermore, because the significance of the discrepancy is dominated by the statistical experimental uncertainties, it cannot easily be explained by an underestimate of the theoretical uncertainties. Additional lattice QCD calculations of are needed to either confirm or reduce the inconsistency. If the disagreement holds up, however, it is evidence for a large new physics contribution to a tree-level Standard Model process at the few percent-level. Therefore, although semileptonic decay provides a theoretically clean method for determining within the framework of the Standard Model, we should keep in mind that new physics could appear in transitions.

Understanding and controlling all sources of systematic uncertainty in lattice QCD calculations of hadronic weak matrix elements is essential in order to allow accurate determinations of Standard Model parameters and reliable searches for new physics. The hadronic amplitudes for , in particular, can be calculated accurately using current lattice QCD methods because the decay process is “gold plated”, i.e., there is only a single stable hadron in both the initial and final states. Lattice calculations with staggered quarks allow for realistic QCD simulations with dynamical quarks as light as , multiple lattice spacings, large physical volumes, and high statistics. The resulting simulations of many light-light and heavy-light meson quantities with dynamical staggered quarks are in excellent numerical agreement with experimental results Davies et al. (2004). This includes both postdictions, such as the pion decay constant Aubin et al. (2004a), and predictions, as in the case of the meson mass Allison et al. (2005). Such successes show that the systematic uncertainties in these lattice QCD calculations are under control, and give confidence that additional calculations using the same methods are reliable.

The publicly available MILC gauge configurations with three flavors of improved staggered quarks Bernard et al. (2001) that have enabled these precise lattice calculations make use of the “fourth-root” procedure for removing the undesired four-fold degeneracy of staggered lattice fermions. Although this procedure has not been rigorously proven correct, Shamir uses plausible assumptions to argue that the continuum limit of the rooted theory is in the same universality class as QCD Shamir (2005, 2007). The rooting procedure leads to violations of unitarity that vanish in the continuum limit; both theoretical arguments Bernard (2006); Bernard et al. (2008a) and numerical simulations Prelovsek (2006); Bernard et al. (2007a); Aubin et al. (2008), however, show that the unitarity-violating lattice artifacts in the pseudo-Goldstone boson sector can be described and hence removed using rooted Staggered Chiral Perturbation Theory (rSPT), the low-energy effective description of the rooted staggered lattice theory Lee and Sharpe (1999); Aubin and Bernard (2003); Sharpe and Van de Water (2005). Given the wealth of numerical and analytical evidence supporting the validity of the rooting procedure, most of which is reviewed in Refs. Dürr (2006); Sharpe (2006); Kronfeld (2007), we work under the plausible assumption that the continuum limit of the rooted staggered theory is QCD. We note, however, that it is important to have crosschecks of lattice calculations of phenomenologically-important quantities using a variety of fermion formulations, since they all have different sources of systematic uncertainty.

Both existing unquenched lattice calculations of the form factor use the MILC configurations. When combined with the Heavy Flavor Averaging Group’s latest determination of the experimental decay rate from ICHEP 2008 Di Lodovico (), they yield the following values for :


where the errors are experimental and theoretical, respectively. Both analyses primarily rely upon data generated at a “coarse” lattice spacing of fm, and use a smaller amount of “fine” data at fm to check the estimate of discretization errors. Neither is therefore able to extrapolate the form factor to the continuum (). The most significant difference in the two calculations is their use of different lattice formulations for the bottom quarks. The HPQCD Collaboration Dalgic et al. (2006) uses nonrelativistic (NRQCD) heavy quarks Lepage et al. (1992), whereas we use relativistic clover quarks with the Fermilab interpretation El-Khadra et al. (1997) via heavy quark effective theory (HQET) Kronfeld (2000); Harada et al. (2002a, b). Both methods work quite well for heavy bottom quarks. The Fermilab treatment, however, has the advantage that it can also be applied to charm quarks; we can therefore use the same method for other semileptonic form factors such as , , and  Aubin et al. (2005); Bernard et al. (2008b). The two unquenched lattice calculations of the form factor, which have largely independent sources of systematic uncertainty, nevertheless lead to consistent values of with similar total errors of .

In this paper we present a new model-independent unquenched lattice QCD calculation of the semileptonic form factor and . Our work builds upon the previous Fermilab-MILC calculation and improves upon it in several ways. We now include data on both the coarse and fine MILC lattices, and can therefore take the limit of our data which is generated at nonzero lattice spacing. We also have additional statistics on a subset of the coarse ensembles. The most important improvements, however, are in the analysis procedures.

We have removed all model-dependent assumptions about the shape in of the form factor from the current analysis. Our result is therefore theoretically cleaner and more reliable than those of previous lattice QCD calculations. The first refinement over previous unquenched lattice form factor calculations is in the treatment of the chiral and continuum extrapolations. We simultaneously extrapolate to physical quark masses and zero lattice spacing and interpolate in the momentum transfer by performing a single fit to our entire data set (all values of , , and ) guided by functional forms derived in heavy-light meson staggered chiral perturbation theory (HMSPT) Aubin and Bernard (2007). We thereby extract the physical form factor in a controlled manner without introducing a particular ansatz for the form factor’s dependence. The second refinement over previous unquenched lattice form factor calculations is in the combination of the lattice form factor result and experimental data for the decay rate to determine the CKM matrix element . We fit our lattice numerical Monte Carlo data and the 12-bin BABAR experimental data Aubert et al. (2007b) together to the model-independent “-expansion” of the form factor given in Ref. Arnesen et al. (2005), in which the form factor is described by a power series in a small quantity with the sum of the squares of the series coefficients bounded by unitarity constraints. We leave the relative normalization factor, , as a free parameter to be determined by the fit, thereby extracting in an optimal, model-independent way. Others have also fit lattice and experimental results together using different, equally-valid, parameterizations Flynn and Nieves (2007); Bourrely et al. (2008). This work, however, is the first to use the full correlation matrices, derived directly from the data, for both the lattice calculation and experimental measurement.

This paper is organized as follows. In Sec. II we describe the details of our numerical simulations. We discuss the gluon, light-quark, and heavy-quark lattice actions, and present the parameters used, such as the quark masses and lattice spacings. We then define the matrix elements needed to calculate the semileptonic form factors and discuss the method for matching the lattice heavy-light current to the continuum. Next we describe our analysis for determining the form factors in Sec. III. This is a three-step procedure. We first fit pion and -meson 2-point correlation functions to extract the meson masses. We then fit the 3-point function, using the masses and amplitudes from the 2-point fits as input, to extract the lattice form factors at each value of the light quark mass and lattice spacing. Finally, we extrapolate the results at unphysical quark masses and nonzero lattice spacing to the physical light quark masses and zero lattice spacing using HMSPT. In Sec. IV we estimate the contributions of the various systematic uncertainties to the form factors, discussing each item in the error budget separately. We then present the final result for with a detailed breakdown of the error by source in each bin. We combine our result for the form factor with experimental data from the BABAR Collaboration to determine the CKM matrix element in Sec. V. We also define the model-independent description of the form factor shape that we use in the fit and discuss alternative parameterizations of the form factor. Finally, in Sec. VI we compare our results with those of previous unquenched lattice calculations. We also compare our determination of with inclusive determinations and to the preferred values from the global CKM unitarity triangle analysis. We conclude by discussing the prospects for improvements in our calculation and its impact on searches for new physics in the quark flavor sector.

Ii Lattice Calculation

In this section we describe the details of our numerical lattice simulations. We first present the actions and parameters used for the light (up, down, strange) and heavy (bottom) quarks in Sec. II.1. We then define the procedure for constructing lattice correlation functions with both staggered light and Wilson heavy quarks in Sec. II.2. Finally, in Sec. II.3, we show how to match the lattice heavy-light vector currents to the continuum with a mostly nonperturbative method, so that lattice perturbation theory is only needed to estimate a small correction.

ii.1 Actions and Parameters

We use the ensembles of lattice gauge fields generated by the MILC Collaboration and described in Ref. Bernard et al. (2001) at two lattice spacings, and  fm, in our numerical lattice simulations of the semileptonic form factor. The ensembles include the effects of three dynamical staggered quarks — two degenerate light quarks with masses ranging from and one heavier quark tuned to within 10–30% of the physical strange quark mass. This allows us to perform a controlled extrapolation to both the continuum and the physical average - quark mass. The physical lattice volumes are all sufficiently large () to ensure that effects due to the finite spatial extent remain small.

For each independent ensemble we compute the light valence quark in the 2-point and 3-point correlation functions only at the same mass, , as the light quark in the sea sector. Thus all of our simulations are at the “full QCD” point. Note, however, that we still have many correlated data points on each ensemble because of the multiple pion energies. Table 1 shows the combinations of lattice spacings, lattice volumes, and quark masses used in our calculation.

(fm) (fm)   Volume   Configs.
2.4 4.1 557  1.476 0.0923  0.09474
2.4 5.8 518  1.476 0.0923  0.09469
2.9 3.8 529  1.72 0.086  0.09372
2.4 3.8 836  1.72 0.086  0.09372
2.4 4.5 592  1.72 0.086  0.09384
2.4 6.2 460  1.72 0.086  0.09368
Table 1: Lattice simulation parameters. The columns from left to right are the approximate lattice spacing in fm, the bare light quark masses , the linear spatial dimension of the lattice in fm, the dimensionless factor (corresponding to the taste-pseudoscalar pion composed of light sea quarks), the dimensions of the lattice in lattice units, the number of configurations used for this analysis, the clover term and bare value used to generate the bottom quark, and the improvement coefficient used to rotate the bottom quark field in the vector current.

For bottom quarks in 2-point and 3-point correlation functions we use the Sheikholeslami-Wohlert (SW) “clover” action Sheikholeslami and Wohlert (1985) with the Fermilab interpretation via HQET El-Khadra et al. (1997); Kronfeld (2000), which is well-suited for heavy quarks, even when . Because the spin-flavor symmetry of heavy quark systems is respected by the lattice regulator, the expansion in of the heavy-quark lattice action has the same form as the expansion of the heavy-quark part of the QCD action. Discretization effects in the lattice heavy-quark action are therefore parameterized order-by-order in the heavy-quark expansion by deviations of effective operator coefficients from their values in continuum QCD. Thus, in principle, the lattice heavy-quark action can be improved to arbitrarily high orders in HQET by tuning a sufficiently large number of parameters in the lattice action. In practice, we tune the hopping parameter, , and the clover coefficient, , of the SW action, to remove discretization effects through next-to-leading order, , in the heavy-quark expansion.

The SW action includes a dimension-five interaction with a coupling that must be adjusted to normalize the heavy quark’s chromomagnetic moment correctly El-Khadra et al. (1997). In our calculation we set the value of , as suggested by tadpole-improved, tree-level perturbation theory Lepage and Mackenzie (1993). We determine the value of either from the plaquette ( fm) or from the Landau link ( fm). The tadpole-improved bare quark mass for SW quarks is given by


such that tuning the parameter to the critical quark hopping parameter leads to a massless pion. Before generating the correlation functions needed for the form factor, we compute the spin-averaged kinetic mass on a subset of the available ensembles in order to tune the bare value for bottom (and hence the corresponding bare quark mass) to its physical value El-Khadra et al. (1997). We then use the tuned value of for the form-factor production runs. Table 1 shows the values of the clover coefficient and tuned used in our calculation.

In order to take advantage of the improved action in the calculation of the form factor, we must also improve the flavor-changing vector current to the same order in the heavy-quark expansion. We remove errors of in the vector current by rotating the heavy-quark field used in the matrix element calculation as


where is the symmetric, nearest-neighbor, covariant difference operator. We set to its tadpole-improved tree-level value of El-Khadra et al. (1997)


The values of the rotation parameter used in our calculation are given in Table 1.

In order to convert dimensionful quantities determined in our lattice simulations into physical units, we need to know the value of the lattice spacing, , which we find by computing a physical quantity that can be compared directly with experiment. We first determine the relative lattice scale by calculating the ratio on each ensemble, where is related to the force between static quarks,  Sommer (1994). These estimates are then smoothed by fitting to a smooth function of the gauge coupling and quark masses. This scale-setting method has the advantage that the ratio can be determined precisely from a fit to the static quark potential Bernard et al. (2000); Aubin et al. (2004b). We convert all of our data from lattice spacing units into units before performing any chiral fits in order to account for slight differences in the value of the lattice spacing between ensembles. In this work we use the value of obtained by combining a recent lattice determination of  Bernard et al. (2007b) with the PDG value of MeV Yao et al. (2006) to convert lattice results from units to physical units.

ii.2 Heavy-light meson correlation functions

The semileptonic form factors parameterize the hadronic matrix element of the quark flavor-changing vector current :


where is the momentum transferred to the outgoing lepton pair. For calculations on the lattice and in HQET, it is more convenient to write the matrix element as El-Khadra et al. (2001)


where is the four-velocity of the -meson, is the component of the pion momentum orthogonal to , and is the energy of the pion in the -meson rest frame (). In this frame the form factors and are directly proportional to the hadronic matrix elements:


We therefore first calculate the hadronic matrix elements in Eqs. (9) and (10) in the rest frame of the -meson to obtain and , and then extract the standard form factors and using the following relations:


These relations automatically satisfy the kinematic constraint .

The 2-point and 3-point correlation functions needed to extract the lattice matrix element for decay are


where and are interpolating operators for the -meson and pion, respectively, and is the heavy-light vector current on the lattice.

In practice, to construct the heavy-light bilinears we must combine a staggered light quark, which is a 1-component spinor, with a 4-component Wilson-type bottom quark; we do so using the method established by Wingate et al. in Ref. Wingate et al. (2003). For the meson we use a mixed-action interpolating operator:


where are spin indices and . The fields and are the 4-component clover quark field and 1-component staggered field, respectively. Based on the transformation properties of under shifts by one lattice spacing, plays the role of a (fermionic) taste index Wingate et al. (2003); Kronfeld (2007). Once is summed over hypercubes in the correlation functions that we compute, also takes on the role of a taste degree of freedom, in the sense of Refs. Gliozzi (1982); Kluberg-Stern et al. (1983). Because the heavy quark field is slowly varying over a hypercube, it does not affect the construction of Refs. Gliozzi (1982); Kluberg-Stern et al. (1983).

For the pion we use the local pseudoscalar interpolating operator,


where . We take the vector current to be


where is the rotated heavy-quark field given in Eq. (5). When forming and , we sum over the taste index. This yields the same correlation functions, with respect to taste, as in Ref. Wingate et al. (2003). Our principal difference with Ref. Wingate et al. (2003) is to use 4-component heavy quarks instead of 2-component non-relativistic quarks, and to derive the correlators in the staggered formalism, without the introduction of naive fermions.

We work in the rest frame of the -meson, so only the pions carry momentum. We compute both the 2-point function and the 3-point function at discrete values of the momenta , and allowed by the finite spatial lattice volume. We use only data through momentum , however, because the statistical errors in the correlators increase significantly with momentum.

We use a local source for the pions throughout the calculation, while we smear the -meson wavefunction in both the 2-point function and the 3-point function :


where is the spatial smearing function. This reduces contamination from heavier excited states and allows a better determination of the desired ground state amplitude. In our study of choices for how to smear the -meson, we found that a wall source, , worked extremely well for suppressing excited state contamination, but at the cost of large statistical errors in the 2-point and 3-point correlation functions. In contrast, use of a 1S wavefunction, , optimized to have good overlap with the charmonium ground state led to smaller statistical errors at the cost of undesirably large excited state contributions to the 3-point function that would make it difficult to extract the ground state amplitude. In order to achieve a balance between small statistical errors and minimal excited state contamination, we tune the coefficient of the exponential in the 1S wavefunction to the smallest value (i.e., the widest smearing) for which the -meson 2-point effective mass is still well-behaved; we find a value of for the coarse ensembles. We note that our determination of the optimal -meson smearing function is consistent with the theoretical expectation that the -meson wavefunction should be wider than the charmonium wavefunction.

For the calculation of the 3-point function, we fix the location of the pion source at and the location of the -meson sink at , and vary the position of the operator over all times in between. If the source-sink separation is too small then the entire time range is contaminated by excited states, but if the source-sink separation is too large then the correlation function becomes extremely noisy. In practice, we set the sink time to on the coarse lattices; we have checked, however, that our results using this choice are consistent with those determined from using and . On the fine lattices we scale the source sink separation by the approximate ratio of the lattice spacings, , and use .

In order to minimize the statistical errors given the available number of configurations in each ensemble, we compute the necessary 2-point and 3-point correlations not only with a source time of , but also with source times of , and ( is the temporal extent of the lattice) and the sink time shifted accordingly. We then average the results from the four source times; this effectively increases our statistics by a factor of four.

ii.3 Heavy-light current renormalization

In order to recover the desired continuum matrix element, the lattice amplitude must be multiplied by the appropriate renormalization factor :


where and are the lattice and continuum vector currents, respectively. This removes the dominant discretization errors from the lattice current operator. In terms of the form factors, Eq. (20) can be rewritten as


where explicit expressions relating and to correlation functions are given in Eqs. (40) and (41).

In this work, we calculate via the mostly nonperturbative method used in the earlier quenched Fermilab calculation El-Khadra et al. (2001). We first rewrite as


The flavor-conserving renormalization factors and account for most of the value of  Harada et al. (2002a). They can be determined from standard heavy-light meson charge normalization conditions:


where the light-light and heavy-heavy lattice vector currents are given by


respectively. In order to reduce the statistical errors in , we compute the lattice matrix element using a clover charm quark as the spectator in the 3-point correlation function. We eliminate contamination from staggered oscillating states in the determination of by using a clover strange quark for the spectator in the 3-point correlation function . Once and have been determined nonperturbatively, the remaining correction factor in Eq. (23), , is expected to be close to unity because most of the radiative corrections, including contributions from tadpole graphs, cancel in the ratio Harada et al. (2002a). We therefore estimate from 1-loop lattice perturbation theory Lepage and Mackenzie (1993).

The matching factor has been calculated by a subset of the present authors, and a separate publication describing the details is in preparation El-Khadra et al. (). The corrections to can be expressed as a perturbative series expansion in powers of the strong coupling:


where is the renormalized coupling constant in the -scheme and is determined from the static quark potential with the same procedure as is used in Ref. Mason et al. (2005). The scale , which should be the size of a typical gluon loop momentum, is computed via an extension of the methods outlined by Brodsky, Lepage, and Mackenzie Lepage and Mackenzie (1993); Brodsky et al. (1983) and Hornbostel, Lepage, and Morningstar Hornbostel et al. (2003). The value of ranges from 2.0–4.5 GeV for the parameters used in our simulations. The 1-loop coefficient, , and higher moments are calculated using automated perturbation theory and numerical integration as described in Refs. Lüscher and Weisz (1986); El-Khadra et al. (2007). We find that the perturbative corrections to matrix elements of the temporal vector current, , are less than a percent, while the corrections to matrix elements of the spatial vector current, , are 3–4%.

Iii Analysis

In this section, we describe the three-step analysis procedure used to calculate the semileptonic form factor, . In the first subsection, Sec. III.1, we describe how we fit the pion and -meson 2-point correlation functions in order to determine the pion energies and -meson mass. We use both of these quantities in the later determination of the lattice form factors and . Next, in Sec. III.2, we construct a useful ratio of the 3-point correlation function to the 2-point functions. We then fit this ratio to a simple plateau ansatz to extract the desired form factors. Finally, in Sec. III.3, we extrapolate the form factors calculated at unphysically heavy quark masses and finite lattice spacing to the physical light quark masses and zero lattice spacing using next-to-leading order (NLO) HMSPT expressions extended with next-to-next-to-leading order (NNLO) analytic terms. (We perform a simultaneous extrapolation in and and interpolation in .) We then take the appropriate linear combination of and to determine the desired form factor, , with statistical errors.

iii.1 Two-point correlator fits

The pion and -meson 2-point correlators obey the following functional forms:


In the above expressions, terms with odd contain the prefactor . This leads to visible oscillations in time in the meson propagators; such behavior arises with staggered quarks because the parity operator is a composition of spatial inversion and translation through one timeslice Golterman and Smit (1984); Kilcup and Sharpe (1987). The contributions of the opposite-parity oscillating states are found to be significant throughout the entire time range and must therefore be included in fits to extract the desired ground state energy.

Because the statistical errors in the pion energies and -meson mass contribute very little to the total statistical error in the form factor, we use a simple procedure to fit the 2-point functions. Although this does not optimize the determinations of and , it is sufficient for the purpose of this analysis. We first select a fit range, , that allows a good correlated, unconstrained fit including only contributions from the ground state and its opposite-parity partner. We then reduce by one timeslice and redo the fit. If the correlated confidence level is too low (), we increase the number of states and try the fit again with the same time range. Otherwise, if the fit is good, we reduce by one more timeslice and repeat the fit. We repeat this procedure until we can no longer get a good fit without using a large number (greater than 4) of states. We note that, by including only as many states as the data can determine, we minimize the possibility of spurious solutions in which the fitter exchanges the ground state with one of the same-parity excited states. We have, however, checked that this method yields the same results within statistical errors as a constrained fit that includes up to three or four pairs of states.

Figure 1 shows examples of both vs. (left plot) and vs. (right plot) on the coarse ensemble, which has the largest light quark mass of the coarse ensembles. The masses are stable as is reduced, and the statistical errors in become smaller as additional timeslices are added to the fit. The statistical errors are determined by performing a separate fit to 500 bootstrap ensembles; each fit uses the full single elimination jackknife correlation matrix which is remade before every fit. The size of the statistical errors does not change when the number of bootstrap ensembles is increased by factors of two or four. We select the time range to use in the analysis based on several criteria: a good correlated confidence level, relatively symmetric upper and lower bootstrap error bars, no - or greater outliers in the bootstrap distribution, and no sign of excited state contamination. The red (filled) data points in Fig. 1 mark the chosen fit ranges for the ensemble in the example plots. Figures 2 and 3 show the corresponding pion and -meson correlator fits, respectively, which go through the data points (shown with jackknife errors) quite well.

Figure 1: Pion mass (left plot) and -meson mass (right plot) versus minimum timeslice in 2-point correlator fit. The red (filled) data points show the fit ranges selected for use in the form factor analysis.
Figure 2: Pion correlator fit corresponding to the red data point in the left-hand graph of Fig. 1. The left plot shows the fit (red line) to the zero-momentum pion propagator on a log scale, while the right plot shows the deviation of the fit from the data point for each timeslice. On both plots the dashed vertical line indicates . Single elimination jackknife statistical errors are shown.
Figure 3: -meson correlator fit corresponding to the red data point in the right-hand graph of Fig. 1. The left plot shows the fit (red line) to the -meson propagator on a log scale, while the right plot shows the deviation of the fit from the data point for each timeslice. On both plots the dashed vertical line indicates . Single elimination jackknife statistical errors are shown.

The gauge configurations have been recorded every six trajectories, and the remaining autocorrelations between consecutive configurations cannot be neglected. We address this by averaging a block of successive configurations together before calculating the correlation matrix and performing the fit. We determine the optimal block size by increasing the number of configurations in a block until the single elimination jackknife statistical error in the correlator data remains constant within errors. This is shown for a representative timeslice of the pion propagator on a coarse ensemble in Fig. 4. We find that it is necessary to use a block size of 5 on the coarse ensembles and 8 on the fine ensembles, and we use these values for the rest of the form factor analysis. We note that the size of the statistical errors that arises from blocking by 5 on the coarse ensemble is consistent with that estimated based on a calculation of the integrated autocorrelation time.

Figure 4: Single-elimination jackknife error versus block size in the zero-momentum pion propagator at . The statistical errors in the errors are calculated with an additional single elimination jackknife loop. The red line is an average of the errors for block sizes 5–8 and is only to make it easier to see that the statistical error plateaus after a block size of 5; it is not used in the form factor analysis.

The pion energy that is extracted from fitting the 2-point function, , should satisfy the dispersion relation in the continuum limit due to the restoration of rotational symmetry. Similarly, the pion amplitude, , should be independent of as . As shown in Fig. 5, our results are consistent with these continuum relations within statistical errors.111As this analysis was being completed we generated data with four times the statistics on the coarse ensemble. In order to make the comparison to the continuum expectation clearer, we use the higher statistics data in Fig. 5. We therefore replace the pion energy by when calculating the lattice form factors and in order to reduce the total statistical uncertainty. The pion amplitude drops out of the form factor calculation, however, because we take suitable ratios of 3-point to 2-point correlators.

Figure 5: Comparison of pion energy (left plot) and amplitude (right plot) with the prediction of the continuum dispersion relation. We also show a power-counting estimate for the size of momentum-dependent discretization errors, which are of , as dashed black lines.

iii.2 Three-point correlator fits

The 3-point correlator obeys the following functional form:




Writing out the first four terms of makes the behavior of the 3-point correlator as a function of both and more transparent:


As in the case of the pion and -meson propagators, the leading contributions from the opposite-parity excited states (the and terms) change sign when ; these produce visible oscillations in the correlation function along the time direction. The subleading contribution from the opposite-parity excited states (the term), however, only changes sign when the source-sink separation is varied, e.g., ; this contribution is not as clearly visible in the data as those that oscillate with the time slice .

The lattice form factors are related to the ground-state amplitude of the 3-point function as follows:


where, as before, and . The pion and -meson energies and amplitudes are known from the 2-point fits described in the previous subsection. Thus, the goal is to determine the 3-point amplitude for along both the spatial and temporal directions.

In principle, the easiest way to determine the coefficient is to divide the 3-point function by the appropriate 2-point functions and fit to a constant (plateau) ansatz in a region of time slices that are sufficiently far from both the pion and -meson sources, such that excited state contamination can be neglected. In practice, however, oscillating excited-state contributions are significant throughout the interval between the pion and -meson, so our raw correlator data cannot be fit to such a simple function. Therefore we construct an average correlator in which the oscillations are reduced before performing any fits. This method for determining the form factors requires knowledge of and ; we use the values determined in the 2-point fits described in the previous subsection and propagate the bootstrap uncertainties in order to properly account for correlations.

The final ratio of correlators used to determine entails several pieces. To begin consider the carefully constructed average of the value of the -meson propagator at time slice with that at :


where . By removing the leading exponential behavior from the correlator before taking the average we suppress the leading oscillating contribution by a factor of the mass-splitting while leaving the desired ground state amplitude unaffected. Note also that, while this procedure affects the size of the excited state amplitudes, it does not alter the functional form of the correlator, nor does it alter the energies in the exponentials. Therefore the average in Eq. (36) is equivalent to using a smeared source that has a smaller coupling to the opposite-parity excited states. This averaging procedure can be iterated in order to make the oscillating terms arbitrarily small. Empirically, we find that two iterations are sufficient for all of our numerical data:


At our various light quark masses and lattice spacings the mass-splittings lie in the range in lattice units; thus use of the iterated average in Eq. (37) reduces the leading oscillating state amplitude by a factor of 50–400 such that it can be safely neglected.

In the case of the 3-point correlation function, we wish to reduce both the oscillating contributions and the less visible non-oscillating contributions that arise from the cross-term between the lowest-lying pion and -meson opposite-parity states. If these contributions are reduced sufficiently, we can safely neglect all of them when extracting the ground-state amplitude . We therefore construct a slightly more sophisticated average which combines the correlator both at consecutive time slices ( and ) and at consecutive source-sink separations ( and ):


This average reduces the unwanted parity states’ contamination significantly. It eliminates both the leading and subleading contributions to the oscillating term, the two lowest-order contributions to the oscillating term, and the contributions to the non-oscillating term. The size of the remaining term is a factor of 7–20 times smaller than in the unsmeared 3-point correlator.

We can now safely ignore contamination from opposite-parity states and determine the lattice form factors in a simple manner. We construct the following ratio of the smeared correlators:


The lattice form factors are then:


We fit and as defined in Eqs. (40)–(41) to a plateau in the region where ordinary excited state contributions can be neglected. Figure 6 shows the determinations of (left plot) and (right plot) for all of the momenta that we use in the chiral extrapolation on the coarse ensemble with . In practice, we fit a range of four time slices, choosing the interval that results in the best correlated confidence level. We have cross-checked the determination of the form factors via Eqs. (40)–(41) against determinations of the form factor that explicitly include excited state dependence in the fit ansatz and find that the results agree within errors. Our preferred method, however, yields the smaller statistical uncertainty in the form factors.

Figure 6: Determination of the form factors (left plot) and (right plot) from plateau fits to the ratios defined in Eqs. (40) and (41). The statistical errors on the data points are from a single-elimination jackknife. The statistical errors in the plateau determination are from separate fits of 500 bootstrap ensembles.

iii.3 Continuum and chiral extrapolation

The quark masses in our numerical lattice simulations are heavier than the physical up and down quark masses. The effects of non-zero lattice spacings in Asqtad simulations are also too large to be neglected. In order to account for these facts, we calculate the desired hadronic matrix elements for multiple values of the light quark masses and lattice spacing, and then extrapolate to the physical quark masses and continuum using functional forms from heavy-light meson staggered chiral perturbation theory (HMSPT) Aubin and Bernard (2007). The HMSPT expressions are derived using the symmetries of the staggered lattice theory, and therefore contain the correct dependence of the form factors on the quark mass and lattice spacing. In the case of the form factors, the HMSPT expressions are also functions of the pion energy (recall that we work in the frame where the -meson is at rest).

HMSPT is a systematic expansion in inverse powers of the heavy quark mass. In the chiral and soft pion limits ( and , the leading-order continuum HMPT expressions for and take the following simple forms:


where ,