and meson spectroscopy
Abstract
Results are presented for the lowlying spectrum of and mesons calculated in lattice QCD using 2+1 flavor CloverWilson configurations made available by the PACSCS collaboration. For the heavy quark, the Fermilab method is employed. The main focus is S and Pwave states of charmed and charmedstrange mesons, where previous lattice QCD results have been mostly from quenched calculations. In addition to the ground states, some excited states are extracted. To check the method, calculations of the charmonium spectrum are also carried out. For charmonium, the lowlying spectrum agrees favorably with experiment. For heavystrange and heavylight systems substantial differences in comparison to experiment values remain for some states.
pacs:
11.15.Ha, 12.38.GcI Introduction
The spectrum of charmedstrange mesons^{1}^{1}1Following Particle Data Group usage Nakamura et al. (2010) mesons are referred to as charmedstrange and only mesons as charmed. contains a number of wellestablished states Nakamura et al. (2010), notably the “Swave” states of quantum numbers (where is the spin and P is parity) and , the and the and the “Pwave” states with quantum numbers (), ( and ) and (). In the heavy quark limit, these states would form three massdegenerate multiplets and are characterized by the total angular momentum of the light quark. Within many quark models (see for example Godfrey and Isgur (1985)), the physical states corresponding to the doublet with , the and were expected to have considerable width and masses above the and thresholds. In experiment, both states are below threshold and narrow.
In addition there are a number of states which have been observed more recently. These are the (observed by both BaBar Aubert et al. (2006, 2009) and Belle Brodzicka et al. (2008)), the (observed by BaBar Aubert et al. (2006, 2009)), the ) (observed by BaBar Aubert et al. (2009)) and an unconfirmed state previously observed by SELEX Evdokimov et al. (2004), the . While the is commonly believed to have quantum numbers , there are several possibilities for the other states which are not ruled out by experiment. The has natural parity and is most often identified with a state Colangelo and De Fazio (2010); Ebert et al. (2010); Chen et al. (2009), while some still argue the possibility of a identification van Beveren and Rupp (2010). The has unnatural parity and is commonly interpreted as either a Colangelo and De Fazio (2010); Ebert et al. (2010); Sun and Liu (2009); Chen et al. (2009) or a state Colangelo and De Fazio (2010).
Until recently, the number of observed charmed meson states was more limited. Again the “Swave” states of quantum numbers () and () as well as the “Pwave” states with quantum numbers (), ( and ) and () are wellestablished Nakamura et al. (2010). In addition, the Particle Data Group Nakamura et al. (2010) lists the seen in Zdecays Abreu et al. (1998) which lacks confirmation. More recently BaBar observed several new charmed mesons del Amo Sanchez et al. (2010) whose quantum numbers are a matter of active discussion Sun et al. (2010); Li et al. (2011); Wang (2011); Zhong (2010); Chen et al. (2011).
Lattice QCD (LQCD) provides the possibility to elucidate the spectrum without resorting to model assumptions. To reach this goal, several systematic sources of uncertainty have to be controlled. Of particular importance is that full dynamical simulations are necessary with up and down quarks light enough so that an extrapolation to the physical point is not needed or can be done in a meaningful way. Also, there should be an extrapolation to the continuum limit. Such calculations have been achieved for the lightquark ground state meson and baryon spectrum (see, for example, Dürr et al. (2008)).
Dealing with heavy charm or bottom quarks on the lattice introduces complications of its own, and so far a similar precision for the spectrum of charmed and charmedstrange states has not been attained. Most previous lattice calculations have been within the quenched approximation or at unphysically large pion masses. Advances in lattice simulations enable us to revisit this problem, removing one deficiency of previous calculations by using dynamical gauge ensembles where pion masses are close to the physical pion mass.
Among the lowlying and mesons accessible to current lattice calculations, the and , whose properties can not be explained within most quark models, are of particular interest. Possible explanations for their unanticipated mass and width have been reviewed in Swanson (2006); Zhu (2008). In early LQCD calculations, the ground states in the and channels were often found to be quite a bit heavier than the experimental resonances (for a short summary of previous results please refer to Section IV). This has further sparked speculations about the nature of these states and provides motivation for this study.
The next section provides some details of our calculational setup. There is also a short description of the procedure used to tune the charm quark mass. To test our setup, the spectrum of charmonium states below multiparticle threshold was calculated and is discussed in Section III. Section IV presents our results and discussion of the charmed and charmedstrange mesons. Finally, in the last section, a summary of our findings is presented. Also included is a short description of the role of scattering thresholds for the calculation. Some preliminary results have already been presented in Mohler and Woloshyn (2010).
Ii Calculational setup
A subset of the dynamical 2+1 flavor configurations generated with the CloverWilson fermion action by the PACSCS collaboration Aoki et al. (2009), spanning seaquark pion masses from 702 MeV down to 156 MeV, are used in this work^{2}^{2}2After work on this project started dynamical simulations with even smaller pion masses, for the first time at Aoki et al. (2010a) and around Dürr et al. (2010a, b) the physical point, were reported.. The number of lattice points is for all ensembles, and using the single lattice spacing of 0.0907(13) fm, as determined in Aoki et al. (2009) this corresponds to a box size of roughly 2.9fm in spatial direction and 5.8fm in time direction. Table 1 shows the parameters for the runs. With exception of the ensemble #5 we use configurations on each ensemble.
Ensemble  #configs  

