A Table of results

Excited and exotic charmonium spectroscopy from lattice QCD


We present a spectrum of highly excited charmonium mesons up to around  GeV calculated using dynamical lattice QCD. Employing novel computational techniques and the variational method with a large basis of carefully constructed operators, we extract and reliably identify the continuum spin of an extensive set of excited states, states with exotic quantum numbers (, , ) and states with high spin. Calculations are performed on two lattice volumes with pion mass and the mass determinations have high statistical precision even for excited states. We discuss the results in light of experimental observations, identify the lightest ‘supermultiplet’ of hybrid mesons and comment on the phenomenological implications of the spectrum of exotic mesons.

a,1]Liuming Liu\noteliuming@maths.tcd.ie a]Graham Moir a,2]Michael Peardon\notemjp@maths.tcd.ie a,3]Sinéad M. Ryan\noteryan@maths.tcd.ie a,4]Christopher E. Thomas\notethomasc@maths.tcd.ie a]Pol Vilaseca \affiliation[a] School of Mathematics, Trinity College, Dublin 2, Ireland b,c]
Jozef J. Dudek b]Robert G. Edwards b]Bálint Joó b]David G. Richards \affiliation[b]Jefferson Laboratory, 12000 Jefferson Avenue, Newport News, VA 23606, USA \affiliation[c]Department of Physics, Old Dominion University, Norfolk, VA 23529, USA \affiliation (for the Hadron Spectrum Collaboration) \preprint TCDMATH 12-04 JLAB-THY-12-1510 \arxivnumber1204.5425

1 Introduction

There has been a resurgence of interest in charmonium spectroscopy over the last decade driven by experimental observations at the -factories, CLEO-c, the Tevatron and BES of rather narrow states close to or above open-charm thresholds. Some had been predicted but hitherto unobserved whereas many were unexpected, for example the enigmatic , and other “X,Y,Z”s. Prior to these discoveries, the relatively small number of observed mesons all fitted into the simple pattern predicted by quark models. The new states have spurred theoretical discussion as to their nature with suggestions including conventional quark-antiquark states appearing at unexpected masses, hybrid mesons in which the gluonic field is excited, molecular mesons and tetraquarks. To date, however, there has been no observation of a state in the charmonium region with manifestly exotic ( and are respectively spin, parity and charge-conjugation quantum numbers); an observation of such an exotic would be a smoking gun for physics beyond quark potential models. Experimental interest continues with BESIII, experiments at the LHC, and the planned PANDA experiment at GSI/FAIR. An extensive review of the experimental and theoretical situation is given in Ref. [1].

From this perspective, the charmonium system is a particularly interesting probe of the strong regime of QCD. Being reasonably non-relativistic it is well described below the open-charm () threshold by modelling it as a quark and anti-quark bound together in a confining potential. It has also been studied using effective field theory approaches and in the framework of QCD sum rules. No real consensus has emerged about the nature of many of the new states.

A vital part of understanding and testing QCD is to make precise predictions of the hadron spectrum it implies, as opposed to the spectrum in a QCD-inspired approach, and to test these predictions against high-quality experimental data. The charmonium system provides an important arena for this and, in particular, discerning the puzzling states within QCD is a crucial aspect of understanding the strong interaction. Lattice QCD provides a method for performing such ab-initio calculations within QCD: the theory is discretised on a finite lattice and then energies and other properties are extracted from numerical computations of Euclidean-time correlation functions.

The precision calculation of the masses of the lower-lying states has long been an important benchmark of lattice calculations; examples of recent studies of the lightest charmonium mesons are given Refs. [2, 3, 4, 5] and reviewed in Ref. [6]. More recently, progress has been made in using lattice QCD to study higher-lying states. A pioneering quenched calculation of the excited charmonium spectrum was presented in Ref. [7] and there have since been some dynamical studies [8, 9]; we comment on these previous studies later when we discuss the results of this work.

In lattice calculations, the theory is discretised on a finite hypercubic grid and the full rotational symmetry of the continuum is reduced to the symmetry group of a cube. States at rest are then not classified according to spin, , but by the irreducible representations, irreps, of this reduced symmetry group; the various components of high spin states are distributed across several irreps. The Hadron Spectrum Collaboration has overcome the challenges this reduced symmetry presents by using a large basis of carefully constructed operators with a combination of techniques and dynamical anisotropic lattices. This methodology has been shown to be very effective at extracting extensive spectra of light isovector and isoscalar mesons [10, 11, 12] and baryons [13, 14, 15], including highly excited states, states with high spin and those with exotic quantum numbers, and in reliably determining the continuum of the extracted states.

In this article we present a dynamical lattice calculation of the excited spectrum of charmonium mesons, the first application to charmonium of these techniques. The results presented here go beyond any previous dynamical calculation of the excited charmonium spectrum. We will show how we can reliably extract an unprecedented number of states, including the spectrum of exotics below around 4.5 GeV, and identify their continuum quantum numbers. In particular, identification of the will be vital in order to compare with experiment and so to disentangle the nature of the many puzzling resonances. Preliminary results from this work have been presented in Ref. [16].

2 The anisotropic lattice action

This study aims to determine the spectrum of charmonium states, including excitations and any states with an intrinsic gluonic component. In Euclidean space-time, this spectrum can be computed by observing the behaviour of correlation functions of appropriately constructed operators, which fall exponentially as the time separation is increased with rate in proportion to the energies of the eigenstates of the Hamiltonian. Since the charmonium spectrum lies above 3 GeV, these correlation functions decrease rapidly at large Euclidean time separations. Correlation functions of creation operators which excite high-lying states also suffer from substantial statistical variance in Monte Carlo determinations. The anisotropic lattice provides a very useful framework to ameliorate these technical problems, as it enables the correlator to be determined with a finer temporal resolution at a manageable computational cost.

The anisotropic lattice offers a further advantage when studying charmonium physics. The charm quark mass, , has proved problematic in simulations in the past. It is too light for non-relativistic approaches to be valid whilst it is also large enough that only for the finest isotropic ensembles currently available. A number of successful approaches have been developed, including the Fermilab and HISQ actions [2, 17]. For this study, we will use a discretisation in which the spatial lattice spacing, , and temporal lattice spacing, , are related via . The fine temporal discretisation will ensure that we carry out our calculation with and the standard relativistic formulation of fermion actions can be used to study states where the typical spatial momentum of a charm quark inside a hadron at rest is small. This should be a reasonable approximation for the spectroscopy measurements we make in the current work. Nevertheless, in the simulations carried out here, is also less than unity. We explore the size of errors in Section 5.3.

2.1 The lattice action

For this study, we make use of the substantial ensembles of anisotropic gauge-field configurations generated by the Hadron Spectrum Collaboration [18, 19]. The gauge action is tree-level Symanzik-improved while in the fermion sector we use the anisotropic Shekholeslami-Wohlert action with tree-level tadpole improvement and three-dimensional stout-link smearing [20] of gauge fields. The fields are drawn from the importance sampling distribution appropriate for studies with two dynamical light quarks and a dynamical strange quark (); more details are given in Refs. [18, 19].

The charm quark action is the same as that used for the light and strange quarks, except that the parameter describing the weighting of spatial and temporal derivatives in the lattice action is chosen to give a relativistic dispersion relation for low-momentum mesons containing charm quarks. This non-perturbative determination [18] is described in more detail in Section 2.2. The charm quark fields are not treated dynamically in the configuration generation stage and so the quenching of these fields introduces a systematic uncertainty in our analysis. We expect the effects on the spectrum of this non-unitary treatment of charm quarks to be fairly small. The bare mass parameter for the valence charm quarks is determined by ensuring the ratio of the meson and the baryon masses takes its experimental value.

