The Search for Beauty-fully Bound Tetraquarks Using Lattice Non-Relativistic QCD

Motivated by multiple phenomenological considerations, we perform the first search for the existence of a $\bar{b}\bar{b}bb$ tetraquark bound state with a mass below the lowest non-interacting bottomonium-pair threshold using the first-principles lattice non-relativistic QCD methodology. We use a full $S$-wave colour/spin basis for the $\bar{b}\bar{b}bb$ operators in the three $0^{++}$, $1^{+-}$ and $2^{++}$ channels. We employ four gluon field ensembles at multiple lattice spacing values ranging from $a = 0.06 - 0.12$ fm, all of which include $u$, $d$, $s$ and $c$ quarks in the sea, and one ensemble which has physical light-quark masses. Additionally, we perform novel exploratory work with the objective of highlighting any signal of a near threshold tetraquark, if it existed, by adding an auxiliary potential into the QCD interactions. With our results we find no evidence of a QCD bound tetraquark below the lowest non-interacting thresholds in the channels studied.


I. INTRODUCTION
Tetraquarks were first considered theoretically decades ago in the context of light-quark physics in order to explain, amongst other experimental features, the a 0 (980) and f 0 (980) broad resonances [1]. 1 More recently, there has been exciting experimental evidence indicating the potential existence of tetraquark candidates amongst the so-called XYZ states -states whose behaviour differs from predictions of the heavy quark-antiquark potential model. The observed XYZ states apparently contain two heavy quarks, (cc) or (bb), and two light quarks [4]. The dynamics of these systems involves both the short distance and long distance behaviour of QCD and hence theoretical predictions are difficult. Consequently, many competing phenomenological models currently exist for these states [5]. Lattice QCD studies of the observed XYZ states are also difficult because these states are high up in the spectrum as well as being in the threshold region for strong decays into two heavy flavour mesons. While there are theoretical arguments that some tetraquark states with doubly heavy flavor (e.g., bbūd, bbūs and bbds) should be bound and stable against all strong decays [6], no general arguments exist for tetraquarks with heavy quark-antiquark content such as QQ qq states.
A tetraquark system composed of four heavy quarks is a much cleaner system to study theoretically as longdistance effects from light-quarks are expected not to be appreciable, as opposed to systems which are a mixture of heavy and light quarks. In the limit of very heavy quarks a chughes@fnal.gov b eichten@fnal.gov c christine.davies@glasgow.ac.uk 1 Recent lattice studies of the scattering amplitude pole do indicate that these states are in fact resonances as opposed to cusp effects, etc. [2,3] perturbative QCD single-gluon exchange will dominate [7] and so the dynamics are relatively simple. This makes these systems particularly useful to study in order to shed light on the aforementioned XYZ states. In fact, there is a multitude of phenomenological models (with a quark mass ranging from the bottom to the very heavy limit) which predict the existence of aQQQQ bound tetraquark [8][9][10][11][12][13][14][15]. However, these are not calculations from firstprinciples and have an unquantifiable systematic error associated with the choice of four-body potential. In reality, the heaviest possible tetraquark system in nature would be abbbb tetraquark. For this, non-perturbative QCD cannot be ignored, making a first-principles lattice QCD study essential. If such a boundbbbb tetraquark did exist, how it would be observed at the LHC has already been addressed [16,17]. Given these pressing theoretical motivations, in this work we perform the first lattice QCD study of thē bbbb system. The sole objective of this exploratory work is to determine if the dynamics of QCD generates enough binding force between thebbbb to produce a tetraquark state with a mass below the lowest noninteracting bottomonium-pair threshold, ensuring it is stable against simple strong decays. Searching for such a boundbbbb tetraquark candidate is particularly wellsuited to the first-principles lattice QCD methodology because this state, if it existed, would be the groundstate in thebbbb system. This means it should be relatively easy to identify. Further,bbbb annihilation effects are strongly suppressed by the heavy quark mass, as in the bottomonium system, and so can be ignored. This paper is organised as follows: in Section II the interpolating operators used in this study are discussed, in Section III the computational methodology is given, while the majority of the results are presented in Section IV. In Section V we explore a novel method of adding an auxiliary potential into QCD with the objective of highlighting a possible tetraquark signal. We then discuss our conclusions in Section VI. I. The colour representations of the different quark combinations. Note that, as described in the text, once the colour representation of the (anti-) diquark is chosen, the Pauli-exclusion principle enforces certain spin combinations in S-wave. Also given are the SU (3) colour contractions needed for thebbbb operators.

II. INFINITE-VOLUME CONTINUUM EIGENSTATES, OPERATORS AND TWO-POINT CORRELATORS
The QCD Fock space contains all colour-singlet singleparticle states such as the conventional mesons |η b (k) , |Υ(k) etc. and, if abbbb bound state also exists, a tetraquark state |T 4b (k) . In addition, there are also the two-particle states which can be labelled by appropriate quantum numbers as |Ptot, J P C ; |krel|, J P1C1 1 , J P2C2 2 , Lrel where Ptot (J P C ) is the total (angular) momentum of the two-particle system, with J PiCi i the quantum numbers of the individual particles and krel (Lrel) the relative (orbital angular) momentum between the two particles.
The sole motivation of this work is to search for a possiblebbbb tetraquark candidate within QCD that couples to a bottomonium-pair and which lies below the lowest threshold. The bottomonium mesons we study can be classified as J P C = 2S+1 L J . As any orbital angular momenta is expected to raise the internal kinetic energy of the state (and hence its rest mass) we focus on two-body S-wave systems (L = 0) with no orbital angular momentum between them (Lrel = 0).
With the quantum numbers of the Υ/η b being J P C = 1 −− /0 −+ , the S-wave 2η b and 2Υ can have (through the addition of angular momenta) a quantum number of 0 ++ , while the Υη b has 1 +− and the 2Υ can also be in a 2 ++ configuration. We now want to construct a full basis of S-wave colour/spin interpolating operators that has overlap with these quantum numbers. To do so, we start by forming all possible colour combinations that the 2b and 2b can be in. These are specified in Table I.
We can construct meson interpolating operators as where Γ M = iγ 5 , γ k projects onto the quantum numbers of the η b and Υ respectively, and G 1 (8) ef g is the colour projection onto the singlet (octet). In addition, it is also possible to construct a (anti-) diquark operator as O3 (6) D (t, x) = G3 (6) where (bĈ) α = C αβbβ is the charge-conjugated field with C = −iγ 0 γ 2 . 2 As the two quarks have the same flavour, the Pauli-exclusion principle applies and the wavefunction has to be completely anti-symmetric. With our choice to focus on S-wave combinations of particles, the spatial wave-function must be symmetric. As the colour (triplet) sextet has a (anti-) symmetric colour wavefunction, this forces the spin-wavefunction to be in a (triplet) singlet with (Γ = γ k ) Γ = iγ 5 . With these building blocks, we can form four classes ofbbbb colour-singlets by contracting the colour factors G in any irreducible representation (irrep) with its conjugate colour factors, i.e., 1 c × 1 c , 8 c × 8 c 3 , 3 c ×3 c and 6 c ×6 c . These SU (3) invariant colour contractions are given in Table I. After doing this, we need to project the operators onto a specific angular momentum J P by using the standard SO(3) Clebsch-Gordan coefficients (using a spherical basis of spin-matrices [18]) as with (P, Q) describing the blocks this configuration is built from, i.e, (η b , η b ), (Υ, Υ), (D, A), etc. We also allow the possibility of the two blocks being separated by a distance r. For the r = 0 case, the operators project onto a definite total angular momentum J. For the r = 0 case, one can Taylor expand O J1,m1 P (t, x + r) around r = 0 to notice that the operator projects onto a superposition of quantum numbers. Consequently, it is possible to utilise this to further search for the lowest ground state of the four quark system. When dealing with the diquark components, to project onto a definite value of chargeconjugation in Eq. (4) one can form the linear combina- In fact, not all of these colour combinations are independent. Fierz relations constrain the number of independent colour-spin operators that are possible. For the local operators in S-wave, the relations between the twomeson and diquark-antidiquark states are given in Table  II.
The simplest quantity that can be calculated on the lattice in order to extract particle masses is the Euclidean two-point correlator. This is defined as where the sum is over all distinct two-particle states X 2 with quantum numbers J P C and Ptot = 0, M S i (M K i ) is the static (kinetic) mass of the particles − as defined in is the reduced mass and Z 2l X2 are non-perturbative coefficients. Energies of states can be extracted using the above functional form once the correlator has been computed. Examining Eq. (5) in the path-integral formalism we can perform the connected Wick contractions 4 so that the correlator can be written as an integral over the gluonfields with the integrand consisting of products of b-quark propagators. For each two-meson type operator, e.g., , as all quarks have the same flavour there are four connected Wick contractions. These are shown diagrammatically in Figure 1.
The first Wick contraction for the two-meson correlator, called Direct1 and shown in Figure 1a, has the expression while the second, called Xchange2, is given by 4 Annihilation diagrams are suppressed by powers of the heavy where the sign function is defined by sgn(X) T = ±1 if X T = ±X. The ± in Eq. (10) corresponds to the 3 or 6 colour representation. It is this prefactor with the signs which enforces the Pauli exclusion principal: the sum cancels for spin combinations that do not make the wavefunction overall anti-symmetric. The spin-triplet/singlet configurations we consider here obey (Cγ k ) T = +(Cγ k ) and (Cγ 5 ) T = −(Cγ 5 ). Diagrammatically the four connected Wick contractions contributing to the diquark correlator are shown in Figure 2.
To calculate the two-point correlators described above within the first-principles Feynman path-integral approach to QCD needs the methodology of lattice QCD. We now discuss our lattice QCD approach.
quark mass [19] and are expected to be negligible.