1  1.52617  0.13700  0.13640  200/200 
2  1.52493  0.13727  0.13640  /200 
3  1.52381  0.13754  0.13640  200/200 
4  1.52327  0.13754  0.13660  /200 
5  1.52326  0.13770  0.13640  200/348 
6  1.52264  0.13781  0.13640  198/198 
ii.1 Source construction
The lowlying spectrum is extracted using the variational method Lüscher and Wolff (1990); Michael (1985). For each combination of quantum numbers (where is the Spin and is parity) a matrix of interpolators is constructed
(1) 
Subsequently, the the generalized eigenvalue problem is solved for each time slice
(2)  
At large time separation only a single mass contributes to each eigenvalue. is in general given by the energy difference between the energy level in consideration and the neighboring level. For more details the reader is referred to the discussion in Blossier et al. (2009). For the determination of the charm quark hopping parameter we will use projections to the lowest few lattice momenta, while masses are obtained from a projection to vanishing spatial momentum. In addition to the eigenvalues, the eigenvectors provide useful information and can serve as a fingerprint for a given state.
The basis is constructed from two types of sources. The first type are simple Jacobismeared Güsken et al. (1989); Best et al. (1997) Gaussianshape sources
(3)  
Here, is the spatial hopping term containing the link variables and the smearing has two parameters and . In addition, we also use derivative sources
(4)  
The parameters for the Gaussianshape sources are and . For the derivative sources the number of iterations is slightly larger while using the same . The full list of interpolating fields used in the simulation can be found in Appendix A.1. Similar source constructions have been used for both quarkonium Liao and Manke (2002); Dudek et al. (2008) and light meson Lacock et al. (1996); Gattringer et al. (2008); Dudek et al. (2009) spectroscopy. To minimize the calculation of expensive propagators, we use both Gaussian and derivative sources for the charm quark but only Gaussian sources for the strange and light quarks. This has an additional implication: The interpolators with derivative sources are constructed from a derivative which is not symmetrically applied to both quarks and these interpolators do not have good charge conjugation quantum number. This is relevant in the case of charmonium, where the same basis is used. For further comments regarding charge conjugation, please refer to Section III.
ii.2 Quark propagators
For the construction described in the previous section we need a single full propagator for the strange (light) quark and 4 full propagators (one Gaussian and three derivatives) for the charm quark. These propagators are calculated for 8 source time slices on each gauge configuration. The source locations are chosen randomly within the time slice. For the calculation of the light and strange quark propagators the dfl_sap_gcr inverter from Lüscher’s DDHMC package Lüscher (2007a, b) is employed. For the charm quark propagators the corresponding inverter without deflation is used. The large number of sources needed for this work makes the use of a deflation inverter extremely efficient, particularly for light quarks on Ensemble 6, with a pion mass of 156MeV.
ii.3 Fitting methodology
This section presents the methodology for fits to both single correlators as well as eigenvalues from the variational method. All of these fits take into account autocorrelation in Euclidean time . The estimate of the covariance matrix Cov employs a singleelimination jackknife
(5) 
where N denotes the number of configurations and the bars denote averages. To obtain a stable estimate, the window for the inversion of the covariance matrix is restricted to the chosen fit range. The inverse covariance matrix is estimated once on the ensemble average and used for each jackknife block. This method has been dubbed “jackknife reuse” in Toussaint and Freeman (2008). To determine the window in Euclidean time for performing fits we calculate effective masses
(6) 
and eigenvector components of the regular eigenvalue problem
(7) 
and choose as our fit range the interval where both the effective masses and the eigenvector components are constant within their statistical uncertainties. Then two parameter fits are carried out within this window. As a cross check, 4 parameter fits with two exponentials and a larger fit range are performed. In some cases, the 4 parameter fits are unstable. Where stable, the 4 parameter fits lead to results that are consistent within errors with the strategy described above. The dependence of the results on the time slice was investigated and, where necessary, the reference time was increased until the results were qualitatively unaffected by a further increase. Appendix A.2 collects various fit results, fit intervals and the associated .
Due to the method employed for the heavy quark, the energy levels extracted should not be identified with the physical energy levels directly. Instead, as final results we always quote energy differences with respect to the spinaveraged 1S states, given by
(8)  
for the case of heavystrange mesons, heavylight mesons and charmonium respectively. As the measurements are on the same gauge ensemble and therefore correlated, taking into account this correlation leads to reduced statistical errors for the mass differences.
ii.4 Heavy quark action and quark mass tuning
To determine the mass parameter for the heavy charm quark we use the Fermilab method ElKhadra et al. (1997) in the form employed by the Fermilab Lattice and MILC collaborations Bernard et al. (2010) for their efforts involving quarkonium Burch et al. (2010). This section describes this method in brief and points out where our method differs from Bernard et al. (2010). In addition to this simplest prescription, an improved Fermilab action has been constructed Oktay and Kronfeld (2008) and preliminary results are encouraging Detar et al. (2010). Within this approach, the charm quark hopping parameter is tuned to the value where the spinaveraged kinetic mass assumes its physical value. In this simplest formulation the heavy quark hopping parameter is set to its tadpole improved Lepage and Mackenzie (1993) value , where is the average link. For simplicity, is estimated to be the fourth root of the average plaquette, unlike in Bernard et al. (2010), where it is given by the Landau link. The lattice dispersion relation takes the general form Bernard et al. (2010)
(9) 
where the momentum and is the spatial extent of the lattice. To determine the kinetic mass the six lowest lattice momenta are used and averaged over all possible vectors . These are denoted symbolically as 000, 100, 110, 111, 200 and 210. Larger momenta are very noisy and do not help to constrain the fits further. As there are not enough points to constrain a fit to the fourparameter form of Equation 9, two simplified methods are considered:

neglect the term with coefficient and fit , and .

