Study of doubly heavy tetraquarks in Lattice QCD

We present results of a lattice calculation of tetraquark states with quark contents $q_1q_2\bar{Q}\bar{Q}, \, q_1,q_2 \subset u,d,s,c$ and $Q \equiv b,c$ in both spin zero ($J=0$) and spin one ($J=1$) sectors. These calculations are performed on three dynamical $N_f = 2 + 1 + 1$ highly improved staggered quark ensembles at lattice spacings of about 0.12, 0.09 and 0.06 fm. We use the overlap action for light to charm quarks while a non-relativistic action with non-perturbatively improved coefficients with terms up to $\mathcal{O}(\alpha_s v^4)$ is employed for the bottom quark. While considering two heavy quarks as charm or bottom, we calculate the energy levels of various four-quark configurations with light quark masses ranging from the physical strange quark mass to that of the corresponding physical pion mass. Results for the spin one states show the presence of ground state energy levels which are below their respective thresholds for all the light flavor combinations with both doubly heavy quarks and particularly for the bottom quarks. Further, we identify a trend that the energy splittings, defined as the energy difference between the ground state energy levels and their respective thresholds, increase with decreasing the light quark masses and are maximum at the physical point for all the spin one states. The rate of increase is however dependent on the light quark configuration of the particular spin one state. We also present a study of hadron mass relations involving tetraquarks, baryons and mesons arising in the limit of infinitely heavy quark and find that these relations are more compatible with the heavy quark limit in the bottom sector but deviate substantially in the charm sector. The ground state spectra of the spin zero tetraquark states with various flavor combinations are seen to lie above their respective thresholds.