III. LATTICE QCD METHODOLOGY
A. Second Generation N f = 2 + 1 + 1 Gluon Ensembles Our lattice calculation uses gauge field configurations generated by the MILC collaboration [20]. For the gauge fields, they used the tadpole-improved Lüscher-Weisz gauge action correct to O(α s a 2 ) [21] and include 2+1+1 flavours in the sea, the up and down quarks (treated as two degenerate light quarks with mass m l ), the strange quark, and the charm quark. The sea-quarks are included using the Highly Improved Staggered Quark formulation [22].
Four ensembles are chosen in this study. As one might expect roughly double the discretisation errors in a 2bb system relative to thebb one, to ensure that the heavy quark potential is accurately represented (where short distance details may be important for a compact tetraquark candidate) we utilise three ensembles that span relatively fine lattice spacings ranging from a = 0.06 − 0.12 fm. Details of the ensembles are given in Table III. Due to the computational expense, most of the ensembles use heavier m l than in the real world. However, to test m l dependence, we use one ensemble (Set 2 in Table III) that has physical am l /am s . Additionally, the ensembles have been fixed to Coulomb gauge to allow non-gauge invariant non-local operators to be used (as constructed in Eq. (4)).

B. b-quarks Using iNRQCD
A non-relativistic effective field theory is appropriate for physical systems where the relative velocity of the constituent particles inside the bound state is much smaller than one (in Planck units). When applied to the strong force, this framework is called non-relativistic QCD (NRQCD). It is well known that b-quarks are very nonrelativistic inside their low-lying bound states, where v 2 ≈ 0.1 [23].
For lattice NRQCD, the continuum NRQCD action is discretised onto the lattice [23] with operators included to a predetermined level in v 2 . Here we use an action accurate through O(v 4 ) with additional spin-dependent terms at O(v 6 ). 5 Operators are also added to correct for discretisation effects. We make the further systematic improvement here, introduced in [24], to include coefficients of O(v 4 ) operators that have been matched to continuum QCD through O(α s ). We call this improved-NRQCD (iNRQCD). This action has already been used to successfully determine bottomonium S, P and D wave mass splittings [24,25], precise hyperfine splittings [26,27], B meson decay constants [28], Υ and Υ leptonic widths [29], B, D meson mass splittings [27] and hindered M1 radiative decays between bottomonium states [30].
The iNRQCD Hamiltonian evolution equations can be 5 The spin-independent O(v 6 ) terms are subleading effects for thē bbbb energies relevant to this study.  [24,31], amq are the sea quark masses, Ns × NT gives the spatial and temporal extent of the lattices in lattice units and n cfg is the number of configurations used for each ensemble. We use 16 time sources on each configuration to increase statistics. Ensembles 1 and 2 are referred to as "coarse", 3 as "fine," and 4 as "superfine".
Set β a (fm) am l ams amc Ns × NT n cfg written as with The parameter n is used to prevent instabilities at large momentum from the kinetic energy operator. A value of n = 4 is chosen for all am b values. We evaluate the propagator using local sources Eq. (11). Here, am b is the bare b quark mass, ∇ is the symmetric lattice derivative, with∇ the improved version, and ∆ (2) , ∆ (4) are the lattice discretisations of Σ i D 2 i , Σ i D 4 i respectively.Ẽ, B are the improved chromoelectric and chromomagnetic fields, details of which can be found in [24]. Each of these fields, as well as the covariant derivatives, must be tadpole-improved. We take the mean trace of the gluon field in Landau gauge, u 0L = 1 3 Tr U µ (x) , as the tadpole parameter, calculated in [24,28]. The matching coefficients c 1 , c 2 , c 4 , c 5 , c 6 in the above Hamiltonian have been computed perturbatively to one-loop [24,32]. c 3 was found to be close to the tree-level value of one nonperturbatively [24] and we set it, as well as the rest of the matching coefficients, to their tree-level values. The quark mass am b is tuned fully nonperturbatively in iN-RQCD [26] using the spin-averaged kinetic mass of the Υ and η b (which is less sensitive to spin-dependent terms in the action). The above Hamiltonian neglects the fourfermion operators which appear at O(α 2 s v 3 ), as well as other operators which appear at higher order in the nonrelativistic expansion. Simple power-counting estimates [24,26] would lead us to expect contributions of order a few percent (or a few MeV) at most to binding energies from these terms. In the case where a tetraquark bound state is observed then we can estimate the systematic effect from neglecting these contributions.
The parameters used in this study are summarised in Table IV. There, c 1 , c 5 and c 6 are the correct values for an O(v 4 ) iNRQCD action [24]. For Set 4, all parameters are those for the O(v 4 ) action. The small changes to these coefficients in going to an O(v 6 ) iNRQCD action (which are similar in magnitude to the two-loop corrections) are not appreciable for the purpose of this work: whether or not a tetraquark candidate exists below the lowest bottomonium-pair threshold. All other parameters listed in Table IV are taken from [26,30].
Within iNRQCD the single-particle energy-eigenstates can be decomposed in the standard non-relativistic expansion as with M S , M K the static and kinetic masses respectively [24]. Because the quark mass term is removed from the iNRQCD Hamiltonian, the static mass is unphysi-cal, differing by a constant shift from the physical (kinetic) mass. This means that only static mass differences determined fully non-perturbatively can be compared to experimental results. However, this constant shift can be calculated in lattice perturbation theory if required [33], or by using an additional experimental input. The kinetic mass does not suffer this problem as it acquires the quark mass contributions from the quark kinetic terms [24]. Also, within iNRQCD the Dirac field Ψ can be written in terms of the quark ψ and anti-quark χ as Ψ = (ψ, χ) T . The propagator is then where G ψ (x|y) is the two-spinor component quark propagator and G χ (x|y) is the two-spinor component antiquark propagator. Taken together, we can now compute the b-quark propagator via the iNRQCD evolution equations on the gluon ensembles listed in Table III. The last piece needed to calculate the two-point correlators, and hence energies, are the discretised finite-volume versions of the interpolating operators.