fit and neglect the term arising from the mismatch of , and
(10)
For both methods, two values of the charm quark hopping parameter close to the physical value are used and the results are interpolated linearly. To determine the energy levels, a matrix of interpolating fields corresponding to interpolators A and B in Appendix A.1 is used in both the and sectors. We choose for these fits. Table 2 shows the results for the spinaveraged ground state from the first method. For illustration an example fit from this method with is shown in Figure 1.
0.8633(5)  0.8931(5)  
0.9337(73)  0.9716(76)  
0.8638(274)  0.8855(284)  
1.0815(86)  1.0878(88)  
2.0315(158)(291)  2.1137(166)(303) 
The three parameter fits are stable and both the rest mass and the kinetic mass are welldetermined. The mass is not used in our analysis but it is encouraging that the result is consistent with the expectation . The error on the kinetic mass in physical units has two contributions and the last row of Table 2 lists both the statistical uncertainty and the uncertainty from setting the lattice scale. This second error from the conversion to physical units is dominant for method 1.
0.8634(5)  0.8932(5)  
1.0889(116)  1.0955(118)  
2.0454(215)(293)  2.1293(227)(305) 
The results from method two can be found in Table 3. This second approach leads to slightly larger values for the kinetic mass . The results of the two methods are however consistent within the statistical uncertainties, which are slightly larger for the second method. With both methods is determined sufficiently well, in the sense that the error due to the lattice scale setting is dominant. The calculation of the heavylight and charmonium spectra in Sections III and IV uses the estimate of from fits with the second method. For this, a plot of the results for the spinaveraged ground state is shown in Figure 2. To account for the slightly unphysical strange quark mass used in the simulation we also calculated the spinaveraged rest mass with a partially quenched strange quark mass corresponding to the estimate of the physical strange quark mass from Aoki et al. (2009). The red curve (squares) in Figure 2 results from shifting our original results to account for this difference. One can then proceed to read off the value with which the simulations to determine the lowlying spectrum are performed.
Having tuned the charm quark hopping parameter using mesons, a stringent test of our setup can be performed by calculating the lowlying charmonium spectrum, where there are several wellestablished states below multiparticlethreshold.
Iii Charmonium results
In recent years many new charmonium resonances have been observed, primarily at the Bfactories. While the agreement between the wellestablished lowlying states and quark model calculations (see for example Godfrey and Isgur (1985)) is quite good, many of the newly observed X, Y and Z states do not fit within a simple quark model picture. For a recent review regarding these charmoniumlike states please refer to Godfrey and Olsen (2008). For an emphasis on recent results see Godfrey (2009); Biassoni (2010); Lange (2010). The charmonium states considered in this section are all well established and considered to be predominantly regular quarkantiquark states with masses well below the and thresholds. This makes the calculation of the lowlying charmonium spectrum a good test of the heavyquark methodology. Where appropriate, we somewhat sloppily refer to those states with their spectroscopic nomenclature and call them S and Pwave states respectively.
As mentioned in Section II.1 some of our interpolators, which were chosen with heavylight mesons in mind, mix different charge conjugations. In the case of the (), the () and the channels () the corresponding admixtures would come from exotic states, which are expected to play a role higher in the spectrum but not in the vicinity of the lowlying states extracted in this section. In the case of the channel(s) the basis is restricted to those interpolators with good charge conjugation and for the interpolators containing the we just extract the ground state, which again should be safe from possible contamination.
Figure 3 shows the chiral behavior of the charmonium hyperfinesplitting compared to the experimental value. The errors in the plot are purely statistical and errors from the tuning of the charm quark mass are nonnegligible^{3}^{3}3Please refer to Section IV for an estimate of this uncertainty in the case of the hyperfinesplitting.. Within our approach, the hyperfinesplitting is expected to show significant discretization errors Burch et al. (2010) which we can not further investigate, as the library of PACSCS lattices used contains only one lattice spacing. This result indicates that discretization effects of a similar magnitude may be expected in other observables. In addition, disconnected diagrams which are neglected in our simulation will have a small effect which may be nonnegligible for the hyperfinesplittings. In Levkova and DeTar (2011) these contributions to the charmonium hyperfinesplitting have been calculated for a reasonable fitting Ansatz. They are found to be negative and therefore can not explain the difference of our results compared to experimental values. Our results also show a small but statistically significant dependence on the strange sea quark mass. Notice however that our charm quark tuning was performed for the ensemble with lightest sea quark mass and . For all other charmonium mass differences the dependence on the strange sea quark mass is very mild and there is no statistically significant difference between the two ensembles at .
Next consider the spintriplet Pwave states, which are the ground states with quantum numbers , and corresponding to the , and . Figure 4 shows the difference between the masses of the individual states and the spinaveraged ground state mass for all ensembles. At the lowest quark mass the ground state Pwaves calculated at finite lattice spacing show good qualitative agreement with experiment, with the ground state showing the largest discrepancy. While our results at the three heaviest and at the lowest quark masses suggest that the chiral extrapolation is within errors well described by the leading order term linear in the pion mass squared, the data from the ensemble at deviates from this behavior. This peculiarity is found in most of the observables we investigate. Note that a similar behavior can be seen in Mahbub et al. (2010) where the first positive parity excitation of the Nucleon, the Roper resonance is calculated on the same gauge configurations. In Mahbub et al. (2010) this is interpreted as “significant chiral curvature”. Our data are not suggestive of this explanation. Without additional simulation data we refrain from interpreting this behavior and will instead quote as our final results those at the lightest pion mass . The effect of the remaining small chiral extrapolation would be negligible compared to current statistical and other systematical uncertainties.
In addition, we also form the combinations (see Bernard et al. (2010) and references therein)
(11)  
(12)  
(13) 
is the Spinaveraged mass of the P wave triplet. is sensitive to the spinorbit interaction while is sensitive to the spinspin interaction Bernard et al. (2010). Figure 5 shows the splitting between the spinaveraged ground state mass and . While we are unable to perform a continuum extrapolation and small discrepancies to experiment are expected, the figure suggests that the spinaveraged quantities are well reproduced within our approach.
The results for and are plotted in Figure 6. The upper panel shows the SpinOrbit splitting of the 1P states, which is slightly underestimated in our simulation, similar to the hyperfinesplitting. The lower panel shows the tensor splitting of the 1P states. Compared to experiment we obtain the correct sign and a similar magnitude with the value off about 30%.
So far the discussion was restricted to properties of ground states. In the and channels our basis allows for a determination of excited energy levels. While these states are more noisy than the ground states and the fit ranges (see tables in the Appendix A.2) are more limited, at least the 2S states are well determined from our charmonium data. Figure 7 shows the splitting between the 2S and 1S spinaveraged states. Again the result agrees qualitatively with experiment suggesting that our approach is well suited to reproduce spinaveraged quantities.
Figure 8 shows the 2S hyperfinesplitting. Unlike the 1S hyperfinesplitting this quantity is not determined very accurately from our data. Just like in experiment the 2S hyperfinesplitting is smaller than the 1S hyperfinesplitting. This again gives us confidence that we capture the important physics well. This concludes the discussion of charmonium results. Final numbers from the ensemble with lightest sea quarks can be found in Section V along with a table that summarizes the results and their respective uncertainties. A table with charmonium energy levels can be found in Appendix A.2.
Iv Heavylight mesons
In the heavy quark limit, the spectrum of S and Pwave heavylight mesons can be classified by the total angular momentum of the light degrees of freedom^{4}^{4}4Note that we use small when referring to the heavy quark limit and when referring to the Spin away from the heavy quark limit.. This classification contains one Swave doublet with containing the and mesons in case of the heavylight states and the and mesons in case of the heavystrange states. The Pwave states fall into two multiplets, one with and one with . Taking a look at the experimental situation one would identify the charmedstrange mesons and as the charmedstrange version of the Pwave heavy quark doublet with and the and the as the charmedstrange version of the multiplet with . While the mass of states corresponding to and states in the heavy mass limit is quite close to the values expected from the quark model Godfrey and Isgur (1985), the states corresponding to the heavystrange doublet are substantially lighter Nakamura et al. (2010), putting them below the and thresholds. Even more puzzling however is a comparison to the spectrum of heavylight mesons, where similar multiplets are expected in the heavy quark limit. In particular, the charmedlight version of the doublet is given by the and the . This corresponds to a masssplitting between these states and the spinaveraged Swave states which is much larger for mesons than for mesons. Such a behavior is not expected in simple potential models. Furthermore, chiral perturbation theory for staticlight mesons suggests the opposite Becirevic et al. (2004).
Before proceeding to report results, let us discuss previous lattice calculations of charmed and charmedstrange mesons focusing in particular on results for the doublet. Lewis and Woloshyn Lewis and Woloshyn (2000) calculated the spectrum of S and Pwave heavylight mesons in lattice NRQCD using quenched gauge configurations. At the time of this study, the states were not yet measured. While their results for the heavier doublet agree within errors with the experimental results, the , and the were predicted at masses much larger than the experimental masses. Another early calculation using lattice NRQCD was performed by Hein et al. Hein et al. (2000). In Boyle (1997, 1998) the UKQCD collaboration reports on results using a relativistic clovertype action based on the Fermilab method ElKhadra et al. (1997). Results from the studies Hein et al. (2000); Boyle (1997, 1998) along with a calculation in the static limit on dynamic gauge configurations have been presented by Bali Bali (2003). All these results share the common feature that the ground state in the channel is estimated at a mass exceeding the mass of the experimental by roughly 130200 MeV. Furthermore, where calculated, the splitting between the Swave states and the Pwave states for heavystrange mesons was determined to be larger or of the same size than for heavylight states.
In addition to the above simulations, the UKQCD collaboration also presented results from both quenched and one dynamical ensemble in Dougall et al. (2003). While they report a somewhat smaller discrepancy with respect to experiment, their scale setting results in a quite large value of for the Sommer scale . Typical values for the Sommer scale from recent dynamical simulations however suggest a significantly smaller value . In particular the scale setting for our simulations, as determined by the PACSCS collaboration, corresponds to Aoki et al. (2009). A smaller value of results in larger mass splittings. The scale setting itself is intimately related to the determination of the physical point, which again emphasizes the need for light dynamical quarks.
More recently, the QCD collaboration presented preliminary results from both quenched Dong and Liu (2008) and dynamical Dong et al. (2009) lattices using relativistic quarks and overlap valence fermions. Their initial results for both quenched and dynamical simulations suggest good agreement between lattice results and the lowlying charmedstrange mesons.
The rest of this section presents results for the spectrum of heavystrange and heavylight mesons. The final numbers and an overview of the results can be found in Section V.
Figure 9 shows the hyperfinesplittings for the (upper panel) and (lower panel) mesons. Similar to Charmonium, the hyperfinesplitting turns out to be somewhat to small, which we attribute partly to nonnegligible discretization effects. This discrepancy is somewhat smaller for heavylight states than for charmonium. Like for charmonium the hyperfinesplitting on the ensemble with is slightly smaller than for the other ensembles. Another nonnegligible source of error is the uncertainty in the determination of . From the tuning runs, this uncertainty is estimated to be roughly 3.2 MeV for mesons and we assume a similar dependence for mesons. In the case of mesons where statistical uncertainties are small, neglecting the data point at our second lowest quark mass, the remaining data could be fit with a good by a fit linear in the pion mass squared. Although the results have much larger statistical errors, chiral effects seem to be more significant for mesons.
Figure 10 shows results for the (left panels) and (right panels) Pwave states. As for charmonium, mass differences relative to the respective spinaveraged Swave state are plotted. The top row shows the ground state in the and () channels. In both cases, the mass splitting gets smaller as the light sea quark masses approach the chiral limit. Unlike suggested by experiment, the splitting for the is smaller than for the . While the discrepancy between our data and the is smaller for almost physical light quarks, calculation at light sea quarks alone can not resolve the apparent puzzle between the experimental results and the lattice data. For the mesons, statistical errors are larger. As previously remarked for the charmonium hyperfinesplitting in Section III, the data calculated on Ensemble 5 (see Table 1) seem to lead to larger mass splittings for all channels, while the rest of the data suggest that the leading order term proportional to the pion mass squared dominates the chiral extrapolation. This is especially evident for the meson. For the final numbers and overview plots in Section V, we therefore quote the numbers obtained on the ensemble with the lightest pion mass. Any sensible functional form for a chiral extrapolation would alter these results by an amount much smaller than the dominant statistical and systematic errors.
The middle panels of Figure 10 show the result for the two lowest states in the and () channels. Both plots are from a matrix of correlators which includes interpolating fields corresponding to both positive and negative charge conjugation in the massdegenerate case. In Bali (2003) the importance of this mixing has been emphasized, but the actual calculation neglected the mixing. In Boyle (1998) the mixing has been considered, but only the ground state could be resolved. Figure 11 illustrates the effect of this mixing. In particular the standard interpolators A and D (see Appendix A.2) both couple strongly to the ground state leading to a much smaller splitting than the mixed basis. The results from the basis in Figure 10 exhibit a mass splitting which is still somewhat small compared to the experimental splitting but enhanced compared to the analysis without mixing. For the the results for the ground state lead to a splitting substantially larger than the , while the difference between the first excited state obtained in the calculation and the is compatible with the nonnegligible discretization errors one expects. For the mesons, the ground state displays a somewhat larger slope in the chiral extrapolation, while the splitting for the first excited state is comparable to the case.
Finally the bottom two panels of Figure 10 show the ground states of the (lefthand side) and (righthand side) with . For some ensembles we encounter a systematic uncertainty (of a size similar to the purely statistical error) with regard to the choice of fit range in this channel. On the ensemble with lowest quark mass, the state is well determined though, and the difference from the experiment value is small. The ground state is somewhat less well determined and the difference from experiment is of a similar magnitude as for its partner state in the heavy quark multiplet.
For the mesons, it was possible to extract some information about excited states in the and channels. The top panel of Figure 12 shows the resulting splitting between the spinaveraged 1S and 2S states. To compare with experiment we also plot the single excited and states in the middle and bottom panels of Figure 12. The displayed errors are statistical only. For both 2S states, the result for the splitting at the lightest quark mass is substantially smaller than on all other ensembles. On ensemble 4 with the results turn out to be somewhat larger than for all other ensembles. From experiment, a excited state, the is known. The result at the lowest quark mass is somewhat larger than experiment, but seems reasonable given the statistical and systematic uncertainties, which are larger for excited states than for the ground states. (Regarding the importance of possible multihadron thresholds, remarks similar to those about the and ground states apply.) To complete the analysis of the 2S states, we also plot the 2S hyperfinesplitting in Figure 13. In general one expects that the hyperfinesplitting decreases for states higher up in the spectrum. In quark models, (see for example Godfrey and Isgur (1985)) a typical splitting for the 2S states would be . For mesons the BaBar collaboration recently identified candidates del Amo Sanchez et al. (2010) for the corresponding states and the splitting between those meson states is of comparable size. While the errors are large compared to the small splitting, we obtain values consistent with these expectations.
V Summary & discussion
We presented results from a dynamical lattice QCD calculation of heavylight mesons. To investigate the effects of light dynamical quarks, configurations generated by the PACSCS collaboration Aoki et al. (2009), with light quarks corresponding to a pion mass as light as 156 MeV, have been used. To test the computational setup, we performed a calculation of the lowlying charmonium spectrum, where scattering thresholds are less important and previous calculations within the same framework exist Bernard et al. (2010). The results from the ensemble with a pion mass close to the physical pion are summarized in Figure 14. The qualitative features of the lowlying charmonium spectrum are well reproduced which gives us confidence in the procedure used to put the heavy charm quark on the lattice. The numerical values for charmonium mass differences and estimates of some of the uncertainties are provided in Table 4. With a calculation at only one lattice spacing it is not possible to make a quantitative estimate of discretization effects. However, the expected truncation errors for the Fermilab action as presented in Figs. 3 and 4 of Oktay and Kronfeld (2008) can provide a rough guide. With the action used here, discretization effects in spinaveraged mass differences of a few percent could be expected. For more sensitive quantities, e.g., the hyperfinesplitting, only simulation can provide a reliable estimate.
Mass difference  This paper [MeV]  Experiment [MeV] 