Theoretical studies of exotic hadrons are not new. Among the exotics, perhaps tetraquarks are the most studied states. Historically, they were introduced by Jaffe [30] as color neutral states of diquarks and antidiquarks 1 in the context of describing light scalar mesons as tetraquarks and later for exotic spectroscopy [31,32]. Subsequently, the diquark picture of tetraquarks was investigated in detail by many authors through various models [18][19][20]27,28]. Phenomenologically, a four-quark state can also be modeled as molecules [33,34], hadroquarkonia [35,36], and also as threshold cusps [37,38], depending on how the four quarks interact mutually.
Though these models are effective with varying degrees in describing these states, it is essential to have a first principles description of these strongly interacting hadrons. Lattice QCD, being a first principles nonperturbative method, ideally provides such an avenue to investigate these states comprehensively. The success of lattice QCD, however, is still limited for these exotic states for multiple reasons. First, almost all such states that are observed lie very close to their threshold energy levels. Though substantial progress has been made for resolving close-by states, it is essential to use novel techniques like distillation [39] that allows for the construction of a large set of operators with the desired overlap onto the ground state which can then be computed using the variational principle [40,41]. Second, to identify a resonance state unambiguously from its noninteracting thresholds, one has to perform the rigorous finite volume analysis [41] of the discrete spectrum on multiple volumes and/or multiple momentum frames. Moreover, these heavy hardons are very much susceptible to discretization error, and a precise statement cannot be made unless one takes a controlled continuum limit of the results obtained at finite lattice spacings. All these issues amount to a very large computationally intensive calculation, which presumably will be carried out in the future but currently is beyond the scope of any lattice group.
Current lattice QCD methods with available computational resources can, however, be a useful tool for studying hadrons which are far below their strong decay thresholds. For example, taking advantage of these methods and available computational resources, one can study the deeply bound multiquark states to investigate whether such state exist in nature. One can employ lattice methodology for a systematic search for these states using various spin and flavor combinations of interpolating operators and then, dialing the quark masses, spanning over a wide range, can study the onset of a stable state with a large binding energy. In fact, it was already speculated several years ago that there may exist deeply bound tetraquark states in the heavy quark limit. Using one pion exchange between the ground state Qq mesons, Manohar and Wise showed that QCD contains stable (under strong interactions) four-quark QQqq hadronic states in the infinite quark mass limit, and for the bottom quark, this binding could well be sufficiently large [42].
The heavy tetraquarks were also studied recently using heavy quark effective theory [24,43], quark models [23,[44][45][46][47][48][49][50], QCD sum rules [51][52][53], and large N c calculations [54][55][56]. 2 The proposed doubly bottom tetraquark state and its isospin cousins are believed to be strong interaction stable states with relatively long lifetimes. Recently, lattice QCD calculations [25,57] and a lattice-QCD-potential based study [58][59][60] also identified a particular exotic flavor-spin combination of two bottom quarks, namely udbb, with a prediction of a deeply bound state which lies below its noninteracting two-meson threshold. It is thus quite crucial to investigate such and similar states using a detailed lattice QCD study by incorporating various heavy and light flavor combinations along with different spin combinations and at multiple lattice spacings.
In this work, we carry out such a calculation in which we use both the charm and bottom as heavy quarks and then vary the light quark masses from the strange quark mass to the corresponding lower pion masses leading to various tetraquark states: q 1 q 2QQ ; q 1 ; q 2 ⊂ u; d; s; c and Q ≡ b, c with both spin 0 (J ¼ 0) and spin 1 (J ¼ 1). These are computed at three lattice spacings of approximately 0.12, 0.09, and 0.06 fm to investigate the discretization effects on these heavy hadrons. We use the relativistic overlap action, for light to charm quarks, while a nonrelativistic action with nonperturbatively improved coefficients with terms up to Oðα s v 4 Þ is employed for the bottom quark. Our results for the spin-1 tetraquarks indicate the presence of energy levels below the respective thresholds for all light flavor combinations with doubly heavy, in particular, for doubly bottom, quarks. The results for spin-0 tetraquarks, which are the flavor symmetric cousin states of the spin-1 counterparts, however, indicate the respective energy levels are above their lowest strong decay two-meson thresholds. In addition to computing the ground state spectra, we also present a lattice study of the hadron mass relations between tetraquarks, heavy baryons, and mesons arising from the heavy quark symmetry. In the future, we will incorporate also the finite volume study so that more quantitative conclusions about the pole structures of these tetraquark states can be made, particularly for the near-threshold states.
The paper is organized as follows. In Sec. II, we elaborate the lattice setup, actions employed, and quark mass combinations that we use for this work. Section III provides details of the tetraquark operators and the flavorspin combinations that we employ in this work. In Sec. IV, with the details of analysis method, we present our results, first for the spin-1 sector followed by the spin-0 sector. Finite volume effects on our results are discussed thereafter. A discussion on the hadron mass relations with the heavy quark symmetry follows afterward. Finally, conclusions from this work are discussed in Sec. V.

II. LATTICE SETUP
We perform this calculation on three dynamical 2þ1þ1 flavors lattice ensembles generated by the MILC Collaboration [61]. These ensembles, with lattice sizes 24 3 × 64, 32 3 × 96 and 48 3 × 144, at gauge couplings 10=g 2 ¼ 6.00, 6.30 and 6.72, respectively, were generated with the highly improved staggered quark (HISQ) action and with the one-loop, tadpole improved Symanzik gauge action with coefficients corrected through Oðα s a 2 ; n f α s a 2 Þ [62]. The masses of strange and charm quarks on these ensembles are set to their physical values, while the light sea quark masses are set such that m s =m l ¼ 5. The lattice spacings as measured using the r 1 parameter for the set of ensembles used here are 0.1207(11), 0.0888 (8), and 0.0582(5) fm, respectively [61]. Further details of these lattice QCD ensembles can be found in Ref. [61].
In the valence sector, for light, strange, and charm quarks, we employ the overlap fermion action [63,64], which has exact chiral symmetry at finite lattice spacings [63][64][65] and is automatically OðmaÞ improved. The numerical implementation of the overlap fermion is carried out following the methods in Refs. [66,67]. A wall source smearing is utilized to calculate the light to charm quark overlap propagators on Coulomb gauge fixed lattices. In Table I, we list the quark masses and corresponding pion masses that we use for this calculation. The strange quark mass is tuned by equating the lattice estimate of thess pseudoscalar meson mass to 688.5 MeV [68][69][70]. We follow the Fermilab prescription of heavy quarks for tuning the charm quark mass [71]. We tune it by equating the spin-averaged kinetic mass of the 1S charmonia [aM kin ð1SÞ ¼ 3 4 aM kin ðJ=ψÞ þ 1 4 aM kin ðη c Þ] to its experimental value, 3068.6 MeV [22]. The tuned bare charm quark masses are found to be 0.528, 0.427, and 0.290 on coarse to fine lattices, respectively, all of which satisfy m c a ≪ 1, ensuring reduced discretization artifacts in this calculation. Details on the charm quark mass tuning can be found in Refs. [69,70]. For the bottom quarks, we employ a nonrelativistic QCD (NRQCD) formulation [72]. In the NRQCD Hamiltonian, we include all the terms up to 1=M 2 0 as well as the leading term of the order of 1=M 3 0 , where M 0 ¼ am b is the bare mass of the bottom quarks in lattice units [73]. The bottom quark propagators are obtained by the usual time evolution of the NRQCD Hamiltonian, H ¼ H 0 þ ΔH, where the interaction term, ΔH, is given by Here, c 1 …c 6 are the improvement coefficients, and for the fine lattice, we use their tree level values, while for the coarser two lattices, we employ their nonperturbative values as estimated by the HPQCD Collaboration [74] on the same set of lattices. To tune the bottom quark mass, we first calculate the kinetic mass of the spin average 1S bottomonia, from the relativistic energy-momentum dispersion relation aM Kin ¼ ððapÞ 2 − ðaΔEÞ 2 Þ=ð2aΔEÞ, and then equate it with its experimental value. Details on the bottom quark mass tuning is given in Ref. [75]. With this setup of light, strange, charm, and bottom quark propagators, we proceed to calculate the tetraquark correlators from the interpolating fields with various flavorspin combinations that we discuss in the next section.

III. FOUR-QUARK INTERPOLATING OPERATORS
In this section, we describe four-quark interpolating fields (operators) that we employ in this work. We construct these operators with two heavy and two light quarks and with the total spin J ¼ 0 and 1. As in Ref. [25], for both spins, we construct two types of operators, with the goal that one overlaps onto a tetraquark state of given quantum numbers and the other one overlaps onto the lowest strong decay two-meson states of the same quantum numbers. The tetraquark-type operators are constructed using the diquark prescription of Jaffe [31,32] where a color neutral hadronic operator is constructed as a product of diquarks and antidiquarks. These diquarks (antidiquarks) can be in thē 3 c ð3 c Þ or 6 c ð6 c Þ of the color SU(3) irreducible representation (irreps). Phenomenologically, the one gluon exchange model [31,32] favors an attractive interaction of two quarks and is in the3 irrep of SU (3). In this work, we construct tetraquark operators with both irreps of SU(3). In the spin J ¼ 1 sector, we use diquarks and antidiquarks with the following configuration: ðQ;QÞ → ð3 c ; 1; F s Þ: ð3Þ The light quark (l 1 ; l 2 ; l 1 ≠ l 2 ) combinations are constructed with color, spin, and flavor degrees of freedom (d.o.f.) as antisymmetric and are restricted within ⊂ ðu; d; s; cÞ. The heavy quark combination (Q,Q) is constructed with color antisymmetric 3 c , forced by (l 1 , l 2 ) being in the3 c , and since flavor is manifestly symmetric, the spin is also symmetric. This combination is restricted to only heavy flavors ⊂ ðc;bÞ with a further restriction of Q ≠ l 1 ≠ l 2 . With these diquarks and antidiquarks, a spin-1 tetraquark-type operator of flavor (l 1 l 2QQ ) is constructed as The label x is a shorthand notation for (⃗ x, t), where ⃗ x is the spatial local site and t is the time slice. We then construct the two-meson-type operators corresponding to each flavor of the (l 1 l 2QQ ) tetraquark operator, T 1 ðxÞ, with the appropriate flavor antisymmetry as The tetraquark operator T 1 ðxÞ is related to the two-meson product M 1 ðxÞM Ã 2 ðxÞ via a Fierz transformation, and the relation is explicitly shown in the Appendix of Ref. [76] with the appropriate change in flavor labels. The various flavor and isospin (I) combinations that we explore for these spin-1 tetraquark-type and two-meson-type operators are tabulated in Table II. For the spin-0 sector, we employ following diquark/ antidiquark configuration in which both diquarks are with spin 0: The combination (l, l) being manifestly flavor symmetric requires the color d.o.f. to be in the 6 c . For the combination (Q,Q), the color d.o.f. is consequently restricted to6 c , while the flavor d.o.f. is manifestly symmetric. In the above expression, for the combination (l, l), we incorporate the flavors (u, s, c), while both c and b are used for Q. A spin-0 tetraquark-type operator of flavor (llQQ) constructed from the product of the aforementioned diquarks and antidiquarks is given by As previously done, we also construct a two-meson-type operator with the same quantum number as that of (llQQ), and it is given by In Table III, we tabulate the spin-0 tetraquark configurations with the possible flavor combinations with the above flavor-spin configurations.
With the operators so constructed, we proceed to compute the correlation matrices of all the possible combinations of these operators for a given spin and flavor and then extract the associated energy states from the generalized eigenvalue solutions. In the next section, we discuss this in detail. II. The tetraquark-type and two-meson-type operators that we study in this work with possible flavor combinations and allowed isospin (I) in the spin-1 sector. The last column shows the range of pion masses that we use for the light quarks on the coarsest lattice spacing.
TABLE III. The tetraquark-type and two-meson-type operators for various flavors of in the spin-0 sector. The range of pion masses used for uubb and uucc states is indicated in the last column. All other states are computed at their physical quark mass.

IV. RESULTS
In this section, first we elaborate the analysis procedure that we utilize to extract the energy levels from the matrix of correlation functions constructed from the interpolating fields mentioned above. Results obtained will be discussed after that.

A. Analysis methods
To evaluate the energy levels corresponding to the operators discussed in Sec. III, we first construct a correlator matrix of these operators and then use the variational method [40,41]. This matrix of correlation functions C ij ðtÞ is given as where the operator O i ð⃗ x; tÞ ∈ fT k ð⃗ x; tÞ; M k ð⃗ x; tÞg is either a tetraquark-type operator or a two-meson-type operator of a particular spin k. For the spin-1 tetraquark states, O i 's correspond to Eqs. (4) and (5), whereas for the spin-0 states, these are from Eqs. (7) and (8). We analyze each spin sector separately. After constructing the correlation matrix, CðtÞ, for a given spin and flavor combination, we solve a generalized eigenvalue problem (GEVP) to obtain the two energy levels [40,41]. The standard methods for GEVP [40,41,77,78] are typically suited for a Hermitian correlator matrix. We note that, since we are using a wall source, the correlator matrix is non-Hermitian. 3 Hence, we employ a variation of the GEVP method, named the eigenvector method, involving eigenvector projection in evaluating the ground state energies [79]. The method involves using the left and right eigenvectors of the correlator matrix to construct the principal correlator as discussed below: (1) Compute left and right eigenvectors of the correlator matrix CðtÞ at chosen time slices (t 1 , t 0 ) as The time slices (t 1 , t 0 ) are chosen such that t 1 =t 0 > 2 and t 1 chosen in the region where the correlator is expected to be dominated by the ground state. (2) The eigenvectors v L;R;n ðt 1 ; t 0 Þ are then used to construct the principal correlator as and the effective masses are then obtained from m n;eff ¼ logðΛ n ðtÞ=Λ n ðt þ δtÞÞ.
For a Hermitian correlator matrix, the left and right eigenvectors will be identical, and hence this method will be the same as standard methods [40,41,77,78]. For a non-Hermitian correlator, the source and sink operators are accordingly rotated by the left and right eigenvectors, respectively. To check the effects of non-Hermiticity, we also solve the GEVP with the standard methods [40,41,77,78]. We find that results obtained with either GEVP methods are consistent with each other, while the results from the eigenvector method is observed to be more stable.
The principal correlators thus obtained correspond to two energy levels, and the ground state energy is computed from the lowest one. On the other hand, we calculate the noninteracting two-meson threshold from the sum of the ground state masses of the two mesons involved. We then compare the lowest energy level obtained from the GEVP solution with the noninteracting two-meson threshold and evaluate the energy splitting between them as where E T k is the ground state energy obtained from the principal correlator of the GEVP, while is the energy of the noninteracting two-meson (M 1 and M 2 ) threshold. The above energy splitting (ΔE k ) can be evaluated directly by fitting the two datasets separately and then computing the difference on each resample. Alternatively, this can also be evaluated by taking the jackknife ratio of the principal correlator [ΛðtÞ] of the GEVP to two-meson correlators, M 1 ðtÞ × M 2 ðtÞ), as A fit to the ratio correlator [Λ 0 ðtÞ] will then directly yield the energy splitting with respect to the relevant threshold. Such a construction offers the advantage of reducing the systematic errors through jackknifing. However, in using such an effective correlator, caution must be exercised as this construction can produce spurious effects since the saturation of the ground states of the numerator and the denominator may not happen at the similar time slices. In this work, in estimating the energy splitting, we utilize both the direct and ratio methods and find consistent results. However, as expected, we find smaller uncertainties in the ratio method. We now present the results obtained through the above-mentioned analysis.
B. Spin one tetraquarks J P = 1 + We begin with presenting data for the spin-1 doubly bottom tetraquark states. As described earlier, we compute a matrix of correlation functions of the tetraquark T 1 ðxÞ and two-meson operators M 1 ðxÞ. The diagonal correlators of this matrix correspond to the same source-sink operators, while the off-diagonal correlators have a tetraquark operator at the source and a two-meson operator at the sink and vice versa. The correlator matrix is non-Hermitian, and as mentioned earlier, in obtaining our final results, we employ the eigenvector method of diagonalization.
As a representative plot on analysis, in Fig. 1, we show the effective mass of the lowest energy level obtained from such a diagonalization along with the effective mass of the noninteracting two-meson threshold correlator for the case of usbb. The data in orange are the effective mass of the noninteracting two-meson correlator, which in this case is obtained from the product of the correlators of the B and B Ã s mesons. 4 The data in green are the effective mass of the lowest eigenvalue (the ground state), which is clearly below the effective mass of the threshold correlator. We also find that the effective mass corresponding to second eigenvalue overlaps with the effective mass of the threshold correlator in its approach to the plateau. However, as expected, it is noisier and needs a bigger basis of operators to extract it reliably. As discussed previously, for each flavor combination, we calculate the energy splitting ΔE 1 directly from Eq. (12) by fitting the individual correlators as well as from the ratio of correlators using Eq. (13).
Following the above procedure, we calculate the energy splittings (ΔE 1 ) for all the doubly bottom tetraquarks with various flavor combinations mentioned in Table II. This is performed on three different lattices (a ∼ 0.12, 0.09, and 0.06 fm), and on each one, we vary the light quark masses over a wide range as listed in Table II. In Fig. 2, we show these results; in the left panel, we plot these energies computed at various pion masses. The results for the flavor combinations, uqbb with q ∈ ðd; s; cÞ, are shown by red, green, and blue colored data, respectively. For a representative plot, we choose to show results at the coarse lattice spacing since here we have the maximum number of pion masses and therefore can show the pion mass dependence of these energy splittings (ΔE 1 ) more prominently. The result for the udbb state exhibits larger uncertainties at lower pion masses due to the presence of two light quarks, while the state usbb allows us to extract results at much lower pion masses. For ucbb, we could extract results even at the physical light quark mass.
It can be noted that for all the flavor combinations there is a trend of increment of ΔE 1 with the lowering of pion masses, and we will discuss the details shortly. The availability of a large number of data points allows us to perform the chiral extrapolation very reliably. At each lattice spacing, we first perform the chiral extrapolation of ΔE 1 and then perform a continuum extrapolation from the results obtained at three lattice spacings. We use the following simple quadratic Ansatz for both chiral and continuum extrapolations: Here, the label k for the spin is kept general since we will also use these Ansätze for both spin sectors. We perform two fittings: one including all data points to show the pion mass dependence over a wide range of pion masses and the other with only the lower few pion masses to perform the chiral extrapolation. The fit results are shown in Table IV, in which in the second column we show the relevant slope parameter labeled as c 1;m π 2 , which is indicative of the pion mass dependence of the energy splitting ΔE 1 . It is instructive to compare c 1;m π 2 parameters for different tetraquark states with different flavor combinations at a given lattice spacing. The fits indicate that the state udbb exhibits the most pronounced trend in the increase of ΔE 1 , followed by the state usbb, while the state ucbb exhibits a very minute variation. The results at the finest lattice spacings do not indicate such a clear trend as we do not have data points at much lighter pion masses at this lattice spacing.
For the second fit, i.e., for the chiral extrapolation, we use the Ansatz in Eq. (14) and employ cuts on the largest pion masses and include data corresponding to as low pion masses as can be afforded by meaningful uncertainties in the extrapolation. The results of the chiral extrapolation are shown in Table IV  shown in the last column. We then use these chirally extrapolated ΔE 1 j m phys π from three different lattice spacings and perform a continuum extrapolation using the Ansatz in Eq. (15). The results of this extrapolation are shown in the right panel of Fig. 2, and the fit results are listed in Table V. The slope parameter c 1;a 2 in this case will be an indicator of the lattice spacing dependence of the particular state. For udbb and usbb, these are consistent with zero, indicating no dependence on lattice spacing. The parameter c 1;a 2 for the ucbb state indicates a mild dependence on the lattice spacing. The state scbb, which is the SU(3) symmetric state of ucbb, requires no chiral extrapolation since all quark masses are at their physical values. The corresponding lattice spacing dependence parameter, c 1;a 2 , as shown in Table V, indicates no dependence on lattice spacing of this state. The continuum extrapolated results ΔE 1 j m phys π a¼0 are shown in Fig. 3.
It can be noted that at the finest lattice spacing the lowest pion mass available is m π ¼ 545 MeV, which may not be low enough for a chiral extrapolation. Because of this, the chirally extrapolated results at this lattice spacing may have a systematic effect arising from the absence of lower pion masses, and that may reflect in the lattice spacing dependence of some of our findings such as for the ucbb state. Hence, we also report our results without including data  We now discuss the results of the spin-1 doubly charm tetraquarks. In Fig. 4, we show those results where the left panel shows the pion mass dependence and the chiral extrapolation on the coarse lattice. The right panel represents results for the continuum extrapolation. The relevant lowest thresholds for the flavor combinations udcc and uscc are the noninteracting D-D Ã and D-D Ã s mesons, respectively. For both cases, we find an energy level below their relevant strong decay thresholds, while the other energy level appears at the threshold. As in the doubly bottom cases, we calculate the energy splittings [ΔE 1 in Eq. (12)] between the lowest energy levels and the threshold states by direct fitting as well as from the ratio of correlators [as in Eq. (13)]. We represent the fitted results for udcc by red data points, while results for uscc are shown by green points. The fitted results for pion mass dependence and chiral extrapolation are shown in Table IV, while the results for continuum extrapolation are shown in Table V. In the case of udcc, similar to udbb, we observe a trend in the increase of ΔE 1 with the lowering of the light quark constituents. This is evident from the fits for the pion mass dependence and is indicated by the c 1;m π parameter on the coarsest two lattice spacings. The finest lattice spacing results do not clearly indicate this trend due to the lack of lower pion masses at that lattice spacing. The pion mass dependence of the energy splitting for uscc, color coded in green, is much flatter in comparison to udcc, and this trend is reflected in the c 1;m π 2 coefficient. The continuum extrapolations for both udcc and uscc indicate no discernible dependence on the lattice spacing.
In column 4 of Table V, we show the continuum extrapolated results for doubly charmed tetraquarks. Column 5 shows the average results obtained the two coarser lattices. Both columns show the presence of energy levels below their respective thresholds for both udcc and uscc. However, they are very close to their respective strong decay thresholds, as was also observed in Ref. [80]. Because of their close proximity to thresholds, a careful finite volume analysis [41] is needed to make conclusive statements about the nature of these states. Though they could be stable under strong interaction, they may not appear as bound states because of threshold effects.