C. Discrete Finite-Volume Operators
Together, the isotropic discretisation and the periodic finite-volume break the infinite-volume continuum SO(3) symmetry of NRQCD to the octahedral group, O h [34]. Thus, while the operators constructed in Sec. II have well-defined J P C quantum numbers associated with SO(3), we need to construct operators which transform within the irreps of the O h symmetry group (relevant for lattice calculation). This can be achieved by the method of subduced representations, where an operator with a specific J P C can be taken to a specific lattice irrep 6 Λ P C by using the subduction coefficients found in Appendix A of [18]. At rest, each of our J P C = 0 ++ and J P C = 1 +− operators trivially subduce into a single lattice irrep labelled by A ++ 1 and T +− 1 respectively. However, the J P C = 2 ++ case is slightly more complicated and subduces into two lattice irreps labelled by T ++ 2 and E ++ (which are three-and two-dimensional). We construct both the T 2 /E operators as √ 2O [2] T2/E = O J=2,m=2 ± O J=2,m=−2 , which are correctly subduced from the J = 2 ++ operators defined in Sec. II. In principle, each lattice irrep allows mixing between different J states, e.g., the A 1 irrep contains not only the J = 0 states but also the J = 4 [34]. However in practice since TABLE V. Thebbbb interpolating operators used in this study. Operators in each column are subduced from the infinite-volume continuum quantum numbers J P C given in the first row. The superscript on each operator denotes the lattice irrep of that operator and the subscript denotes the building blocks of the operator, as explained in the text. We generate each operator with three different spatial configurations as shown in Eq. (4): where the building blocks are separated by a distance rx = 0, 1 or 2 lattice units in the x-direction.
we are only looking for the ground state of thebbbb correlators we are not sensitive to these mixing effects.
A complete list ofbbbb interpolating operators used to produce the correlator data herein is given in Table  V. In fact, this is an over-constrained set due to the Fierz identities (shown in Table II) which relate the twomeson and diquark-antidiquark correlators. We include this over constrained system and ensure that we reproduce the Fierz relations to numerical precision, performing a non-trivial check on our data. Additionally, we also reproduce the relations between the 8 c × 8 c colour combination and the others [36] on a subset of the data.
It has been found [37] that separating each hadron within the two-hadron interpolating operator by a specific distance r can significantly aid in the extraction of the (ground) state energy. In this direction, we use three different spatial configurations of thebbbb correlators where the individual building blocks are separated by a distance of r x = 0, 1 or 2 lattice units in the xdirection 7 .
Finally, the subduced lattice interpolating operators are defined in terms of the Dirac fields as in Eq. (1), and the correlators are defined analogously, as in Eq. (10). We can then use the decomposition of K −1 (x|y) given in Eq. (15) with suitable boundary conditions to write the correlator in terms of the iNRQCD quark propagator G ψ (x|y). Due to our use of iNRQCD, there are no backward propagating valence anti-quarks in our calculation.
Consequently, the appreciable finite-temporal effects seen in relativistic two-meson lattice correlators 8 [38] do not arise in our calculation, simplifying the analysis. With this methodology, it is now possible to compute the lowest energy levels associated with thebbbb system.
In order to determine if there is an energy eigenstate below the 2η b threshold, we first need to find the noninteracting thresholds on each ensemble listed in Table  III. This is achieved by computing the bottomonium η b and Υ two-point correlators as described above, then fitting them to the functional form given in Eq. (6) to extract the single particle energies. As in the range of studies listed in Sec. III B, amongst others, we utilise the well-established Bayesian fitting methodology [39] in this work and refer the reader to [24] for technical details. Although we fit the correlator data in order to extract fit parameters, in the following we display effective mass plots so the reader can visualise the data. The single particle effective mass is constructed as As can be inferred from Eq. (17), the greater the mass gap E 1 − E 0 or the larger the ground state overlap Z 0 then the quicker aE eff J P C converges to a plateau, which gives aE J P C . The η b and Υ effective masses are shown in Figure 3, where the returned fitted ground-state energy from the correlator fits is also shown overlaid in black. The large difference in energy between the η b (Υ ) and the ground state η b (Υ) means that the effective mass plots fall rapidly to a plateau given by the ground state energy, indicating that the fit to the correlator data will extract the ground state energy precisely.
Also evident from the effective mass plot is the constant signal-to-noise in the η b data, as might be expected from straightforward application of the Parisi-Lepage arguments [40,41] for noise growth in a system where all the quarks are the same and no 0 ++ bound tetraquark exists. The argument specifies that the noise of the two-point correlator should behave like The effective mass plot for the η b and Υ on the "superfine" ensemble (Set 4 listed in Table III). The effective mass plots on the other ensembles are qualitatively identical.
is the lowest energy eigenstate of the bottomonium operator O J P C constructed to have the quantum numbers J P C and E GS is the lowest energy eigenstate of the mean squared correlator which controls the noise. Thus, it would be surprising if a tetraquark candidate did exist below the 2η b threshold from the lattice perspective alone as then E GS < 2E η b and the noise of the η b data would grow exponentially. However, the lattice calculation still needs to be performed for a conclusive statement to be made about the existence of this tetraquark candidate since the Parisi-Lepage arguments do not allow for raw crossed Wick contractions that would contribute to either the full twomeson or tetraquark correlator. Lattice correlators are affected by both the discrete nature of the space-time lattice and separately by its finite volume. Each has a separate but calculable effect on the extracted lattice energies. Corrections in energies due to discretisation effects are proportional to a k , where k depends on the level of improvement. Here systematic discretisation errors are reduced to α 2 s a 2 by the improvements made, as for those studies listed in Sec. III B, and we expect this to be small enough to have little impact. We can assess this from our results with different values of the lattice spacing.
Finite-volume effects for single-particle energies (arising from the lightest particle in the sea propagating around the spatial boundaries) are known to behave like exp(−aM π N s ) [42] and are not appreciable for the ensembles used here which have aM π N s ≈ 4 [43]. In fact, N s on Set 1 and 2 differ by a factor of two, giving a basic test of volume-dependence. However, finite-volume interactions can shift a two-particle energy by an amount that depends on the infinite-volume scattering matrix. Further, these shifts are non-trivial to parameterise (see for example [44][45][46][47][48][49][50]). As the specific purpose of this study is to search for a hypothetical tetraquark bound state be- Assuming a tetraquark exists 100 MeV below the 2η b threshold, different normally distributed couplings inbbbb correlator mock-data produces different effective mass curves on the "superfine" ensemble (Set 4 listed in Table III) as described in the text. (colour online) low the lowest threshold, we will not attempt to quantify these finite-volume interactions. Eq. (7) describes the non-relativistic two-particle contribution to the correlator after t = 1 fm (as shown in Appendix A). In this case, because of the additional 1/t 3 2 time-dependence, the effective mass formula for these contributions differs from their single-particle counterparts. Removing the leading order time-dependence yields an effective mass defined as For the 0 ++ , 1 +− and 2 ++ operators that are constructed, if a stable tetraquark exists below the 2η b , η b +Υ or 2Υ thresholds then it would show up as the ground state of each correlator and hence also in the effective masses. Otherwise, each threshold would be the lowest energy eigenstate. Higher energy states will also appear in each correlator. For example, the 2Υ and η b η b in the 0 ++ case, the Υη b in the 1 +− and the ΥΥ and χ b0 χ b2 (1P ) in the 2 ++ . Of these, if no tetraquark state is present, only the 2Υ in the 0 ++ might be noticeable while studying the ground state, as it is O(100) MeV above the 2η b . All other excited states have similar energy differences to those appearing in the η b /Υ effective masses shown in Figure 3, which rapidly falls to the ground state.
It is helpful at this stage to generate mock correlator data to illustrate how we might expect thebbbb correlator results to behave in the presence or absence of a low-lying tetraquark state (neglecting two-particle finite-volume effects). Using the extracted lattice η b and Υ masses, we can compute the non-interacting 2η b and 2Υ thresholds on our ensembles. Further, for a fixed value of Z 0 X2 , we can also compute their leading order two-particle contribution to the correlator from Eq. (7). Additionally we can infer the values of the non-interacting η b η b and ΥΥ masses on our ensembles using the experimental PDG [4] values as input in order to include their two-particle contributions in the mock correlator data also. Further, in the 0 ++ channel, if we assume a tetraquark boundstate exists 100 MeV below the 2η b threshold 9 , for a fixed value of the tetraquarks non-perturbative overlap, Z 4b , this hypothetical state's contribution to the correlator is given by Eq. (6). Then, given the effective mass formula defined in Eq. (16), for each different choice of the non-perturbative coefficients we can generate a separate effective mass curve. In practice, we choose different values of the coefficients from a normal distribution with zero mean and unit variance. Figure 4 shows such a plot for the "superfine" ensemble (Set 4 in Table III) where the solid blue curves represent different values of the normally distributed coefficients. As the non-perturbative coefficients show up through the ratio Z 0 X2 /Z 4b in the effective mass formula, analogously to Eq. (17), once the energy difference E 1 − E 0 is set the effective mass is only sensitive to the relative size of the of the tetraquark overlap to that of the lowest threshold (once the contribution of excited states has become negligible). With this knowledge, the lower red dotted curve gives mock-data in the situation where there is only a new state present in the correlator data (Z 0 X2 = 0), the middle red dashed curve indicates the case when the tetraquark and two-meson states have the same value of coefficient while the upper dot-dashed red curve is the mock-data in the case where there is no new state in the data (Z 4b = 0). The blue curves below the middle red one have a larger overlap onto the new state while those above have an increasingly vanishing one. With our overconstrained colour-spin basis of S-wave operators at least one operator should have an appreciable overlap onto a tetraquark state below threshold (if it exists) and, as illustrated by the mock-data, give a easy/clean signal similar to the middle red dashed curve where the effective mass drops below the 2η b threshold. Even though Figure  4 is mock-data, the upper red curve with no tetraquark state present looks very like the real data shown in Figures 5 and 6.
One important point to note is the additional factor of 1/t 3 2 appearing in the two-particle contribution in Eq. (7) relative to the single-particle one in Eq. (6). This factor suppresses the two-particle contribution relative to the single-particle state, e.g., t = 100 gives a suppression of the two-particle state of (0.01) 1.5 . This is one reason why the middle red dashed curve has a particularly rapid fall to the new state which lies only 100 MeV below the 2η b (compared to Figure 5d where the 2Υ is O(100) MeV above the 2η b ). Overall this effect would produce an enhancement of the stable tetraquark state if it exists. Further, while we used the "superfine" ensemble for 9 The smallest binding of a below thresholdbbbb tetraquark from the phenomenological studies (which are shown in Figure 17) was 108 MeV [9].     illustrative purposes, the other ensembles with larger lattice spacings have similar features but with less temporal resolution due to the larger numerical value of the lattice spacing.
If one is searching for a bound state below threshold then both single-and two-particle contributions appear in the correlator and one may use the single particle effective mass formula given in Eq. (18) in order to highlight the bound ground state mass (as done in the mock data above). However if no bound state exists in the data, which we do not know a priori, then not removing the 1/t butions has an important effect: an additional factor of 1.5 log(1 + 1/t) is introduced into the effective mass formula in Eq. (18), giving a contamination that vanishes slowly as t → ∞. This would produce a confusing picture of what is actually contributing. The two-particle effective mass formula Eq. (19) removes this contribution.
In the results to be reported now, we will overlay both single-and two-particle effective masses on the same plot for the reader's convenience.
We generate thebbbb correlator data for the operators given in Table V using the ensembles listed in Table III and fit this data simultaneously with the bottomonium meson data so to include correlations between data sets. All thebbbb data within a specific irrep and those which are unrelated by a Fierz relation 10 are fit using Eq. (7) for the two-particle contributions and Eq. (6) for a hypothetical tetraquark state below threshold. The mean of the prior energy of the 2η b , η b + Υ and 2Υ thresholds are roughly estimated based on the effective masses and then given a suitably wide prior width of 100 MeV while a tetraquark state prior energy is taken to be 250(100) MeV below each threshold. As can be seen in Figure 5, since the data plateaus to the non-interacting 2η b threshold, no energy eigenstate is found below this threshold and variations of the tetraquark prior energy are insignificant. Similar behaviour is seen with the other quantum numbers.
Again, while we fit the correlator data in order to extract particle energies, so that the reader can visualise this data we display effective mass plots on the different ensembles. The "superfine" ensemble (Set 4) 0 ++ two-meson effective masses are shown in Figures 5, while the 0 ++ diquark-antidiquark are given in Figure 6, the "physical coarse" (Set 2) 1 +− two-meson and diquarkantidiquark in Figure 7 and the "fine" (Set 3) 2 ++ twomeson and diquark-antidiquark correlators subduced into the T 2 lattice irrep are shown in Figure 8. Each plot has the fitted ground state energy overlaid in black for comparison. The 2 ++ subduced into the E irrep is similar to the T 2 case. Further, the behaviour of the lattice data on all ensembles is qualitatively similar to those shown. The extracted ground state energies in each channel are given in Table VI and a comparison of the energies is shown in graphical form in Figure 9.
It should also be noted that the numerical value of the effective mass (shift) plateau in two-hadron correlators has been shown, in certain cases, to be sensitive to the choice of interpolating operators [51,52]. There, the authors found that when the noise growth in the correlator data restricts the study of effective masses to a maximum propagation time of approximately 2 fm, fake plateaus can appear. These fake plateaus can be a consequence of different choices of source and sink (smeared) 10 Simultaneously fitting data sets related by a Fierz identity would mean the correlation matrix would have a zero eigenvalue and thus not be invertable for use in a least-squares minimisation. aE eff D3×3→3×3 (t)  6. Thebbbb effective masses for the 0 ++ diquark-antidiquark3 ×3 and 6 ×6 correlators. E eff and E eff,t are the single-and two-particle effective masses defined in Eq. (18) and Eq. (19) respectively. The diquarks are separated by a distance rx in the x-direction when constructing the diquarkantidiquark interpolating operator as given in Eq. (4). (color online) operators: this can cause the a negative sign in the Z 1 term in the effective mass formula Eq. (17) and a dip below threshold can appear for a short time range which can be misinterpreted as a bound state. Wall-sources were shown to be particularly prone to this behaviour of producing a "false dip" and obtaining an appreciably different "plateau" than a Gaussian source. Here we only use local quark sources. In addition, the elastic scattering states can also have a dependence on the choice of operator, which can cause a slow decay to the ground state and mimic a slowly varying effective mass that can be mistaken for a plateau over a short time range. As noted in these studies, a necessary check for a real effective mass plateau when using different source and sink operators is the convergence of all data to a single plateau at times larger than approximately 2 fm. As we separate the operators by r x = 0, 1 and 2 lattice units (as described in Section III C) and propagate to t > 8 fm, this is a consistency check we satisfy. t (fm) a0.12fm Ns48 Thebbbb effective masses for the 1 +− Υη b and diquark-antidiquark3 ×3 correlators. E eff and E eff,t are the single-and two-particle effective masses defined in Eq. (18) and Eq. (19) respectively. (color online) A few notable features of thebbbb effective mass plots are evident. First and foremost, no value of the effective mass is observed below the lowest non-interacting bottomonium-pair threshold in any channel, in line with what one would expect if no stable tetraquark candidate existed below threshold. Indeed, thebbbb effective mass plots are strikingly similar to the upper dot-dashed curve in the mock data in Figure 4 where no bound tetraquark state is present. Additionally, the 2η b → 2η b effective mass shown in Figure 5a plateaus very early due to the larger overlap onto the 2η b threshold, while the 0 ++ 2Υ → 2Υ effective mass shown in Figure 5d falls more slowly to the 2η b threshold due to the larger overlap onto the nearby 2Υ threshold. The cross correlators 2η b → 2Υ show how the different operators converge to a single plateau at a time greater than t ≈ 4 fm, a necessity for a true plateau as discussed above. The local diquark-antidiquark 0 ++ correlator data is a linear combination of the 2η b and 2Υ two-meson data, related by the Fierz identities given in Table II, and the effective masses shown in Figure 6 reflects this. It is empirically observed that separating the diquark from the   anti-diquark by too large a distance r x results in larger noise due to the separation of colour sources. The 1 +− Υη b and diquark-antidiquark effective masses are shown in Figure 7, where the noise starts to increase after t ≈ 7 fm due to the Parisi-Lepage argument mentioned above with the noise being set by at least the 2η b threshold. Based on this, one would also expect the signal-to-noise to be worse for the 2 ++ data, which is also evident from the correlator data subduced into the T 2 irrep shown in Fig. 8. As Set 2 has physical m l /m s corresponding to a pion mass of O(131) MeV, while the other ensembles have nonphysical m l /m s corresponding to pion masses of O(300) MeV [53], no sensitivity to light sea-quarks is observed. This would be expected from the smallness of the Van-der-Waals potential generated by the two-pion exchange between two 1S bottomonium mesons [54]. As can be seen in Figures 8 and 9, the 2 ++ ground state obtained from the lattice is slightly higher than that of the non-interacting threshold. However, this is the state which has the largest signal-to-noise, restricting the data to shorter time regions and it is possible that we are sen-  sitive to the same aforementioned issue of a slowly varying fake plateau. Alternatively, this positive shift in the two-particle energy could potentially indicate appreciable infinite-volume continuum scattering arising from finitevolume interactions [44], but quantifying these phaseshifts is outside the remit of this study. Regardless, these effects do not indicate that a bound tetraquark state exists in this channel. For illustration purposes, we also show the effective masses of the individual Wick contractions contributing to the 2η b → 2η b correlator in Figure 10. As is evident, in each individual Wick contraction the effective mass drops below the 2η b threshold but then rises slowly to threshold. However, importantly, when all Wick contractions are added together to yield the full correlator (shown Figure 5a) the effective mass falls rapidly to threshold from above. This behaviour will be discussed further in Sec. VI.
After analysing all our data, as hinted by the effective mass plots, there is no indication of a bound tetraquark state below the non-interacting thresholds on any ensemble, as shown in Fig. 9. We see no evidence of any change in the ground-state energies (with respect to the thresholds) as we vary lattice spacing or sea quark masses.
Searching for a new tetraquark candidate has at least a two-dimensional parameter space: the hypothetical state would have an energy and also an overlap onto a specific operator. Using the lattice data presented here, we   can determine a relationship between these parameters. Assuming that a tetraquark does exist below the lowest bottomonium-pair threshold in our data, at a certain time t * the correlator can be modeled with a two-state ansatz. Specifically for the 0 ++ channel, given that the tetraquark has an energy E 4b and an overlap Z 4b onto a particular operator, at a large enough time the only other appreciable contribution will come from the higher 2η b threshold which has an overlap Z 2η b with the same operator. In this case the correlator is given by perturbative overlaps yields the constraint with a∆E = aE 2η b − aE 4b > 0. As can be seen, if the tetraquark is not observed by a time t * then the overlap onto this new state must be (at least) exponentially suppressed with the binding of the tetraquark state, e.g, with −∆E. This point illustrates that if a tetraquark did exist with E 4b < E 2η b then it is possible that it was not observed in our data because Z 4b ≈ 0 within statistical precision. In this scenario, we can use the constraint Eq. (21) to estimate an upper bound on the magnitude of the overlaps given that no clear evidence of the tetraquark is observed within our statistical precision. The needed inputs for the constraint include the value of t * where the two-state ansatz is valid, aE eff (t * ) from correlator data constructed with a specific operator, as well as aE 2η b = 2aM η b . For the local O (η b ,η b ) operator on the a = 0.06 fm ensemble, by examining Figure 5a, a choice of t * = 143 ensures that the two-state ansatz is valid (given the long plateau at the 2η b threshold). Here, aE eff (t * = 143) = 0.87634 (61) can also be precisely obtained. The value of aM η b , given in Table VI, is found from the η b -meson data. Then, using this data in the constraint, for a certain choice of E 4b , a numerical value of the ratio of overlaps is found such that it is consistent with zero within its small 1σ statistical error. Any value of the ratio of overlaps larger than this 1σ error is inconsistent, at this level of confidence, with our data observing a tetraquark at this value of E 4b . We use this model to estimate how small the hypothetical tetraquark overlap would need to be so that the tetraquark was not observed within our statistical precision. A 1σ, 3σ and 5σ exclusion plot of the parameter space is given in Figure 11. As the input data into the constraint has a long propagation time past t * > 8 fm and a statistically precise value of aE eff (t * ) which does not fall below the threshold, a significant amount of parameter space is excluded. It should be understood that this figure is only valid for a particular operator in a certain channel. The given 0 ++ channel in Figure 11 excludes the largest amount of parameter space as it is the most statistically precise. Also Figure 11 is constructed from data on the "superfine" ensemble alone, where discretisation effects are smallest and would not change the quantitative behaviour significantly.
To conclude this section, we find no evidence of a stable tetraquark candidate below the non-interacting thresholds by studying a full S-wave colour-spin basis of QCD operators. In the next section we will perform an exploratory and complementary study of an alternative approach so to ensure the robustness of our conclusions.