1S hyperfine  
1P spinorbit  
1P tensor  
2S hyperfine 
Figure 15 shows the results for heavylight and heavystrange mesons as calculated on the ensemble with the lightest quark masses. In both cases differences between the state of interest and the spinaveraged ground state are shown. Table 5 provides numerical values for these splittings and for the 1S and 2S (for the ) hyperfinesplittings. Even at pion masses close to the physical pion mass, significant differences in comparison to experiment remain for the case of mesons. This is especially evident for the states corresponding to a doublet in the heavy quark limit. For the mesons, errors are larger and the difference from experimental values is more pronounced for the states. The results for the doublet on the other hand are qualitatively different from the heavystrange case.
Mass difference  This paper [MeV]  Experiment [MeV] 

  
  
1S hyperfine  
2S hyperfine    
1S hyperfine 
While light dynamical quarks improve the agreement with experiment for heavystrange and heavylight mesons, considerable discrepancies remain in the case of the and channels. There are a number of possible reasons for this. For one, future investigations will have to include a full continuum extrapolation. At our lattice spacing, discretization errors are nonnegligible. However, in the light of previous results for the charmonium spectrum, it seems unlikely that discretization errors alone are to blame for the observed behavior. Another possible reason for discrepancies is the lack of an extrapolation to the infinite volume limit. While the effects of finite volume are small at large sea quark masses, there might be a concern that this is not the case for the ensemble with the lightest sea quark where However, finite volume effects have been studied for heavy mesons in Colangelo et al. (2010) and for the ground state D meson are found to be smaller than 1% for the pion mass and lattice volume used here. Overall the following picture emerges: Given the limitations of the current analysis the and multiplets are in reasonable agreement with experiment, while the large masses of the heavystrange doublet resulting from our calculation are hard to explain by these limitations. In particular, we consider it unlikely that a continuum extrapolation will change these results qualitatively.
In the case of the heavystrange states, the and thresholds may play an important role. Figure 16 shows the lowest energy levels for the (left panel) and (right panel) ground states. The physical and thresholds are indicated with a red plus. The black plus signs show the central value of these scattering states on the ensembles with the three lightest sea quark masses. While the mesons have been determined in this paper, the kaon mass has been taken from Aoki et al. (2009). The leading discrepancy between the scattering level at the physical point and the scattering level on the lattice results from the unphysically heavy strange quark mass used in our simulations. An important observation is that the overlap of regular quarkantiquark interpolating fields on multihadron states (which has to be present in dynamical simulations) is most likely small and current state of the art spectrum calculations Dudek et al. (2010); Engel et al. (2010) do not see most of the possible scattering states with a simple basis. (For states in a Pwave, the contribution from scattering states is expected to be volume suppressed. However, for many channels even the expected Swave levels are absent.) Comparing the eigenvectors at our lightest and thirdlightest pion mass suggests that we observe the same state on all ensembles and that the expected Swave scattering level is also absent from our data. The obvious solution is to include multihadron states in the basis for the variational method, which is a challenging task left for future calculations. Including multihadron states in the basis would also enable one to treat these mesons correctly as resonances. The mass and width of hadronic resonances below the inelastic threshold can be determined within Lüscher’s finite volume framework Lüscher (1986, 1991). First attempts to include scattering states for the calculation of the resonance mass and width have been published in Aoki et al. (2007); Feng et al. (2010); Aoki et al. (2010b); Frison et al. (2010) and similar studies are currently underway for other lightquark resonances. For charmonium excitations, multiparticle states have been included in Bali and Ehmann (2009). For heavylight mesons first steps in this direction have been made in Liu et al. (2008).
Finally, there is also the possibility that the states in question are indeed not conventional quarkantiquark states, but rather tetraquark states or other states with which our basis of interpolating fields has poor overlap. For the case of the as a tetraquark, this has been investigated in lattice QCD recently Gong et al. (2011) and no lowlying tetraquark states have been identified.
Acknowledgements.
We thank the PACSCS collaboration for access to their gauge configurations and Martin Lüscher for making his DDHMC software available. The calculations were performed on computing clusters at TRIUMF and York University. We thank Sonia Bacca and Randy Lewis for making these resources available. D.M. would like to thank Carleton DeTar, Georg Engel and Sasa Prelovsek for helpful discussions. This work is supported in part by the Natural Sciences and Engineering Research Council of Canada (NSERC).Appendix A Tables
a.1 Interpolating fields
Table 6 collects the interpolators used for both charmonium and heavylight mesons. For the details of the source construction, please refer to Section II.1.
Lattice irrep  Quantum numbers in irrep  Interpolator label  Operator 