To quote energies and lengths in physical units we follow Ref. [19] and consider the ratio of the mass of the baryon measured on these ensembles,  [14], to the experimental mass,  MeV [21]. Setting the scale in this way, we find the parameters in the lattice action used in this study correspond to a spatial lattice spacing of  fm and a temporal lattice spacing approximately 3.5 times smaller, GeV. The lattices employed here have . Two volumes are used, and , corresponding to spatial extents of fm and fm respectively. The ensembles used in this study are summarised in Table 1 and full details of the lattice actions and parameters are given in Refs. [18, 19].

Volume MeV
396 96 128 64
396 552 32 162
Table 1: The gauge-field ensembles used in this work. and are respectively the number of gauge-field configurations and time-sources per configuration used; refers to the number of eigenvectors used in the distillation method [22].

2.2 Determining the action parameters for the charm quark

On the anisotropic lattice, the weights of any action terms that represent temporal and spatial derivatives must be determined in order that a physical observable takes its experimental value when calculated in the simulation. This is of particular importance for dynamical simulations to ensure the correct continuum limit [23]. For these ensembles where the target ratio of scales is , this determination for the gauge fields along with the dynamical light and strange quarks is presented in Ref. [18]. This ratio has been measured in the light meson sector from the pion dispersion relation, giving  [24].

As part of this study, we determined the dispersion relation for the meson on the ensemble. The dispersion relation for meson is


and since all quark fields making up the meson obey periodic boundary conditions for the finite spatial extent of the lattice, the momentum values are quantized, so with . Since we plan to follow up this work with a study of the scattering of mesons including charm quarks and hope to investigate the physics of states above the open-charm threshold, it is important to check that the dispersion relation is correct and that is consistent with that in the light quark sector.

In the left panel of Fig. 1, the dependence of the meson energy with increasing quantized momenta, labelled by , is plotted, along with the best straight-line fits. The triangular symbols show the dispersion relation when the bare parameter, , in the charm quark action is the same as that used for the light sea quarks,  [18]. From this data, the measured ratio of scales is , a discrepancy from the target value of about .

To correct for this discrepancy, simulations were performed to determine the charm-quark action parameters that recover a consistent dispersion relation. Increasing the bare parameter to , gives the second dispersion relation denoted by squares in the left panel of Fig. 1. The fit now yields a measured value of , consistent with the target value. Note that the statistical uncertainty on this determination is below the percent-level. We employ this in the charm quark action in all subsequent computations.


*[width=0.48]PLOTS/dispersion_etac.pdf \includegraphics*[width=0.48]PLOTS/dispersion_heavylight.pdf

Figure 1: Left panel: the dispersion relation of the meson measured on 96 configurations for two choices of the bare parameter, , in the charm quark action. As described in the text, is the parameter in the light sea-quark action and is determined by ensuring the correct dispersion relation for the meson. Right panel: the dispersion relation for the heavy-light meson measured on 39 configurations with in the charm quark action. Both on volumes.

As mentioned above, our future plans include a study of scattering of mesons, so it is important to check whether the self-consistent action parameters presented above also result in a relativistic dispersion relation for the meson. The light quark action parameters are left unchanged from those in the sea. This dispersion relation is presented in the right panel of Fig. 1 and the resulting slope yields a determination of the ratio of scales as , just different from the target value. The meson mass extracted from this analysis is MeV, but note that the simulated light quark mass is greater than its physical value,  MeV.

We also note that the data in each of the figures agree with the relativistic dispersion relation up to . All five momenta are used in the linear fits which have .

3 Operator construction and correlator analysis

In lattice calculations, meson masses are extracted from two-point correlation functions,


where the interpolating operators create the state of interest at and annihilate the state at a later Euclidean time . Inserting a complete set of eigenstates of the Hamiltonian we obtain the spectral decomposition,


where the vacuum-state matrix elements,


are referred to as overlaps, and the sum is over a discrete set of states because the calculation is performed in a finite volume. It is essential that the operators overlap well with the states under consideration. A well-established way to achieve this is to apply a smearing to the quark fields in the operators to emphasise the relevant low-energy modes.

We use the distillation technique, described in Ref. [22], which implements a smearing and provides an efficient method to calculate correlation functions with large bases of operators. Furthermore, it allows us to project onto definite momentum at both the source and the sink with reasonable computational cost. Distillation defines a smearing function,


where in this work and are the eigenvalues and eigenvectors of the gauge-covariant lattice Laplacian, , evaluated on the background of the spatial gauge-fields of time-slice and sorted by eigenvalue. Here is some smearing profile and in this work we set . The distillation operator, , is a projection operator into the subspace spanned by these eigenmodes and so . Smeared quark fields are constructed by applying this distillation operator to each quark field appearing in the interpolating operators, . As described in Ref. [22], the outer product structure of allows a correlation function to be factorised. The computational cost is manageable because inversions are performed in the smaller space of distillation vectors and spin, with dimensionality , where is the number of components of a lattice Dirac spinor. Further, this factorisation allows the efficient computation of correlation functions with a large basis of operators, as well as three-point functions and multi-hadron two-point correlation functions. All the quark fields referred to subsequently are smeared using the distillation operator and we denote them by .

The simplest interpolating operators are smeared local fermion bilinears of the form , where is any product of Dirac gamma matrices and for clarity we have suppressed spin, flavour and colour indices. However, such operators do not give access to or even to the complete set of quantum numbers for . In addition, as will be discussed later, we require a large basis of operators in each quantum number channel in order to reliably and precisely extract the energies of excited states.

To give us a large basis of operators with a range of spatial structures in each channel, we use the derivative-based construction for operators described in Refs. [10, 11]: gauge-covariant spatial derivatives are combined with a gamma matrix within a fermion bilinear. The operator is of the general form


where is a lattice-discretised gauge-covariant derivative. As discussed in Ref. [11], operators with definite can be constructed using these “forward-backward” derivatives. The derivatives and vector-like gamma matrices are first expressed in the circular basis so that they transform as . For creation operators the basis is


Using the standard Clebsch-Gordan coefficients, we can combine any number of derivatives with gamma matrices to construct operators with the desired quantum numbers, , where is the component. The choice of and the way in which the derivatives have been coupled determine the parity, , and charge conjugation, , quantum numbers. The notation we use follows that of Ref. [11]: an operator contains a gamma matrix, , with the naming scheme given in Table 2, and derivatives coupled to spin , overall coupled to spin .

Table 2: Gamma matrix naming scheme.
Table 3: Continuum spins, , subduced into lattice irreps, , where dim is the dimension of the irrep.

In lattice QCD calculations the full three-dimensional rotational symmetry of an infinite volume continuum is reduced to the symmetry group of a cube. Instead of the infinite number of irreducible representations labelled by spin, , we instead have a finite number of lattice irreducible representations (irreps); the five single-cover irreps for states at rest are , , , and . Note that parity and any flavour quantum numbers such as charge conjugation are still good quantum numbers on the lattice. The distribution of the various components, , of a spin- state across the lattice irreps is known as subduction and in Table 3 we show the pattern of these subductions for . Because a state with is split across many lattice irreps and each lattice irrep contains an infinite tower of , it is not straightforward to identify the of states extracted in lattice computations. Refs. [10, 11] presented a method to address this problem and we discuss this further in Section 4.



Figure 2: Principal correlator fits according to Eq. \eqrefequ:fitprincorr for a selection of low-lying states in the irrep on the volume. Plotted are data and the fits for showing the central values and one sigma statistical uncertainties; the grey points are not included in the fits. From left to right the states have , , and . In all cases, for the used we obtain reasonable fits with .

Refs. [10, 11] also discussed how to construct operators which transform in a definite lattice irrep, , and row, , from the continuum operators, ,


where are subduction coefficients. In this work we include operators built from all possible combinations of gamma matrices with up to three derivatives and then subduced into lattice irreps. Table 4 lists the numbers of operators in each irrep.

