Excited and exotic charmonium, $D_s$ and $D$ meson spectra for two light quark masses from lattice QCD

We present highly-excited charmonium, $D_s$ and $D$ meson spectra from dynamical lattice QCD calculations with light quarks corresponding to $M_{\pi} \sim 240$ MeV and compare these to previous results with $M_{\pi} \sim 400$ MeV. Utilising the distillation framework, large bases of carefully constructed interpolating operators and a variational procedure, we extract and reliably identify the continuum spin of an extensive set of excited mesons. These include states with exotic quantum numbers which, along with a number with non-exotic quantum numbers, we identify as having excited gluonic degrees of freedom and interpret as hybrid mesons. Comparing the spectra at the two different $M_\pi$, we find only a mild light-quark mass dependence and no change in the overall pattern of states.


Introduction
The experimental status of the charm sector of Quantum Chromodynamics (QCD) has changed dramatically over the last decade [1]. The discovery of a plethora of unexpected charmonium-like states, commonly known as "X, Y, Z's", has highlighted the need for a more complete theoretical understanding of the spectrum. Many different interpretations have been put forward: some are suggested to be hybrid mesons (a quark-antiquark pair with excited gluonic degrees of freedom) and others two quarks and two antiquarks in a tightly-bound configuration (tetra-quark), a molecular-like combination of two mesons, or a charmonium-like core surrounded by light degrees of freedom (hadro-quarkonium). There are similar puzzles in the open-charm sector (D and D s mesons) where the measured masses and widths of the low-lying D * s0 (2317) ± and D s1 (2460) ± states are significantly smaller and narrower than expected from quark models. For some recent reviews see Refs. [2][3][4][5][6].
In principle these states can be understood within Quantum Chromodynamics (QCD) using lattice QCD, a non-perturbative, ab initio formulation of the theory. Spurred on by the experimental situation, there have been many lattice QCD calculations of hidden and open-charm mesons. The majority have focused on lowest-lying states below threshold, achieving unprecedented precision with the various systematic effects under control (some recent examples can be found in Refs. [7][8][9][10][11]). On the other hand, there have been a number of investigations of excited charmonia and open-charm mesons [12][13][14][15][16][17][18], all of which have some systematic uncertainties not fully accounted for and extract a more limited set of states than we consider here.
In a previous lattice QCD study, the Hadron Spectrum Collaboration used large bases of interpolating operators with various structures to robustly extract many excited and high-spin states and, crucially, to identify their continuum quantum numbers. Highlights included the presence of states with exotic quantum numbers (i.e. those forbidden with solely a quark-antiquark pair) and the identification of "supermultiplets" of hybrid mesons. However, these calculations were performed with unphysically-heavy light quarks corresponding to M π ∼ 400 MeV. The results provided useful benchmarks for other approaches such as nonrelativistic effective field theories, for example see Ref. [19].
The current work extends these earlier investigations by performing similar calculations with light-quark masses significantly closer to their physical values, corresponding to M π ∼ 240 MeV. The spectra at the two light quark masses are compared, focusing on the overall qualitative picture and, in particular, whether changes in the pattern of states with exotic quantum numbers or other hybrid mesons are observed. This allows us to explore the light-quark mass dependence of excited heavy quarkonia which has been suggested to be significant [20].
In this study the unstable nature of states above threshold is not considered -a point discussed in [21][22][23] -and so the spectra should only be considered a guide to the pattern of resonances. In the charm sector, we have already addressed this limitation for a variety of states appearing as bound-states and resonances in coupled-channel Dπ, Dη and D sK scattering [24] for M π ∼ 400 MeV and investigations of various other channels involving charm quarks are underway. This paper lays the foundation for extending those studies to M π ∼ 240 MeV, where the additional light-quark mass, closer to the physical value, will enable us to study the evolution with light-quark mass of hidden and open-charm bound-states and resonances.
The remainder of the manuscript is organised as follows. In Section 2 we describe the lattice ensembles used in this study, provide some details on the tuning of the anisotropy and charm-quark mass, and give a brief overview of the analysis of two-point correlation functions. In Section 3 we present and interpret the charmonium, D s and D meson spectra from the calculations with M π ∼ 240 MeV. In Section 4 we compare these spectra to those from earlier computations with M π ∼ 400 MeV and we present a summary in Section 5.