C. Spin zero tetraquarks J P = 0 +
In the spin-0 sector, we compute the energy levels of the tetraquark states with various flavor combinations that are listed in Table III. These tetraquark states are flavor symmetric cousins of those listed in Table II. As in the case of the spin-1 sector, we compute a matrix of correlation functions consisting tetraquark-type, T 0 ðxÞ, and two-meson-type, M 0 ðxÞ, interpolating fields and employ the eigenvector method of diagonalization in obtaining our final results.
We shall begin by discussing the spin-0 doubly charmed and doubly bottom tetraquark states with I ¼ 1. The effective masses of the principal correlators, obtained from GEVP analysis, for the flavor combination uubb are shown in Fig. 5. This representative figure is obtained on the fine lattice and at the pion mass m π ¼ 688 MeV. The relevant strong decay threshold in this case is the two noninteracting B mesons. The effective mass of the product correlator of two B mesons is represented by the orange data. The effective mass of the lowest eigenvalue, shown in green, is seen to coincide with the threshold correlator. This behavior is in contrast when compared with its flavor antisymmetric partner udbb where there is a clear indication of the ground state level being below the relevant threshold. The energy splitting [ΔE 0 in Eq. (12)] of the tetraquark state uubb is shown in the left panel of Fig. 6 by red colored data points, and results are obtained at various pion masses (on the coarser lattice) to explore the pion mass dependence. We note that the determination of these energy splittings is significantly noisier in comparison to the energy splittings in spin-1 udbb state with the same statistics. This limits us to using much lighter pion masses for uubb. Furthermore, this also forces us to use the entire dataset for exploring both the pion mass dependence as well as the chiral extrapolation. We perform a chiral extrapolation with the Ansatz in Eq. (14) at each lattice spacing, and the results are listed in Table VI. The fits for the parameter c 0;m π 2 indicate a dependence on pion mass for a ¼ 0.1207 fm, and no dependence is seen for the other two lattice spacings, since c 0;m π 2 is consistent with zero. It can be noted that this behavior again is in contrast with the pion mass dependence of the udbb state where a nontrivial dependence was clearly identified. After the chiral extrapolation, we perform the continuum extrapolation using the Ansatz in Eq. (15), and fits are shown in Table VII. The slope parameter c 0;a 2 for the state uubb is consistent with zero, indicating no dependence on the lattice spacing. The physical and continuum extrapolated result for uubb clearly indicates that there is no energy level below its lowest strong decay threshold with any statistical significance and is consistent with zero.
The green data points in Fig. 6 show the results for ΔE 0 (on a ¼ 0.1207 fm lattice) for the spin-0 doubly charmed tetraquarks uucc. In this case, the GEVP solutions also display similar qualitative features, as the corresponding doubly bottom states where the ground state coincides with the threshold and a well-separated second state lies above that. Here, the threshold is that of the two noninteracting D mesons. As in the previous case, we use the entire dataset for the pion mass dependence as well as chiral extrapolation. The chiral extrapolation fits at each lattice spacing shown in Table VI indicate no dependence on the pion mass since the parameter c 0;m π 2 is found to be consistent with zero. The continuum extrapolation for this case, color coded in green, is shown in the right panel of Fig. 6, which indicates a mild dependence on the lattice spacing. The physical and continuum extrapolated results (ΔE 1 j m phys π a¼0 ) are shown in the fifth column of Table VII, and all are found to lie above the respective threshold states. As in the spin-1 case, we have also calculated the average values of these energy splittings from the results obtained on two coarse lattices and show that in the last column of Table VII. With our available quark propagators, we are also able to study I ¼ 0, J ¼ 0 tetraquark states, ssbb, sscc, and ccbb, where the strange, charm, and bottom quark masses are tuned to their physical values. Energy levels obtained for these states will thus be at the physical points, and there is no need for any chiral extrapolation. The thresholds for these states are the noninteracting B s B s , D s D s , and B c B c , respectively. These require only a continuum extrapolation, which is shown in the two panels of Fig. 7, and the fitted results are shown in Table VII. The estimates of the energy splitting ΔE 0 for the state ssbb (color coded in red) show no lattice spacing dependence, and the final result is consistent with zero, indicating the absence of any bound state. For the state sscc, we also find similar results, and the continuum extrapolated result lies above its respective threshold, which is most likely to be a scattering state. Results for the state ccbb indicate a mild lattice spacing dependence, and the continuum result is also most likely be a scattering state. In conclusion, our analysis on the I ¼ 0, spin-0, tetraquarks with flavor combinations ssbb, sscc, and ccbb suggest the absence of any bound state, and the observed energy levels correspond to the scattering states. Recently, a potential based lattice QCD study in Ref. [59] for doubly bottom spin-0 states also concluded the same.

