The semileptonic form factor from threeflavor lattice QCD: A modelindependent determination of
Abstract
We calculate the form factor for meson semileptonic decay in unquenched lattice QCD with 2+1 flavors of light sea quarks. We use Asqtadimproved 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 heavylight 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, modelindependent manner. We find a value of .
pacs:
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 heavytolight quarkflavor 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 CabibboKobayashiMaskawa (CKM) matrix element with all sources of systematic uncertainty under control. In the Standard Model, the differential decay rate for this process is
(1) 
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 matterantimatter asymmetry of the universe. Thus we know that it is incomplete, and expect new physics to affect the quarkflavor 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 treelevel 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 heavyquark flavor sector such as asymmetries in Lin et al. (2008), constraints on from neutral meson mixing and 1loop penguininduced 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 treelevel Standard Model process at the few percentlevel. 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 lightlight and heavylight 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 “fourthroot” procedure for removing the undesired fourfold 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 unitarityviolating lattice artifacts in the pseudoGoldstone boson sector can be described and hence removed using rooted Staggered Chiral Perturbation Theory (rSPT), the lowenergy 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 phenomenologicallyimportant 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 :
(2)  
(3) 
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 ElKhadra 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 modelindependent unquenched lattice QCD calculation of the semileptonic form factor and . Our work builds upon the previous FermilabMILC 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 modeldependent 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 heavylight 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 12bin BABAR experimental data Aubert et al. (2007b) together to the modelindependent “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, modelindependent way. Others have also fit lattice and experimental results together using different, equallyvalid, 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, lightquark, and heavyquark 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 heavylight current to the continuum. Next we describe our analysis for determining the form factors in Sec. III. This is a threestep procedure. We first fit pion and meson 2point correlation functions to extract the meson masses. We then fit the 3point function, using the masses and amplitudes from the 2point 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 modelindependent 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 heavylight 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 2point and 3point 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  
For bottom quarks in 2point and 3point correlation functions we use the SheikholeslamiWohlert (SW) “clover” action Sheikholeslami and Wohlert (1985) with the Fermilab interpretation via HQET ElKhadra et al. (1997); Kronfeld (2000), which is wellsuited for heavy quarks, even when . Because the spinflavor symmetry of heavy quark systems is respected by the lattice regulator, the expansion in of the heavyquark lattice action has the same form as the expansion of the heavyquark part of the QCD action. Discretization effects in the lattice heavyquark action are therefore parameterized orderbyorder in the heavyquark expansion by deviations of effective operator coefficients from their values in continuum QCD. Thus, in principle, the lattice heavyquark 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 nexttoleading order, , in the heavyquark expansion.
The SW action includes a dimensionfive interaction with a coupling that must be adjusted to normalize the heavy quark’s chromomagnetic moment correctly ElKhadra et al. (1997). In our calculation we set the value of , as suggested by tadpoleimproved, treelevel perturbation theory Lepage and Mackenzie (1993). We determine the value of either from the plaquette ( fm) or from the Landau link ( fm). The tadpoleimproved bare quark mass for SW quarks is given by
(4) 
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 spinaveraged 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 ElKhadra et al. (1997). We then use the tuned value of for the formfactor 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 flavorchanging vector current to the same order in the heavyquark expansion. We remove errors of in the vector current by rotating the heavyquark field used in the matrix element calculation as
(5) 
where is the symmetric, nearestneighbor, covariant difference operator. We set to its tadpoleimproved treelevel value of ElKhadra et al. (1997)
(6) 
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 scalesetting 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 Heavylight meson correlation functions
The semileptonic form factors parameterize the hadronic matrix element of the quark flavorchanging vector current :
(7) 
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 ElKhadra et al. (2001)
(8) 
where is the fourvelocity 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:
(9)  
(10) 
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:
(11)  
(12) 
These relations automatically satisfy the kinematic constraint .
The 2point and 3point correlation functions needed to extract the lattice matrix element for decay are
(13)  
(14)  
(15) 
where and are interpolating operators for the meson and pion, respectively, and is the heavylight vector current on the lattice.
In practice, to construct the heavylight bilinears we must combine a staggered light quark, which is a 1component spinor, with a 4component Wilsontype 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 mixedaction interpolating operator:
(16) 
where are spin indices and . The fields and are the 4component clover quark field and 1component 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); KlubergStern et al. (1983). Because the heavy quark field is slowly varying over a hypercube, it does not affect the construction of Refs. Gliozzi (1982); KlubergStern et al. (1983).
For the pion we use the local pseudoscalar interpolating operator,
(17) 
where . We take the vector current to be
(18) 
where is the rotated heavyquark 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 4component heavy quarks instead of 2component nonrelativistic 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 2point function and the 3point 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 2point function and the 3point function :
(19) 
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 2point and 3point 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 3point 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 2point effective mass is still wellbehaved; 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 3point 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 sourcesink separation is too small then the entire time range is contaminated by excited states, but if the sourcesink 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 2point and 3point 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 Heavylight current renormalization
In order to recover the desired continuum matrix element, the lattice amplitude must be multiplied by the appropriate renormalization factor :
(20) 
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
(21)  
(22) 
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 ElKhadra et al. (2001). We first rewrite as
(23) 
The flavorconserving renormalization factors and account for most of the value of Harada et al. (2002a). They can be determined from standard heavylight meson charge normalization conditions:
(24)  
(25) 
where the lightlight and heavyheavy lattice vector currents are given by
(26)  
(27) 
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 3point correlation function. We eliminate contamination from staggered oscillating states in the determination of by using a clover strange quark for the spectator in the 3point 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 1loop 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 ElKhadra et al. (). The corrections to can be expressed as a perturbative series expansion in powers of the strong coupling:
(28) 
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 1loop coefficient, , and higher moments are calculated using automated perturbation theory and numerical integration as described in Refs. Lüscher and Weisz (1986); ElKhadra 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 threestep 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 2point 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 3point correlation function to the 2point 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 nexttoleading order (NLO) HMSPT expressions extended with nexttonexttoleading 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 Twopoint correlator fits
The pion and meson 2point correlators obey the following functional forms:
(29)  
(30) 
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 oppositeparity 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 2point 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 oppositeparity 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 sameparity 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.






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.
The pion energy that is extracted from fitting the 2point 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.^{1}^{1}1As 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 3point to 2point correlators.