V. NRQCD WITH A HARMONIC OSCILLATOR POTENTIAL
A stable tetraquark state in the 0 ++ , 1 +− and 2 ++ channels would overlap with the full basis of S-wave colour-spin operators utilised above. In this section, we go one additional step by exploring an alternative approach in order to further investigate the possibility of a tetraquark state. Adding a central confining potential to the quark interactions can produce a more deeply bound tetraquark relative to the threshold as the strength of this interaction is increased (as we will see below). Furthermore, an appropriate choice of additional interaction can reduce the fiducial volume of the lattice and thus thin the allowed discrete momenta states of the two meson degrees of freedom. Adding an external attractive scalar central-potential to the QCD interactions yields these desired effects.
The harmonic oscillator potential is a particularly suitable choice of scalar interaction between quarks. For a particle of mass m at position x away from the centre x 0 the potential is just κr 2 /2 ≡ κ|x − x 0 | 2 /2. Defining ω = (κ/m), the ground state energy and wavefunction are E 0 = 3ω/2 and ψ(r) = C exp (−mωr 2 /2). Additionally, the separability of the combined QCD and harmonic oscillator potential into total and relative coordinates ensures that solutions of multiquark systems can be split into two parts with the total coordinate piece analytically solvable. This follows from the nature of the harmonic oscillator potential 11 . Thus we expect that for  small values of ω the ground state for 2η b mesons is approximately 2M η b (ω) + 3ω. The 3ω term comes from two colour-singlet η b mesons in the harmonic oscillator potential and the mass of the η b is shifted slightly from the value at ω = 0 because of the additional harmonic oscillator interaction combined with the QCD interactions that bind the two heavy quarks into a η b 12 . However, for a compact tetraquark state the mass would be M (bbbb)(ω) + 3ω/2, as there is only one colour-singlet state in the central harmonic oscillator potential (3ω/2) and if the tetraquark state is also a compact state (on the scale of η b and much less than the effective lattice volume) then its mass will also receive only a modest positive correction due to ω.
Hence if there were a tetraquark state near threshold then this additional interaction could drive it further below threshold, giving a much cleaner and distinct signal for the tetraquark candidate in our calculation. As the potential model framework describing the η b has been largely successful, we can use this framework as a general guide for the exploratory non-perturbative lattice calculation when including the harmonic oscillator potential. Of course, we are mainly interested in QCD (and not the harmonic oscillator) and, as such, if we do find a stable tetraquark state when including the harmonic oscillator potential then we must take the ω → 0 limit. Thus, the objective of this section is to determine if a stable tetraquark state exists when the quarks are exposed to an auxiliary potential (which could push the tetraquark increasingly lower than the threshold) and if it does, will it survive the QCD limit. The harmonic oscillator potential is defined as 13 where the quarks are pulled towards the fixed point x 0 with a strength ω. We choose x 0 to be the same position as the source. Intuitively, as the quarks/η b 's start to propagate further away from the source, the harmonic oscillator potential pulls them closer together. In turn, this restricts the quarks/η b 's to be in a certain volume. First, it is necessary to determine which volumes the quarks/η b 's are confined into by the addition of the harmonic oscillator potential. The root mean square distance, r rms , gives an indication of this. We determine the r rms of the η b based on solutions of the Schrodinger equation using a Cornell potential 14 . The Matthieu equation can describe the behaviour of two free η b 's in a harmonic oscillator potential on a periodic box, and the solutions of which can yield r rms for the 2η b state. The results from such a calculation are plotted in Figure 12.
The harmonic oscillator is implemented through a minor modification of the NRQCD evolution equations via where e −aH is the purely NRQCD evolution equation defined in Eq. (12). This implementation was chosen so that the evolution equation is still time-reversal symmetric. Here, l is a stability parameter akin to n in Eq. (12) which is used to prevent possible numerical instabilities [23]. Values of l = 13 and 10 were chosen for the calculations on Set 1 and 3 respectively. Following these details, we are now able to present results from the nonperturbative lattice calculations.