D. Finite volume effects
For all the spin-1 tetraquark states with various flavor combinations listed in Table II, we have found the energy levels below their respective strong decay thresholds. In some cases, the energy splittings (ΔE 1 ) between the ground state and the threshold state are very large, while for others, they are close to and below their respective thresholds. However, all these energy levels are obtained within a single volume of about 3 fm. It is thus necessary to estimate the finite volume effects on these energy differences and obtain their infinite volume estimates, which can then be interpreted as the binding energies of the corresponding bound states. However, repeating these calculations on multiple lattice volumes is computationally very expensive and so is beyond the scope of this work.
However, it is possible to identify a few states for which the finite volume corrections will be suppressed, i.e., could be very small. The estimation of ΔE 1 on a single large enough volume for such a case, in fact, would be close to its binding energy (B ∞ ). As demonstrated in Refs. [81][82][83], the finite volume corrections Δ FV to energy levels corresponding to an infinite volume bound state with energy E ∞ scale as where E FV is the energy level computed a cubic lattice, k ∞ is the binding momentum of the infinite volume state, and (m 1 , m 2 ) are the masses of the two noninteracting particles with the threshold energy m 1 þ m 2 . It should be noted from the above expression that the finite volume effects are suppressed by the threshold mass (m 1 þ m 2 ) and that this suppression is significantly enhanced for the cases in which the threshold states are heavy mesons, such as those we are studying here. In addition to that, if ΔE is also large, then the finite volume corrections will further be suppressed since it also enters in the exponential. Therefore, in the doubly bottom sector, tetraquark states with the flavor combinations, udbb and usbb, for which the ΔE values are found to be more than 150 and 100 MeV, respectively, will have small finite volume corrections. For these cases, it is quite natural to expect that the energy splitting ΔE will be closer to their infinite volume binding energy. Therefore, these states will be stable under strong interactions. However, for the cases, particularly for the doubly charmed tetraquarks, which are below but closer to their thresholds (i.e., ΔE values are closer to zero), it will be difficult to get any qualitative estimate for their finite volume corrections.
In those cases, one needs to perform a detail finite volume study [41] to make any conclusive statement about their infinite volume pole structures.