12 6 13 5
4 6 5 5
18 26 22 22
18 18 22 14
14 12 17 9
Table 4: The number of operators used in each lattice irrep, . All combinations of gamma matrices and up to three derivatives are included.

We use the variational method [25, 26] to extract spectral information from the two-point correlation functions with the particular implementation described in detail in Ref. [11]. In brief, for each lattice irrep we form a matrix of correlators, , where and label operators in our basis for that irrep. The best extraction (in the variational sense) of the energies and overlaps then follows from solving a generalised eigenvalue problem,


where an appropriate reference time-slice is chosen as described in Refs. [11, 7], the energies are determined by fitting the dependence of the eigenvalues, , called the principal correlators, on , and the overlaps are related to the eigenvectors, .

We fit the principal correlators using the fit form,


where the fit parameters are , and . This form allows for a second exponential and empirically we find that the contribution of this second exponential decreases rapidly as we increase . Some example fits to the principal correlators are shown in Fig. 2. These are plotted with the dominant time-dependence due to state divided out and so we would see a horizontal line at 1 in the case where a single exponential dominates the fit.

The combination of the variational method, our operator constructions, the distillation technique, and these anisotropic lattice ensembles has been shown to be very effective in studies of excited light isovector [10, 11] and isoscalar [12] mesons and baryons [14, 15]. In the following section we review our method for identifying the continuum spin of extracted states and present its first application to charmonium spectroscopy.

4 Determining the spin of a state

In principle, the spin of a single-hadron state can be determined by extracting the spectrum at various lattices spacings and then extrapolating to the continuum limit. There, where full rotational symmetry is restored, energy degeneracies between different lattice irreps should emerge. For example, according to the pattern of subductions (Section 3), a spin-2 state would be expected to appear as degenerate energies within the and irreps. However, there are difficulties with this procedure. First, this requires high precision calculations with successively finer lattice spacings and so with increasing computation cost; this is not practical with currently available resources. Second, the continuum spectrum can exhibit degeneracies and this is particularly true in the charmonium system where hyperfine splittings are small. The question then arises as to how to identify, without infinite statistics, which degeneracies are due to the lattice discretisation and which are physical degeneracies.



Figure 3: The normalised correlation matrix, () with in the irrep on the volume. The operators are ordered such that those subduced from spin 1 appear first followed by spin 3 and then spin 4. The correlation matrix is observed to be approximately block diagonal.

An alternative approach, proposed in Refs. [10, 11], is to consider the overlaps, Eq. \eqrefequ:Z, using this additional information to determine the spin of extracted states at a single lattice spacing. Here we will only review the important points. The subduced operators constructed in Section 3 respect the symmetries of the lattice. However, we find that each operator, carries a “memory” of the continuum spin, , from which it is subduced. If our lattice is reasonably close to restoring continuum symmetry at the hadronic scale, we can expect that an operator subduced from spin will only overlap strongly onto those states of continuum spin  [27].

\includegraphics[height=0.25]PLOTS/Histo_T2mm.pdf \includegraphics[height=0.25]PLOTS/Histo_Emm.pdf \includegraphics[height=0.25]PLOTS/Histo_A2mm.pdf \includegraphics[height=0.25]PLOTS/Histo_A1mm.pdf


Figure 4: The overlaps, , of a selection of operators onto some of the lower-lying states in each of the lattice irreps, , on the volume; states are labelled by their mass in lattice units, . The -values are normalized so that the largest value for that operator across all states is equal to 1. The lighter area at the top of each bar represents the one sigma statistical uncertainty on either side of the mean. Different color shading represents different operators.

This is apparent in both the correlation matrix and the overlaps of individual states. Fig. 3 shows the correlation matrix for at time-slice 5 where the operators are ordered according to the spin from which they are subduced. It can be seen that the correlation matrix is approximately block diagonal. In Fig. 4 we show the overlaps of a selection of operators onto states in the irreps . The operators are colored black, are red, are green, are blue and are gold. Every state has a clear preference to overlap only onto those operators subduced from a single continuum spin. For clarity we only show a subset of the operators but this pattern is observed for our full operator basis.

To be more quantitative, we can compare the overlaps obtained in different irreps. The continuum operators have definite spin so that and therefore we have . The value of is common to different lattice irreps: because the subduction coefficients form an orthogonal matrix, the spectral decomposition will contain terms proportional to ; these do not depend on the subduced into, up to discretisation effects. This suggests that we can identify the spin of a state by comparing the -values obtained independently in different irreps. For example, a state created by a operator will have the same -value in each of the and irreps.

In Fig. 5 we show the -values for states conjectured to have , 3 and 4 in the irreps . It can be seen that there is generally good agreement between the -values extracted in different irreps although there are some deviations from exact equality. The various discretisation effects have been discussed previously [11].

\includegraphics[height=0.36]PLOTS/Z_values_2mm.pdf \includegraphics[height=0.36]PLOTS/Z_values_3mm.pdf \includegraphics[height=0.36]PLOTS/Z_values_4mm.pdf
Figure 5: A selection of -values for states conjectured to have in irreps on the volume. From left to right, the operators in the left panel are , , , , , ; the operators in the middle panel are , , , , ; the operator in the right panel is .

For the states which have components distributed across different irreps, the question then arises as to what we should quote as our estimate of the mass of the state. The mass obtained from fits to principal correlators in each irrep can be different due to discretisation effects and fitting fluctuations. To minimize the fluctuations caused by changes in the fitting of the principal correlators we perform a joint fit to the principal correlators [11]. There is a common mass, , in the fit but we allow the second exponential in the principal correlator [Eq. \eqrefequ:fitprincorr] in each irrep to be different and so the fit parameters are , and . Our final results are determined from such joint fits and we find they are generally reasonable, with . In Fig. 6, as an example, we present the joint fit of the lightest state in the four irreps , , and .



Figure 6: A simultaneous fit to the four principal correlators of the lightest state on the volume using a common mass, . Plotted are data and the fit; the grey points are not included in the fit.

In summary, we have demonstrated that the -values of our carefully-constructed subduced operators can be used to identify the continuum spin of the extracted states. This was shown in the light meson sector [11] and here we have tested the method in an explicit calculation of the charmonium spectrum. When we present our final spin-assigned spectra, for states we use the results of joint fits to the principal correlators.

5 Stability of the extracted spectrum

In this section we present the results of some systematic tests of the stability of the extracted spectrum. We consider variations of the reference time-slice, , used in the variational method, the number of distillation vectors and changes to the spatial clover coefficient, .

5.1 Variational analysis and

In our implementation of the variational method the choice of a suitable is guided by how well the correlator is reconstructed, as described in Ref. [11]. In this study we find that typically is appropriate, depending on the irrep being considered.

Fig. 7 shows the variation with of the lightest ten extracted masses in the irrep on the ensemble. It can be seen that for large enough but not too large , the extracted spectrum is stable under variations in . We suggest that this stability is in a large part due to the inclusion of a second exponential term in Eq. \eqrefequ:fitprincorr. This term can ‘mop up’ other states which leak into the principal correlator because our operator basis is finite. We find that typically the mass in this second exponential, , is larger than that of all of the extracted states. As increases, the contribution of this second exponential falls rapidly due to it having a smaller amplitude, , and a larger .



Figure 7: Extracted mass spectrum as a function of for the lower-lying states in the irrep on the volume; the top left pane shows the lightest ten states and the other panes show the lightest four states in more detail. Horizontal bands show the masses extracted at with the width from the mean. For large enough but not too large the spectrum is seen to be stable under changes in .

In Figs. 8 and 9 we show the extracted overlaps as a function of for a selection of the states in the irrep. We observe that the overlaps show somewhat more sensitivity to than the masses. This might be expected because the -values are related to the eigenvectors of the generalised eigenvalue problem; these are likely to be strongly effected by a ‘wrong’ orthogonality when the number of contributing states is larger than the number of operators in the basis.