A. Numerical Results
All correlator data from Set 1 and 3 discussed in Section IV was generated again with the inclusion of the harmonic oscillator potential at the four different ω values given above. The harmonic oscillator alters both the 13 In a periodic box of length L the harmonic oscillator is also periodic. 14 The potential form is − 4αs 3r + r b 2 with αs = 0.36, b = 2.34 GeV and reduced mass 2.59 GeV ('hhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhh',) single-and two-particle contributions to the correlator so that they become (as derived in Appendix B) dependent on ω as First, Figure 13 shows the effective masses, as defined in Eq. (16), of the η b when including the harmonic oscillator. Also overlaid are the effective masses when removing the 1 + e −2ωt dependence to enable a better comparison with the data when no harmonic oscillator is included. As can be seen, the dip in the harmonic oscillator effective masses is from this additional time dependence. Physically, this can be understood to be due to the bquarks travelling non-relativistically and so it takes time for the harmonic oscillator to have an effect. The η b correlator data when including the harmonic oscillator potential is fit to the functional form given by Eq. (24) in order to extract the lowest energy eigenstate M (ω) η b +3ω/2 from the asymptotic behaviour. We show the fitted result overlaid on the effective mass plot in Figure 13. As before, the long plateau indicates that the ground state will be extracted accurately. To compare to the potential model predictions, we subtract the η b mass with no harmonic oscillator included (M (ω = 0) η b ) and then plot the energy differences against ω, as shown in Figure 14. Good qualitative agreement between the lattice results and the potential model predictions is observed.
For thebbbb system, we show the 0 ++ effective masses on Set 3 in Figure 15. It is evident that the 0 ++ and the η b data contains more noise when a harmonic oscillator potential is included. While fitting the data to the form in Eq. (25) can be performed, it is not necessary as the purpose of this exploratory work is to determine if a stable tetraquark exists when ω = 0. As can be seen, there is no fall below the 2η b threshold for any value of ω. Similar behaviour is seen with the data on Set 1.
We show the effective masses for the individual Direct1 and Xchange2 Wick contractions of the 2η b → 2η b correlator in Figure 16. As before, the effective masses of the individual Wick contractions drop below the 2η b threshold, even though importantly, when added together to yield the full correlator shown in Figure 15a the effective mass is always above threshold. As this was also seen in the pure NRQCD data shown in Sec. IV, it may be a problematic feature of models that utilise a phenomenologically motivated four-body potential for this system.
To conclude this section, despite adding an auxiliary potential into the QCD interactions that should push a near threshold tetraquark candidate increasingly lower we find no indication of any state below the 2η b threshold. The conclusions of this section then agree with those of Sec. IV.