iii.2 Threepoint correlator fits
The 3point correlator obeys the following functional form:
(31) 
where
(32) 
Writing out the first four terms of makes the behavior of the 3point correlator as a function of both and more transparent:
(33)  
As in the case of the pion and meson propagators, the leading contributions from the oppositeparity 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 oppositeparity excited states (the term), however, only changes sign when the sourcesink 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 groundstate amplitude of the 3point function as follows:
(34)  
(35) 
where, as before, and . The pion and meson energies and amplitudes are known from the 2point fits described in the previous subsection. Thus, the goal is to determine the 3point amplitude for along both the spatial and temporal directions.
In principle, the easiest way to determine the coefficient is to divide the 3point function by the appropriate 2point 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 excitedstate 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 2point 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 :
(36)  
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 masssplitting 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 oppositeparity 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:
(37)  
At our various light quark masses and lattice spacings the masssplittings 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 3point correlation function, we wish to reduce both the oscillating contributions and the less visible nonoscillating contributions that arise from the crossterm between the lowestlying pion and meson oppositeparity states. If these contributions are reduced sufficiently, we can safely neglect all of them when extracting the groundstate amplitude . We therefore construct a slightly more sophisticated average which combines the correlator both at consecutive time slices ( and ) and at consecutive sourcesink separations ( and ):
(38)  
This average reduces the unwanted parity states’ contamination significantly. It eliminates both the leading and subleading contributions to the oscillating term, the two lowestorder contributions to the oscillating term, and the contributions to the nonoscillating term. The size of the remaining term is a factor of 7–20 times smaller than in the unsmeared 3point correlator.
We can now safely ignore contamination from oppositeparity states and determine the lattice form factors in a simple manner. We construct the following ratio of the smeared correlators:
(39) 
The lattice form factors are then:
(40)  
(41) 
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 crosschecked 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.


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 nonzero 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 heavylight 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 leadingorder continuum HMPT expressions for and take the following simple forms:
(42)  
(43) 
where ,