Figure 8: Extracted overlaps, , with operator as a function of for the lightest three states in the irrep on the volume. The -values for the 1st and 2nd excited states have been arbitrarily scaled by factors of 3 and 1.25 respectively to fit on the plot. Coloured bands show fits to an exponential plus a constant or a constant over an appropriate range of . The -values are seen to plateau for sufficiently large but can show significant curvature at small .


Figure 9: Extracted overlaps, , with operator as a function of for the lightest state in the irrep on the volume. The horizontal band shows a fit to a constant over an appropriate range of . The -value is seen to plateau for large but shows significant curvature at small .

In summary, we conclude that our implementation of the variational method gives a reliable extraction of the spectrum if is chosen appropriately. The masses extracted show very little variation with whereas the -values can show a more significant dependence on . This was also observed in a study of the spectrum of light isovector mesons [11].

5.2 Varying the number of distillation vectors

As described in Section 3, in the distillation method there is a choice of how many eigenvectors of the Laplacian, , to include in Eq. \eqrefequ:distillation. Using fewer distillation vectors reduces the computational cost but using too few results in over smearing and limits our ability to extract high-spin and excited states [11]; the unsmeared limit corresponds to where 3 corresponds to the number of colours. To get the same smearing operator on a larger volume requires to be scaled by a factor equal to the ratio of the spatial volumes [22]. Therefore, to optimise the efficiency of these calculations, it is important to use the smallest number of distillation vectors for which the states of interest can be extracted reliably.



Figure 10: Extracted spectrum for the irrep on the volume as a function of the number of distillation vectors, , using lower statistics than the main results. Horizontal bands correspond to the masses extracted with and give the 1 sigma range on either side of the mean. The spectrum is observed to be stable for .

In Fig. 10 we show the low-lying states in the irrep on the volume as a function of the number of distillation vectors used. These correlators were calculated using lower statistics, 90 gauge-field configurations and 6 time-sources per configuration, compared to the main results and so the extracted spectrum at is not identical to that in the full results presented later. The spin of some states for could not be reliably identified and these are not included in the figure. The spectrum is reasonably stable for , although the statistical precision reduces slightly with compared to . The quality of the spectrum degrades for fewer vectors and this is particularly apparent for the excited and high-spin states. A possible explanation for this behaviour was given in Ref. [11].

In summary, the number of distillation vectors must be sufficiently large in order to reliably extract high-spin and excited states, and once the number of vectors is large enough the spectrum is stable under changes in the number. These results are consistent with those of Ref. [11]. The observation that vectors are sufficient on the volume suggests that vectors is reasonable for the volume.

5.3 effects and the hyperfine splitting

An accurate determination of the hyperfine splitting between the ground-state pseudoscalar and vector mesons made from two heavy quarks has been a long-standing problem for lattice investigations. The determination of this energy is known to be delicately sensitive to artefacts arising from the discretisation of the Dirac operator on the lattice. The result for the splitting determined on our ensemble using tree-level improvement is 80(2) MeV, consistent with that obtained on the ensemble and substantially undershooting the physical value of 117 MeV. This determination neglects the effects of the disconnected Wick graphs contributing to the correlation function of charmonium mesons. These diagrams are not expected to change the hyperfine splitting substantially [28] due to OZI suppression and their absence is unlikely to explain in full the discrepancy noted above. We have made some first tests of the effect of disconnected diagrams in our calculation and come to the same preliminary conclusion. One other source of systematic uncertainty is the use of up and down quarks heavier than their physical values in our lattice action. Other lattice calculations [8] suggest this will also result in only a small effect which will not explain the discrepancy we see.



Figure 11: The mass shifts, measured on the ensemble, for the lowest-lying states in a selected set of charmonium lattice irreps as , the size of the chromomagnetic clover term, is varied from tadpole-improved tree-level value of (green) to (cyan). Masses presented here are measured relative to . Experimental data is shown as solid lines; note that experimentally the has a significant hadronic width.

The lattice action for charm quark fields used in this study is -improved at tree level in tadpole-improved perturbation theory. It would be expected that a non-perturbative determination [29] of action parameters would yield a larger value of , the coefficient of the spatial Pauli-like term, , appearing in the lattice action, than the value we use. Studies show that increasing the value of yields a larger hyperfine splitting at finite lattice spacing. To investigate the effect of a larger value of , a study of some charmonium states using was carried out on the volume. To be clear, the choice of is an ad hoc one, chosen from the consideration that a non-perturbative determination is expected to be larger than the tadpole-improved tree-level value which in our simulations is .

The lowest lying , , and states were computed, corresponding to the lightest S and P-wave states in a quark model and the lightest exotic, and Fig. 11 shows these masses with the -mass subtracted. Choosing gives a mass difference between the and of 114(2) MeV, very close to the experimental value and also leads to significantly better agreement in the P-wave system. The spin-two state is split on the lattice across two irreps, and . For both values of , these two states are degenerate within the statistical precision of this study. This is still consistent with the explanation that the hyperfine splitting is underestimated due to lattice artefacts. Since all dimension-five operators consistent with lattice symmetries do not break the continuum rotation group, any differences between these lattice states should only arise at in a Symanzik expansion and so would be expected to be small.

The lightest charmonium exotic, the state, shows mild dependence on . This suggests the energies of such states in our study should be close to their continuum values. Note that the results presented later in Section 6 use the tree-level tadpole-improved value of , not the boosted version of this test. This computation does however enable us to assign an approximate scale of 40 MeV to the size of the systematic uncertainty arising from the leading-order correction and observe that it does not increase for the higher-lying states in the spectrum where this has been tested. However, a more detailed investigation of these discretisation effects has not yet been performed and would be required for more complete control of this systematic uncertainty.

6 Results

In this section we present the results from our lattice calculations. We first give the results by lattice irrep and compare the spectra obtained on the two lattice volumes. We then show the final results with states labelled by their continuum quantum numbers, . Following on from this, in Section 7 we interpret the results and discuss the hybrid phenomenology suggested by them.



Figure 12: Extracted spectra by lattice irrep, , on the (lighter shading) and (darker shading) volumes. The vertical size of each box gives the one sigma statistical uncertainty on either side of the mean and the colour coding indicates the continuum spin.

6.1 Results by lattice irrep and volume comparison

In Fig. 12 we show the results of the variational analysis for each lattice irrep, , with the spectra obtained from the two volumes side by side. The one sigma statistical uncertainty on either side of the mean, as determined from the variational analysis, is indicated by the vertical size of the boxes. The colour coding gives the continuum spin determined by considering the operator overlaps as described in Section 4. States identified as are coloured black, are red, are green, are blue and are gold.



Figure 13: The extracted spin-identified charmonium spectrum labelled by ; the (open boxes) and (filled boxes) volumes agree well. The vertical size of the boxes represents the one sigma statistical uncertainty on either side of the mean. The dashed lines indicate the lowest non-interacting and levels using the and masses measured on the ensemble.

By comparing the overlaps between different irreps, we are able to match levels which correspond to the different components of a state with definite continuum spin subduced into lattice irreps; generally we see no significant discrepancy between the masses determined in different irreps. The dense spectrum of states at and above would appear to be impossible to disentangle using solely the extracted masses. An advantage of using anisotropic lattices is that all masses we consider have , which is not easily achieved for charmonium on isotropic lattices. The spatial lattice spacing is  fm and this is fine enough for an effective restoration of rotational symmetry at the hadronic scale, shown by the fact that we can successfully determine the continuum spin of extracted states.