VI. DISCUSSION AND CONCLUSIONS
In this work we have studied the low-lying spectrum of thebbbb system using the first-principles lattice nonrelativistic QCD methodology in order to search for a stable tetraquark state below the lowest non-interacting bottomonium-pair threshold in three different channels: the 0 ++ which couples to the 2η b and 2Υ, the 1 +− which couples to Υη b and the 2 ++ which couples to 2Υ.
In Section III we describe our numerical methodology. Four gluon ensembles were employed with lattice spacings ranging from a = 0.06 − 0.12 fm, and one ensemble which has physical light-quark masses. All ensembles have u, d, s and c quarks in the sea.
In Sec. IV we presented the majority of the results in this work. Here, we determined the lowest energy eigenstate of thebbbb system with the quantum numbers 0 ++ , 1 +− and 2 ++ using an over-constrained Swave colour/spin basis (arising from Fierz relations between the diquark-antidiquark and two-meson systems as shown in Table II). We did not observe any state below the lowest non-interacting bottomonium-pair threshold in any channel, as can be seen in Figures 5, 6, 7 and 8, and a summary of our results from this section is given in Figure 9.
In Sec. V, to ensure the robustness of our conclusions, we performed an exploratory calculation of a novel method which added an auxiliary scalar potential into the QCD interactions with the objective of pushing a near threshold tetraquark increasingly lower than the threshold. This would give a more distinct and cleaner signal for its presence in our calculation. The harmonic oscillator was found to be a suitable central scalar potential. For the η b -meson with this potential, we first verified agreement between the non-perturbative lattice calculations and a potential model (as shown in Figure 14) and then used this potential model as a general guide to choose multiple appropriate values of the potential strength. Despite studying thebbbb system with this additional scalar potential on the lattice, no indication of the QCD tetraquark was observed.
This work is the only first-principles study of the lowlyingbbbb spectrum in the literature. However, there are others which utilise different methodologies. For example, [8] predicts the tetraquark mass by solving the twoparticle Schrodinger equation with a phenomenologically motivated non-confining potential between the point-like diquark and anti-diquark, and finds a 0 ++ , 1 +− and 2 ++ tetraquark to be bound by 44, 51 and 5 MeV respectively 15 . The authors of [14] used a diquark model including a confining linear potential, but neglected spin effects, and found a 0 ++ tetraquark to be bound by 48 MeV. However, it has also been found that the root-meansquare distance between the diquark and anti-diquark inside the tetraquark in this model is similar in magnitude to the distance between the quarks inside each diquark [9]. Consequently, such an approach is internally inconsistent. More recently, [10] used a Hamiltonian including only spin-spin interactions mediated by a one-gluon exchange. Here, all other effects such as the chromoelectric interactions, colour confinement and the b-quark mass need to be set separately. The authors set these additional contributions in two ways: by estimating the effects using an effective heavy quark [10] or by using the experimental meson masses as input. In this way, the authors find that the tetraquark state could be either below the 2η b threshold or lie inbetween the 2η b and 2Υ thresholds (where in both cases the thresholds were determined using the experimental meson masses). In [55] the author also uses a model including only chromomagnetic interactions and finds an unbound tetraquark. Within the QCD sum-rules framework, [11] finds a tetraquark candidate approximately 300 MeV below the experimental 2η b threshold while [12] finds a tetraquark lying inbetween the experimental 2η b and 2Υ thresholds. Using phenomenological arguments, [13] also finds a tetraquark candidate lying inbetween the thresholds. Indeed, in the limit of very heavy quarks where the force proceeds through one-gluon exchange containing only colour-Coulomb contributions (safely neglecting spin and long distance effects), the authors of [7] used a variational methodology to determine that a bound tetraquark exists for a QQQ Q system when m Q /m Q < ∼ 0.15 (where both m Q and m Q are heavy relative to Λ QCD ). However, if m Q /m Q is varied and the tetraquark becomes unbound, then as the free two-meson eigenstate becomes the ground state of the system this numerical methodology has an increasingly slow convergence to a solution [7] (being numerically ill-posed) due to a redundant degree of freedom in the minimisation procedure. Indeed [7] indicates that for m Q /m Q = 1 the solution is unstable, a hint that no bound tetraquark exists for all identical quarks in the very heavy mass limit. The authors of [14] assume that the b-quark is sufficiently heavy so to use the one-gluon exchange with only colour-Coulomb contributions, and by also neglecting the mixing between different colour-components of the 2 × 2 potential matrix, finds a tetraquark bound by 78 (20) MeV (by using the experimental η b mass to determine the threshold). In an orthogonal direction to the above work, the authors of [15] only include a linear string contribution in the one-gluon exchange (neglecting spin effects and the appreciable short-distance Coulomb contributions), and without mixing between the different colour components of the potential matrix, find a bound tetraquark when m Q /m Q = 1. However, in subsequent work [56,57], by modelling the aforementioned mixing they concluded that no bound tetraquark exists. Perhaps the most sophisticated non first-principles methodology used to study the four-bodybbbb tetraquark is the diffusion Monte Carlo method utilised by [9]. Here, one determines the ground state of a phenomenologically motivated Hamiltonian by solving the Schrodinger equation and examining the stability of n e −(En−E0)t Ψ n (x) to determine E 0 , where E n (Ψ n ) is the n-th energyeigenstate (eigenfunction). The authors include both the colour-Coulomb and linear contributions in the gluon exchange but neglect the mixing between the different     colour components in the potential matrix, and find a stable tetraquark candidate 108 MeV below the 2η b threshold (determined from the experimental meson mass).
Consequently, there is no study in the literature which is not from first-principles that includes all the appreciable effects relevant for the bbbb system: treating the bottom-quarks as fundamental particles, including both short and long distance effects in the gluon exchange and including the mixing between the different colour components in the 2 × 2 potential matrix. It should be emphasised however that these studies, unlike ours, are not from first-principles and thus have an unquantifiable systematic error associated with the choice of four-body potential. To emphasise this further, thinking of each Wick contraction (shown diagrammatically in Figure 1) as a different potential contributing to the QCD dynamics, then only studying a subset of these interactions can change the energies of states. This can lead to the misindentification of a new state below threshold. For example, the effective masses of the individual Wick contractions contributing to the 2η b → 2η b correlator are shown in Figure 10. As is evident there, in each individual Wick contraction the effective mass drops below the 2η b threshold but rises slowly to threshold even though when all Wick contractions are added together to yield the full correlator (shown Figure 5a) the effective mass falls rapidly to threshold from above. This behaviour is even more pronounced in the data with the additional scalar potential, shown in Figure 16, possibly indicating that this may be a problematic feature of models that utilise a phenomenologically motivated four-body potential: a subset of the interactions show    behaviour that may be misinterpreted as a bound state below threshold, while when all interactions are included no bound state is seen. Particularly, the slow rise to threshold from below could make the diffusion Monte Carlo method practically difficult due to the slowly varying stability condition combined with the fact that a long evolution time (greater than 8 fm) is necessary.
In conclusion: we find no evidence of abbbb tetraquark with a mass below the lowest non-interacting bottomonium-pair thresholds in the 0 ++ , 1 +− or 2 ++ channels. We give a constraint in Eq. (21) that future phenomenological models must satisfy if such QCD states are postulated. For the 0 ++ channel, we use this constraint to estimate how small the non-perturbative overlap of the hypothetical tetraquark (onto a particular operator) would need to be, relative to the 2η b , so that it was not observed within our statistical precision. A ject to the value of its non-perturbative overlap as given in Figure 11. In this comparison, we take our ground state energy obtained on the "superfine" ensemble (Set 4 in Table III) as a representative because it has the smallest discretisaion effects and the statistical error encompasses the results on the other ensembles as shown in Figure 9. The y-axis labels results from different phenomenological models [7][8][9][10][11][12][13]  1σ, 3σ and 5σ exclusion plot of the parameter space is shown in Figure 11, and discussed in Sec. IV. As we have propagation times longer than 8 fm and statistically precise data, we can exclude all but the most finely-tuned parameter space. Our lattice results then rule out the phenomenological models discussed above that predict a tetraquark below the lowest bottomonium-pair thresholds which have a value of non-perturbative overlap that is excluded by Fig. 11. A comparison of these results with ours is shown in Figure 17.
Further studies of possible heavy tetraquark channels that include orbital angular momentum either between the mesons in the tetraquark or between the quarks in the meson could be performed with the methodology used here 16 . Similarly one could also study whether stablē cccc,bcbc orbbcc tetraquarks exist or not. Additionally, two-hadron systems receive a finite-volume energy shift which depends on the infinite-volume scattering amplitude which is non-trivial to parameterise. Here we do not calculate these finite-volume energy shifts. Doing so in a more extended study would allow statements to be made about the existence of resonant tetraquark states above the lowest thresholds, that likely do exist in nature. Quantifying these shifts would be an exciting avenue for future work.
Finally, recent work based on heavy-quark symmetry [6] and phenomenological arguments [58] indicates that a J P = 1 +bb ud tetraquark will be stable in QCD. In fact, by extracting a potential from the lattice in the static heavy-quark limit and solving the Schrodinger equation [59][60][61] also finds binding in this channel. Initial lattice calculations hint that such a state exists [62] but calculations are difficult because of a signal-to-noise problem for heavy-light states [63]. Lattice QCD calculations in this direction are essential for a conclusive first-principles statement to be made and to give further motivation for a targeted experimental search for these tetraquark configurations of nature.
The unaveraged correlator data, that have been analysed to produce the results in this work, are publicly available in a SQLite database from any of the authors upon request 17 .