, ,  A  
B  
C  
D  
, ,  A  
B  
C  
, , ,  A  
B  
C  
D  
E  
F  
, , ,  A  
B  
C  
D  
E  
F  
, , ,  A  
B  
, , ,  A  
B  
, ,  A  
B  
, ,  A  
B 
a.2 Results
In this appendix, the combination of interpolating fields used for the fits, the respective fit intervals and results for the energy levels are tabulated in Tables 7, 8 and 9. The ensembles are labeled by their number in Table 1. We also provide the resulting of the correlated fits and the reference time slice whenever a matrix of interpolators was used. The states are labeled by their quantum numbers. Excited states are listed right below the respective ground states. For a full list of interpolating fields please refer to Section A.1.
State  Interpolators  Ensemble  Energy  

ground state  A,B,C  1  7  26  1.22924(23)  1.22  4 
2  7  27  1.22098(21)  2.03  4  
3  7  27  1.21369(20)  0.75  4  
4  7  27  1.21017(24)  0.71  4  
5  7  27  1.21268(15)  0.71  4  
6  7  27  1.20765(21)  1.82  4  
first excitation  A,B,C  1  7  11  1.5595(142)  1.21  4 
2  7  15  1.5341(136)  0.53  4  
3  7  14  1.5395(97)  0.48  4  
4  7  12  1.5247(122)  2.12  4  
5  6  14  1.5385(71)  1.78  4  
6  6  13  1.5151(82)  0.59  4  
ground state  A,B,C,D  1  7  26  1.27877(37)  1.01  4 
2  7  27  1.26860(33)  2.08  4  
3  7  27  1.25988(31)  1.39  4  
4  7  27  1.25499(34)  0.99  4  
5  7  27  1.25904(22)  0.93  4  
6  7  27  1.25258(35)  1.45  4  
first excitation  A,B,C,D  1  7  15  1.5876(110)  0.85  4 
2  7  15  1.5824(110)  0.83  4  
3  7  13  1.5535(95)  0.09  4  
4  7  12  1.5467(99)  3.16  4  
5  7  13  1.5618(71)  2.35  4  
6  7  13  1.5373(86)  0.83  4  
ground state  A,B,C  1  7  24  1.4597(29)  0.46  3 
2  7  22  1.4468(27)  1.83  3  
3  7  23  1.4256(29)  0.57  3  
4  7  23  1.4223(29)  0.62  3  
5  7  23  1.4289(18)  0.33  3  
6  7  24  1.4165(28)  0.70  3  
ground state  A  1  7  24  1.4937(30)  0.73   
2  7  25  1.4756(30)  1.41    
3  7  24  1.4584(27)  0.72    
4  7  25  1.4538(26)  0.75    
5  7  24  1.4639(22)  0.69    
6  7  26  1.4475(27)  0.77    
ground state  A,B in T2 irrep  1  7  22  1.5095(34)  0.97  3 
2  7  22  1.4961(33)  0.74  3  
3  7  21  1.4808(29)  0.86  3  
4  7  26  1.4735(27)  0.96  3  
5  7  24  1.4856(22)  0.91  3  
6  7  24  1.4661(26)  0.88  3  
ground state  D  1  7  21  1.4932(31)  0.80   
2  7  22  1.4826(39)  0.98    
3  7  24  1.4670(31)  0.39    
4  7  21  1.4627(34)  1.19    
5  7  24  1.4711(27)  0.97    
6  7  25  1.4571(33)  1.38   
State  Interpolators  Ensemble  Energy  