In Fig. 13 we show the spin-identified spectrum labeled by continuum with the results from the two volumes presented side by side; only the well-determined low-lying states are included. The dashed lines indicate the lowest thresholds for open charm decay, the non-interacting and levels with the and masses measured on the ensemble. We observe that there is no significant difference between the two volumes at the level of the statistical uncertainties, even above the open-charm thresholds, evidence that we are not seeing multi-hadron states in our extracted spectra. We used higher statistics for the calculation compared to the calculation (see Section 2.1). On the volume we can reliably extract and identify the spin of more states and determine their masses with higher statistical precision, and we will therefore quote these as our final results.

In summary, we are able to interpret our spectra in terms of single mesons of definite continuum spin subduced into various lattice irreps. This, along with the lack of any significant volume dependence, is evidence that we are not observing multi-hadron states in our extracted spectra; the distribution of such states across irreps would have a different pattern. Further evidence is provided by observing that we do not extract energy levels where we would expect non-interacting two-meson levels to appear. These points have been discussed in detail in Ref. [11].

6.2 Final spin-identified spectrum



Figure 14: Summary of the charmonium spectrum up to around labelled by . The red and green boxes are the masses calculated on the volume; black lines are experimental values from the PDG [21]. We show the calculated (experimental) masses with the calculated (experimental) mass subtracted. The vertical size of the boxes represents the one sigma statistical uncertainty on either side of the mean. The dashed lines indicate the lowest non-interacting and levels using the and masses measured on the ensemble (fine green dashing) and using the experimental masses (coarse grey dashing).

Our final results, the well-determined states on the volume, are shown in Fig. 14 along with experimental masses taken from PDG summary tables [21]; these results are also tabulated in Appendix A. The is not shown because its ( or ) has not been determined experimentally. In the plot, the calculated (experimental) mass of the has been subtracted from the calculated (experimental) masses in order to reduce the systematic error from the tuning of the bare charm quark mass (see Section 2). The dashed lines indicate the lowest non-interacting and levels using the and masses measured on the ensemble (fine green dashing) and using the experimental masses (coarse grey dashing). We note that higher up in the spectrum the states become less well determined. This is in part because, although we have a relatively large number of operators in each channel, the basis is still restricted in size. In order to better determine these states we would need to include operators with different radial structures and, in order to determine states with , different angular structures, for example higher numbers of derivatives.

In the following section we interpret our results and compare with the experimental situation, but before that we discuss some other lattice calculations. Ref. [7] presented a pioneering quenched lattice QCD calculation of the excited charmonium spectrum which was interpreted further in Ref. [30]. In comparison, apart from being a dynamical calculation, the current work uses a larger basis of operators and these operators have been constructed so as to facilitate the identification of the continuum spin of the extracted states. This means that we can precisely extract a much larger number of excited states and states with high spin, and reliably identify their continuum .

Results from an calculation (dynamical up and down quarks but quenched strange quarks) of excited charmonia are presented in Ref. [9] using isotropic lattices with a few different lattice spacings and . They also consider mixing with light mesons in the pseudoscalar channel and mixing with some low-lying multi-mesons states in a few channels. In addition, preliminary results from dynamical () calculations of excited charmonium spectra have been presented in Ref. [8]. A range of pion masses and lattice spacings are considered and chiral and continuum extrapolations are attempted. However, both these references use a smaller operator basis compared to ours making a robust identification of the continuum spin of the extracted states difficult. Therefore a limited number of states can be reliably extracted and in particular the lightest exotic hybrid () is difficult to unambiguously disentangled from a non-exotic spin four state.

7 Interpretation of the results

In this section we give some interpretation of our results, compare with the experimental situation and then discuss hybrid meson phenomenology in more detail. Many of the states with non-exotic in Fig. 14 appear to follow the pattern predicted by quark potential models, where is the spin of the quark-antiquark pair, is the relative orbital angular momentum, is the total spin of the meson and is the radial quantum number. The quantum numbers of the members of various quark-antiquark supermultiplets are listed in Table 5. The assignments for our extracted states are determined by considering the operator-state overlaps [30].

Table 5: Supermultiplets for quark-antiquark pairs with spin and relative orbital angular momentum . Also shown are some hybrid supermultiplets, discussed in Section 7.2, where are the quantum numbers of the gluonic excitation; exotic are shown in bold.

In the left panel of Fig. 14 we show the negative parity states. We observe the ground-state S-wave pair and at 700 MeV the first excitation; a second excitation is observed at 1150 MeV. There is a complete D-wave set [, ] just above the threshold. In the region above 1200 MeV there is a complete set of D-wave excitations and parts of a G-wave indicated by the presence of spin-4 states. In addition, there are three states at 1200 to 1300 MeV with which do not appear to fit with the pattern expected by quark models. We note that these three states all have a relatively large overlap onto operators proportional to the gluonic field-strength tensor, something not observed for the states that fit into quark-antiquark supermultiplets. For example, this , the sixth state in the irrep in Fig. 4, is the only state shown which has a large overlap onto the operator, an operator proportional to the chromomagnetic part of the field-strength tensor. Following Ref. [11] we suggest that these ‘excess’ states can be identified as non-exotic hybrid mesons and we return to this in Section 7.2.

In the middle panel we show the positive parity states. Below threshold there is a clear P-wave set [, ]. In the region of 1000 MeV there is a P-wave excitation and, slightly higher, an F-wave [, ]. The band of states around 1400 MeV probably contains part of the second excitation of the P-wave set and several non-exotic hybrids which lie above the lightest negative-parity hybrids; we will discuss these hybrid candidates in Section. 7.2.

The states with exotic are presented in the right panel. These cannot consist solely of a quark-antiquark pair – more degrees of freedom are necessary. In general, there are several possible interpretations such as hybrid mesons where the gluonic field is excited, molecular mesons or tetraquarks. The exotic states in our spectrum have significant overlap onto operators with non-trivial gluonic structure and this suggests that we can identify them as hybrid mesons. Further, the lack of evidence for multi-hadron states in any channel, as discussed in Section 6, suggests that these exotics are not multi-hadron states. The is the lightest with 1200 MeV and nearly degenerate with the three non-exotic hybrid candidates [, ] observed in the negative parity sector. Higher in the spectrum there is a pair of states, , at 1400 MeV and a second state slightly higher again. In Section 7.2 we discuss the hybrid meson phenomenology suggested by our spectra.


[width=0.48]PLOTS/Z_values_Pwave.pdf \includegraphics[width=0.48]PLOTS/Z_values_Dwave.pdf

Figure 15: Overlaps, , for the lightest P-wave (left panel) and D-wave (right panel) supermultiplets with operators and respectively.

In order to further test the assignment of states into the lightest P and D-wave supermultiplets, we follow Ref. [31] and compare operator-state overlaps for different states within a supermultiplet. Consider the operators with and with , where and .1 These operators have the structure of a quark-antiquark pair in a gauge-covariant version of a P-wave with () or (). The operators with and with have the structure of a quark-antiquark pair in D-wave with or .

In Fig. 15 we show the overlaps of these operators with states suggested to be members of the lightest P-wave (left panel) and D-wave (right panel) supermultiplets, and observe that these overlaps are significant. More generally, we find that the overlap onto non-relativistic operator combinations [31, 30], e.g. and , are much larger than the complementary combinations, e.g. and , which would be zero in the non-relativistic limit. This might be expected since the charmonium system is fairly non-relativistic.

For states within a given supermultiplet, it is expected that the -values for each of these operators, projected into the relevant lattice irreps, will be similar because they correspond to the same underlying spatial wavefunction with differing internal spin and angular momentum couplings. This is seen to be the case in Fig. 15, supporting our identification of these states as members of the P and D-wave supermultiplets. We note that only qualitative conclusions can be drawn from the -values for lattice-regularised matrix elements; these require renormalisation to be compared with continuum matrix elements and this renormalisation can mix operators.

In summary, our extracted spectrum features states that follow the pattern expected by quark models along with others, having both exotic and non-exotic quantum numbers, that do not fit into this pattern. We have argued that these extra states can be identified as hybrid mesons and discuss their phenomenology below after we compare our results with the experimental situation.