Calculation Details
In this study we use an anisotropic lattice formulation where the temporal lattice spacing, a t , is smaller than the spatial lattice spacing, a s ≈ 0.12 fm, with an anisotropy ξ ≡ a s /a t ≈ 3.5. The gauge sector is described by a tree-level Symanzik-improved anisotropic action, Lattice Volume M π (MeV) N cfgs N tsrcs for cc, cs, cl N vecs 24 3 × 128  391  553  32, 16, 16  162  32 3 × 256  236  484  1, 1, 2  384   Table 1. The lattice gauge field ensembles used. The volume is given as (L/a s ) 3 × (T /a t ) where L and T are respectively the spatial and temporal extents of the lattice. The number of gauge field configurations used, N cfgs , and the number of perambulator time-sources used per configuration, N tsrcs , are shown along with the number of eigenvectors used in the distillation framework [48], N vecs .
while the fermionic sector uses a tadpole-improved anisotropic Sheikholeslami-Wohlert (clover) action with stout-smeared gauge fields [49] and N f = 2 + 1 flavours of dynamical quarks. For both ensembles the heavier dynamical quark is tuned to approximate the physical strange quark, but the ensembles differ in light quark mass giving the two different pion masses. Table 1 summarises these lattice ensembles -full details are given in Refs. [50,51]. We use the same relativistic action for the charm quark as for the light and strange quarks (with tadpole-improved tree-level clover coefficients). The charm-quark mass and anisotropy parameters are tuned to reproduce the physical η c mass and a relativistic dispersion relation -this process was described for the M π ∼ 400 MeV ensemble in Ref. [22]. Throughout this work we do not correct experimental data for electromagnetic effects. For the M π ∼ 240 MeV ensemble, the momentum dependence of the η c energy after tuning is shown in Figure 1. The momentum is quantised by the periodic boundary conditions on the cubic spatial volume, p = 2π L n, where n = (n x , n y , n z ) and n i ∈ {0, 1, 2, . . . , L/a s − 1}. A reasonable fit to the dispersion relation, is obtained giving ξ ηc = 3.456(4), in agreement with the anisotropy measured from the pion dispersion relation on this ensemble, ξ π = 3.453(6) [52]. The fit gives M ηc = 2945 (17) MeV compared to the experimental value 2983.6(6) MeV [1], and so we estimate that the systematic uncertainty from tuning the charm-quark mass is of order 1%. Figure 1 also shows the momentum dependence of the D meson energy; a fit to Eq. 2.1 gives ξ D = 3.443(7), in reasonable agreement with ξ π and ξ ηc . To give results in physical units, we set the scale via a −1 t = M phys Ω /(a t M Ω ) using the Ω baryon mass measured on this ensemble, a t M Ω = 0.2789 (16) [52], leading to a −1 t = 5997 MeV. When quoting masses we reduce the already small systematic uncertainty from tuning the charm-quark mass by subtracting M ηc ( 1 2 M ηc ) from the mass of charmonia (open-charm mesons), rendering it negligible compared to other systematic uncertainties.
The aim of this work is to study how the spectra change as we vary the light-quark mass and only statistical uncertainties are given in the spectra we present in the following sections. While a full error budget is beyond the scope of this work, the uncertainties arising from working at a finite lattice spacing and in a finite volume were discussed in Ref. [22], where they were estimated to be small and have no overall qualitative effect on the spectrum. The uncertainty arising from the ambiguity in how to set the scale can be estimated by choosing a different reference observable. For example, setting the scale on the M π ∼ 240 MeV ensemble using the h c -η c mass splitting gives a −1 t = 5960 MeV, 0.6% lower than from using M Ω . On the other hand, using the η c (2S) -η c (1S) mass splitting gives a −1 t = 5787 MeV, 4% lower than when using the Ω baryon mass. Another source of systematic uncertainty comes from ignoring the unstable nature of states above threshold (see Refs. [21][22][23]). Although this is difficult to estimate, for a narrow resonance a conservative approach is to consider the uncertainty to be of the order of the width [21].