ground state  A,B,C  1  6  26  0.86273(33)  1.42  3 
2  6  26  0.85158(37)  1.04  3  
3  6  24  0.84000(36)  1.16  3  
4  7  25  0.82848(40)  1.37  3  
5  6  24  0.83929(26)  1.22  3  
6  7  25  0.83149(34)  1.80  3  
first excitation  A,B,C  1  6  11  1.2423(175)  1.69  3 
2  6  11  1.2361(205)  1.74  3  
3  6  11  1.2406(157)  1.19  3  
4  6  12  1.2418(168)  1.65  3  
5  6  12  1.2294(117)  1.32  3  
6  6  13  1.1782(123)  1.24  3  
ground state  A,B,C,D  1  6  28  0.93297(59)  0.73  3 
2  6  26  0.91864(63)  1.51  3  
3  6  25  0.90429(60)  1.57  3  
4  7  27  0.89015(69)  1.20  3  
5  6  24  0.90429(43)  0.68  3  
6  7  30  0.89268(59)  1.88  3  
first excitation  A,B,C,D  1  6  13  1.2772(135)  0.63  3 
2  6  10  1.2732(150)  0.91  3  
3  6  12  1.2600(126)  0.95  4  
4  6  12  1.2794(117)  2.07  3  
5  6  12  1.2545(88)  1.03  3  
6  6  13  1.2113(96)  0.56  4  
ground state  A,B,C  1  6  21  1.1057(35)  0.98  3 
2  6  21  1.0784(40)  0.59  3  
3  6  20  1.0509(43)  0.92  3  
4  7  13  1.0399(44)  1.28  3  
5  6  24  1.0570(25)  0.68  3  
6  6  20  1.0342(37)  1.54  3  
ground state  A,C,D,F  1  7  13  1.1574(35)  0.12  4 
2  7  15  1.1319(35)  1.22  4  
3  7  15  1.1073(34)  1.01  4  
4  7  13  1.0907(33)  2.63  4  
5  7  15  1.1130(23)  1.04  4  
6  7  14  1.0888(31)  0.79  4  
first excitation  A,C,D,F  1  7  13  1.1654(54)  0.82  4 
2  7  15  1.1458(55)  1.02  4  
3  7  15  1.1273(49)  0.27  4  
4  7  13  1.1155(51)  2.92  4  
5  7  15  1.1320(35)  1.22  4  
6  7  14  1.1047(44)  1.01  4  
ground state  A,B in T2 irrep  1  7  20  1.1810(52)  0.76  3 
2  7  17  1.1662(60)  1.11  3  
3  7  18  1.1491(52)  0.72  3  
4  7  18  1.1353(51)  1.47  3  
5  7  19  1.1530(36)  1.41  3  
6  7  20  1.1241(44)  1.29  3 
State  Interpolators  Ensemble  Energy  