7.1 Comparison with experimental results

It can be seen from Fig. 14 that the pattern of our extracted masses is in general agreement with the rather limited number of experimentally observed states. We note that states above threshold can have large hadronic widths and a conservative approach is to only consider our mass values accurate up to the hadronic width [11]. Discrepancies could be due to discretisation effects arising from the finite lattice spacing, systematic uncertainties within the variational method and because we have unphysically heavy up and down quarks. We do not expect the unphysical up and down quarks masses to introduce a large error in the determination of lower-lying states [8]. However, a more careful treatment of these effects is needed in the vicinity of open-charm thresholds [32]. In addition, we have used the tree level value for the spatial clover coefficient and, as discussed in Section 5.3, this means that we underestimate the S-wave hyperfine splitting. Further sources of systematic uncertainty arise from not including disconnected contributions, though these are expected to give only a small contribution in charmonium. We have shown in Section 6.1 that we do not see any finite volume effects at our level of statistical precision.

In the channel we identify the with our excited and the with our ground state . For these extracted states, we find that the admixture of in the dominantly state is small, as is the admixture of in the state. We could not extract a significant signal for the mixing angle.

Higher up in this channel there are two states in the PDG review, and , whose robustness is not beyond question. These states appear in the mass region where we have the second excitation of the with a relatively large uncertainty on the mass. In light of the observations in Section 5.2, a larger number of distillation vectors may be required in order to extract this state more precisely. Alternatively, we may need to add operators with more radial structures to our basis, for example by including operators with additional spatial derivatives or different smearing profiles.

The at 1300 MeV with is one of the states that appears to be supernumerary compared to the pattern expected by quark models. A possible interpretation is as a non-exotic hybrid meson [33, 34, 35]. The mass of the hybrid candidate in our calculation agrees well with the mass of the , supporting the hybrid interpretation. However, solely comparing masses does not allow strong conclusions to be drawn. We find conventional charmonia in the same mass region and from our calculation we also can not exclude the possibility that the is a multi-meson state or a tetraquark.

We find that the hybrid candidate has a significant overlap onto an operator which has the structure of a colour-octet quark-antiquark pair with in S-wave coupled to the gluonic field, . This is in contrast to the conventional mesons which have and is interesting because the decays to  [21]. The folklore was that such decays should conserve the spin of the heavy quarks and so, because the has , this would imply that the also has , at odds with our hybrid. However, the observation by CLEO [36] of a cross section for (the has ) comparable to that of at a center of mass energy of 4170 MeV suggests that this folklore may not be correct. Even so, the being produced in annihilation suggests that it has some admixture. In order to test the interpretation of the and draw firmer conclusions about the nature of these vector mesons, the unsmeared decay constants could be calculated for all the vector states we extract; this is beyond the scope of the present work. In Section 7.2 we discuss the more general hybrid phenomenology suggested by our results.

The at 890 MeV, very close to threshold, is the ‘unexpected’ state that has been confirmed by the largest number of experiments and perhaps is one of the most enigmatic; its structure remains unclear and is the subject of much discussion. Possible values for the are and . In our results, the D-wave state has around 30 MeV below the , while the first excitation of the P-wave is around 110 MeV above the . Another possible interpretation for the is as a molecular state of and mesons. As we have argued, we do not appear to see multi-hadron states in our extracted spectrum and so we would not expect to observe such a molecular state in this calculation. To pin down the nature of the we plan to supplement our basis with operators constructed to overlap well with multi-meson states [37, 24].

7.2 Exotic mesons and hybrid phenomenology



Figure 16: Charmonium spectrum up to around showing only channels in which we identify candidates for hybrid mesons. Red (dark blue) boxes are states suggested to be members of the lightest (first excited) hybrid supermultiplet as described in the text and green boxes are other states, all calculated on the volume. As in Fig. 14, black lines are experimental values and the dashed lines indicate the lowest non-interacting and levels.

In Fig. 16 we show the charmonium spectrum for the subset of channels in which, by considering operator-state overlaps, we identify candidate hybrid mesons. A state is suggested to be dominantly hybrid in character if it has a relatively large overlap onto an operator proportional to the commutator of two covariant derivatives, the field-strength tensor. We note that within QCD non-exotic hybrids can mix with conventional charmonia. We find that the lightest exotic meson has and is nearly degenerate with the three states observed in the negative parity sector suggested to be non-exotic hybrids, . Higher in mass there is a pair of states, , and a second state slightly above this. Not shown on the figures, an excited appears at around 4.6 GeV, there is an exotic state at around 4.8 GeV and the lightest exotic does not appear until above 5 GeV.

The observation that there are four hybrid candidates nearly degenerate with , coloured red in Fig. 16, is interesting. This is the pattern of states predicted to form the lightest hybrid supermultiplet in the bag model [38, 39] and the P-wave quasiparticle gluon approach [40], or more generally where a quark-antiquark pair in S-wave is coupled to a chromomagnetic gluonic excitation as shown Table 5. This is not the pattern expected in the flux-tube model [41] or with an S-wave quasigluon. In addition, the observation of two states, with one only slightly heavier than the other, appears to rule out the flux-tube model which does not predict two such states so close in mass. The pattern of of the lightest hybrids is the same as that observed in light meson sector [11, 31]. They appear at a mass scale of above the lightest conventional charmonia. This suggests that the energy difference between the first gluonic excitation and the ground state in charmonium is comparable to that in the light meson [31] and baryon [15] sectors.

To explore this hypothesis of a lightest hybrid multiplet further, we follow Ref. [31] and consider in more detail operator-state overlaps. The operators with and with are discussed in that reference. These operators have the structure of colour-octet quark-antiquark pair in S-wave with () or (), coupled to a non-trivial chromomagnetic gluonic field with where , and refer to the quantum numbers of gluonic excitation. Fig. 17 shows that the four states suggested to form the lightest hybrid supermultiplet have considerable overlap onto operators with this structure.



Figure 17: Overlaps, , for the proposed lightest hybrid supermultiplet with operators .

For states within a given supermultiplet, it is expected that the -values for each of these operators, projected into the relevant lattice irreps, will be similar as discussed above. The relevant overlaps presented in Fig. 17 suggest that the four hybrid candidates have a common structure, supporting our identification of them as the members of the lightest hybrid supermultiplet consisting of S-wave coupled to a gluonic excitation. We note again that only qualitative conclusions can be drawn from the -values for lattice-regularised matrix elements.

A heavier hybrid supermultiplet composed of a P-wave colour-octet quark-antiquark pair coupled to a gluonic field with would contain states with , , , , , , as shown in Table 5. In Ref. [31] this was identified as the first excited hybrid supermultiplet in the light meson sector. The exotic state and two states are observed in our spectrum with small mass splittings and these have relatively large overlap onto operators with the structure of a P-wave coupled to a gluonic . In the non-exotic positive-parity sector, we have evidence from similar operator overlaps that in the region around 1400 to 1500 MeV there are , , and hybrids as well as candidates for three hybrids. We therefore observe candidates for all expected members of this excited hybrid supermultiplet and these are coloured dark blue in Fig. 16. To pin these states down more precisely we would need to add operators with more derivatives and spatial structures to our basis; this is beyond the scope of the current work.

In summary, we have identified possible lightest and first excited hybrid supermultiplets and suggest that these can be interpreted in terms of pairs in, respectively, S-wave and P-wave, coupled to a chromomagnetic gluonic excitation with . The appearing at a much higher mass scale (above 5 GeV) may suggest a different type of gluonic excitation, for example a S-wave quasigluon, or have some other origin. We note that the inclusion of multi-hadron operators in our basis could modify the interpretation of states above open-charm threshold.

8 Summary

