A lattice investigation of exotic tetraquark channels

We perform an $n_f=2+1$ lattice study of a number of channels where past claims exist in the literature for the existence of strong-interaction-stable light-heavy tetraquarks. We find no evidence for any such deeply-bound states, beyond the $J^P=1^+$, $I=0$ $ud\bar{b}\bar{b}$ and $I=1/2$ $ls\bar{b}\bar{b}$ states already identified in earlier lattice studies. We also describe a number of systematic improvements to our previous lattice studies, including working with larger $m_\pi L$ to better suppress possible finite volume effects, employing extended sinks to better control excited-state contamination, and expanding the number of operators used in the GEVP analyses. Our results also allow us to rule out several phenomenological models which predict significant tetraquark binding in channels where no such binding is found.


I. INTRODUCTION
Results from multiple recent lattice studies [1][2][3][4][5] now rather firmly establish the existence of an exotic, doubly bottom, I = 0, J P = 1 + udbb tetraquark state, bound with respect to BB * , and hence strong-interaction stable. Though results for the binding energy vary, all correspond to masses below BB threshold, making the state stable with respect to, not only strong-interaction, but also electromagnetic, decays. The results of Refs. [3,4] also predict a strong-interaction-stable flavour3 F , I = 1/2, J P = 1 + strange partner. Refs. [6,7], using heavy quark symmetry arguments, supplemented by either phenomenological input [6] or phenomenological plus lattice input [7] for leading finite heavy mass corrections, concur on the strong-interaction stability of both the I = 0 and I = 1/2 doubly bottom states [8]. The binding in these channels appears to be driven by a combination of the attractive short-distance colour Coulomb interaction between twob antiquarks in a colour 3 c and the spin-dependent interaction of light quarks in the3 F , spin 0, colour3 c "good light diquark" configuration, whose attractive character is known phenomenologically from observed heavy baryon splittings [3,9]. Neither of these interactions is accessible in a state consisting of two well-separated heavy mesons.
While the existence of this3 F of J P = 1 + strong-interaction-stable tetraquark states seems well established theoretically, detecting these states, and confirming this prediction, remains experimentally challenging. Production rates are likely to be very low, given that two bb pairs must first be produced. While the doubly charmed Ξ cc baryon has now been observed by LHCb [10,11], and double-Υ production has been reported by CMS [12], to date, not even bottom-charm, let alone doubly bottom baryon states have been detected.
The recently proposed inclusive search strategy, in which a B c meson originating from a displaced vertex would signal the presence of doubly bottom hadron production at the LHC [13], though clearly of interest, would still leave open the question of whether such a signal contained a doubly bottom tetraquark component, or was due entirely to doubly bottom baryon production.
In this paper we use lattice QCD calculations to investigate whether strong-interactionstable (hereafter "bound") tetraquarks exist in other channels likely to be more amenable to experimental detection. A number of such channels have been investigated in the literature using various approaches, including QCD-inspired models  and QCD sum rules [43][44][45][46][47][48][49][50][51][52][53][54][55], and we will study a range of channels where claims for the existence of such bound tetraquark states have been made. Our investigations will also serve to test those aspects of the models not constrained by fits to the ordinary meson and baryon spectrum as well as approximations employed in implementing the QCD sum rule framework in studies of tetraquark channels.
In the case of model studies, the parameters of the models (most typically constituent quark models) are fixed from earlier fits to the ordinary meson and baryon spectrum.
Interactions of a constituent quark pair in a colour3 c or a constituent quark and antiquark in a colour 1 c configuration are, thus, phenomenologically constrained. In tetraquark (and other multi-quark) channels, however, additional colour configurations, which have totally unconstrained constituent interactions, are also present. In those channels already identified by the lattice as supporting deeply bound, doubly bottom tetraquark states, where the physics of the binding is believed to be dominated by a combination of the colour Coulomb attraction in the anti-diquark colour 3 c and the diquark attraction in the good-light-diquark 3 c configuration, these dominant interactions are phenomenologically constrained, and one would, therefore, expect the models to successfully predict tetraquark binding in these channels. This is, indeed, typically the case. However, even in these channels, differences in the assumed form of the model interactions, which, by construction, produce no differences in the ordinary hadron spectrum (since the model parameters have been fit to ensure the spectrum is reproduced) can produce significant differences in an exotic channel. An example is provided by the comparison of results for predicted binding in the I = 0, J P = 1 + udbb channel produced by a range of dynamical chiral [21,23,25,26,29,30,36,41] and non-chiral [16, 17, 19, 20, 22, 24, 26, 27, 29-32, 37, 40, 42] quark models. The latter typically employ a purely one-gluon-exchange (OGE) form for the spin-spin interaction, while spinspin interactions for the former are produced by a combination of effective Goldstone-boson exchange, acting between the light constituents only, and nominal OGE acting between all constituents. Even if the tetraquark state is entirely dominated by the 3 c anti-diquark,3 c good-light-diquark configuration, this configuration includes both 1 c and 8 cb -light quark pairs. The fits of the models to the ordinary meson and baryon spectrum place no constraints on these residual 8 c heavy-light interactions. In addition, the colour dependence of the OGE interaction, and the absence of Goldstone-boson-exchange contributions to the spindependent heavy-light interactions in the chiral quark model framework mean one must expect significant differences in the residual heavy-light interactions, and hence in the predictions for the tetraquark binding, in the two classes of models. This expectation is also borne out in the literature, where the dynamical non-chiral models of Refs. [16, 17, 19, 20, 22, 24, 27, 29-31, 37, 40, 42] produce binding energies between 54 and 160 MeV, in contrast to the results of the dynamical chiral models of Refs. [21,23,29,30,36,41], which lie between 214 and 497 MeV.
Different systematic questions arise for predictions generated using QCD sum rules, where the underlying dispersive representations are rigorously valid but practical applications require approximations on both the spectral and OPE sides. All tetraquark studies we are aware of work with Borel-transformed sum rules and the SVZ spectral ansatz, which consists of a single low-lying "pole" (characterized by its mass, M , and coupling, f , to the relevant interpolating current) and a continuum (approximated using the operator product expansion (OPE)) for s above a "continuum threshold" s 0 [56,57]. The OPE side is taken as input, and the sum rules used to fix M , f , and s 0 . The OPE, however, involves not just known quantities, like quark masses, α s and qq , but also unknown higher dimension condensates. Contributions up to dimensions D = 8 or 10 are retained in recent tetraquark analyses. As stressed in the earliest of the sum rule studies of doubly heavy tetraquarks [43], estimating such higher D contributions (typically in terms of products of known lower-dimension condensates) requires use of the factorization approximation. The accuracy of this approximation is not known in general, but for the I = 1 vector and axial-vector current-current two-point functions, extractions of the two relevant D = 6 four-quark condensates from finite-energy sum rule analyses of OPAL [58] and ALEPH [59] hadronic τ decay data found it to be off by as much as a factor of ∼ 6 [60,61]. On the spectral side of the sum rules, the single-narrow-resonanceplus-continuum SVZ ansatz may not be suitable for all channels. It would, for example, represent a poor choice for the I = 0, J P = 0 − channel, which has not one, but three narrow, low-lying η resonances. Its use in exotic channels, where limited prior knowledge of the qualitative features of the spectral distribution is likely, thus has the potential to produce difficult-to-assess systematic uncertainties. One way to investigate this issue is to consider sum rules for different interpolating operators with the same quantum numbers. The associated spectral functions will all receive contributions from all states with these quantum numbers, with only the contribution strengths operator dependent. Finding compatible pole masses and similar continuum thresholds from all the sum rules would increase confidence in the reliability of the SVZ ansatz. Finding incompatible pole masses would, in contrast, signal either that the approximations used on the OPE sides were unreliable, or that the channel has more than one narrow state, raising questions about the suitability of the use of the SVZ ansatze.
The lattice approach is ideally suited to investigating channels with a deeply bound tetraquark ground state, and to testing model and sum rule predictions for the existence of such states. Provided one employs interpolating operators with reasonable ground state overlap, the sizeable gap to the lowest meson-meson threshold increases the likelihood the ground state will dominate the corresponding two-point function at moderate Euclidean times, before the signal is lost in noise. Judicious source and sink choices (discussed in more detail in Sec. IV below) are also relevant to improving the ground state signal. Weakly bound states represent more of a challenge because finite volume (FV) effects have the potential to produce a FV analogue of what will become a meson-meson scattering state in the continuum that is shifted below the continuum meson-meson threshold. FV effects are expected to be small for bound states, but not necessarily negligible for continuum scattering states. They are typically handled either by extrapolations using simulations at multiple volumes, or by working with ensembles where the product m π L, with L the lattice length, is large enough to suppress the dominant "round-the-world" FV effects. A general rule of thumb is that this requires m π L > 4. The fact that the good-light-diquark configuration is believed to play an important role in binding for the tetraquarks so far identified and that, phenomenologically, the associated contributions to binding are expected to grow with decreasing light quark mass, also puts a premium on working with m π as close to the physical point as possible, subject to maintaining a detectable signal and keeping FV effects under control.
For the 32 3 × 64 ensembles used in our previous studies of the J P = 1 + ,3 F udbb and sbb ( = u, d) channels [3] and J P = 1 + , I = 0 udcb channel [9] the m π L values were 6.1 for m π 415 MeV and 4.4 for m π 299 MeV, but only 2.4 for m π 164 MeV. For the current study, which focuses on searching for channels with bound, non-molecular (i.e., non-weakly-bound) tetraquark ground states, we have thus generated a new ensemble with the same lattice spacing but a larger volume (48 3 × 64) and a pion mass, m π 192 MeV, still close to physical, but with m π L = 4.2 > 4. With only this single m π , we will be unable to identify channels in which a shallow bound state might exist at physical m π , but not at the slightly higher m π of our simulation. Channels identified as having a moderate-to-deeply bound tetraquark ground state for m π 192 MeV, will, however, certainly also have such a ground state at slightly lower, physical m π .
We now turn to the channels to be investigated in this paper. The likelihood of binding is increased by restricting our attention to those channels where no relative spatial excitations are required, and where light quark pairs have access to the favourable3 F , spin 0, colour 3 c , good-light-diquark configuration. The antiquark pair must then have access to the 3 c configuration and, with no relative spatial excitation, have spin 1 if the two antiquarks are the same. Both spin 0 and 1 are possible if the antiquarks are different. These considerations lead us to investigate channels with either I = 0 or 1/2, J P = 1 + and two identical antiquarks, or I = 0 or 1/2, J P = 0 + or 1 + and two non-identical antiquarks.
In what follows, we will first briefly revisit our previous study [3] of the doubly bottom, J P = 1 + ,3 F channels already well established to support bound tetraquark states. Our goal here is not to re-establish this fact, but rather to illustrate an important improvement to the methods employed in Ref. [3]. In that analysis (as well as in our analysis of the I = 0, binding is already clearly established, the ultimate goal is to carry out simulations at not just the larger volume and single m π of the current study, but also, at this same volume, for multiple m π , in order to provide updated results for the binding, extrapolated to physical m π . The results of this extended update, which is in progress, will be reported elsewhere.
The rest of this paper is organized as follows. In Sec. II, we list the additional channels to be studied and review existing predictions for tetraquark binding in each. Sec. III provides a list of the operators used in these studies and Sec. IV a discussion of our gauge-fixed wall sources, and the extended, Box-Sink construction, introduced to improve our ground state signals. Sec. V provides details of our lattice setup and calculations, and Sec. VI the results of these calculations. Our final conclusions, together with a discussion of the implications of our results, are given in Sec. VII.

PREDICTIONS THEREIN
In this section, we specify the additional channels to be considered in this study, providing a review of existing tetraquark binding predictions for, and the reasons for an interest in, each such choice. The predictions are to be tested against results from the lattice in Sec. VII.
The existence of sizeable B c data sets already establishes the existence of significant simultaneous production of bb and cc pairs. A bound tetraquark in one of the two mixed-heavy charm-bottom channels is thus expected to be much more amenable to experimental detection than either of the J P = 1 + ,3 F doubly bottom analogues discussed above. Replacing one of the twob antiquarks by ac, however, is expected to both reduce the Coulomb attraction of the 3 c heavy antiquark pair and increase the residual spin-dependent heavy-light interactions, which heavy baryon spectrum splittings suggest is likely to further reduce binding. The bottom-charm states (if any) are thus qualitatively expected to be more weakly bound than their doubly bottom counterparts [9].
The mixed-heavy I = 0, J P = 1 + udcb channel was investigated previously on the lattice [9]. Evidence of a possible bound ground state with binding between ∼ 15 and 61 MeV was found, though with rather short, late-time plateaus. The heavy quark symmetry arguments of Refs. [6,7], which produce binding compatible with that found from the lattice for the doubly bottom I = 0, J P = 1 + ground state, in contrast, predict no binding. A number of non-chiral models [17,19,20,34,37,42] also find an either unbound or only weakly bound ground state, with binding energies between −1 MeV (unbound) and 20 MeV (depending on the details of the potential used) and 23 MeV reported in Refs. [20] and [34], Non-chiral models [18-20, 33, 34, 37, 42] predict an either unbound or only weakly bound ground state, with Refs. [20,33] and [34] quoting binding energies of between −11 MeV (unbound) and 13 MeV (depending on the details of the potential used), 11 ± 13 MeV, and 23 MeV, respectively. Significantly larger binding is found in most recent chiral model studies, with Refs. [36,38,41], for example, reporting bindings of 136 ± 12 MeV, 196 MeV and 178 MeV, respectively. Two recent QCD sum rule studies differ in their pole mass results, with Ref. [45], which includes OPE contributions up to D = 8, finding 7.14 ± 0.10 GeV, while Ref. [48], which includes contributions up to D = 10, finds 6660 ± 150 MeV. The latter corresponds to a very deep, 485 ± 150 MeV, binding relative to DB threshold.
The J P = 0 + and 1 + scb channels have yet to be studied on the lattice. Analogous J P = 1 + , sb b states, with variableb mass as low as ∼ 0.6 m b , have, however, been considered, for m π = 299 MeV only, in Ref. [9]. With bothb andb treated using NRQCD, this study could not be extended tob masses as low as m c ; the observed variableb mass dependence, however, was argued to make a bound scb state extremely unlikely in the J P = 1 + channel at this pion mass. Both the heavy quark symmetry arguments of Refs. [6,7] and non-chiral model of Ref. [42] predict no tetraquark bound state in either this channel or the J P = 0 + analogue. In contrast, the QCD sum rule study of Ref. [54], which includes contributions up to D = 8, finds binding energies of 200 ± 130 and 180 ± 110 MeV for the J P = 0 + and 1 + ground states, respectively. Even deeper J P = 0 + bindings, exceeding 400 MeV, are reported in the QCD sum rule studies of Refs. [49] and [53], which include OPE contributions up to The reported observation by the D0 collaboration [62,63] of a narrow state, X(5568), decaying to B s π ± , and hence having four distinct flavours, prompted speculation that the related singly-heavy udsb and udsc channels might support bound tetraquark states. The initial argument was based on the expectation that the putative X(5568) should have a U -spin-interchanged SU (3) udsb partner with similar mass, coupled with the observation that the threshold, 5773 MeV, for the lowest-lying, two-meson udsb state, BK, lies well above the reported X(5568) mass. While this initial motivation is weakened by the fact that searches by the LHCb [64], CMS [65], CDF [66] and ATLAS [67] collaborations failed to confirm the D0 result, a number of model and QCD sum rule studies exist predicting bound tetraquark ground states in the I = 0, J P = 0 + and 1 + , udsb and/or udsc channels. The quark colour delocalization screening model, for example, predicts J P = 0 + and 1 + udsb ground states bound by 74 and 58 MeV relative to the respective two-meson, BK and B * K, thresholds [39].
Another SU (3) chiral quark model study, that of Ref. [35], obtains similar bindings, 70 and 68 MeV, for the J P = 0 + and 1 + bottom-strange states. Somewhat smaller bindings, 19 and 16 MeV, respectively, are found in the alternate chiral quark model study of Ref. [41]. A bound tetraquark with even lower mass, 5380 ± 170 MeV (corresponding to a binding of 394 ± 170 MeV relative to BK threshold), is also found for the J P = 0 + bottom-strange channel in the QCD sum rule study Ref. [51]. In fact, the only bottom-strange sector investigation we are aware of which does not produce a bound tetraquark state is the non-chiral-model, [17]. This study also found no bound tetraquark in the I = 0, J P = 1 + udsc channel. An absence of binding was also found for the related I = 0, J P = 0 + udsc channel in the sum rule study of Ref. [51]. To our knowledge, the only other udsc channel study is that of Ref. [41], which finds I = 0, J P = 0 + and 1 + states bound by 15 and 9 MeV, respectively. Any bound state found in these channels would be of considerable phenomenological interest, since it would involve light degrees of freedom in the novel strangeness +1, colour 3 c configuration about which no phenomenological information is currently known.
While it might appear natural to include the I = 0, J P = 1 + udcc and I = 1/2, J P = 1 + scc channels in our study, we have not done so, since evidence already exists that bound doubly charmed tetraquark states, even if they exist in these channels, will be at most weakly bound. This evidence comes from both the lattice and the heavy-quark-symmetry arguments of Refs. [6,7], whose results are compatible with those of lattice studies for the doubly bottom bound states. The latter predicts no binding in either of the doubly charmed J P = 1 + channels. On the lattice, in the I = 0, J P = 1 + udcc channel, (i) the results of Ref. [68] clearly establish the absence of binding at m π = 391 MeV, (ii) Ref. [4] (which reaches a lower, 257 MeV, pion mass for the coarsest of the three lattices studied) finds an extrapolated, physical-m π binding of 23±11 MeV, sufficiently small for the authors to comment that further FV studies are required to reach a firm conclusion concerning binding, and (iii) results for the b mass dependence of the binding energies for I = 0, J P = 1 + udb b states at m π = 299 MeV, reported in Ref. [9], strongly disfavour the possibility that binding survives to physical m π and b masses as low as m c . Only one lattice result, whose binding, 8 ± 8 MeV, is compatible with a two-meson-threshold ground state, exists for the I = 1/2, J P = 1 + scc channel [4]. QCD sum rule studies also find no evidence for binding in either of the non-strange [43,44,46,55], or strange [44,46,55] doubly charmed J P = 1 + channels. Non-chiral models whose I = 0, J P = 1 + udbb bindings are compatible with those found on the lattice, similarly, predict either unbound [17,22,27,29,33,36,37], or only weakly bound [24,26,30], doubly charmed J P = 1 + ground states. While many chiral quark model studies predict significant binding in the I = 0, J P = 1 + udcc channel [21,23,26,29,30,38,41], the models underlying these predictions typically also over-bind the I = 0, J P = 1 + , udbb ground state relative to results known from the lattice. We conclude that the absence of deeply bound tetraquark states in the two doubly charmed channels is already established, and hence have not included these channels in the current study. Additional, larger volume ensembles would be required to reliably investigate such channels, where at-best-weak binding makes FV studies mandatory.
In future, we plan to generate additional ensembles with similarly low m π , but larger volume than considered here, and will revisit the doubly charmed J P = 1 + channels when these become available.
Two final channels are studied in this paper. These are the triply-heavy J P = 1 + ucbb and scbb channels, considered previously only in the lattice study of Ref. [4], which found ground state masses compatible with the corresponding lowest two-meson thresholds. These channels are, of course, not among those where a bound tetraquark state, if it existed, would be more amenable to current experimental detection than the theoretically established doubly bottom states. They are included here only to provide an additional, independent check of the results of Ref. [4].
It is also worth emphasizing that chiral and non-chiral quark models, all of which admit parametrizations which successfully reproduce the ordinary meson and baryon spectra, make qualitatively different predictions for the tetraquark channels discussed above. The non-chiral models, once one moves beyond the doubly bottom sector, predict only a few, usually atmost-weakly-bound tetraquark candidates. The majority of chiral model studies, in contrast, predict a larger number of much more deeply bound states. The lattice studies below should thus allow us to rule out at least one of these classes of models as providing an acceptable phenomenological representation of QCD in the low-energy regime.

III. OPERATORS
In this paper we work with the following set of operators having couplings to states with two quarks and two antiquarks with quark flavours ψ, φ, θ, and ω, and having either "meson-meson" or "diquark-antidiquark" spin-colour-flavour structure: The Dirac structures Γ 1 and Γ 2 relevant to the channels we consider are specified below.
The short-hand terminologies "meson-meson" and "diquark-antidiquark" are meant only to compactly characterize the spin-colour-flavour structures, and should not be over-interpreted; there is nothing, for example, to prevent one of the "meson-meson" operators from coupling to a state with a spatially extended diquark-antidiquark structure or one of the "diquarkantidiquark" operators from coupling to a meson-meson scattering state. An illustration is provided by the deeply bound ground state in the doubly bottom, I = 0, J P = 1 + channel, which was found to couple strongly to both the "meson-meson" and "diquark-antidiquark" operators considered in Ref. [3].
The lattice explicitly breaks rotational symmetry, breaking the continuum rotation group down to the little group O h . The operators listed above lie in either the A 1 or T 1 irrep, which map directly to the continuum J = 0 and J = 1 quantum numbers respectively. For the operators that are produced from the product  out the scalar and vector components respectively [5].
The operator combinations which in principle couple to states in the J P = 0 + and 1 + channels considered below are listed in Table I. While the diquark-antidiquark operators, D and E, were included in initial explorations of the udsc and udsb channels, they were found to couple only weakly to the corresponding ground states and hence not included in the final analyses of these channels. They were also omitted from the final version of the analyses of the udcb channels where using the set of meson-meson operators was found to give more statistically precise results and produce no significant change to the energies of the lowest-lying states.
Throughout this work we will obtain our ground states from a generalised eigenvalue problem (GEVP) approach [70][71][72]. We use the solutions of the GEVP (sample by sample) to extract the "optimised correlator" where V is the matrix with columns made of the eigenvector solutions of τ (the "diagonalization time" [73]) in Eq. (2) is to be chosen large enough to produce improved projection onto the ground state. We then perform a correlated single-exponential fit to the resulting optimised correlators to obtain our final levels. In the rest of this work we will quote results obtained using t 0 /a = 2 and τ /a = 4, values for which the ground state masses are found to display good stability upon variation of these parameters. Larger values of t 0 and τ tend to make the solution statistically less precise, and eventually unstable.

IV. WALL SOURCES AND THE BOX-SINK CONSTRUCTION
Throughout this work we will use Coulomb-gauge-fixed wall sources [74,75] (fixed to high-precision using the FACG algorithm of [76]). These sources have some benefits, as well as some peculiarities, which are discussed below.
A gauge-fixed wall source is a time-slice of point sources on a Coulomb gauge-fixed background. Here we assume the gauge condition has been applied over the whole volume.
The elements of the source do not need to be stochastically drawn as the gauge fixing condition connects colors appropriately; this means that these sources can be used for baryons [77][78][79].
Gauge-fixed wall source results can be considerably more statistically precise, especially at lighter pion masses, than point sources at the same cost due to the volume-sampling. Like other wall sources, these automatically zero-momentum project, and if momentum is needed it must be put into the source explicitly as a momentum source or introduced via partial twisting [80,81], as was done for the NRQCD tuning (see App.A).
Use of a guage-fixed wall source is much like that of a point source, though with some potential complications for channels like those studied here where two-meson thresholds exist.
When combined with a local sink, one typically finds a negative sign for the amplitude of (at least) the first excited state, and hence a "Wall-Local" correlator whose effective mass approaches the ground state mass from below. If one cannot measure the signal to large enough t to ensure the ground state plateau has been reached, the ground state mass may then be under-estimated, and lead one to either conclude that a bound state exists when in fact one does not or over-estimate the binding in the case one does exist. This issue does not arise for point or stochastic wall sources, where the effective masses approach the ground-state mass from above and so offer a rigorous upper bound on that mass.
Often it useful to compute "Wall-Wall" correlators [82], defined by contracting objects in the usual way, with propagators that are summed over the spatial sites of the sink time-slice, The Another would be to simply take advantage of the gauge condition and construct propagators by summing over points lying inside a sphere around each reference sink point x, With this construction, varying R 2 between 0 and its maximum value, 3(L/2) 2 , produces a continuous interpolation between the Wall-Local and Wall-Wall cases. We refer to sinks   Values were chosen so the corresponding effective mass plateaus begin reasonably early.
One should bear in mind that, in contrast to what happens in the Wall-Wall case, Wall-Box and Wall-Local correlation matrices are not symmetric [83]. The GEVP method, however, doesn't necessarily need this symmetry to obtain the underlying real eigenvalues. Careful monitoring of the imaginary parts of the eigenvalues is needed and in practice the lack of symmetry was not found to be an issue.
An illustration of the utility of the Box-Sink construction is provided, for the case of the B c meson channel, in Fig. 1. We see that intermediate values of R 2 , indeed, exist which extend the ground state plateau to significantly lower t than is the case for the Wall-Local or Wall-Wall combinations. Results for all Box-Sink radii must, of course, approach the same ground state mass at sufficiently large t. The point here is that the improvement in the ground state plateau occurs for small enough R 2 that the onset of the significant increase in noise that occurs as R 2 grows towards the Wall-Wall limit is avoided in the optimal R 2 plateau region. The figure also illustrates that a simple Wall-Local correlator can have significant excited state contamination and only approach the plateau at very large times.
This is less of a problem in single-meson channels, like the B c channel, where one can follow the signal out to large enough t, but becomes an issue in channels with larger signal-to-noise problems, such as the tetraquark channels we study here, where we are typically only able to follow our signal out to t/a ≈ 20 due to an exponential signal to noise problem. Such issues have been discussed in the literature for a long time in the context of the nucleon and ρ meson, see, e.g., Refs. [77,78,84].
A further illustration of the utility of the Box-Sink construction, now for tetraquark signals, is provided by the effective mass plots of the optimised ground-state correlators for the J P = 1 + , I = 0 udbb and I = 1/2 sbb channels at various Box-Sink radii R 2 , shown in Fig. 2. While the errors grow, as anticipated, as R 2 increases, the ground state signals plateau much earlier, and in the region of much lower noise, than was the case for the corresponding Wall-Local plateaus shown in Ref. [3]. This is especially true of the optimised R 2 region, R 2 36 or 49. While the signals degrade quickly for larger t/a (and larger R 2 ), it is clear that a "window of consistency" exists in the region of the Box-Sink plateaus, for both channels. It is also clear that, while the Wall-Local (R 2 = 0) results reach the same plateaus, they do so only at much later t/a, rather close to the region where the signals are lost.
The fact that enhanced excited-state contamination is seen in the Wall-Local results of both the B c and tetraquark channels suggests this behavior may be common for Wall-Local correlators. Fig. 2 illustrates, in more detail, the potential danger this creates. In the absence of the longer, more reliable Box-Sink plateaus, it would be easy, with only the Wall-Local results, to mistakenly identify the start of the Wall-Local plateaus as occurring around t/a = 14, leading to an overestimate of the binding in both channels. The utility of the Box-Sink construction in avoiding this problem is clear. Given the results shown in Fig. 2, we expect our earlier Wall-Local-correlator-based estimates for the physical point bindings in these channels [3] to be reduced once the ongoing improved, multiple-m π Box-Sink analysis is completed. For now, the results for these channels, at the single m π of the current study, serve to establish the Box-Sink construction as an improved method for identifying channels supporting deeply-bound tetraquark states and quantifying the binding therein, and motivate our use of it in the studies of the other channels, reported below.

V. LATTICE DETAILS
For all channels investigated in this work we have used a single n f = 2 + 1, 48 3 × 64 Wilson-Clover ensemble generated as an extension to the PACS-CS ensembles [85] used in our previous works [3,9]. As this ensemble is generated with the same β and C SW , the lattice-spacing (for which we use a −1 = 2.194(10) GeV [86]) and the charm-quark action parameters should be the same as for those previous ensembles.
This ensemble was generated using openQCD-1. 6 [87] and all propagator inversions were performed using a version of this code (modified [88] for the charm quarks) as well. For the light and strange quarks we use the DFL+SAP+GCR solver algorithm [89][90][91], and solve for multiple wall source positions as the same deflation subspace can be used. For the charm quarks, we just use the SAP+GCR solver as low-mode deflation is practically ineffective for very heavy quarks and there is a significant overhead in generating the deflation subspace compared to the solve time. For the charm-quark action we use the relativistic heavy quark action of Refs. [92,93] with the same parameters as listed in [94]. The new ensemble has some attractive features. In particular, its pion mass, ≈ 192 MeV , is reasonably close to physical while still maintaining m π L > 4. As for the previous PACS-CS ensembles, we have a strange sea quark heavier than physical and thus have to use a partiallyquenched strange (using the κ s -value of [95]) to get a respectable (507 MeV) kaon mass. The Tsukuba tuning also gives close to physical D and D * meson masses. The lattice-valued masses we obtain are presented in Table VI in App. C. Note that, for NRQCD, where there is an additive mass renormalisation, ∆, only mass differences have physical meaning. We will present results in this lattice-valued form for the majority of this work.
As we use NRQCD for the b quarks, channels with one or more b quarks do not have a rigorous continuum limit. We do, however, include terms that should cancel up to O(a 4 ) discretisation effects in our NRQCD Hamiltonian [96]. The c quarks have been tuned to have a flat dispersion relation and so are expected to have a very mild approach to the continuum.
The light-quarks are O(a)-improved non-perturbatively with the C SW term. In all, we expect discretisation effects to be relatively small.  In this section, we present the measured spectra (GEVP eigenvalues) and associated two-meson thresholds for the channels discussed above. The eigenvalues are obtained from The results for the J P = 0 + and 1 + scb channels are shown in Fig. 3b. Here we see the ground states in both channels pushed marginally below the corresponding lowest two-meson thresholds, DB s and DB * s , respectively. With no other states in the vicinity, however, we suspect this reflects either a residual systematic (from either the GEVP procedure or our fitting to the optimised ground-state correlators) or a small (< 10 MeV) FV shift, and that the ground states in these channels correspond to the lowest two-meson thresholds. With only a single volume, the possibility of weakly bound ground states cannot, however, be unambiguously ruled out. Interestingly, the three highest excited states in the J P = 1 + channel all lie just below higher-lying two-meson thresholds, with splittings consistent with those between these thresholds. This "pushing down" of the eigenvalues is of the order of a few MeV at most. This pattern is not, however, observed for the J P = 0 + channel, where the third excited state lies much higher than the third excited two-meson threshold, D * B * s . once more correspond to the two-meson threshold states. Results for the analogous I = 0, J P = 0 + and 1 + udsc channels are shown in Fig. 4b. The ground states lie slightly above the lowest two-meson thresholds, with, once more, no sign of an additional, bound state below threshold in either channel.
Finally, in Fig. 5 we show our results for the triply-heavy J P = 1 + , ucbb and scbb channels.
The ground states are again consistent with the corresponding two-meson thresholds (BB * c and B s B * c , respectively), confirming the earlier results for these channels reported in Ref. [4].

VII. DISCUSSION AND CONCLUSIONS
In this work we have used lattice simulations to investigate a number of exotic tetraquark channels in which model and/or QCD sum rule studies and/or heavy-quark symmetry arguments predict bound, strong-interaction-stable tetraquark states to exist. We have also revisited our earlier lattice studies of the doubly bottom, J P = 1 + ,3 F udbb and sbb and bottom-charm, J P = 1 + , I = 0 udcb channels [3,9], introducing a number of significant improvements over those earlier works. In particular, we have (i) generated a new ensemble with close-to-physical ( 192 MeV) pion mass and significantly larger volume (m π L = 4.2 > 4) to better control possible finite volume effects; (ii) added additional operators to improve our ability to resolve the tower of states contributing to each correlator, and (iii) made major improvements to our ground-state signals via the Box-Sink construction.
We have also, as described in more detail in App. A, re-tuned our NRQCD bare mass to a much more precise and reliable value, using partially-twisted boundary conditions. All of these represent dramatic improvements over our previous work.
The only channels in which we find evidence for deeply bound, strong-interaction-stable tetraquark states are the doubly bottom J P = 1 + ,3 F , udbb and sbb channels already identified as supporting such bound states in earlier lattice studies, including our own. Our new results for these channels, albeit so far only at a single m π 192 MeV, display significantly improved ground-state plateaus, largely due to the improved, Box-Sink construction. With the bindings at this m π reduced compared to those implied for the same m π by the fits to the m π dependences reported in Ref. [3], we expect to find similarly reduced physical-m π results once we complete our ongoing multiple-m π analysis of these channels and are able to perform the physical-point extrapolation. With the binding in these channels expected to increase with decreasing m π , the final physical-point bindings will be deeper than those of the current study. The current, single-m π results, thus provide further confirmation of the conclusions reached in previous lattice studies regarding the existence of a3 F of strong-interaction-stable, J P = 1 + , doubly bottom tetraquark states. Once the multiple-m π analysis has been completed, the improved ground state signals produced by the Box-Sink construction should also significantly improve the reliability of the extrapolated physical-m π results.
In contrast to the doubly bottom,3 F , J P = 1 + channels, no bound, non-molecular ground states are found in any of the other ten channels considered here. States with bindings > ∼ 10 − 20 MeV, as predicted by many model and QCD sum rule studies, should have been easily detectable in these analyses.
Among the channels where current results rule out the possibility of a bound, nonmolecular ground state is the I = 0, J P = 1 + udcb channel. This conclusion stands in contrast to that reached in our earlier study of this channel [9], where indications for possible modest tetraquark binding were found, albeit based on rather short, late-time plateau signals, and with worries about possible FV effects for the ensemble with m π = 164 MeV, where m π L = 2.4. The desirability of studying this channel at larger volumes, in order to more strongly test the bound state interpretation of the ground state results, was stressed already in Ref. [9], and the current study provides precisely such a larger-volume extension, one which, moreover, has the advantage of the significantly improved Box-Sink analysis methodology.
The indications of binding seen in the earlier study do not survive the improved, largervolume, Box-Sink analysis: the ground-state mass is found to lie, clearly slightly above, rather than below, DB * threshold. The very good associated ground-state effective-mass plateau is shown in the right-hand panel of Fig. 8 in App. D. We conclude that there is, in fact, no non-molecular bound tetraquark state in this channel, and that the late-time behavior of the earlier Wall-Local results was likely affected by the late-Wall-Local plateauing problem discussed in detail for the B c and doubly bottom3 F , J P = 1 + channels in Sec. IV above.
Turning to a comparison to the results of other approaches, we note first that the predictions of the heavy-quark-symmetry approach [6,7] for the binding energies in the I = 0, J P = 1 + udbb and J P = 1 + sbb channels, as well as for an absence of binding in the J P = 0 + and 1 + , I = 0 udcb and J P = 0 + and 1 + scb channels, are all compatible with current lattice results.
We turn next to tests of QCD sum rule and model predictions. quarks [23,25,29,31,35,36,38], or "SU(2)" form, where it is restricted further to u and d quarks only [30,35,41]. A final ChQM implementation, used to study the I = 0, J P = 0 + and 1 + udsb channels [39], is the quark delocalization color screening model (QDCSM), which retains pseudoscalar exchange between light quarks and OGE between all quarks, but includes no scalar meson exchange, instead employing a modified confinement form that allows two-center cluster configurations with quark delocalization between the clusters. The scalar-nonet implementation is ruled out by its 32 MeV prediction for the I = 0, J P = 1 + udbb binding, which lies well below the range, ∼ 100 − 180 MeV, allowed by current lattice results. The SU (2) implementation is also ruled out, this time by significant over-bindings relative to lattice results: 318 [41] or 404 MeV [30] for the I = 0, J P = 1 + udbb channel and 178 and 199 MeV [41] for the I = 0, J P = 0 + and 1 + udcb channels (which our results show to be at most weakly unbound). The SU (3) implementation is similarly ruled out, with predicted I = 0, J P = 1+ udbb bindings of 341 [23], 214 [29], 217 [31], 278 [36] or 359 [38] MeV, I = 0, J P = 0 + and 1 + udcb bindings of 136 ± 12 and 171 ± 12 MeV [36] or 178 and 217 MeV [38], and I = 0, J P = 0 + and 1 + udsb bindings of 70 and 68 MeV [35] (our lattice results again show the ground states to be at most weakly bound in these channels). Finally, the QCDSM implementation predicts bindings of 74 and 58 MeV for the I = 0, J P = 0 + and 1 + udsb channels, and hence is also ruled out by the absence of such deep binding in our lattice results. We conclude that all implementations of the ChQM approach are ruled out by current lattice results, and hence that the ChQM framework does not provide a reliable phenomenological representation of QCD in the low-energy regime.
The results of non-chiral models fare much better when tested against lattice results.
Multiple non-relativistic quark model studies of tetraquark channels exist using the AL1 [20] and Bhaduri (or BCN) [97] potentials. These potentials are characterized by pairwise linear confinement, colour Coulomb and regularized OGE hyperfine interactions. The AL1 model is a variation of the BCN model in which the strict OGE relation between the strengths of the colour Coulomb and hyperfine interactions has been relaxed. Ref. [20] also investigates three additional models (the AL2, AP1 and AP2 models), characterized by alterations to the radial dependence of the AL1 confinement potential and/or the regularization of the AL1 OGE potentials. Predictions for binding in the I = 0, J P = 1 + udbb; J P = 1 + sbb; I = 0, J P = 0 + and 1 + , udcb; and I = 0, J P = 1 + udcc channels exist for all these models, as well as for the Godfrey-Isgur-Capstick (GIC) model [98,99], a "relativized" variant of this class of model. BCN results may be found in Refs. [17,19,20,24,26,29,30], AL1 results in Refs. [20,24,34,40], AL2, AP1 and AP2 results in Ref. [20], and GIC results in Ref. [42].
Predictions for binding in the I = 0, J P = 1 + udbb; J P = 1 + sbb; I = 0, J P = 0 + udcb; and I = 0, J P = 1 + udcb channels, obtained from the non-relativistic BCN, AL1, AL2, AP1 and AP2 models, range from 131 to 167 MeV [19,20,24,26,29], 40 to 61 MeV [20], 0 to 23 MeV [20,34] and 0 to 23 MeV [20,34], respectively. Taking the ∼ 20 − 35 MeV variations with model choice for these channels seen in Ref. [20] as a measure of the uncertainties associated with the specifics of the implementations of the interactions in this class of model, one sees that the predictions are compatible within errors with current lattice results. Future lattice studies, using larger-volume ensembles to better quantify possible residual FV effects, will be useful for sharpening these tests for the J P = 0 + and 1 + , I = 0 udcb channels, where the results of this paper show no evidence for binding while the most recent AL1 model study [34] finds modest 23 MeV bindings for both. The relativized GIC model, in contrast, predicts only 54 MeV binding in the I = 0, J P = 1 + udbb channel and no binding in the J P = 1 + sbb channel [42], and hence is ruled out by lattice results.
Another model approach ruled out by lattice results, specifically those for the I = 0 udcb channels, is that of the LAMP model. In this approach, an ansatz for a relativistic confinement potential is employed in the multiquark sector. Taking the BD interaction studied in Ref. [32] to be specific, linearly rising mean-field confinement potentials for the B and D are merged into a one-body, two-well, W-shaped confinement potential for the BD system by truncating the individual B and D wells at the midpoint between their centers. Such a potential allows for the delocalization of an initially single-well light-quark wavefunction into the second of the two potential wells. The optimal well-center separation and degree of delocalization are determined by a variational minimization. The result, in which the effect of the OGE hyperfine interaction has been neglected, is a BD binding of 155 MeV, clearly incompatible with the lattice results above. This rules out the LAMP ansatz for the effective confining interaction.
Turning to tests of QCD sum rule results, we see that recent predictions for very deep tetraquark binding (by 569 ± 260 MeV [48] in the I = 0, J P = 1 + udbb channel, 477 ± 250 MeV [53] in the J P = 1 + sbb channel, 485 ± 150 MeV [48] in the I = 0, J P = 0 + udcb channel, 407 ± 160 MeV [49], 467 ± 150 MeV [53] and 200 ± 130 MeV [54] in the J P = 0 + scb channel, 180 ± 110 MeV [54] in the J P = 1 + scb channel, and 394 ± 170 MeV [52] in the I = 0, J P = 0 + udsb channel) are not compatible, within quoted uncertainties, with current lattice results. Other sum sum rule predictions (400 ± 300 MeV [43,44], 80 ± 80 MeV [46] and 240 ± 150 MeV [55] in the I = 0, J P = 1 + udbb channel, 290 ± 300 MeV [44], 140 ± 80 MeV [46] and 180 ± 160 MeV [55] in the J P = 1 + sbb channel, 5 ± 100 and 60 ± 100 MeV [45] in the I = 0, J P = 0 + and J P = 1 + udcb channels), are compatible with lattice results within errors, though, with the exception of the I = 0, J P = 1 + udbb prediction of Ref. [ MeV [48], for binding in the I = 0, J P = 1 + udbb channel obtained from analyses employing the same interpolating operator, essentially identical OPE input, and both truncating the OPE at the same dimension, D = 10. The predictions, 142 ± 80 MeV [46] and 477 ± 250 MeV [53], for binding in the J P = 1 + sbb channel also come from analyses using the same interpolating operator, essentially identical OPE input and truncation of the OPE at the same dimension, D = 10. A similar discrepancy exists in the I = 0, J P = 0 + udcb channel where very different ground state masses, 7140 ± 100 MeV [45] and 6660 ± 150 MeV [48], are predicted by analyses using the same interpolating operator and similar OPE input, differing only in the dimension at which the OPE was truncated (D = 8 in Ref. [45] and D = 10 in Ref. [48]). These observation make clear that significant, yet-to-be-fully-quantified systematic uncertainties exist for applications of the conventional implementation of the QCD sum rule framework to the multiquark sector. The targets provided by current (and future) lattice results for tetraquark binding (or lack thereof) in the channels investigated in this paper should help in further investigations of these issues, and have the potential to aid in the development of improved sum rule implementation strategies.
We close with a reminder of improvements/extensions to the present study either already in progress or planned for the near future. The first is the completion of the extension of the m π 192 MeV Box-Sink analysis reported here to additional 48 3 × 64 ensembles having both lighter and heavier m π . This will allow us to update and sharpen our previous physical-point extrapolations for the bindings of the doubly bottom, J P = 1 + ,3 F states, and further test non-chiral model predictions. This work is ongoing. The next step after that is to generate a new set of ensembles with even larger volume covering a similar range of m π . This will allow us to carry out quantitative FV studies, and extend the current analysis to channels, such as the doubly charmed channels, where shallow binding remains a possibility and experimental detection would be less challenging than in the doubly bottom sector.
A tilde indicates that an improved version has been used (see e.g. [102] for the tadpoleimproved derivative terms) and M 0 is a bare parameter that we tune to recover the spinaveraged kinetic mass of the connected η b and Υ. We use the fourth root of the plaquette for the tadpole improvement and a value of n = 4 for the stability parameter.
For the tuning we use Coulomb gauge-fixed wall-sources and we implement the momenta by solving one of the b propagators on a partially-twisted boundary [80,81]. This is practically implemented by applying a U (1) phase to the gauge field [103], which yields p µ = 2π Lµ θ µ . The benefit of this approach is that one can stay within the same irrep for arbitrary momenta, smoothly interpolating discretisation and finite-volume effects.
We choose to twist equally in all directions to minimise the rotation-breaking hypercubic artifacts and the effects introduced by the twisting. The momenta we use is listed in Table V and the effective masses for these twists is shown in Fig. 6.
where aM 1 is the mass of the ground state and aM K is the kinetic mass. This enables us to describe our data simply by five fit parameters.
We note that, as in Ref. [104], the kinetic mass ordering is inverted compared to nature.
This situation can be improved by the inclusion of yet-higher-order (1/M 3 ) terms in the NRQCD Hamiltonian [105], specifically the parameter c 7 . However, we follow the approach of [104] and tune to the spin average of these kinetic masses instead. Using HPQCD's estimate [104] for the spin-averaged mass, 9.445(2) GeV, and the value a −1 = 2.194 GeV [86] for our lattice spacing, we find the physical b quark corresponds to a bare mass aM 0 = 1.864. (Right) Stability and error reduction as the number of momenta is increased. This tuning is obtained by measuring the kinetic masses for several bare masses and fitting their dependence linearly, as is shown in Fig. 7.
When we compute our heavy-light mesons we pack the backward and forward propagating NRQCD propagators into one as the NRQCD propagator has only 2 Dirac indices whereas the full light, strange, and charm propagators have four, We also apply anti-periodic boundary conditions in time to the NRQCD propagators.

Appendix B: Meson masses
As we use Coulomb gauge-fixed wall sources our amplitudes are partially related to physical matrix elements, but are polluted by the wall source [77]. If we consider for example the pseudoscalar J P = 0 − meson channel, we can use any of the following 8 correlation functions at large asymptotic times, Here we use P to denote the operator with a γ 5 -insertion and A to denote the temporal axial current operator. Using these sources means, for example, that C W L P A is not equal to C W L AP . One can, however, clearly fit all 8 correlators simultaneously with 5 fit parameters. One can similarly analyse correlators involving the V i = γ i and T it = γ i γ t currents for J P = 1 − states or I and γ t for J P = 0 + states. The same fit form can be used when the sink is non-local too, although neither matrix element will correspond to a physical one. We use these fits to determine the masses of mesons not containing a b quark. This fit is not possible for mesons containing b quarks as C AA = C P P at our level of heavy-light approximation.
In this work we use simple, local operators of the form to determine the masses of ordinary meson states. The results, in lattice units, and with the additive offset in the case of mesons containing a b quark, are given in Tab. VI. Where there are several Γ insertions listed, a simultaneous fit of the form of Eq. B1 was used, otherwise a correlated single-cosh fit was performed on the one channel. As some of the particles listed in the table are resonances, the ground-state masses obtained from the associated local-operator correlators should not be directly compared to the masses of these resonances in nature. To simplify notation, we have, nonetheless, used the name of the lowest-lying resonance in the channel to label the channel itself in these cases, trusting this will not lead to any confusion.
Throughout this work we use the shorthand name B I for the scalar, J P = 0 + B meson.       If we consider our correlators C(t) to behave like single exponentials at large t (a reasonable expectation for heavy states) we can define an effective mass via am eff (t) = tanh −1 C(t − 1) − C(t + 1) C(t − 1) + C(t + 1) .
This has the desirable property that tanh −1 (x) is defined for all x in [−1, 1], and provides a definition that is, empirically, better behaved for very noisy data than the more common log or cosh −1 definitions, log C(t) C(t+1) or cosh −1 C(t+1)+C(t−1)

2C(t)
, where fluctuations can lead to negative values of the ratios appearing in the arguments of the log or cosh −1 , causing these definitions to break down.