ground state  A,B,C  1  6  26  0.84022(37)  1.52  3 
3  7  23  0.79580(61)  0.79  3  
5  7  26  0.78798(82)  1.35  3  
6  7  24  0.77646(119)  1.23  3  
ground state  A,B,C,D  1  5  26  0.91237(60)  0.89  3 
3  7  26  0.86327(99)  0.92  3  
5  7  22  0.85776(122)  1.41  3  
6  7  24  0.83656(189)  1.10  3  
ground state  A,B,C  1  6  15  1.0911(50)  0.28  3 
3  6  11  0.9956(80)  0.38  3  
5  6  17  1.0150(66)  1.04  3  
6  6  18  0.9442(79)  0.68  3  
ground state  A,D  1  7  14  1.1394(47)  1.09  4 
3  7  15  1.0631(48)  0.64  4  
5  7  12  1.0601(58)  2.07  4  
6  7  15  1.0050(63)  0.84  4  
first excitation  A,D  1  7  14  1.1538(65)  1.00  4 
3  7  15  1.0874(69)  0.77  4  
5  7  12  1.0883(78)  0.65  4  
6  7  15  1.0629(93)  0.79  4  
ground state  A,B in T2 irrep  1  7  18  1.1869(60)  0.85  3 
6  7  18  1.1153(74)  1.48  3  
3  7  15  1.1166(80)  1.02  3  
5  7  18  1.0868(95)  1.10  3 
References
 Nakamura et al. (2010) K. Nakamura et al. (Particle Data Group), Journal of Physics G: Nuclear and Particle Physics 37, 075021 (2010), URL http://stacks.iop.org/09543899/37/i=7A/a=075021.
 Godfrey and Isgur (1985) S. Godfrey and N. Isgur, Phys. Rev. D32, 189 (1985).
 Aubert et al. (2006) B. Aubert et al. (BABAR), Phys. Rev. Lett. 97, 222001 (2006), eprint arXiv:hepex/0607082.
 Aubert et al. (2009) B. Aubert et al. (BABAR), Phys. Rev. D80, 092003 (2009), eprint arXiv:0908.0806.
 Brodzicka et al. (2008) J. Brodzicka et al. (Belle), Phys. Rev. Lett. 100, 092001 (2008), eprint arXiv:0707.3491.
 Evdokimov et al. (2004) A. V. Evdokimov et al. (SELEX), Phys. Rev. Lett. 93, 242001 (2004), eprint arXiv:hepex/0406045.
 Colangelo and De Fazio (2010) P. Colangelo and F. De Fazio, Phys. Rev. D81, 094001 (2010), eprint arXiv:1001.1089.
 Ebert et al. (2010) D. Ebert, R. N. Faustov, and V. O. Galkin, Eur. Phys. J. C66, 197 (2010), eprint arXiv:0910.5612.
 Chen et al. (2009) B. Chen, D.X. Wang, and A. Zhang, Phys. Rev. D80, 071502 (2009), eprint arXiv:0908.3261.
 van Beveren and Rupp (2010) E. van Beveren and G. Rupp, Phys. Rev. D81, 118101 (2010), eprint arXiv:0908.1142.
 Sun and Liu (2009) Z.F. Sun and X. Liu, Phys. Rev. D80, 074037 (2009), eprint arXiv:0909.1658.
 Abreu et al. (1998) P. Abreu et al. (DELPHI), Phys. Lett. B426, 231 (1998).
 del Amo Sanchez et al. (2010) P. del Amo Sanchez et al. (The BABAR), Phys. Rev. D82, 111101 (2010), eprint arXiv:1009.2076.
 Sun et al. (2010) Z.F. Sun, J.S. Yu, X. Liu, and T. Matsuki, Phys. Rev. D82, 111501 (2010), eprint arXiv:1008.3120.
 Li et al. (2011) D.M. Li, P.F. Ji, and B. Ma, Eur. Phys. J. C71, 1582 (2011), eprint arXiv:1011.1548.
 Wang (2011) Z.G. Wang, Phys. Rev. D83, 014009 (2011), eprint arXiv:1009.3605.
 Zhong (2010) X.H. Zhong, Phys. Rev. D82, 114014 (2010), eprint arXiv:1009.0359.
 Chen et al. (2011) B. Chen, L. Yuan, and A. Zhang Phys. Rev. D 83, 114025 (2011), eprint arXiv:1102.4142.
 Dürr et al. (2008) S. Dürr et al., Science 322, 1224 (2008), eprint arXiv:0906.3599.
 Swanson (2006) E. S. Swanson, Phys. Rept. 429, 243 (2006), eprint arXiv:hepph/0601110.
 Zhu (2008) S.L. Zhu, Int. J. Mod. Phys. E17, 283 (2008), eprint arXiv:hepph/0703225.
 Mohler and Woloshyn (2010) D. Mohler and R. M. Woloshyn, PoS LATTICE2010, 116 (2010), eprint arXiv:1010.2786.
 Aoki et al. (2009) S. Aoki et al. (PACSCS), Phys. Rev. D79, 034503 (2009), eprint arXiv:0807.1661.
 Aoki et al. (2010a) S. Aoki et al. (PACSCS), Phys. Rev. D81, 074503 (2010a), eprint arXiv:0911.2561.
 Dürr et al. (2010a) S. Dürr et al. (2010a), eprint arXiv:1011.2403.
 Dürr et al. (2010b) S. Dürr et al. (2010b), eprint arXiv:1011.2711.
 Lüscher and Wolff (1990) M. Lüscher and U. Wolff, Nucl. Phys. B339, 222 (1990).
 Michael (1985) C. Michael, Nucl. Phys. B259, 58 (1985).
 Blossier et al. (2009) B. Blossier, M. Della Morte, G. von Hippel, T. Mendes, and R. Sommer, JHEP 04, 094 (2009), eprint arXiv:0902.1265.
 Güsken et al. (1989) S. Güsken et al., Phys. Lett. B227, 266 (1989).
 Best et al. (1997) C. Best et al., Phys. Rev. D56, 2743 (1997), eprint arXiv:heplat/9703014.
 Liao and Manke (2002) X. Liao and T. Manke (2002), eprint arXiv:heplat/0210030.
 Dudek et al. (2008) J. J. Dudek, R. G. Edwards, N. Mathur, and D. G. Richards, Phys. Rev. D77, 034501 (2008), eprint arXiv:0707.4162.
 Lacock et al. (1996) P. Lacock, C. Michael, P. Boyle, and P. Rowland (UKQCD), Phys. Rev. D54, 6997 (1996), eprint arXiv:heplat/9605025.
 Gattringer et al. (2008) C. Gattringer, L. Y. Glozman, C. B. Lang, D. Mohler, and S. Prelovsek, Phys. Rev. D78, 034501 (2008), eprint arXiv:0802.2020.
 Dudek et al. (2009) J. J. Dudek, R. G. Edwards, M. J. Peardon, D. G. Richards, and C. E. Thomas, Phys. Rev. Lett. 103, 262001 (2009), eprint arXiv:0909.0200.
 Lüscher (2007a) M. Lüscher, JHEP 07, 081 (2007a), eprint arXiv:0706.2298.
 Lüscher (2007b) M. Lüscher, JHEP 12, 011 (2007b), eprint arXiv:0710.5417.
 Toussaint and Freeman (2008) D. Toussaint and W. Freeman (2008), eprint arXiv:0808.2211.
 ElKhadra et al. (1997) A. X. ElKhadra, A. S. Kronfeld, and P. B. Mackenzie, Phys. Rev. D55, 3933 (1997), eprint arXiv:heplat/9604004.
 Bernard et al. (2010) C. Bernard et al. (2010), eprint arXiv:1003.1937.
 Burch et al. (2010) T. Burch et al., Phys. Rev. D81, 034508 (2010), eprint arXiv:0912.2701.
 Oktay and Kronfeld (2008) M. B. Oktay and A. S. Kronfeld, Phys. Rev. D78, 014504 (2008), eprint arXiv:0803.0523.
 Detar et al. (2010) C. Detar, A. S. Kronfeld, and M. B. Oktay, PoS LATTICE2010, 234 (2010), eprint arXiv:1011.5189.
 Lepage and Mackenzie (1993) G. P. Lepage and P. B. Mackenzie, Phys. Rev. D48, 2250 (1993), eprint arXiv:heplat/9209022.
 Godfrey and Olsen (2008) S. Godfrey and S. L. Olsen, Ann. Rev. Nucl. Part. Sci. 58, 51 (2008), eprint arXiv:0801.3867.
 Godfrey (2009) S. Godfrey (2009), eprint arXiv:0910.3409.
 Biassoni (2010) P. Biassoni (BaBar) (2010), eprint arXiv:1009.2627.
 Lange (2010) J. S. Lange (Belle) (2010), eprint arXiv:1010.2331.
 Levkova and DeTar (2011) L. Levkova and C. DeTar, Phys.Rev. D83, 074504 (2011), eprint arXiv:1012.1837.
 Mahbub et al. (2010) M. S. Mahbub, W. Kamleh, D. B. Leinweber, P. J. Moran, and A. G. Williams (2010), eprint arXiv:1011.5724.
 Becirevic et al. (2004) D. Becirevic, S. Fajfer, and S. Prelovsek, Phys. Lett. B599, 55 (2004), eprint arXiv:hepph/0406296.
 Lewis and Woloshyn (2000) R. Lewis and R. M. Woloshyn, Phys. Rev. D62, 114507 (2000), eprint arXiv:heplat/0003011.
 Hein et al. (2000) J. Hein et al., Phys. Rev. D62, 074503 (2000), eprint arXiv:hepph/0003130.
 Boyle (1997) P. Boyle (UKQCD), Nucl. Phys. Proc. Suppl. 53, 398 (1997).
 Boyle (1998) P. Boyle (UKQCD), Nucl. Phys. Proc. Suppl. 63, 314 (1998), eprint arXiv:heplat/9710036.
 Bali (2003) G. S. Bali, Phys. Rev. D68, 071501 (2003), eprint arXiv:hepph/0305209.
 Dougall et al. (2003) A. Dougall, R. D. Kenway, C. M. Maynard, and C. McNeile (UKQCD), Phys. Lett. B569, 41 (2003), eprint arXiv:heplat/0307001.
 Dong and Liu (2008) S. J. Dong and K. F. Liu (chi QCD), PoS LATTICE2008, 117 (2008), eprint arXiv:0810.2993.
 Dong et al. (2009) S. J. Dong et al. (QCD), PoS LAT2009, 090 (2009), eprint arXiv:0911.0868.
 Colangelo et al. (2010) G. Colangelo, A. Fuhrer, and S. Lanz, Phys.Rev. D82, 034506 (2010), eprint arXiv:1005.1485.
 Dudek et al. (2010) J. J. Dudek, R. G. Edwards, M. J. Peardon, D. G. Richards, and C. E. Thomas, Phys. Rev. D82, 034508 (2010), eprint arXiv:1004.4930.
 Engel et al. (2010) G. P. Engel, C. B. Lang, M. Limmer, D. Mohler, and A. Schafer (BGR [BernGrazRegensburg]), Phys. Rev. D82, 034505 (2010), eprint arXiv:1005.1748.
 Lüscher (1986) M. Lüscher, Commun. Math. Phys. 105, 153 (1986).
 Lüscher (1991) M. Lüscher, Nucl. Phys. B354, 531 (1991).
 Aoki et al. (2007) S. Aoki et al. (CPPACS), Phys. Rev. D76, 094506 (2007), eprint arXiv:0708.3705.
 Feng et al. (2010) X. Feng, K. Jansen, and D. B. Renner (2010), eprint arXiv:1011.5288.
 Aoki et al. (2010b) P.. S. Aoki et al. (CS), PoS LATTICE2010, 108 (2010b), eprint arXiv:1011.1063.
 Frison et al. (2010) J. Frison et al. (BudapestMarseilleWuppertal), PoS LATTICE2010, 139 (2010), eprint arXiv:1011.3413.
 Bali and Ehmann (2009) G. Bali and C. Ehmann, PoS LAT2009, 113 (2009), eprint arXiv:0911.1238.
 Liu et al. (2008) L. Liu, H.W. Lin, and K. Orginos, PoS LATTICE2008, 112 (2008), eprint arXiv:0810.5412.
 Gong et al. (2011) M. Gong et al. (2011), eprint arXiv:arXiv:1103.0589.