Calculation of spectra
We follow the methodology presented in Refs. [21][22][23] to compute the spectra. In brief, meson masses and other spectral information are obtained from the analysis of the time dependence of two-point Euclidean correlation functions, is the creation operator [annihilation operator] and t is the time separation. When computing charmonium correlators, disconnected Wick diagrams, where the charm quark and antiquark annihilate, are not included -these are OZI suppressed and so are expected to only give a small contribution in charmonium. There are no such disconnected contributions to the open-charm meson correlators considered here.
The hypercubic lattice has a reduced symmetry compared to an infinite volume continuum so states at rest are labelled by the irreducible representations (irreps), Λ, of the octahedral group, O h , rather than spin [53]. A method to ameliorate this issue and determine the continuum spin, J, of extracted states is given in Refs. [21][22][23] which also contain demonstrations of its efficacy. Parity, P , and any relevant flavour quantum numbers, e.g. charge-conjugation, C, are still good quantum numbers in our lattice formulation.
In each quantum-number channel, the distillation technique [48] is used to compute correlation functions involving a large basis of derivative-based fermion-bilinear interpolating operators [21]. 1 The resulting matrices of correlation functions, C ij (t), are analysed using a variational procedure [54][55][56] as described in Ref. [21]. This amounts to solving a generalised eigenvalue problem, where t 0 is a carefully chosen reference time-slice. For sufficiently large times, the eigenvalues, λ n (t, t 0 ), known as principal correlators, are proportional to e −Mn(t−t 0 ) where M n is the energy of the n th state. Energies are extracted from a fit to the form, ( , where the fit parameters are M n , A n and M n . The second exponential proves useful in stabilising the fit because it 'mops up' excited state contamination. The eigenvectors, v n j , are related to the operator-state overlaps (or matrix elements), Z i |0 , and contain information on the structure of a state -they are used in our method for determining the continuum spin. Figure 2 shows a selection of principal correlators from charmonium correlation functions in the Λ P C = T −− 1 irrep on the M π ∼ 240 MeV ensemble. The leading time dependence, e −Mn(t−t 0 ) , has been divided out yielding a plateau when a single exponential dominates. Beneath each principal correlator we show the overlap, Z, of each operator onto that state and below that, for comparison, the overlaps for the corresponding state on the M π ∼ 400 MeV ensemble. The operators were constructed to have definite J P C in the continuum: red bars correspond to J = 1, blue to J = 3 and yellow to J = 4. It is clear that each state is dominated by operators from a given J, demonstrating that the spin-identification methodology [21] can be used -this pattern is repeated for each of the spectra we determine. The darker shade of red and lighter shade of blue represent operators that are proportional to the spatial part of the field strength tensor, F ij . We identify a state as hybrid, i.e. a meson with excited gluonic degrees of freedom [21], when overlaps from these operators onto a given state are large compared to their overlaps onto other states 2 .

Charmonium and Open-Charm Spectra
In this section we present the spectra, labelled by J P (C) , computed on the M π ∼ 240 MeV ensemble. Results for charmonium are described first followed by those for D s and D mesons.

Charmonium
The charmonium spectrum computed on the M π ∼ 240 MeV ensemble is shown in Figure 3 and the results are tabulated in Appendix A. For flavour singlets such as charmonium, charge-conjugation, C, and parity, P , are both good quantum numbers and so states are Red parts of the curves show the time regions used in the fits; blue points were not included in the fits. Middle row: the operator-state overlaps, Z, for the state above, normalised so that the largest value for an operator across all states is equal to unity. Colour coding is described in the text and the error bars indicate the one sigma statistical uncertainty. Bottom row: overlaps for the corresponding state on the M π ∼ 400 MeV ensemble. labelled by J P C . As discussed above, masses are presented after subtracting the η c mass to reduce the systematic uncertainty arising from tuning the charm quark mass. Dashed lines indicate the location of some thresholds for strong decay: η c ππ (the lowest threshold if the charm quark and antiquark do not annihilate), DD and DD * . Since the resonant nature of states above threshold is not investigated in this work, a conservative approach is to only consider the mass values accurate up to the order of the hadronic width [21].
As found in Ref. [22], many of the states with non-exotic J P C follow the n 2S+1 L J pattern predicted by quark potential models, where J is the total spin of the meson with  [1]. As discussed in the text, we show the calculated (experimental) masses with the calculated (experimental) η c mass subtracted. The vertical size of the boxes represents the one-sigma statistical (or experimental) uncertainty on either side of the mean. Red and blue boxes correspond to states identified as hybrid mesons grouped into, respectively, the lightest and first-excited supermultiplet, as described in the text. Dashed lines show the location of some of the lower thresholds for strong decay using computed (coarse green dashing) and experimental (fine grey dashing) masses. relative orbital angular momentum L, quark-antiquark spin S and radial quantum number n. We find all states up to J = 4 expected by such models. Figure 3 also shows the states (coloured red and blue) that do not fit the n 2S+1 L J pattern. Four of these have exotic J P C quantum numbers, 0 +− , 1 −+ , 2 +− , and we find that they, as well as the excess states with non-exotic quantum numbers, have relatively large overlaps onto operators that are proportional to the spatial components of the field strength tensor, F ij (i.e. operators that have a non-trivial gluonic structure), something not seen for the other states in the spectrum. Furthermore, on removing operators proportional to F ij from the variational basis we generally observe a reduction in the quality of the signal for these states. We therefore follow Refs. [21,22] and interpret these excess states as hybrid mesons.
As noted in Section 2, these calculations are performed at a single spatial lattice spacing. On the 400 MeV ensemble we estimated a scale of 40 MeV for the discretisation uncertainty arising from O(a s ) corrections to charmonia [22]. Since the 240 MeV ensemble has the same spatial lattice spacing, we expect the 40 MeV scale to also be a reasonable estimate for the discretisation uncertainty here.

D s and D mesons
For flavoured mesons, such as D s and D, charge conjugation is no longer a good quantum number and states are labelled only by J P . Figures 4 and 5 show the D s and D meson spectra respectively; these results are tabulated in Appendix A. Masses are presented with half the mass of the η c subtracted in order to reduce the systematic uncertainty arising from tuning the charm quark mass. Dashed lines indicate some of the lower strong-decay thresholds (DK for the D s spectrum and Dπ and D * π for the D meson spectrum).
As for charmonium, the D s and D spectra can be interpreted in terms of a n 2S+1 L J pattern and we identify complete S, P, D and F -wave multiplets. Within the negative parity sector of both spectra, there are four states, highlighted in red, that do not appear to fit this pattern. Due to their relatively large overlap with operators featuring a nontrivial gluonic structure, these are identified as the members of the lightest hybrid meson supermultiplet in each flavour sector. The pattern is again consistent with a 1 +− gluonic excitation coupled to an S-wave quark-antiquark pair and they appear at an energy ∼ 1.2 -1.3 GeV above the lightest conventional multiplet. However, unlike in charmonium, the first excited hybrid supermultiplet is not robustly determined for open-charm mesons and is not shown here.

Comparison of the spectra at two light quark masses
The principal difference between the spectra presented in Refs. [22,23] and this work is the light quark mass, corresponding to M π ∼ 400 MeV in those references and M π ∼ 240 MeV here. Figures 6, 7 and 8 show comparisons of the charmonia, D s and D spectra at the two light quark masses -it can be seen that, in general, we observe only a mild light quark mass dependence throughout the entire spectra, with no change in the overall pattern of states. The systematic uncertainties were discussed in Section 2.
We note in passing that we achieve a greater statistical precision on the M π ∼ 400 MeV ensemble due to the larger number of time-sources used (see Table 1). In the discussion  that follows some notable features in each spectrum are highlighted and in Section 4.4 we Figure 6. Charmonium spectrum, labelled by J P C , with M π ∼ 240 MeV (left column for each J P C ) compared to the spectrum with M π ∼ 400 MeV from Ref. [22] (right column for each J P C ). As in earlier figures, red and blue boxes highlight states identified as constituents of, respectively, the lightest and first-excited supermultiplet of hybrid mesons. Dashed lines show some of the lower thresholds using computed masses for M π ∼ 240 MeV (coarse dashing) and M π ∼ 400 MeV (fine dashing): green is η c ππ, red is DD and blue is DD * .
investigate the mixing between spin-triplet and spin-singlet open-charm mesons.

Charmonium
In charmonium the light quark dependence enters through the sea quark content in the dynamical gauge field ensembles. As shown in Figure 6, for the low-lying states the masses are generally consistent between the two ensembles within statistical uncertainties. An exception is the hyperfine splitting, M J/ψ − M ηc , where we find a small but statistically significant increase when the light quark mass is decreased.
A second notable feature is that the masses of states higher up in the spectrum are generally larger on the M π ∼ 240 MeV ensemble. This is particularly the case for the hybrids, implying a small increase in their mass as M π is reduced; as a consequence the splitting between the hybrids and low-lying conventional mesons increases, albeit in a rather mild fashion. However, it is important to note that at higher energies the statistical uncertainties are larger and neglecting the unstable nature of states may be more important. We emphasise that the overall pattern of hybrid mesons is unaffected by decreasing the light quark mass.

D s mesons
As for charmonium, and shown in Fig. 7, only mild dependence on the light quark mass is observed throughout the D s meson spectrum. The largest change in the low-lying states is for the lightest 0 + (our candidate for the D * s0 (2317)). However, this state is expected to be heavily influenced by the nearby DK threshold to which it can couple in S wave, and interestingly, it has decreased just enough to remain below the threshold, in agreement with the experimental situation.
Once again we observe a tendency for the hybrid states, coloured red in Fig. 7, to increase in mass, and hence the splitting between the hybrids and the lowest conventional D s mesons to increase, as M π is reduced. However, there is no change to their overall pattern. Figure 8 shows that, in the D meson spectrum, the light quark mass dependence is also relatively mild and there is no change to the pattern of states. As expected, the masses generally decrease with decreasing pion mass -a D meson contains a valence light quark unlike a charmonium or D s meson. The most significant differences are observed for the lightest 0 + and 1 + states and this may be because they are strongly influenced by nearby thresholds that can couple in S wave, namely Dπ and D * π respectively. However, the mass of the second-lightest 1 + , which is also in the vicinity of the D * π threshold, does not change significantly. This may be because the mass difference between the charm quark and the light quark is large enough for the expectations of the heavy-quark limit to be a reasonable guide. In this limit, one of the 1 + states can decay to D * π only in S wave, whereas the other can decay to D * π only in D wave [63]; the latter would be expected to be influenced less by the position of the D * π threshold.

D mesons
At higher energies in the spectrum, there are generally only small or statistically insignificant mass shifts while, as for charmonia and D s mesons, there is a general trend for the hybrid mesons to become heavier as the light-quark mass decreases. This change is somewhat less clear in the D meson spectrum because of the opposing trend for mesons to become lighter as the light-quark mass decreases.

Mixing of spin-triplet and spin-singlet open-charm mesons
As discussed in Section 3.2, charge conjugation is not a good quantum number for opencharm mesons and, consequently, quark model spin-singlet ( 1 L J=L ) and spin-triplet ( 3 L J=L ) states with the same J = L can mix. Quantifying this mixing at different light quark masses can provide an insight into the flavour symmetry breaking. Using a two-state hypothesis and assuming energy-independent mixing we can determine the mixing angle defined in Eq. (6.1) of [64] from ratios of operator overlaps (interpreted non-relativistically) as described in that reference.
In Table 2, we show the mixing angles for the lightest pairs of P -wave (J P = 1 + ), D-wave (J P = 2 − ) and J P = 1 − hybrid states extracted using three different operators for the two different ensembles. The variation between mixing angles determined using the three different operators gives an estimate of the size of the systematic uncertainties as discussed in Ref. [64]. The 1 + mixing angle from the ρ − ρ 2 operator in the charmlight sector is closer to the heavy-quark limit value on the M π ∼ 240 MeV ensemble, but the analogous angle in the charm-strange sector does not differ significantly between the ensembles. For both charm-light and charm-strange mesons, the 2 − mixing angle is closer to the heavy-quark limit value for the lighter pion mass whereas the 1 − hybrid mixing angle shows no significant difference between the two ensembles.
In all cases on both ensembles, the determined mixing angles lie between the flavoursymmetry limit (0 • or 90 • ) and the heavy-quark limit values. This is expected since the charm quark, although much heavier than the light and strange quarks, is not heavy enough for the heavy-quark limit to apply strictly.
Heavy-quark limit  Table 2. Absolute value of the mixing angles for the lightest pairs of 1 + , 2 − and hybrid 1 − states in the charm-strange (c-s) and charm-light (c-l) sectors on the two ensembles. The mixing angles expected in the heavy-quark limit are also shown. In the M π ∼ 400 MeV case highlighted by the , we have subtracted the angle given in Ref. [64] from 90 • so that the mass ordering of the states is consistent between the two ensembles.

Summary and Outlook
We have presented spectra of excited hidden and open-charm mesons obtained from dynamical lattice QCD calculations with a pion mass of approximately 240 MeV. The use of distillation in combination with large variational bases of interpolating operators allows us to extract highly excited mesons, while the spin identification scheme has allowed a robust identification of the J P (C) of states as high as spin four, including states with exotic quantum numbers. The majority of mesons we extract can be interpreted in terms of the n 2S+1 L J pattern expected from quark potential models. However, excess states, with both exotic and non-exotic quantum numbers, that do not fit this pattern are also determined. By examining the operator overlaps we identify these as hybrid mesons, i.e. having excited gluonic degrees of freedom. The supermultiplets of hybrid mesons follow a pattern consistent with a quark-antiquark combination in S or P -wave coupled to a 1 +− gluonic excitation. The pattern and energy scale of hybrids are the same as that found in the light meson and baryon sectors [21,[57][58][59][60], studies of charmed baryons [61,62] and in our earlier work on charmonia and open-charm mesons [22,23].
Comparing the spectra to those from a similar lattice calculation with a pion mass of approximately 400 MeV, we find that the overall qualitative features are the same and, even in the case of charm mesons with a valence light quark, we find only small quantitative differences. The hybrid mesons appear to show a mild increase in mass as the pion mass is decreased but the pattern of states and supermultiplet structure is unchanged.
We also compared the spin-singlet -spin-triplet mixing angles for the lightest pairs of charm-strange and charm-light P -wave (J P = 1 + ), D-wave (J P = 2 − ) and hybrid (J P = 1 − ) states between the two lattice ensembles. Using a non-relativistic interpretation of operator overlaps, our results suggest that the mixing angles for the charm-light 1 + and the charm-light and charm-strange 2 − states become closer to those expected in the heavyquark limit as the pion mass is reduced. Conversely, we find no significant difference in the hybrid 1 − mixing angles between the two ensembles.
As discussed earlier, a limitation of these calculations is that we have not accounted for the unstable nature of states above threshold. This issue has already been addressed for a variety of mesons appearing as bound states and resonances in coupled-channel Dπ, Dη and D sK scattering [24]. The work presented here lays the foundation for extending these scattering calculations to pion masses closer to the physical value, as well as to other scattering channels involving hidden and open-charm mesons.
Illinois. Blue Waters is a joint effort of the University of Illinois at Urbana-Champaign and its National Center for Supercomputing Applications. This work is also part of the PRAC "Lattice QCD on Blue Waters". This research used resources of the National Energy Research Scientific Computing Center (NERSC), a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DEAC02-05CH11231. The authors acknowledge the Texas Advanced Computing Center (TACC) at The University of Texas at Austin for providing HPC resources that have contributed to the research results reported within this paper.
Gauge configurations were generated using resources awarded from the U.S. Department of Energy INCITE program at the Oak Ridge Leadership Computing Facility, the NERSC, the NSF Teragrid at the TACC and the Pittsburgh Supercomputer Center, as well as at Jefferson Lab.

A Tables of Results
In Tables 3, 4 and 5 we present numerical values for, respectively, the charmonium, D s and D meson masses obtained for M π ∼ 240 MeV. Masses are given in MeV with either the mass of the η c subtracted (charmonium) or half the mass of the η c subtracted (open-charm mesons) in order to minimise the systematic uncertainty in tuning the charm quark mass. In all cases the quoted error corresponds to the (one-sigma) statistical uncertainty. As discussed earlier, above the lowest multi-hadron threshold in each channel states can decay strongly into lighter hadrons and, aside from any other systematic uncertainties, we only expect the masses to be correct up to around the width of the state [21].  (7) 865 (7) 1316 (17) 1345(27) 1427(17) (10) 1138 (17) 1569 (26) 1783 (29) (16) 1623(26) 1665 (29) 1925(36) 3 + 1724 (21) 1743(21) 4 + 1804(22) Table 5. Summary of the D meson spectrum presented in Figure 5. Masses are shown with M ηc /2 subtracted. Quoted uncertainties are statistical only.