Using distillation and the variational method with a large basis of carefully constructed operators, we have successfully computed an extensive spectrum of charmonium mesons with high statistical precision and reliably identified the continuum of the extracted states. The large number of extracted states up to , across all combinations, goes beyond any other dynamical calculation. The spin identification methodology has allowed us to identify the continuum spin of the extracted states and this will be crucial in order to compare the higher-lying states against observations once this dense spectrum has been more fully mapped out by experiments.

We find no statistically significant volume-dependence and also that the spectrum is stable with respect to changes in the details of the variational analysis and variation of the number of distillation vectors. Lattice discretisation effects were investigated by changing the spatial clover coefficient; however, a more complete determination of this systematic uncertainty would require additional lattice spacings. The results presented here are for unphysically-heavy up and down quarks corresponding to a pion mass of at a single lattice spacing; lattice ensembles with lighter up and down quarks and further lattice spacings will be considered in future work. In our calculations we include only the connected quark line contributions to charmonium correlators. The disconnected contributions are expected to be small in charmonium since they are OZI suppressed. Nevertheless, in the distillation framework we can efficiently compute these contributions and determine their effect on the spectrum. This has already been done in the isoscalar light meson spectrum [12] and by applying the same methodology to charmonium we can determine the mixings between , and basis states; these calculations are ongoing [16].

For the first time, we have computed the dynamical spectrum of charmonia with exotic quantum numbers (, , ) up to 4.5 GeV. A determination of these exotic states is particularly interesting because it points to the presence of explicit gluonic degrees of freedom. The spectrum of non-exotic states has many features in common with the pattern expected by quark models along with some states which do not appear to fit into such a classification. We have suggested that these extra states can be interpreted as non-exotic hybrids and have identified the lightest hybrid supermultiplet consisting of states with , as well as an excited hybrid supermultiplet. This is consistent with the pattern observed in lattice calculations in the light meson sector [31] and can be interpreted as a colour-octet quark-antiquark pair coupled to a chromomagnetic gluonic excitation. We find that the lightest gluonic excitation has an energy scale of , comparable to that found in the light meson [31] and baryon [15] sectors, suggesting common physics. Our results allow an interpretation of the as a non-exotic vector hybrid meson, but based solely on comparing masses we can not draw a definitive conclusion; we have suggested further work to probe its structure.

We have presented arguments that, as in the study of light isovector mesons [11], we see no clear evidence for multi-hadron states. To study such states we plan to enlarge the basis of operators to include those with a larger number of fermion fields. In Refs. [37, 24] multi-hadron constructions have been presented and the method of Lüscher and its extensions [42, 43, 44, 45] was used to compute phase shifts from the denser spectrum of energy levels extracted. The efficacy of these constructions was shown in application to scattering in isospin-2. By applying this methodology to the charmonium sector we can determine the spectrum of multi-hadron states and then use Lüscher’s method and its extensions to compute phase shifts, and so determine the mass and width of resonances.

As well as the renaissance in charmonium spectroscopy, in recent years there has also been renewed interest in the spectroscopy of open-charm and mesons with unexpected states being observed. The methodology presented in this article can be applied to the study of these open-charm states and this work is in progress.

An important goal of the Hadron Spectrum Collaboration is a full understanding of charmonium spectroscopy including properties of resonances. Experiments such as BESIII, PANDA and those at the LHC will collect vast numbers of states to explore the charm spectrum above and below threshold, including states with intrinsic gluonic excitations. Our successful application of novel methods including distillation and spin-identification to determine the single-particle spectrum of charmonium up to GeV is a timely contribution to this effort.


We thank our colleagues within the Hadron Spectrum Collaboration. Chroma [46] and QUDA [47, 48] were used to perform this work on the Lonsdale cluster maintained by the Trinity Centre for High Performance Computing funded through grants from Science Foundation Ireland (SFI), at the SFI/HEA Irish Centre for High-End Computing (ICHEC), and at Jefferson Laboratory under the USQCD Initiative and the LQCD ARRA project. Gauge configurations were generated using resources awarded from the U.S. Department of Energy INCITE program at the Oak Ridge Leadership Computing Facility at Oak Ridge National Laboratory, the NSF Teragrid at the Texas Advanced Computer Center and the Pittsburgh Supercomputer Center, as well as at Jefferson Lab. This research was supported by the European Union under Grant Agreement number 238353 (ITN STRONGnet) and by the Science Foundation Ireland under Grant Nos. RFP-PHY-3201 and RFP-PHY-3218. CET acknowledges support from a Marie Curie International Incoming Fellowship, PIIF-GA-2010-273320, within the 7th European Community Framework Programme. JJD, RGE, BJ and DGR acknowledge support from U.S. Department of Energy contract DE-AC05-06OR23177, under which Jefferson Science Associates, LLC, manages and operates Jefferson Laboratory. JJD also acknowledges the support of the Jeffress Memorial Fund and the U.S. Department of Energy Early Career award contract DE-SC0006765.

Appendix A Table of results

In Table 6 we tabulate the masses presented in Fig. 14. The calculated mass of the has been subtracted in order to reduce the systematic error from the tuning of the bare charm quark mass (see Section 2).