E. Heavy quark effective theory and hadron mass relations
The Heavy quark effective theory (HQET) is a very useful tool and is often utilized to understand various properties of heavy hadrons including their energy spectra. Using heavy quark symmetries, one can also obtain mass relations between heavy flavored hadrons such as those mentioned in Ref. [24]. Using such symmetry relations, Ref. [24] predicted masses and binding energies of various tetraquarks states including some of those studied in this work. Although such relations are valid in the infinite quark mass limit, they are used at the bottom and even at the charm quark masses. It will therefore be interesting to investigate these relations by a first principles nonperturbative method, such as lattice QCD, with a goal to validate these relations at a given quark mass and access their deviation, if any, from the heavy quark limit. The availability of data on the ground state masses on mesons, baryons, and tetraquarks obtained from this calculation, both at the charm and the bottom quark masses, provides such an opportunity to systematically investigate these relations. Below, we elaborate that.
The work in Ref. [24] states the following relation among the hadrons with heavy quarks, where Q i and q k denote heavy and light quarks, respectively. Here, we use the same notation as in Ref. [24]. The braces f…g and ½… imply the symmetrization and antisymmetrization, respectively, with respect to the flavor d.o.f. In this notation, (fQ i Q j g½q kql ) 5 represents a tetraquark operator with the flavor symmetries indicated by the braces, while (fQ i Q j gq y ), (Q x ½q k q l ) and (Q xqy ) represent a heavy-heavy-light baryon, heavy-light-light baryon, and heavy-light meson, respectively. It should be noted that Ref. [24] provides four such relations depending on the combination of flavor symmetrization/antisymmetrization, and the one shown here corresponds to our operator construction. The relation in Eq. (17) can then be employed  The tetraquark operator used in this work is a complex conjugate of this operator.
to predict the masses of the tetraquark states by substituting the relevant masses of heavy baryons and mesons. In Ref. [24], this was calculated by using the spin average masses of the charmonia, bottomonia, and heavy baryons by inserting their experimental or quark model values.
Here, we aim to study this relation at both the charm and the bottom quark masses. We do not consider the spinaverage mass and instead use the spin-1=2 states for baryons and pseudoscalar mass for the heavy-light meson. Any deviation from the equality in Eq. (17) would be maximum in this choice. In doing so, we will be able to estimate an upper bound of the deviation from the heavy quark limit which originates from all ð1=m Q Þ n corrections. In evaluating Eq. (17), we find it to be convenient 6 to redefine the relation as a ratio which for the charm and bottom quarks is given by In the limit of infinitely heavy quarks, the ratio R Q will be unity. In computing these ratios (R c=b ), we first evaluate the jackknife ratios of the following correlators, which directly provide the difference of masses as shown above. R c=b are then evaluated from the fits to these ratio correlators. In addition, we also fit the individual masses of tetraquarks, mesons, and baryons and calculate R c=b from Eq. (18). We find consistent results with both methods, and the evaluation with Eq. (18) provides improved uncertainties. As we have access to a large number of light quark masses, while keeping the heavy quark mass at the charm and bottom quark, we vary the light quark mass and calculate R c=b for each case. In Fig. 8, we show these results at several pion masses for the coarser lattice (a ∼ 0.12 fm) using the entire dataset in fitting. This is done for other lattice spacings as well. The results clearly indicate a wide separation of ratios between the charm and bottom quarks; while R b is closer to the heavy quark limit of unity, R c deviates from it substantially. After repeating this calculation on the other two lattices, we perform a simplistic chiral and continuum extrapolation according to the Ansatzes in Eqs. (14) and (15). The fit results are shown in Tables VIII and IX at three lattice spacings. For both ratios, R b and R c , we do not observe any appreciable dependence on the pion mass as indicated by the parameter c π 2 in Table VIII. In addition, the continuum extrapolation fit in Table IX does not indicate any lattice spacing   (22) dependence for the bottom and charm quarks. The continuum extrapolated results are listed in the last column of Table VIII; we find R b ¼ 0.837ð38Þ and R c ¼ 0.602ð22Þ. These results clearly indicate that there is a substantial deviation from the heavy quark limit at the charm quark mass, implying there might be large contributions from ð1=m Q Þ n corrections. However, results at the bottom quark mass are much closer to the heavy quark limit. Our results indicate that, as far as the heavy quark symmetry relations such as that is shown in Eq. (17) are considered, the charm quark mass is not heavy enough for the equality, and one certainly needs to incorporate appropriate leading order 1=m Q and then higher order correction terms. However, one can of course use these relations for bottom quarks with higher order 1=m Q corrections.

