Tetraquark operators in lattice QCD and exotic flavour states in the charm sector

We present a general class of operators resembling compact tetraquarks which have a range of colour-flavour-spin structures, transform irreducibly under the symmetries of the lattice and respect other relevant symmetries. These constructions are demonstrated in lattice QCD calculations with light quarks corresponding to $m_\pi =$ 391 MeV. Using the distillation framework, correlation functions involving large bases of meson-meson and tetraquark operators are computed in the isospin-1 hidden-charm and doubly-charmed sectors, and finite-volume spectra are extracted with the variational method. We find the spectra are insensitive to the addition of tetraquark operators to the bases of meson-meson operators. For the first time, through using diverse bases of meson-meson operators, the multiple energy levels associated with meson-meson levels which would be degenerate in the non-interacting limit are extracted reliably. The number of energy levels in each spectrum is found to be equal to the number of expected non-interacting meson-meson levels in the energy region considered and the majority of energies lie close to the non-interacting levels. Therefore, there is no strong indication for any bound state or narrow resonance in the channels we study.


Introduction
There are an abundance of experimentally-observed mesons containing one or more heavy quarks [1] and these provide a window on a rich variety of strong-interaction physics. In particular, many of them, the so-called 'X, Y, Z's', are not compatible with quark model expectations. Clear examples of this incompatibility are the charged charmonium-like Z + c (3900) and Z + c (4430) which cannot be solely cc and must contain at least one additional quark-antiquark pair. One possible explanation for such exotic states is that they are tetraquarks, compact bound states of four quarks. Others, for example Ref. [2], suggest that compact tetraquarks are not required to explain the observed spectrum. Recent reviews of some of the X, Y, Z's, with interpretations such as compact tetraquarks, molecular mesons, hybrid mesons and threshold cusps, can be found in Ref. [3][4][5][6][7][8]. As well as hidden-charm ccqq configurations, doubly-charmed ccqq tetraquarks have been hypothesised [9][10][11][12], but there are currently no experimental candidates for these.
Quantum Chromodynamics (QCD) is the fundamental theory of the strong interaction and, in principle, should predict whether four-quark states exist and whether these are consistent with the expectations of tetraquark models or other interpretations. The only ab-initio framework for performing systematically-improvable calculations at the hadronic scale is lattice QCD: spacetime is discretised on a finite four-dimensional Euclidean lattice and Monte Carlo techniques are used to compute correlation functions from which observables can be extracted. The discrete spectrum of finite-volume energy eigenstates is obtained from calculations of two-point correlation functions involving interpolating operators which have the required quantum numbers. Hadron-hadron scattering amplitudes, and hence the properties of resonances and other scattering phenomena, can be calculated via the Lüscher formalism [13,14] which relates finite-volume spectra to infinite-volume scattering amplitudes. There is currently no extension of this formalism to three or more hadron scattering channels that is practical to use in calculations but this is an active area where progress is being made. A more in-depth review of the Lüscher formalism and a discussion on its applications and extensions can be found in Ref. [15].
The Hadron Spectrum Collaboration has developed a range of interpolating operators resembling quark-antiquark [16,17] and meson-meson [18,19] structures which transform irreducibly under the symmetries of the lattice and efficiently interpolate the states of interest. These operators have proven very successful in recent computations of finitevolume spectra which are then used to determine scattering amplitudes [18][19][20][21][22][23][24][25][26]. As has been emphasised in studies such as those in Refs. [19,22], not including a sufficiently diverse set of relevant operators in the calculations could lead to an unreliable determination of finite-volume spectra and, in turn, incorrect scattering amplitudes. Hence, it is desirable to consider operators with other potentially-relevant colour-flavour-spatial-spin structures, resembling compact tetraquarks, and investigate whether their inclusion has any impact on the extracted spectra. The main goal of this work is to develop a very general class of operators with compact tetraquark structures, which transform irreducibly under the symmetries of the lattice and which respect other relevant symmetries. We will test these constructions in lattice QCD computations of spectra in hidden-charm and doubly-charmed channels. These include isospin-1 J P G = 0 +− , 1 ++ , 1 +− hidden-charm spectra 1 which are relevant for exotic charged charmonium-like states and where the lightest tetraquark multiplet is expected to appear [27], and isospin-0 J P = 0 + , 1 + , 2 + and isospin- 1 2 strange J P = 0 + , 1 + exotic doubly-charmed spectra.
There have been a number of recent lattice QCD studies of tetraquarks containing one or more heavy quarks. Computations have not found any clear indication for the presence of hidden or open-charm tetraquarks [28][29][30][31]. Other recent lattice QCD calculations relevant for the channels we study can be found in Refs. [32][33][34]. In the bottom sector, there is some evidence supporting the existence of a doubly-bottom (I)J P = (0)1 + tetraquark where finite-volume spectrum calculations find an energy level below the relevant mesonmeson thresholds [35]. In addition, a number of computations of the potential between two static quarks in the presence of two light quarks [36][37][38][39][40][41][42] have found evidence for a bound state [36,38,[40][41][42]. We discuss these studies further in the context of our results in Section 6.
The structure of the rest of this paper is as follows. We begin in Section 2 by describing the construction of a general class of tetraquark operators which transform irreducibly under the symmetries of the lattice. In Section 3, the methodology for calculating the finitevolume spectrum with large bases of operators in the distillation framework is presented. The resulting spectra in the hidden-charm and doubly-charmed sectors are presented in Section 4. Some systematic effects and the stability of the extracted spectra are investigated in Section 5. We discuss the results in light of phenomenological and other lattice QCD studies in Section 6 before giving a summary in Section 7. Appendices present some additional properties of diquark and tetraquark operators, give quark model interpretations of the diquark structures and list the operators used to calculate the finite-volume spectra.

Tetraquark operator construction
To construct interpolating operators which resemble a compact tetraquark, we combine a diquark operator with an anti-diquark operator. The diquark operator is built from two quark fields coupled together to obtain appropriate colour, flavour and spin quantum numbers and, analogously, the anti-diquark operator is built from two antiquarks. The diquark and anti-diquark are then combined to form a colour singlet with the desired flavour and spin. These constructions provide, with no loss in generality, a convenient way to build a diverse class of tetraquark operators which have the required quantum numbers and respect appropriate symmetries. In this section we present an overview of these operators: we begin by describing the flavour and colour structures before presenting expressions for the diquark, anti-diquark and tetraquark operators. Further details and a discussion of some additional properties can be found in Appendix A. Further modeldependent understanding on how the different diquark configurations interpolate different states can be found in Appendices B and C.
The diquark operator is constructed by coupling two quark fields together to definite colour, flavour and continuum spin. In colour space, the quarks belong in the fundamental representation of SU(3) C and so the diquark is in either the antisymmetric3 representation or the symmetric 6 representation. In flavour space, we use SU(3) F constructions to form a convenient basis of operators, but this does not imply any assumption of SU(3) F symmetry in the theory -as long as a sufficient basis is used, an arbitrary flavour combination can be constructed from a linear combination of these operators. The up, down and strange (u, d and s) quarks belong in the fundamental representation of SU(3) F and the charm (c) quark is placed in a singlet. The quarks are coupled together to obtain the desired flavour irrep out of 1, 3,3 and 6. For example, coupling two u, d, s quarks as 3 ⊗ 3 →3, this irrep gives a component with flavour quantum numbers (isospin, strangeness) = (0, 0) with flavour structure 1 √ 2 (ud − du), and a ( 1 2 , −1) multiplet with flavour structure 1 √ 2 (us − su) and 1 √ 2 (ds − sd). Alternatively, coupling a c quark with a u, d, s quark as 1 ⊗ 3 → 3 gives a (0, −1) component with flavour structure cs, and a ( 1 2 , 0) multiplet with flavour structure cu and cd. If both quarks are in the same flavour representation, Fermi symmetry requires that the overall operator is antisymmetric under the interchange of the quarks and this constrains the allowed diquark configurations.
The diquark operator in colour irrep R (row r), flavour irrep F (row f ) and continuum (2.1) where spinor indices have been suppressed, q( x, t) is a quark field smeared with the distillation operator as discussed in Section 3.2, D a , d a ; D b , d b |D, d are the SU(2) or SU(3) Clebsch-Gordan coefficients that couple the irreps D a ⊗D a → D with d a , d b and d the irrep rows, C is the charge conjugation matrix such that γ 0 = Cγ T 0 C, and Γ is a Dirac gamma matrix which determines J, m and other properties of the diquark operator as shown in Table 3 in Appendix A. Choosing an appropriate Γ gives access to spins up to J = 1 -in order to access higher spins or excitations, this operator can be generalised by including gauge-covariant derivatives in a similar way to the fermion-bilinear operator constructions discussed in Ref. [16].
In the anti-diquark operator, the antiquarks belong in the anti-fundamental representation of SU(3) C and therefore couple to colour irrep 3 or6. The up, down and strange antiquarks belong in the3 irrep of SU(3) F and the charm antiquark is in the singlet. Possible flavour irreps for the anti-diquark are therefore 1, 3,3 and6. If both antiquarks are in the same flavour irrep, Fermi symmetry again constrains the allowed configurations. The anti-diquark operator is defined in an analogous way to the diquark operator as, is a (smeared) antiquark field. Here the charge conjugation matrix comes after the Dirac gamma matrix so that under the charge conjugation operator Cδ RF , and this ensures a convenient definition of tetraquark operators with definite G-parity as discussed later.
Tetraquark operators are formed by coupling a diquark operator and an anti-diquark operator to a colour singlet with definite flavour and spin. The only possible diquark and anti-diquark colour combinations which give a colour singlet are3 ⊗ 3 and 6 ⊗6 and this restricts the possible diquark-anti-diquark configurations. The flavour quantum numbers of the tetraquark operator are obtained by coupling the appropriate flavour irreps of the diquark and anti-diquark and then choosing the desired row. By projecting onto zero momentum, tetraquark operators have definite parity and, in channels where G-parity is a good quantum number, operators with definite G-parity can be constructed. The tetraquark operator, projected onto momentum p, is, For the remainder of this study, we only consider p = 0 and so this operator has definite parity P which is determined by the gamma matrices Γ 1 and Γ 2 as described in Appendix A. In channels where G-parity is a good quantum number, a tetraquark operator with definite G-parity is given by, The G-parity of this operator is G =Gξ J ξ 1 ξ 3 where ξ J , ξ 1 , ξ 3 are phases arising from the exchange symmetry of the Clebsch-Gordan coefficients in Equation (2.3) as described explicitly in Appendix A.
This tetraquark operator has definite spin J in the continuum but the lattice discretisation breaks rotational symmetry and so J is no longer a good quantum number. With a cubic lattice discretisation and volume, the symmetry group is reduced to the octahedral group O h for states at rest [43] and broken further to the little group for states with non-zero momentum [44]. The distribution, or subduction, of J ≤ 4 into the irreps of O h is tabulated in Table 1. We construct lattice tetraquark operators which transform in lattice irrep Λ (row µ) by subducing the continuum operators as described in Ref. [16], where S are subduction coefficients. The generalisation to p = 0 involves the construction of helicity operators and then the subduction to irreps of the little group of p as discussed in Ref. [17]. We will use these lattice tetraquark operators to calculate correlation functions in lattice QCD from which the finite-volume spectrum can be extracted.

Calculation of the spectrum
To determine the spectrum in each quantum-number channel, we calculate a matrix of two-point correlation functions using a basis of interpolating operators with appropriate quantum numbers, between a creation operator O † j (0) at the source with Euclidean time 0 and an annihilation operator O i (t) at the sink with Euclidean time t. Inserting a complete set of energy eigenstates into this expression gives, where |n is an energy eigenstate with energy E n and the operator-state matrix elements, Z n i ≡ n|O † i (0)|0 , are also referred to as overlaps. Note that in a finite volume, the set of energy eigenstates is discrete. The spectrum can be extracted by utilising the variational method [45][46][47]: a generalised eigenvalue problem C ij (t)v n j = λ n (t, t 0 )C ij (t 0 )v n j is solved for some appropriate choice of t 0 , the eigenvalues λ n , known as principal correlators, are related to E n and the eigenvectors v n i are related to the overlaps. In our implementation of the variational method described in Refs. [16,48], we fit the principal correlators to the function, where the fit parameters are E n , E n , and A n . The second exponential is used to account for possible contamination due to excited states. The eigenvectors can be used to construct optimised operators, [18], the optimal linear combination of the operators that interpolates state n. As will be discussed later, these optimised operators are useful for the construction of operators resembling pairs of mesons.
In principle, any operator can interpolate every state with the same quantum numbers according to Equation (3.2) and one can use a basis containing the tetraquark operators described above to fully calculate the finite-volume spectrum. However, in practice, previous studies such as Refs. [19,22] have highlighted that a sufficiently diverse set of interpolating operators must be used if finite-volume spectrum is to be extracted reliably. Even in the absence of interactions, the spectra we study will contain meson-meson-like states or admixtures of such states (and other multi-hadron combinations at higher energies). It has been show in previous work [18][19][20][21][22][23][24][25][26] that such states can be efficiently interpolated by including meson-meson-like operators. Therefore, to efficiently and reliably extract the finite-volume spectra, our operator bases will contain operators of meson-meson and tetraquark structure.

Meson-meson operators
In this section we briefly review how operators with a meson-meson-like structure can be constructed from the product of two single-meson-like operators -further details are given in Refs. [18,19].
Following Refs. [16,17], fermion-bilinear operators of continuum spin J and momentum p are constructed as, where we have suppressed colour, flavour and spinor indices for clarity. The quark and antiquark fields are distillation-smeared as discussed below, the quark and antiquark flavour representations are chosen and coupled to give the desired flavour quantum numbers, and [Γ ← → D . . . ← → D ] J,m consists of a Dirac gamma matrix Γ and gauge-covariant derivatives ← → D coupled together to give spin J. In the case when p = 0, m refers to the J z component, whilst when p = 0 we construct helicity operators using Wigner-D matrices as described in Ref. [17] and m then refers to the helicity. Since we are working in a finite cubic spatial volume of extent L with periodic boundary conditions, the momentum is quantised to p = 2π L (n x , n y , n z ) where (n x , n y , n z ) is a triplet of integers -we use [n x n y n z ] as a shorthand notation to denote p. As for tetraquark operators, lattice operators which transform in lattice irrep Λ (row µ) are constructed by subducing the continuum operators to obtain O . We refer to these as single-meson operators. Meson-meson operators [18,19] are built from products of two single-meson operators, (3.5) where we have restricted to overall zero momentum, Ω M i Λ i ,µ i is an optimised operator for interpolating meson M i transforming in the lattice irrep Λ i (row µ i ) and the sum runs over the lattice irrep rows and all momentum directionsq related by an allowed lattice rotation to couple Λ 1 ( q) ⊗ Λ 2 (− q) → Λ P ( p = 0) using generalised Clebsch-Gordan coefficients C. The construction of these meson-meson operators follows the methodology given in Refs [18,19] and a more detailed discussion of operators containing mesons with non-zero spin will be presented in a forthcoming publication. Analogous constructions can be used for meson-meson operators with overall non-zero momentum, but in this study we only calculate spectra at overall zero momentum.
A guide to which meson-meson operators should be included in the basis is given by the non-interacting meson-meson energy levels in the energy regions we consider. These are calculated from the relativistic dispersion relation E = m 2 1 + p 2 1 + m 2 2 + p 2 2 for stable single-mesons. In cases when a single-meson has non-zero spin, there can be multiple ways to couple the orbital and spin angular momenta together to a given meson-meson J P which subduce into the same irrep, leading to degenerate levels in the non-interacting limit. For example, a pseudoscalar and vector can be coupled to J P = 1 + in either s-wave or d-wave. To see how this manifests on the lattice, consider the pseudoscalar with p = [100] and the vector with p = [−100] coupled to the lattice irrep Λ P = T + 1 . The mesons transform in the irreps of the little group Dic 4 : the pseudoscalar subduces into the A 2 irrep while the helicity-0 and helicity-1 components of the vector subduce into respectively the A 1 and E 2 irreps. It is possible to obtain Λ P = T + 1 from both A 2 ⊗ A 1 and A 2 ⊗ E 2 [49] and two different linear combinations of these would correspond to s and d wave in the continuum and infinite volume limit. In general, we must include a sufficient number of relevant meson-meson operators that are capable of extracting and disentangling these multiple energy levels. The comparison of spectra calculated with different operator bases in Section 5 demonstrates the importance of including a sufficient basis of meson-meson operators.

Calculation of correlation functions
We choose a basis of meson-meson operators as described above and tetraquark operators as described later in Section 4 and compute the two-point correlation functions using the distillation framework [50]. The combination of distillation and the techniques described here have been demonstrated in, for example, Refs. [16, 18, 19, 21-25, 51, 52]. In brief, the distillation operator on timeslice t which acts in 3-space ( x and y) and colour space (r and s) is defined as, ( xr, ys; t) = Nvecs n=1 ξ n ( xr; t)ξ † n ( ys; t), where in the implementation used here ξ n are the lowest N vecs eigenvectors of the gauge-covariant Laplacian. The quark and antiquark fields in interpolating operators are smeared by the distillation operator, q → q, which removes high-frequency modes and increases the overlap onto lower-lying states. Distillation also allows for a factorisation of the two-point correlation functions where M is the Dirac matrix, and elementals which describe the operators with various structures projected onto definite momentum, and this enables the efficient computation of the correlation function matrix for a large basis of operators. Meson elementals Φ αβ n 1 n 2 ( p, t) are presented in Ref. [19] and tetraquark elementals are given by, (3.6) where n i index distillation vectors, Greek letters label the Dirac spinor indices and C pqrs are combinations of SU(3) C Clebsch-Gordan coefficients that couple the colour representations 3 ⊗3 ⊗ 3 ⊗ 3 → 1 as in the tetraquark operators in Section 2. Meson elementals are matrices with (4N vecs ) 2 independent components and tetraquark elementals are rank-4 with (4N vecs ) 4 independent components. This means that the cost of calculations involving tetraquark operators (multiplying and tracing perambulators and elementals) increases rapidly when the number of vectors is increased. Therefore, if the calculation is to be feasible, the number of vectors must not be too large.
To keep the cost of contractions reasonable by using a relatively small number of vectors for tetraquark operators, whilst maintaining a larger number of vectors for other operators, we introduce a second distillation operator,˜ ( xr, ys; t) = Ñ vecs n=1 ξñ( xr; t)ξ † n ( ys; t), composed of the lowestÑ vecs vectors whereÑ vecs < N vecs . Quark/antiquark fields in tetraquark operators are smeared with˜ whereas those in other operators are smeared with . As an example, consider a meson-meson operator given by O ∼ (c Γ u)(d Γ c) and a tetraquark operator, T ∼ (˜ c) T CΓ(˜ u) (c˜ )ΓC(d˜ ) T , where we are suppressing various indices and factors which are not relevant for this discussion. One of the connected contributions to the correlation function between these two operators is, schematically,  Figure 1. A schematic representation of the types of Wick contractions required to compute the two-point correlation function matrices in this study. We use Φ (grey) to depict the singlemeson elementals, Ψ (red) to depict the tetraquark elementals and the lines joining them to depict perambulators.
square matrices and Ψ(0) is of rank-4 with (4Ñ vecs ) 4 components. The viability of having a lower number of distillation vectors for tetraquark operators and some tests of varying the number of vectors are discussed in Section 5.3. Although not utilised in this study, another possible use of employing more than one distillation operator is to increase the number of operators in the variational basis through including operators with different smearings. To further reduce the computation time, we only calculate one half of the off-diagonal elements in the matrix of two-point correlation functions between a tetraquark operator and meson-meson operator, and then obtain the other half using the Hermiticity of the correlation matrix. In addition, we neglect contributions where a charm quark and antiquark annihilate: these are expected to be small due to OZI suppression and this has been found to be the case empirically in lattice calculations [53]. The elements of the two-point correlation function matrix that we compute are shown in Figure 1 where we show a schematic representation of the types of Wick contractions required.

Results
As a first application of these tetraquark operator constructions, we perform calculations on an anisotropic lattice of volume (L/a s ) 3 × (T /a t ) = 16 3 × 128 where L is the spatial extent of the lattice, T is the temporal extent, a s ≈ 0.12 fm is the spatial lattice spacing and a t is the temporal lattice spacing such that the anisotropy ξ = as at ≈ 3.5. We use 478 configurations generated from a tree level Symanzik improved gauge action and a Clover fermion action with N f = 2 + 1 flavours of dynamical quarks. The mass parameter of the two degenerate light quarks is such that m π = 391 MeV while the strange quark is tuned so that its mass approximates the physical value [54,55]. The quenched Clover charm quark mass parameter is tuned to reproduce the physical η c meson mass [52]. When quoting results in physical units, we set the scale using the mass of the Ω baryon from the measured value on this lattice, a t m latt. Ω = 0.2951 (22) [56], and the experimental mass volumes, the anisotropy was measured to be ξ π = 3.444(6) from the dispersion relation of the pion [18] and ξ D = 3.454(6) from the D [57]. Using only the 16 3 volume, we find ξ ηc = 3.484(2) from the η c . For the purposes of this study, where the anisotropy is only used to compute the location of non-interacting meson-meson energy levels, we will use the value of ξ π .
To indicate the location of non-interacting meson-meson energy levels on plots, for stable mesons we use the relativistic dispersion relation giving, E = m 2 1 + p 2 1 + m 2 2 + p 2 2 , as discussed in Section 3.1. The masses of relevant stable mesons on this lattice ensemble are given in Table 2 and the variationally-optimised operators for these mesons, Ω M Λ,µ , are constructed from linear combinations of single-meson operators as discussed above. For the ρ meson, which is unstable on this lattice, we compute 'non-interacting M -ρ energy levels', where M is a stable meson, using the relativistic dispersion relation for M and the finite-volume ρ energy levels obtained on this ensemble as given in Table 2, The optimised ρ operators, Ω ρ Λ,µ , are linear combinations of meson-meson and single-meson operators [19]. It should be emphasised that in this study the only uses of these non-interacting energy levels are to show their location on plots and as an indication for which meson-meson operators should be included in the operator basis.
The meson-meson and tetraquark operators used to calculate the spectrum in each channel are listed in Appendix D. For this lattice volume, we useÑ vecs = 24 for tetraquark operators and N vecs = 64 for other operators unless stated otherwise. The choice of mesonmeson operators has already been discussed in Section 3.1. For the tetraquark operators, ideally all relevant operators of the form described in Section 2 would be included in the operator basis, but the computational cost can then become too high for the calculations to be practical for the purposes of this study. In the doubly-charmed isospin-0 Λ P = A + 1 , E + , T + 2 channels we are able to include all the tetraquark operator constructions allowed by Fermi symmetry. For the remaining channels, we use a subset of tetraquark operators that are expected to overlap onto lower-lying states. Appendices B and C describe how a nonrelativistic model can provide a guide to which diquark/anti-diquark configurations are expected to overlap most efficiently onto lower-lying tetraquarks. In short, tetraquark operators containing a diquark (anti-diquark) with a γ 5 or γ 0 γ 5 gamma matrix structure and in colour irrep3 (3) Table 2. Ground state masses of stable mesons (left) and the energy of the lowest-lying finitevolume energy level for lattice irrep Λ and momentum relevant for the ρ meson denoted by ρ p Λ (right) as measured on our ensemble [18,19,52,57]. Only the statistical uncertainty is quoted.

1089(5)
found that the γ 5 and γ 0 γ 5 structures do not overlap onto the energy eigenstates in sufficiently distinct ways (the correlation matrix contains approximately linearly-dependent rows/columns). Therefore, instead of having redundant operators, we included a selection of other tetraquark operators to give more diverse bases.
We now present computed spectra for a range of channels, beginning with a detailed discussion of the Λ P G = T ++ 1 irrep in the isospin-1 hidden-charm sector before presenting other isospin-1 hidden-charm results and then moving to the doubly-charmed sector.  isospin-1 hidden-charm channel. The first three operators are tetraquark operators and the remaining are meson-meson operators ordered as in Table 6.

Isospin-1 hidden-charm sector
As an illustration of the results, we first discuss some features of the spectrum computed in the Λ P G = T ++ 1 irrep 3 of the isospin-1 hidden-charm sector (with flavour content ccll where the light quark and antiquark are coupled to isospin-1). The basis of operators used is given in Table 6 of Appendix D. Note that, because we do not include contributions arising from a charm quark and antiquark annihilating, our operator bases do not contain any single-meson operators.
The diagonal elements of the matrix of correlators for the three tetraquark operators are shown in Figure 2 -signals are seen to be precise and significantly non-zero. In Figure 3 we present the two-point correlator matrix on timeslice 3. This shows that some of the off-diagonal elements between tetraquark and meson-meson operators are non-zero.
After applying the variational method, principal correlators for the lowest six energy levels are shown in Figure 4. We fit these to Equation (3.3) and in each case find a reasonable description having χ 2 /N d.o.f ∼ 1 -the resulting spectrum is given in the figure. It can be seen that the number of energy levels in the computed spectrum is equal to the number of non-interacting meson-meson levels expected in the energy region considered and they all lie close to the non-interacting levels. As discussed in Section 3. Energy (MeV) Figure 4. The central plot shows the spectrum in the hidden-charm isospin-1 Λ P G = T ++ 1 channel calculated using the basis of meson-meson and tetraquark operators given in Table 6 of Appendix D. Boxes give the computed energies with their vertical extent representing the one-sigma statistical uncertainty on either side of the mean and, solely as a visual aid, they are coloured according to their dominant meson-meson operator overlap. Horizontal lines denote the non-interacting mesonmeson energy levels with an adjacent number indicating the degeneracy if it is larger than one. The corresponding principal correlators are shown on the left ordered by increasing energy from bottom to top: the data (points) and fits (curves) for t 0 = 9 are plotted as λ n (t, t 0 )e En(t−t0) showing the central values and one sigma statistical uncertainties; in each case the fit is reasonable with The histograms on the right show the operator-state overlaps, Z n i = n|O † i |0 , for each energy level. The operators are given in the legend and the overlaps are normalised so that the largest value for one given operator across all energy levels is equal to one. the figure, some of the non-interacting meson-meson energy levels are degenerate. Because our basis of operators has sufficiently different structures, we are able to cleanly extract nearly-degenerate energy levels.
Normalised operator-state overlaps are also shown in Figure 4 and we see that every energy level has a dominant overlap onto one meson-meson operator. Additionally, the third and fourth levels have dominant overlaps onto two linearly independent J/ψπ operators -this is not surprising since around this energy there are two degenerate non-interacting levels. We cannot draw strong quantitative conclusions about the tetraquark operator overlaps because the absolute normalisations are somewhat arbitrary and renormalisation factors would be needed to relate the overlaps to physical quantities, but we do see that most states have some overlap onto one or more tetraquark operators.
For comparison, Figure 5 (left panel) shows the Λ P G = T ++ 1 spectrum calculated with the full basis of meson-meson and tetraquark operators, with only meson-meson operators and with only tetraquark operators. No significant deviations are observed between the spectrum computed using the full basis and that computed using the basis of only mesonmeson operators. If only tetraquark operators are used, some poorly determined energy levels are found but the spectrum is not reliably extracted and this suggests that these tetraquark operators alone do not constitute a sufficient basis of operators. Moving to other channels in the isospin-1 hidden-charm sector, extracted spectra for the Λ P G = A +− 1 , T +− 1 irreps 4 are shown in Figure 5. In general, a similar pattern of features is seen as was found for Λ P G = T ++ 1 : there are no significant deviations between the spectra calculated using the full basis and using only meson-meson operators, the spectrum is not reliably determined if only tetraquark operators are used, and with a full basis of operators the number of energy levels is equal to the number of non-interacting meson-meson levels expected and they lie close to the non-interacting levels. Furthermore, the operator-state overlaps follow the same qualitative pattern as shown in Figure 4.
From previous studies, when a narrow resonance is present in elastic scattering, an 'extra' finite-volume energy level is observed in that energy region but no evidence for such an extra level is seen in our spectra. The results suggest that there are only weak hadron-hadron interactions and no strong indications of a bound state or narrow resonance in these channels. However, the situation is not as straightforward when one considers coupled-channel scattering or broad resonances [23,24]. To draw rigorous conclusions and determine whether there are bound states or resonances present, a Lüscher analysis, where the finite-volume spectra are related to the scattering amplitudes, is necessary. To reliably constrain the scattering amplitudes, this would require calculations at non-zero momentum and/or different volumes which is beyond the scope of this first study. An important conclusion is that the addition of a class of operators resembling compact tetraquarks has little consequence on the finite-volume spectrum and, in turn, the scattering amplitudes.

Doubly-charmed sector
Turning to the doubly-charmed sector, Figure 6 shows spectra for flavour content ccll in the Λ P = T + 1 , E + , T + 2 isospin-0 channels 5 and Figure 7 shows spectra for flavour ccls with isospin-1 2 in the irreps Λ P = A + 1 , T + 1 . It can be again seen that there are no significant deviations between the spectra including and excluding tetraquark operators, and the spectra can not be reliably extracted using only tetraquark operators. Using the full basis of operators (see Table 7), the number of energy levels in each spectrum is equal to the number of expected non-interacting meson-meson energy levels in the relevant energy region. Because the basis of operators used has sufficiently diverse structures, we are able to extract many nearly-degenerate energy levels. In addition, we find that every energy level has a dominant meson-meson operator overlap. As in the results of the hidden-charm sector, we emphasise that the addition of a class of operators resembling compact tetraquarks does not significantly alter the finite-volume spectrum extracted.
The lowest-lying DD and D * D * levels in s-wave are forbidden in the J P = 0 + , 2 + isospin-0 doubly-charmed channels: the flavour wavefunction is antisymmetric in isospin-0 whilst the spin and spatial wavefunctions are symmetric, giving an overall antisymmetric wavefunction which is forbidden by Bose symmetry. These channels are particularly appealing to look for a tetraquark because if a low-lying state exists, it would lie far below the allowed non-interacting meson-meson energy levels and would be easily identified. Additionally, a low-lying J P = 2 + stable tetraquark would subduce into both of the irreps 4 The lowest spin in each of these irreps is respectively J P G = 0 +− , 1 +− . 5 The lowest spin J P = 1 + appears in T + 1 and the lowest spin J P = 2 + appears in E + , T Λ P = E + , T + 2 and so appear with little ambiguity -no such energy levels are seen in Figure 6. A J P = 0 + tetraquark would appear in the Λ P = A + 1 irrep and, although a plot is not shown, we calculated the spectrum in this channel with the operators given in Table 7. The first allowed non-interacting meson-meson energy level is D D(2S), where D(2S) is the first radial excitation of the D meson. 6 We do not find any energy levels below the DD(2S) threshold at ∼ 4500 MeV.
We now draw particular attention to the spectrum in the Λ P = T + 1 isospin-0 channel where the non-interacting DD * and D * D * levels can have degeneracy two. We have reliably extracted two energy levels (the third and fourth) that have dominant overlap onto the two relevant DD * operators and two energy levels (fifth and sixth) that have dominant overlap onto the two relevant D * D * operators. It can be seen that each pair of energy levels is nondegenerate which suggests there is some interaction. In order to quantify this, a further analysis requiring computations on different volumes and overall non-zero momentum is needed to relate the finite-volume spectrum to the scattering amplitudes via the Lüscher formalism. It is also important to stress that a reliable determination of the coupled s and 6 We use a single-meson operator with structure [γ 5 ← → D ← → D ], where the two derivatives are coupled to J = 0, for the D(2S) rather than an optimised operator. d-wave scattering amplitudes in this channel depends on our ability to robustly extract these multiple energy levels.

Systematics and stability of the extracted spectra
Before discussing the results further, we consider some systematic effects which may have an impact on them and present some tests of varying the operator basis and the number of distillation vectors.

Systematic uncertainties
As a first application of these tetraquark operator constructions, we have performed calculations on a relatively small lattice volume, with spatial extent L ∼ 2 fm, and this may be too small to distinguish the spatial structures of the extended meson-meson and compact tetraquark. The tetraquark operator can be Fierz rearranged as a linear combination of meson-meson operators multiplied by a factor of 1/L 3 [31], suppressing the overlap of a possible tetraquark state onto the meson-meson operators. Further calculations, beyond the scope of this study, would be required to give some indication on how the results vary with the volume. As this is a first demonstration, we have performed calculations with unphysicallyheavy light quarks, corresponding to m π = 391 MeV, and the presence/absence of tetraquarks may depend on the mass of the light quarks. Ultimately, calculations with light-quark masses approaching their physical values are required for comparison with experiment. On the other hand, studying how the spectra change as the quark masses are varied would give insight into the relevant QCD interactions and could be compared with expectations in different models.
Other possible systematic uncertainties include discretisation effects and the tuning of the charm quark mass. These issues were addressed in Ref. [52] and we do not repeat the discussion here.

Varying the operator basis
Some tests of how varying the operator bases affects the results have already been presented in Section 4. In summary, it was found that there were no significant changes in the lowlying spectra when only meson-meson operators were used compared to using the full basis of tetraquark and meson-meson operators. However, reliable spectra could not be extracted if only tetraquark operators were used.
As an illustration of what could happen if a sufficiently diverse set of meson-meson operators is not used, we show in Figure 8 spectra in the doubly-charmed isospin-0 Λ P = T + 1 channel computed using different operator bases. Note that degenerate meson-meson energy levels would be present here in the non-interacting limit as discussed in Section 3.1. Column A shows the spectrum computed using the full basis of meson-meson and tetraquark operators (see Table 7) and we see that the number of energy levels is equal to the number of expected non-interacting meson-meson energy levels in the energy region considered. We make the same conclusion for column B which shows the spectrum calculated using only meson-meson operators. Column C shows the spectrum using only meson-meson operators without the D

[100]
operator which is relevant for the DD * level at ≈ 4100 MeV -it is seen that now one fewer energy level is extracted and the second DD * level moves slightly higher in energy which is as expected when not enough operators are used [19]. The right column D shows the spectrum calculated with the operators as in C supplemented with the tetraquark operators -an additional level is found compared to C high up in the spectrum. This demonstrates the necessity of accounting for all the relevant meson-meson energy levels in the energy region being considered and using a sufficient basis of operators of different structures. Otherwise there is the danger that this level could be mistakenly taken as a signal for the presence of a tetraquark.

Varying the number of distillation eigenvectors
If the number of distillation eigenvectors used for the tetraquark operators is too low, the operator may not efficiently interpolate states of interest as it may be too smeared and so no longer resemble a compact tetraquark. However, as discussed in Section 3.2, the computational cost involving tetraquark operators scales much more strongly than mesonmeson operators with the number of distillation eigenvectors and therefore the number used can not be too large if the calculations are to be feasible. In this section, we test how sensitive the results are to varying number of distillation eigenvectors.
The spectrum in the doubly-charmed isospin-0 Λ P = T + 1 channel is shown in Figure 9 using different numbers of distillation vectors for tetraquark operators,Ñ vecs = 16, 24, 32, with both the full basis of meson-meson and tetraquark operators (see Table 7) and only tetraquark operators. It can be seen that the results are not sensitive to the number of distillation eigenvectors used.
We also computed all the spectra using N vecs =Ñ vecs = 24, i.e. the same number of distillation vectors for both meson-meson and tetraquark operators. The results were found to be consistent with the spectra presented in Section 4 (which usedÑ vecs = 24 for tetraquark operators and N vecs = 64 for meson-meson operators). This shows that the results are also not sensitive to the number of distillation vectors used for the meson-meson operators.
In summary, these tests suggest that the results are not very sensitive to the number of distillation vectors being used. In addition, a recent study in Ref. [58] demonstrated that a small number of distillation vectors is sufficient to extract finite-volume spectra as long as one does not consider higher momenta, higher spin or highly excited states. Because we are considering overall zero momentum and relatively low-lying states, this gives further support to our conclusion.  Figure 6 but showing the Λ P = T + 1 isospin-0 ccll spectrum calculated using tetraquark operators with different numbers of distillation vectors,Ñ vecs = 16, 24, 32. The left three columns are using the full basis of meson-meson and tetraquark operators, and the right three columns are using only tetraquark operators.

Discussion and comparison with previous studies
In this section we discuss the results in the context of expectations from phenomenological models and compare with previous lattice calculations. In a simple one-gluon-exchange model of a diquark as described in Appendix B, the two quarks interact via a colourcolour spin-spin interaction and the most attractive diquark and anti-diquark configurations have (colour irrep, spin) = (3, 0) and (3, 0) respectively. Therefore, the most favourable tetraquark has J P = 0 + which subduces into Λ P = A + 1 . However, a large quark mass suppresses the spin-spin interactions and so some models expect spin-1 diquark configurations involving heavy quarks to occur and form a tetraquark multiplet [27]. This multiplet contains tetraquarks with J P = 1 + and 2 + which subduce into the T + 1 irrep and the E + , T + 2 irreps respectively. Besides the quark mass, this one-gluon-exchange interaction does not depend on the flavours of the quarks. However, when the flavour irreps of the two quarks (antiquarks) are the same, Fermi symmetry requires the overall diquark (anti-diquark) configurations to be antisymmetric and this restricts the allowed structures.
In the hidden-charm isospin-1 sector, there are no identical quarks/antiquarks and so no constraints from symmetry on the allowed configurations. Models [27,59] suggest that the lightest tetraquark multiplet has J P G = 0 +− , 1 ++ , 1 +− , 2 +− and we have performed a thorough investigation in all these channels except for J P G = 2 +− . In the Λ P G = A +− 1 channel, expected to be the most attractive, and the Λ P G = T ++ 1 and T +− 1 channels, there are no hints of a narrow state or any significant interactions in the computed spectra. No experimental candidate has been observed with J P G = 0 +− , nor is there currently any charged charmonium-like candidate with undetermined J P G that is light enough to be identified as the lowest-lying J P G = 0 +− tetraquark. The observed Z + c (3900) has J P G = 1 ++ [1] and has been suggested to be a candidate for a tetraquark. That we see no sign of it is consistent with previous lattice QCD calculations presented in Ref. [29] which also calculated the finite volume spectrum using meson-meson and tetraquark operators. Our results are also consistent with other lattice QCD calculations [32][33][34] which do not find evidence of a bound state or narrow resonance in this channel. There is currently no wellestablished experimental candidate with J P G = 1 +− , but if the X(3872) is a tetraquark its isospin-1 partner would appear in this channel. Again, that we see no clear signal for a state here is consistent with previous lattice QCD calculations [31]. That study also found that the spectrum in this channel was insensitive to the addition of tetraquark-like operators to the operator basis.
In the doubly-charmed sector, possible diquark configurations are further constrained by Fermi statistics. The (3, 0) and (6, 1) cc diquarks are forbidden as they are symmetric under the interchange of quarks and only the (3, 1) and (6, 0) diquarks are allowed. In the one-gluon exchange model given in Appendix B, the colour-colour spin-spin interaction is repulsive for these allowed configurations and is least repulsive for (3, 1). The attractive (3, 0)qq anti-diquark configurations are required to be antisymmetric in flavour in F = 3 and the most attractive configuration has isospin, I = 0. Therefore, the most favourable tetraquark has (I)J P = (0)1 + . Other attractive configurations include (I)J P = (0)0 + , (0)2 + containing a (6, 1) anti-diquark and (I)J P = ( 1 2 )0 + , ( 1 2 )1 + from picking the I = 1 2 components of the (3, 0) anti-diquark. However, no signs of these are seen in any of the computed spectra in the many doubly-charmed channels we studied. That we find no significant deviation between the spectra including and excluding tetraquark operators is consistent with the results presented in Ref. [30] which computed the spectrum in the (I)J P = (0)1 + channel. That study used meson-meson and tetraquark operators but, because the operator basis was more restricted than ours, was unable to extract all of the multiple levels which correspond to degenerate meson-meson levels in the non-interacting limit. Computations presented in Ref. [28] find an attractive interaction in the (I)J P = (0)1 + channel using a less direct approach in which lattice QCD computations are used to extract a potential which is then used to determine scattering amplitudes. They do not find a bound state or resonance for a range of light quark masses corresponding to m π = 410 − 700 MeV and conclude that this attractive interaction gets stronger with decreasing pion mass, further motivating studies of how the results vary as the light quark mass decreases towards the physical point.
In one-gluon exchange models, the colour-colour spin-spin interaction is always repulsive for the cc diquark, but the repulsion is suppressed by the quark mass which suggests that doubly-bottomed tetraquarks may be more favourable than doubly-charmed tetraquarks. This is supported by lattice QCD calculations of finite volume-spectra us-ing bases of meson-meson and tetraquark-like operators which suggest the existence of a (I)J P = (0)1 + doubly-bottomed tetraquark [35]. Further support comes from lattice calculations of the potential between two static bottom quarks in the presence of two light antiquarks [36,[38][39][40]42]. This potential is found to lead to a bound state with (I)J P = (0)1 + . Our doubly-charmed (I)Λ P = (0)T + 1 spectrum is not inconsistent with there being an attractive interaction although there were no obvious signs of a bound state in this channel. This is also consistent with recent phenomological studies [60][61][62] which suggest the doubly-bottom tetraquark is bound and the doubly-charmed tetraquark is unbound. Further calculations using bottom quarks and a Lüscher analysis would be of interest. Computations involving the bottom quark with the fermion action used in this study are not straightforward since discretisation effects would be large. It is possible to implement the bottom quark with alternative actions such as Non-Relativistic QCD but this is beyond the scope of this study.
Overall, our study has improved on previous lattice QCD investigations of tetraquarks in two ways. The first is that we use a diverse set of tetraquark and meson-meson operators so that we can reliably obtain a large number of energy levels in each channel and, for the first time in a lattice QCD calculation, robustly extract the multiple energy levels associated with meson-meson energy levels which are degenerate in the non-interacting limit. 7 This is important for future spectrum calculations involving the scattering of mesons with nonzero spin. The second is that we have computed spectra in a large number of channels proposed to contain the hypothetical lightest tetraquark multiplet -some of these channels have not been studied before.

Summary
We have described the construction of a general class of operators resembling compact tetraquarks which have a range of different diquark-anti-diquark structures, transform irreducibly under the symmetries of the lattice and respect other relevant symmetries. As a first demonstration, these operators have been used in conjunction with meson-meson operators to compute correlation functions in the isospin-1 hidden-charm and doubly-charmed sectors using the distillation framework. Finite-volume spectra were extracted by analysing the correlation functions with the variational method. It was found that the addition of tetraquark operators to a basis of meson-meson operators did not significantly affect the finite-volume spectrum and subsequently, would not affect the scattering amplitudes. Because a diverse set of operators was used, for the first time we were able to reliably extract the multiple energy levels associated with degenerate non-interacting meson-meson levels. In all channels, we find that the number of energy levels is equal to the number of noninteracting meson-meson levels expected in the energy region considered and the majority of energies were at most slightly shifted from the non-interacting levels. Hence, there are no strong indications that there are any bound states or narrow resonances present.
This study sets out the groundwork and technology for future work. Calculations with larger lattice volumes and/or at non-zero overall momentum would be necessary to reliably determine scattering amplitudes via the Lüscher method and so rigorously discern the bound state and resonance content in the various channels. In addition, there is strong motivation to study how the results change when the light quark mass is varied. Calculations with lower quark masses would require large volumes. As discussed in Ref. [50], to maintain a given smearing radius in the distillation approach, the number of distillation vectors used must scale with the volume. Because tetraquark elementals are of rank-4, increasing the number of distillation vectors vastly increases the computational cost and so, to make the calculations feasible, an extension of the distillation framework would be required, for example, a stochastic version [63] or an alternative basis of vectors. Calculations with more distillation vectors would also enable the dependence of the results on the degree of tetraquark-operator smearing to be investigated further. Whilst the extracted spectra in the channels studied did not significantly change upon the addition of tetraquark operators to the operator basis, there are many other channels where tetraquarks have been suggested to exist and more detailed lattice QCD investigations are of interest, for example, in the isospin-0 hidden-charm sector, the open-charm sector, the bottom sector and the light scalar mesons.
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 Oak Ridge National Lab, the NSF Teragrid at the Texas Advanced Computer Center and the Pittsburgh Supercomputer Center, as well as at Jefferson Lab.

A Diquark and tetraquark operators
In this appendix we present some additional details of the diquark and tetraquark operators. Consider a diquark operator, δ( x, t) = CGs q T ( x, t)(CΓ)q( x, t), where CGs refers to the Clebsch-Gordan coefficients in Equation (2.1) and various indices have been suppressed. Under proper Lorentz transformations, this operator transforms in the same way as the analogous fermion bilinear and the continuum spin, J, of the diquark for different choices of Γ is given in Table 3. Under a parity transformation, the operator transforms to Pδ( x, t)P −1 = η P q T (− x, t)CΓq(− x, t), where η P = ±1 is a parity factor that depends on the gamma matrix Γ as given in the table. The analogous anti-diquark operator has the same transformation properties. The parity of the tetraquark operator is the product of the parity factors of the diquark and anti-diquark operators.
Taking the Hermitian conjugate of the diquark, we obtain δ † = h Γ s C s Fδ , where the symmetry of the Dirac gamma matrix h Γ = ±1 is shown in Table 3 and s = ξ 1 ξ 3 are phases arising from the exchange symmetry of the colour (s C ) and flavour (s F ) Clebsch-Gordan coefficients of the diquark operator. The phase ξ 1 = ±1 arises from reversing the order of the SU(3) irreps, and the phase ξ 3 = ±1 arises from complex conjugating the irreps, We use the phase conventions of Refs. [67,68]. For tetraquark operators with G-parity symmetry as in Equation (2.4), the G-parity is given by G =Gξ J ξ 1 ξ 3 where ξ J = (−1) J 1 +J 2 −J is the phase arising when the arguments of the SU (2) Clebsch-Gordan coefficients are interchanged and ξ 1 and ξ 3 here arise from the exchange symmetry of the SU (3) F Clebsch-Gordan coefficients of the tetraquark operator. When the flavour irreps of the quarks in the diquark are identical, the overall colourflavour-spin coupling in the diquark must be antisymmetric due to Fermi statistics. To see this schematically, consider the diquark δ = C ab q a q b with overall coupling coefficients C. If C ab is symmetric, the sum would be exactly zero because q a q b is antisymmetric. The symmetry arising from spin (CΓ) αβ = s Γ (CΓ) βα is given in Table 3 and the symmetries arising from colour and flavour are discussed above. Table 3. For different Dirac gamma matrices, we show the notation Γ used to denote the gamma matrix, the continuum spin J, the parity factor η P , the hermiticity factor h Γ and the spin coupling symmetry s Γ .

B One-gluon exchange model
In a simple one-gluon exchange model of a diquark [69], the two quarks interact via a colour-colour spin-spin interaction term, where A 12 is a model-dependent mass term that behaves like 1/m 1 m 2 in the heavy quark limit, λ are the Gell-Mann matrices that span the Lie algebra of SU(3) C and S is the spin of the quark. The relative factors that arise for various colour irreps R and spin S are given in Table 4. It can be seen that the most attractive diquark is the (R, S) = (3, 0) configuration. Similarly, the most attractive anti-diquark is the (3, 0) configuration. Hence a scalar J P = 0 + tetraquark is expected to be the most favourable. Whilst other configurations are less favourable, this one-gluon exchange interaction is suppressed by the masses of the quarks such that in the heavy quark limit, a rich spectrum of tetraquark states with J P = 0 + , 1 + , 2 + is expected to be observed in models such as Ref. [27].
In the case when the flavour irreps of the quarks within the diquark are identical, Fermi symmetry constrains the number of possible configurations. If the flavour irrep is antisymmetric, then the only allowed diquarks are the attractive configurations, (3,0) and (6,1). On the other hand, when the flavour irrep is symmetric the only allowed diquarks are the repulsive configurations, (3,1) and (6,0). A consequence of this is that the doubly-charmed cc diquark is always repulsive with the least repulsive diquark being (3,1). However, the repulsive interaction is suppressed by the quark mass and so it is expected that such tetraquarks may exist in the heavy quark limit [70].

C Non-relativistic quark model
In a non-relativistic quark model, diquark states at rest with orbital angular momentum L and spin angular momentum S coupled to total angular momentum J can be constructed as,  where b † α ( q) is a creation operator for a quark of momentum q and J z component α, and f nL (| q|) is a model-dependent wavefunction that is determined by some interaction potential and is specified by L and the principal quantum number n. Annihilating this state with the field expansion of the diquark operator, we obtain, where u is a Dirac spinor. Expanding u in the non-relativistic limit where | q| is much smaller than the mass of the quark, we find to leading order for Γ = γ 5 ,  Table 5. 1 γ 5 γ 0 γ 5 γ 0 γ i γ i γ 0 γ 5 γ i [γ i , γ j ] qq( 2S+1 L J ) 3 P 0 1 S 0 1 S 0 -3 S 1 3 S 1 3 P 1 1 P 1

D Operator lists
The interpolating operators used to calculate the spectra are listed in Table 6 for the hidden-charm sector and Table 7 for the doubly-charmed sector.  Table 6. The interpolating operators used to calculate the spectra in the isospin-1 hidden-charm sector. For the tetraquark operators, we use the notation δ Γ1