/ MeV
0 663(3) 1143(13) 1211(13)
80.2(1) 698(6) 840(3) 1154(28) 1301(14) 1339(38)
860(3) 1334(17) 1350(17)
859(5) 1333(18)
867(3) 1269(26) 1392(12)
461.6(7) 972(9) 1361(46) 1488(30)
534(1) 1006(9) 1360(38) 1462(51) 1493(19) 1513(39)
521.6(9) 1002(10) 1415(14) 1484(48)
554(1) 1041(12) 1112(8) 1508(21)
1142(6) 1564(22)
1411(40) 1525(18)
Table 6: Summary of the charmonium spectrum calculated on the ensemble, as presented in Fig. 14, with statistical uncertainties shown. The scale has been set using the -baryon mass as described in Section 2; we note that MeV in these simulations. The calculated mass has been subtracted in order to reduce the systematic error from the tuning of the bare charm quark mass. The masses of states with are from joint fits to principal correlators as described in Section 4.


  1. these are for creation operators and we correct a typo appearing in the definition in Ref. [31]


  1. N. Brambilla et al., Heavy quarkonium: progress, puzzles, and opportunities, Eur. Phys. J. C71 (2011) 1534, [arXiv:1010.5827].
  2. HPQCD Collaboration, E. Follana et al., Highly Improved Staggered Quarks on the Lattice, with Applications to Charm Physics, Phys. Rev. D75 (2007) 054502, [hep-lat/0610092].
  3. T. Burch et al., Quarkonium mass splittings in three-flavor lattice QCD, Phys. Rev. D81 (2010) 034508, [arXiv:0912.2701].
  4. Namekawa, Y. for the PACS-CS Collaboration, Charm quark system on the physical point in 2+1 flavor lattice QCD, PoS(Lattice 2011) 132, [arXiv:1111.0142].
  5. D. Mohler and R. M. Woloshyn, and meson spectroscopy, Phys. Rev. D84 (2011) 054505, [arXiv:1103.5506].
  6. C. DeTar, Charmonium Spectroscopy from Lattice QCD, arXiv:1101.0212.
  7. J. J. Dudek, R. G. Edwards, N. Mathur, and D. G. Richards, Charmonium excited state spectrum in lattice QCD, Phys. Rev. D77 (2008) 034501, [arXiv:0707.4162].
  8. G. Bali et al., Spectra of heavy-light and heavy-heavy mesons containing charm quarks, including higher spin states for , PoS(Lattice 2011) 135, [arXiv:1108.6147].
  9. G. S. Bali, S. Collins, and C. Ehmann, Charmonium spectroscopy and mixing with light quark and open charm states from lattice QCD, Phys. Rev. D84 (2011) 094506, [arXiv:1110.2381].
  10. J. J. Dudek, R. G. Edwards, M. J. Peardon, D. G. Richards, and C. E. Thomas, Highly excited and exotic meson spectrum from dynamical lattice QCD, Phys. Rev. Lett. 103 (2009) 262001, [arXiv:0909.0200].
  11. J. J. Dudek, R. G. Edwards, M. J. Peardon, D. G. Richards, and C. E. Thomas, Toward the excited meson spectrum of dynamical QCD, Phys. Rev. D82 (2010) 034508, [arXiv:1004.4930].
  12. J. J. Dudek et al., Isoscalar meson spectroscopy from lattice QCD, Phys. Rev. D83 (2011) 111502, [arXiv:1102.4299].
  13. J. Bulava et al., Nucleon, and excited states in lattice QCD, Phys. Rev. D82 (2010) 014507, [arXiv:1004.5072].
  14. R. G. Edwards, J. J. Dudek, D. G. Richards, and S. J. Wallace, Excited state baryon spectroscopy from lattice QCD, Phys. Rev. D84 (2011) 074508, [arXiv:1104.5152].
  15. J. J. Dudek and R. G. Edwards, Hybrid Baryons in QCD, Phys.Rev. D85 (2012) 054016, [arXiv:1201.2349].
  16. L. Liu, S. M. Ryan, M. Peardon, G. Moir, and P. Vilaseca, Charmonium spectroscopy from an anisotropic lattice study, PoS(Lattice 2011) 140, [arXiv:1112.1358].
  17. A. X. El-Khadra, A. S. Kronfeld, and P. B. Mackenzie, Massive Fermions in Lattice Gauge Theory, Phys. Rev. D55 (1997) 3933–3957, [hep-lat/9604004].
  18. R. G. Edwards, B. Joo, and H.-W. Lin, Tuning for Three-flavors of Anisotropic Clover Fermions with Stout-link Smearing, Phys. Rev. D78 (2008) 054501, [arXiv:0803.3960].
  19. H.-W. Lin et al., First results from 2+1 dynamical quark flavors on an anisotropic lattice: light-hadron spectroscopy and setting the strange-quark mass, Phys. Rev. D79 (2009) 034502.
  20. C. Morningstar and M. J. Peardon, Analytic smearing of SU(3) link variables in lattice QCD, Phys. Rev. D69 (2004) 054501, [hep-lat/0311018].
  21. Particle Data Group Collaboration, K. Nakamura et al., Review of particle physics, J.Phys.G G37 (2010) 075021. and 2011 partial update for the 2012 edition.
  22. Hadron Spectrum Collaboration, M. Peardon et al., A novel quark-field creation operator construction for hadronic physics in lattice QCD, Phys. Rev. D80 (2009) 054506, [arXiv:0905.2160].
  23. R. Morrin, A. Ó Cais, M. Peardon, S. M. Ryan, and J.-I. Skullerud, Dynamical QCD simulations on anisotropic lattices, Phys. Rev. D74 (2006) 014505, [hep-lat/0604021].
  24. J. J. Dudek, R. G. Edwards, and C. E. Thomas, S and D-wave phase shifts in isospin-2 pi pi scattering from lattice QCD, arXiv:1203.6041.
  25. C. Michael, Adjoint Sources in Lattice Gauge Theory, Nucl. Phys. B259 (1985) 58.
  26. M. Luscher and U. Wolff, How to calculate the elastic scattering matrix in two- dimensional quantum field theories by numerical simulation, Nucl. Phys. B339 (1990) 222–252.
  27. Z. Davoudi and M. J. Savage, Restoration of Rotational Symmetry in the Continuum Limit of Lattice Field Theories, arXiv:1204.4146.
  28. L. Levkova and C. DeTar, Charm annihilation effects on the hyperfine splitting in charmonium, Phys. Rev. D83 (2011) 074504, [arXiv:1012.1837].
  29. M. Luscher, S. Sint, R. Sommer, P. Weisz, and U. Wolff, Non-perturbative O(a) improvement of lattice QCD, Nucl. Phys. B491 (1997) 323–343, [hep-lat/9609035].
  30. J. J. Dudek and E. Rrapaj, Charmonium in lattice QCD and the non-relativistic quark- model, Phys. Rev. D78 (2008) 094504, [arXiv:0809.2582].
  31. J. J. Dudek, The lightest hybrid meson supermultiplet in QCD, Phys. Rev. D84 (2011) 074023, [arXiv:1106.5515].
  32. F.-K. Guo and U.-G. Meissner, Light quark mass dependence in heavy quarkonium physics, arXiv:1203.1116.
  33. S.-L. Zhu, The possible interpretations of Y(4260), Phys. Lett. B625 (2005) 212, [hep-ph/0507025].
  34. F. E. Close and P. R. Page, Gluonic charmonium resonances at BaBar and Belle?, Phys. Lett. B628 (2005) 215–222, [hep-ph/0507199].
  35. E. Kou and O. Pene, Suppressed decay into open charm for the Y(4260) being an hybrid, Phys. Lett. B631 (2005) 164–169, [hep-ph/0507119].
  36. CLEO Collaboration, T. K. Pedlar et al., Observation of the using collisions above DDbar threshold, Phys. Rev. Lett. 107 (2011) 041803, [arXiv:1104.2025].
  37. J. J. Dudek, R. G. Edwards, M. J. Peardon, D. G. Richards, and C. E. Thomas, The phase-shift of isospin-2 pi-pi scattering from lattice QCD, Phys. Rev. D83 (2011) 071504, [arXiv:1011.6352].
  38. T. Barnes, F. E. Close, F. de Viron, and J. Weyers, Q anti-Q G Hermaphrodite Mesons in the MIT Bag Model, Nucl. Phys. B224 (1983) 241.
  39. M. S. Chanowitz and S. R. Sharpe, Hybrids: Mixed States of Quarks and Gluons, Nucl. Phys. B222 (1983) 211.
  40. P. Guo, A. P. Szczepaniak, G. Galata, A. Vassallo, and E. Santopinto, Heavy quarkonium hybrids from Coulomb gauge QCD, Phys. Rev. D78 (2008) 056003, [arXiv:0807.2721].
  41. N. Isgur and J. E. Paton, A Flux Tube Model for Hadrons in QCD, Phys. Rev. D31 (1985) 2910.
  42. M. Luscher, Signatures of unstable particles in finite volume, Nucl. Phys. B364 (1991) 237–254.
  43. K. Rummukainen and S. A. Gottlieb, Resonance scattering phase shifts on a nonrest frame lattice, Nucl.Phys. B450 (1995) 397–436, [hep-lat/9503028].
  44. C. Kim, C. Sachrajda, and S. R. Sharpe, Finite-volume effects for two-hadron states in moving frames, Nucl.Phys. B727 (2005) 218–243, [hep-lat/0507006].
  45. N. H. Christ, C. Kim, and T. Yamazaki, Finite volume corrections to the two-particle decay of states with non-zero momentum, Phys.Rev. D72 (2005) 114506, [hep-lat/0507009].
  46. SciDAC Collaboration, R. G. Edwards and B. Joo, The chroma software system for lattice qcd, Nucl. Phys. B. Proc. Suppl. 140 (2005) 832, [hep-lat/0409003].
  47. M. A. Clark, R. Babich, K. Barros, R. C. Brower, and C. Rebbi, Solving Lattice QCD systems of equations using mixed precision solvers on GPUs, Comput. Phys. Commun. 181 (2010) 1517–1528, [arXiv:0911.3191].
  48. R. Babich, M. A. Clark, and B. Joo, Parallelizing the QUDA Library for Multi-GPU Calculations in Lattice Quantum Chromodynamics, in International Conference for High Performance Computing, Networking, Storage and Analysis (SC), pp. 1–11, 2010. arXiv:1011.0024.