V. DISCUSSION AND CONCLUSIONS
Recently, there has been tremendous amount of activities in studying multiquark states both theoretically and experimentally. In particular, heavy tetraquarks are being investigated at various laboratories as well as studied theoretically through different models and by lattice QCD calculations. In this work, using lattice QCD, we have performed a detailed study on the doubly heavy tetraquark states with quark contents q 1 q 2QQ ; q 1 ; q 2 ⊂ u; d; s; c, and Q ≡ b, c, in both spin-0 (J ¼ 0) and spin-1 (J ¼ 1) sectors. Not only do we study udbb and usbb, as was studied in Ref. [25], but we also explore ucbb, udcc, and uscc states and additionally include the spin-0 sector of doubly heavy tetraquarks. In doing so, we have presented a systematic dependence of the ground state spectra of such states on their light quark constituents over a wide range of quark masses starting from the quark mass corresponding to the physical pion mass to the strange quark mass. Since all these hadrons involve heavy quarks, naturally, like any heavy flavored hadrons, they are susceptible to heavy quark discretization effects in a lattice calculation. To check the lattice spacing dependence, we have obtained results at three lattice spacings, the finest one being at 0.0582 fm. At a given lattice spacing, we perform a chiral extrapolation using several quark masses and then perform a continuum extrapolation to get the final results. For all the states in the spin-1 sector, we observe the presence of energy levels below their respective two-meson thresholds, the deepest one being for the doubly bottom tetraquark, udbb. Furthermore, for various flavor combinations of the tetraquark states, we find that there is a clear trend of increase in the energy splitting (ΔE) as the light quark masses of such states are decreased and it becomes maximum at the physical quark mass. This energy splitting in the infinite volume limit of such a state can be interpreted as its binding energy. This trend was first indicated in the lattice calculation in Ref. [25] for the states udbb and usbb.
Here, we confirm that over a wide range of quark masses. Additionally, we find that such a trend holds for all the spin-1 states considered here, including the doubly charm tetraquark states. For the doubly charmed tetraquark states, udcc and uscc, we also find that the ground states are below their respective thresholds. However, they are quite close to their thresholds, which was also observed in Ref. [80]. Though they could be stable under strong interactions, one needs to carry out finite volume analysis to establish their bound state properties, if there are any. We would also like to point out that most of these states, except ucbb, show either no discernible dependence or very mild dependence on lattice spacing. However, this will be clear when in a future study we include much lower pion masses on the fine lattice. Our final results for doubly heavy spin-1 tetraquarks states from this calculation are summarized in Table X. Our estimates for the udbb and usbb are in agreement with those of Ref. [25] at a lattice spacing (approximately 0.09 fm) where both of our data points are available.
We also provide a comparison of global results of spin-1 doubly heavy tetraquark states with various flavors and show that in Fig. 9. The results from Refs. [23,24] are based on HQET and the potential model, respectively, while the rest are lattice calculations. All results agree with the existence of deeply bound spin-1 tetraquark states, udbb and usbb, which are stable under strong interactions. Our results for the doubly bottom states agree well with those from the HQET predictions [24] as well as that of the result in Ref [25] at similar lattice spacings (approximately 0.09 fm). Reference [25] used N f ¼ 2 þ 1 PACS-CS gauge field configurations and Coulomb gauge fixed wall sources with clover action in the valence sector. The results were extracted at a single lattice spacing (a ∼ 0.09 fm) at three pion masses, and a chiral extrapolation with m 2 π was performed to obtain the final result. The results from Ref. [60] were obtained from the potential based lattice QCD study in which potentials of two B mesons were computed in the static approximation for various spinisospin combinations. These were then fitted to a phenomenologically motivated Ansatz and were further used to solve a Schrödinger equation to determine a bound state. These calculations were performed at three pion masses ranging from m π ∼ 340-650 MeV, and the final results were obtained after chiral extrapolation. Reference [80] used an anisotropic N f ¼ 2 þ 1 clover action, and results were obtained at a single lattice spacing (a t ∼ 0.0035 fm with anisotropy 3.5) and at a single pion mass (m π ¼ 391 MeV). For the doubly charm states, our results are in disagreement with those from the HQET results [24]. As we showed earlier, this discrepancy is due to the deviation of HQET relations at the charm quark mass. Inspired by the results in the spin-1 sector, we also explore the spin-0 tetraquark states with doubly bottom as well as with doubly charm quarks. Here, we have computed flavor symmetric uubb and uucc states and also explored the pion mass dependence by dialing the light quark mass. To check the lattice spacing dependence of the observed results, we perform the calculation on three different lattice spacings. In addition, we have also computed the following flavor symmetric states, namely, ssbb, sscc, and ccbb at the physical strange, charm, and bottom quark masses. For the doubly bottom state uubb, we find that the energy splittings (ΔE 0 ) are generally noisy and do not clearly exhibit a trend of increase in ΔE 0 as the pion mass is lowered. Contrary to the results of its flavor antisymmetric cousin udbb, the ground state energy of uubb coincides with its threshold at lower pion masses with no clear indication of any level below the threshold. For the doubly charm state, uucc, the extracted energy levels clearly lie above their respective thresholds with no discernible dependence on pion mass, again contrary to the results of its flavor antisymmetric cousin udcc. In performing the continuum extrapolation, no lattice spacing dependence is observed for the uubb state, while the uucc exhibits a mild dependence on the lattice spacing. The flavor symmetric states ssbb, sscc, and ccbb exhibit similar qualitative features in that all the energy levels are found to be above their respective thresholds and no significant lattice spacing dependence is observed in the continuum extrapolation. Our final results for the spin-0 sector are shown in Table XI. In conclusion, the states in the spin-0 sector do not indicate energy levels below their thresholds, suggesting it is very unlikely that there exists any doubly heavy bound tetraquark states with spin 0.
The availability of energy values of spin-1 tetraquark states for a large number of light quark masses provide us an opportunity to investigate the mass relations [Eq. (17)] between different heavy flavored hadrons due to the heavy quark symmetry, as mentioned in Ref. [24]. For this, we redefine the relation as a ratio (R) between different hadron masses [Eq. (18)] where a value of unity justifies the validity of such a mass relation, and any deviation from unity indicates the amount of breaking of the heavy quark symmetry at a given heavy quark mass. We find that for bottom quarks R b ¼ 0.837ð38Þ, indicating that the bottom quark is very close to the heavy quark limit. On the contrary, at the charm quark mass, we find R c ¼ 0.602ð22Þ, which substantially deviates from the heavy quark limit. This clearly suggests that the charm quark is not heavy enough to impose heavy quark symmetry relations among hadron masses such as in Eq. (17); i.e., as far those mass relations are concerned, one needs to be careful while treating the charm quark within HQET.
The tetraquark states studied in this work are computed in a single volume. In order to make conclusive statements about their scattering amplitudes and complex poles, one needs to carry out similar studies on multiple volumes followed by a finite volume analysis [41]. Such an analysis will especially be useful for the states which are close to their thresholds. However, a comprehensive finite volume analysis for a calculation that is reported here requires significantly large computational resources. Currently, that is beyond the scope of this work, but we intend to pursue such finite volume analysis in the near future. However, it is worth noting that the finite volume corrections for many heavy tetraquarks, particularly for which ΔE values are large, will be substantially suppressed. This is because, as has been pointed out before [81][82][83], such corrections to the observed energy splitting are suppressed not only because of its large value but also for the large masses of the threshold states, which are two heavy mesons in these FIG. 9. Comparison of global results on the spin-1 doubly bottom and charm tetraquark states with various flavor combinations. ΔE is the energy difference between the ground state and the lowest strong decay threshold. Various flavor combinations represented on the horizontal axis are color coded as blue, green, red, magenta, and grey for the states udbb, usbb, ucbb, udcc, and uscc, respectively. cases. It is thus expected that such tetraquark states will be stable under strong interactions. Other errors related to our calculations, namely, unphysical sea quark mass, quark mass tuning, scale setting, mixed action effects, and excited state contamination, together will be much smaller compared to the statistical error [84], and the conclusion reached here will be unaffected by those. It will therefore be very useful to search experimentally spin-1 doubly heavy tetraquarks particularly with two bottom quarks, such as udbb. However, it is very unlikely that there exists any doubly heavy bound tetraquark state with spin 0.

ACKNOWLEDGMENTS
We are thankful to the MILC Collaboration and in particular to S. Gottlieb for providing us with the HISQ lattices. We would like to thank R. V. Gavai for discussions, particularly on HQET relations, and M. Hansen for the discussions on finite volume corrections. We also thank R. Lewis

APPENDIX
In Table XII, we show the energy splittings, ΔE, defined as the difference between the threshold energy and the ground state energy levels, of tetraquark states with various flavor-spin combinations as studied in this work.