ACKNOWLEDGMENTS
We would like to thank William Bardeen and Zhen Liu for the many insightful discussions on tetraquarks, as well as Gavin Cheung who in addition gave guidance on the tetraquark implementation. We are also grateful to the MILC collaboration for the use of their gauge configurations. This manuscript has been authored by Fermi Research Alliance, LLC under Contract No. DE-AC02-07CH11359 with the U. S. Department of Energy, Office of Science, Office of High Energy Physics. The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this manuscript, or allow others to do so, for United States Government purposes. 16 It should be noted that although we focused on S-wave combinations of quarks, the channels we study also exclude certain combinations of orbital angular momentum from producing a bound tetraquark. For example, the 0 ++ overlaps with 2Υ in an orbital angular momentum D-wave configuration. If this state produced a low-lying bound tetraquark it would also show up in our calculation. 17 The unaveraged correlator data is too large to be hosted for free on a remote server. all other back-to-back states separated by multiples of this value. Consequently, due to the bottomonium mass being large compared to the smallest allowed momentum, adjacent back-to-back states are sufficiently close in energy that fitting the momentum states as a discrete sum would require a vast set of correlators projected onto each separate |k| 2 /2µ r (with the methodology used in [38]). Practically, this would be overly computationally expensive and instead, the fact that the states with k = 0 are related by the dispersion relation (and are not independent as the sum would assume) should be included. This can be achieved by first expanding the nonperturbative coefficient Z X 2 (k) as a polynomial in |k| 2 /µ 2 r , as dictated by rotational symmetry and by ensuring the Taylor coefficients have the same dimension, then keeping all terms needed to a certain precision. After this the correlator can be written as where going from the first to the second line we have replaced the finite sum by an integral. Taking the limits of the integral to ±∞ and performing the integrals over k analytically yields the fit function given in Eq. (7). Once it is shown that it is possible to replace the finite sum by the indefinite integral within our statistical precision then it is valid to use the above fit function with our data. To do so, using spherical coordinates in Eq. (A6), we define the quantities that we need to compare as The integrands of both are shown diagrammatically in Figure 18, where it is observed that due to the Gaussian time-dependence the peaks of the integrand move towards the origin with larger t. As such, one objective is to choose a large enought such that a sufficient majority of the integrand is within the maximum momentum π/a. We can replace the discrete finite-volume fit function with it's infinite-volume counterpart if the relative difference between them is less than our statistical errors.
where l max is the maximum number of moments to be included in the fit function, the inequality in the second line holds as the moment integrands are positive (shown diagrammatically in Figure 18), in the third line the Cauchy inequality has been used, and in the fourth line it is assumed that the leading moment gives the dominant contribution (Z 2,(l) ≤ Z 2,(0) ). Studying Eq. (A11) instead of Eq. (A9) is a conservative option. Each part of the first term in Eq. (A11) represents how similar I (l) (t) and D (l) (t) need to be in order to be considered equivalent within statistical precision. This is shown in Figure 19. For a particular l max , the second term represents when the higher moments look like noise within statistical precision, also shown in Figure 19. Each Figure was generated with the coarse ensemble parameters (listed as Set 1 in Table III) as this ensemble has the largest lattice spacing (and hence smallest π/a value − the upper limit on the integral of interest) and also the smallest N s (the number of discrete momenta used in the finite-volume sum). As such, the other ensembles will give better approximations and studying Set 1 is conservative. Overlaid on each plot is the smallest relative statistical error from the data on any ensemble. Due to the constant signal-to-noise ratio, the number of configurations and the size of the lattice spacing, the smallest FIG. 19. The difference between the discrete finite-volume and infinite-volume continuum moments (upper) and which moments can be neglected compared to our statistical precision (lower) as discussed in the text. (color online) statistical error was the 2η b correlator on the fine ensemble. Only examining situations below this curve is the most conservative option for all data generated. As can be observed in Figure 19, the discrete finite-volume sums are well represented by the indefinite integrals. Additionally, in order to neglect the higher moments within our statistical precision, a choice oft = 1 fm and l max = 2 is sufficient.
Expanding the finite-volume two-particle energy nonrelativistically in Eq. (A4) neglected a possible finitevolume energy shift which depends on the infinite-volume scattering amplitude. In the small scattering-length limit, the energy shift is known to be volume suppressed [40]. The two-particle systems under study are in this limit as the η b and Υ are compact due to the heavyquark mass, with a size of 0.2 − 0.3fm. As such, the low-momentum energy shifts are not expected to be appreciable given the large volume ensembles we employ. Energy shifts in higher momentum states from Eq. (A4) are exponentially suppressed due to the Gaussian integral in Eq. (A6). Consequently, these too are not appreciable and no large influence of finite-volume energy shifts are seen (c.f., the effective mass figures in Sec. IV). Quantifying these finite-volume scattering shifts is outside the remit of this study. Regardless, the scattering shifts would be positive and push the finite-volume two-particle energy higher and not contribute to a misidentification of a bound tetraquark below the non-interacting threshold.
Finally, the equal mass two-particle propagator starting at a common origin x 0 = (0, 0) and ending at a time t with P tot = 0, in the presence of an external harmonic oscillator potential, is found from using Eq. (B7) in Eq. (B2), to givẽ ∆(t) = mω 2π 1 sinh(2ωt) 3/2 e −2mt . (B9) By comparing (B9) to (B3), and noting that m = 2µ r , we see that the free two-meson harmonic oscillator propagator reduces to the non harmonic oscillator case as ω → 0.