The Search for Beauty-fully Bound Tetraquarks Using Lattice Non-Relativistic QCD
FFERMILAB-PUB-17-381-T
The Search for Beauty-fully Bound Tetraquarks UsingLattice Non-Relativistic QCD
Ciaran Hughes,
1, a
Estia Eichten,
1, b and Christine T. H. Davies
2, c Fermi National Accelerator Laboratory, Batavia, Illinois, 60510, USA SUPA, School of Physics and Astronomy, University of Glasgow, Glasgow, G12 8QQ, UK
Motivated by multiple phenomenological considerations, we perform the first search for the exis-tence of a ¯ b ¯ bbb 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 ¯ b ¯ bbb operators in the three 0 ++ , 1 + − and 2 ++ channels. We em-ploy four gluon field ensembles at multiple lattice spacing values ranging from a = 0 . − .
12 fm,all of which include u , d , s and c quarks in the sea, and one ensemble which has physical light-quarkmasses. Additionally, we perform novel exploratory work with the objective of highlighting anysignal of a near threshold tetraquark, if it existed, by adding an auxiliary potential into the QCDinteractions. With our results we find no evidence of a QCD bound tetraquark below the lowestnon-interacting thresholds in the channels studied. PACS numbers: 12.38.Gc, 13.20.Gd, 13.40.Hq, 14.40.Pq
I. INTRODUCTION
Tetraquarks were first considered theoretically decadesago in the context of light-quark physics in order to ex-plain, amongst other experimental features, the a (980)and f (980) broad resonances [1]. More recently, therehas been exciting experimental evidence indicating thepotential existence of tetraquark candidates amongstthe so-called XYZ states - states whose behaviour dif-fers from predictions of the heavy quark-antiquark po-tential model. The observed XYZ states apparentlycontain two heavy quarks, ( c ¯ c ) or ( b ¯ b ), and two lightquarks [4]. The dynamics of these systems involves boththe short distance and long distance behaviour of QCDand hence theoretical predictions are difficult. Conse-quently, many competing phenomenological models cur-rently exist for these states [5]. Lattice QCD studies ofthe observed XYZ states are also difficult because thesestates are high up in the spectrum as well as being inthe threshold region for strong decays into two heavyflavour mesons. While there are theoretical argumentsthat some tetraquark states with doubly heavy flavor(e.g., bb ¯ u ¯ d , bb ¯ u ¯ s and bb ¯ d ¯ s ) should be bound and stableagainst all strong decays [6], no general arguments existfor tetraquarks with heavy quark-antiquark content suchas Q ¯ Q (cid:48) q ¯ q (cid:48) states.A tetraquark system composed of four heavy quarksis a much cleaner system to study theoretically as long-distance effects from light-quarks are expected not to beappreciable, as opposed to systems which are a mixture ofheavy and light quarks. In the limit of very heavy quarks a [email protected] b [email protected] c [email protected] Recent lattice studies of the scattering amplitude pole do indicatethat 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 makesthese systems particularly useful to study in order to shedlight on the aforementioned XYZ states. In fact, thereis a multitude of phenomenological models (with a quarkmass ranging from the bottom to the very heavy limit)which predict the existence of a ¯ Q ¯ QQQ bound tetraquark[8–15]. However, these are not calculations from first-principles and have an unquantifiable systematic errorassociated with the choice of four-body potential. In re-ality, the heaviest possible tetraquark system in naturewould be a ¯ b ¯ bbb tetraquark. For this, non-perturbativeQCD cannot be ignored, making a first-principles latticeQCD study essential. If such a bound ¯ b ¯ bbb tetraquark didexist, how it would be observed at the LHC has alreadybeen addressed [16, 17].Given these pressing theoretical motivations, in thiswork we perform the first lattice QCD study of the¯ b ¯ bbb system. The sole objective of this exploratorywork is to determine if the dynamics of QCD gener-ates enough binding force between the ¯ b ¯ bbb to producea tetraquark state with a mass below the lowest non-interacting bottomonium-pair threshold, ensuring it isstable against simple strong decays. Searching for sucha bound ¯ b ¯ bbb tetraquark candidate is particularly well-suited to the first-principles lattice QCD methodologybecause this state, if it existed, would be the ground-state in the ¯ b ¯ bbb system. This means it should be rela-tively easy to identify. Further, ¯ b ¯ bbb annihilation effectsare strongly suppressed by the heavy quark mass, as inthe bottomonium system, and so can be ignored.This paper is organised as follows: in Section II theinterpolating operators used in this study are discussed,in Section III the computational methodology is given,while the majority of the results are presented in SectionIV. In Section V we explore a novel method of adding anauxiliary potential into QCD with the objective of high-lighting a possible tetraquark signal. We then discuss ourconclusions in Section VI. a r X i v : . [ h e p - l a t ] J u l TABLE I. The colour representations of the different quarkcombinations. Note that, as described in the text, once thecolour representation of the (anti-) diquark is chosen, thePauli-exclusion principle enforces certain spin combinations in S -wave. Also given are the SU (3) colour contractions neededfor the ¯ b ¯ bbb operators. b ¯ b ¯ bb ¯ b ¯ b bb Colour Irrep 3 c ¯3 c c , c c , ¯6 c ¯3 c , c G efg G ef (cid:48) g (cid:48) δ fg δ f (cid:48) g (cid:48) G efg G ef (cid:48) g (cid:48) δ fg (cid:48) δ f (cid:48) g − δ fg δ f (cid:48) g (cid:48) / G efg G ef (cid:48) g (cid:48) ( δ ff (cid:48) δ gg (cid:48) − δ fg (cid:48) δ gf (cid:48) ) / G efg G ef (cid:48) g (cid:48) ( δ ff (cid:48) δ gg (cid:48) + δ fg (cid:48) δ gf (cid:48) ) / II. INFINITE-VOLUME CONTINUUMEIGENSTATES, OPERATORS AND TWO-POINTCORRELATORS
The QCD Fock space contains all colour-singlet single-particle states such as the conventional mesons | η b ( k ) (cid:105) , | Υ( k ) (cid:105) etc. and, if a ¯ b ¯ bbb bound state also exists, atetraquark state | T b ( k ) (cid:105) . In addition, there are also thetwo-particle states which can be labelled by appropriatequantum numbers as | P tot , J P C ; | k rel | , J P C , J P C , L rel (cid:105) where P tot ( J P C ) is the total (angular) momentum of thetwo-particle system, with J P i C i i the quantum numbers ofthe individual particles and k rel ( L rel ) the relative (orbitalangular) momentum between the two particles.The sole motivation of this work is to search for a pos-sible ¯ b ¯ bbb tetraquark candidate within QCD that couplesto a bottomonium-pair and which lies below the lowestthreshold. The bottomonium mesons we study can beclassified as J P C = S +1 L J . As any orbital angular mo-menta is expected to raise the internal kinetic energy ofthe state (and hence its rest mass) we focus on two-body S -wave systems ( L = 0) with no orbital angular momen-tum between them ( L rel = 0).With the quantum numbers of the Υ /η b being J P C =1 −− / − + , the S -wave 2 η b and 2Υ can have (through theaddition 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 basisof S -wave colour/spin interpolating operators that hasoverlap with these quantum numbers. To do so, we startby forming all possible colour combinations that the 2¯ b and 2 b can be in. These are specified in Table I.We can construct meson interpolating operators as O M ( t, x ) = G efg ¯ b f Γ M b g ( t, x ) (1)where Γ M = iγ , γ k projects onto the quantum numbersof the η b and Υ respectively, and G efg is the colour pro-jection onto the singlet (octet). In addition, it is also possible to construct a (anti-) diquark operator as O ¯3(6) D ( t, x ) = G ¯3(6) efg ¯ b ˆ Cf Γ D b g ( t, x ) (2) O A ( t, x ) = G efg ¯ b f Γ A b ˆ Cg ( t, x ) (3)where ( b ˆ C ) α = C αβ ¯ b β is the charge-conjugated field with C = − iγ γ . As the two quarks have the same flavour,the Pauli-exclusion principle applies and the wavefunc-tion has to be completely anti-symmetric. With ourchoice to focus on S -wave combinations of particles, thespatial wave-function must be symmetric. As the colour(triplet) sextet has a (anti-) symmetric colour wavefunc-tion, this forces the spin-wavefunction to be in a (triplet)singlet with (Γ = γ k ) Γ = iγ .With these building blocks, we can form four classesof ¯ b ¯ bbb colour-singlets by contracting the colour factors G in any irreducible representation (irrep) with its con-jugate colour factors, i.e., 1 c × c , 8 c × c , 3 c × ¯3 c and6 c × ¯6 c . These SU (3) invariant colour contractions aregiven in Table I. After doing this, we need to project theoperators onto a specific angular momentum J P by usingthe standard SO (3) Clebsch-Gordan coefficients (using aspherical basis of spin-matrices [18]) as O J,m ( P,Q ) ( t, x ) = (cid:88) m ,m (cid:104) J, m | J , m , J , m (cid:105)× O J ,m P ( t, x + r ) O J ,m Q ( t, x ) (4)with ( P, Q ) describing the blocks this configuration isbuilt from, i.e, ( η b , η b ), (Υ , Υ), (
D, A ), etc. We also allowthe possibility of the two blocks being separated by a dis-tance r . For the r = case, the operators project onto adefinite total angular momentum J . For the r (cid:54) = case,one can Taylor expand O J ,m P ( t, x + r ) around r = to notice that the operator projects onto a superposi-tion of quantum numbers. Consequently, it is possible toutilise this to further search for the lowest ground stateof the four quark system. When dealing with the diquarkcomponents, to project onto a definite value of charge-conjugation in Eq. (4) one can form the linear combina-tion O ,m D O ,m A ± O ,m D O ,m A .In fact, not all of these colour combinations are inde-pendent. Fierz relations constrain the number of inde-pendent colour-spin operators that are possible. For thelocal operators in S -wave, the relations between the two-meson and diquark-antidiquark states are given in TableII.The simplest quantity that can be calculated on thelattice in order to extract particle masses is the Euclideantwo-point correlator. This is defined as C J PC i,j ( t, P tot = ) = (cid:90) d x (cid:104)O J,m i i ( t, x ) O J,m j j (0 , ) † (cid:105) (5) γ = diag(1 , −
1) in the convention used by NRQCD. Since the 8 c irrep of SU (3) is the adjoint, it is similar to itsconjugate ¯8 c . TABLE II. Fierz relations in the ¯ b ¯ bbb system relating the two-meson and the diquark-antidiquark bilinears. J PC Diquark-AntiDiquark Two-Meson0 ++ ¯3 c × c − |
0; ΥΥ (cid:105) + √ | η b η b (cid:105) ++ c × ¯6 c √ |
0; ΥΥ (cid:105) + | η b η b (cid:105) + − ¯3 c × c √ ( |
1; Υ η b (cid:105) + | η b Υ (cid:105) )2 ++ ¯3 c × c |
2; ΥΥ (cid:105) where we choose to project to zero spatial-momentumand i, j label potentially different operators at the sourceand sink with the same J P C , e.g., i = ( η b , η b ), j = (Υ , Υ).The single-particle contributions to the correlator are de-termined by inserting a complete set of single-particlestates in the Hilbert-space formalism into Eq. (5) to yield C J PC i,j ( t, P tot = 0) = (cid:88) n Z in Z j, ∗ n e − E n t (6)with Z in = (cid:104) |O J,m i i | n (cid:105) the non-perturbative overlap ofthe operator to the eigenstate | n (cid:105) and E n | n (cid:105) = H | n (cid:105) theenergy eigenvalue. Note that all states | n (cid:105) with the samequantum numbers contribute to this correlator, e.g., forthe bottomonium 0 − + pseudoscalar correlator the | η b (cid:105) aswell as all radial excitations contribute. The two-particle contributions to the correlator areslightly more complicated. In this case, as derived inAppendix A, the non-relativistic two-particle states givea contribution to the correlator that is C J PC i,j ( t, P tot = 0) = (cid:16) µ r πt (cid:17) (cid:88) X e − ( M S + M S ) t × (cid:26) Z X + Z X tµ r ) + Z X tµ r ) + · · · (cid:27) (7)where the sum is over all distinct two-particle states X with quantum numbers J P C and P tot = 0, M Si ( M Ki ) isthe static (kinetic) mass of the particles − as defined inEq. (14) − , µ r = M K M K / ( M K + M K ) is the reducedmass and Z lX are non-perturbative coefficients.Energies of states can be extracted using the abovefunctional form once the correlator has been computed.Examining Eq. (5) in the path-integral formalism we canperform the connected Wick contractions so that thecorrelator can be written as an integral over the gluon-fields with the integrand consisting of products of b -quarkpropagators. For each two-meson type operator, e.g., O c η b O c η b , as all quarks have the same flavour there arefour connected Wick contractions. These are shown dia-grammatically in Figure 1.The first Wick contraction for the two-meson corre-lator, called Direct1 and shown in Figure 1a, has theexpression G RAde G RAd (cid:48) e (cid:48) G RBgh G RBg (cid:48) h (cid:48) (cid:88) x Tr (cid:104) Γ M K − ( t, x ; 0 , z ) e (cid:48) g Γ † M K − ( t, x ; 0 , z ) † hd (cid:48) (cid:105) Tr (cid:104) Γ M K − ( t, x (cid:48) ; 0 , z (cid:48) ) eg (cid:48) Γ † M K − ( t, x (cid:48) ; 0 , z (cid:48) ) † h (cid:48) d (cid:105) (8)while the second, called Xchange2, is given by −G RAde G RAd (cid:48) e (cid:48) G RBgh G RBg (cid:48) h (cid:48) (cid:88) x Tr (cid:104) Γ M K − ( t, x ; 0 , z ) e (cid:48) g Γ † M K − ( t, x (cid:48) ; 0 , z ) † hd Γ M K − ( t, x (cid:48) ; 0 , z (cid:48) ) eg (cid:48) Γ † M K − ( t, x ; 0 , z (cid:48) ) † h (cid:48) d (cid:48) (cid:105) (9)where x (cid:48) = x + r . The other diagrams, Direct3 andXchange4, have similar expressions. For the diquark-antidiquark type operators, e.g., O c D O ¯6 c A , there are also four Wick contractions which can be combined into oneexpression as C J PC ( t, P tot = 0) = (cid:2) ± sgn ( C Γ D ) T ± sgn ( C Γ A ) T + sgn ( C Γ D ) T sgn ( C Γ A ) T (cid:3) G RAde G RAd (cid:48) e (cid:48) G RBgh G RBg (cid:48) h (cid:48) × (cid:88) x Tr (cid:104) C Γ D K − ( t, x ; 0 , z ) eg (cid:48) Γ † D CK − ( t, x ; 0 , z ) Tdh (cid:48) (cid:105) Tr (cid:104) Γ A CK − ( t, x (cid:48) ; 0 , z (cid:48) ) ∗ e (cid:48) g C Γ † A K − ( t, x (cid:48) ; 0 , z (cid:48) ) † hd (cid:48) (cid:105) (10) Annihilation diagrams are suppressed by powers of the heavy (a) Direct1 (b) Xchange2(c) Direct3 (d) Xchange4
FIG. 1. There are four connected Wick contractions for the two-meson type correlator when the quarks have the same flavour.The grey region represents a colour neutral meson, the blue line a quark and the red line an antiquark. We call these the (a)Direct1 contraction where each meson propagates to itself. (b) Xchange2 where an anti-quark is exchanged between the mesonpair, (c) Direct3 where each meson propagates to the other, and (d) Xchange4 where a quark is exchanged between the mesonpair. (colour online) where the sign function is defined by sgn ( X ) T = ± X T = ± X . The ± in Eq. (10) corresponds to the 3 or 6colour representation. It is this prefactor with the signswhich enforces the Pauli exclusion principal: the sumcancels for spin combinations that do not make the wave-function overall anti-symmetric. The spin-triplet/singletconfigurations we consider here obey ( Cγ k ) T = +( Cγ k )and ( Cγ ) T = − ( Cγ ). Diagrammatically the four con-nected Wick contractions contributing to the diquarkcorrelator are shown in Figure 2.To calculate the two-point correlators described abovewithin the first-principles Feynman path-integral ap-proach 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 METHODOLOGYA. Second Generation N f = 2 + 1 + 1 GluonEnsembles
Our lattice calculation uses gauge field configurationsgenerated by the MILC collaboration [20]. For the gaugefields, they used the tadpole-improved L¨uscher-Weiszgauge action correct to O ( α s a ) [21] and include 2+1+1flavours in the sea, the up and down quarks (treated astwo degenerate light quarks with mass m l ), the strangequark, and the charm quark. The sea-quarks are includedusing the Highly Improved Staggered Quark formulation[22].Four ensembles are chosen in this study. As onemight expect roughly double the discretisation errors ina 2¯ bb system relative to the ¯ bb one, to ensure that theheavy quark potential is accurately represented (whereshort distance details may be important for a com-pact tetraquark candidate) we utilise three ensembles (a) (b)(c) (d) FIG. 2. There are four connected Wick contractions for the diquark-antidiquark type correlator when the quarks have the sameflavour. The blue shaded region represents a diquark, the red shaded region the antidiquark, a blue line a quark and the redline an antiquark. The uncrossing of the lines in Figure 2b to produce Figure 2a gives a ± , as discussed in the text, whichenforces the Pauli-exclusion principle. (colour online) that span relatively fine lattice spacings ranging from a = 0 . − .
12 fm. Details of the ensembles are given inTable III. Due to the computational expense, most of theensembles use heavier m l than in the real world. How-ever, to test m l dependence, we use one ensemble (Set 2in Table III) that has physical am l /am s . Additionally,the ensembles have been fixed to Coulomb gauge to al-low 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 appropri-ate for physical systems where the relative velocity ofthe constituent particles inside the bound state is muchsmaller than one (in Planck units). When applied tothe strong force, this framework is called non-relativisticQCD (NRQCD). It is well known that b -quarks are verynonrelativistic inside their low-lying bound states, where v ≈ . v . Here we use an action accu-rate through O ( v ) with additional spin-dependent termsat O ( v ). Operators are also added to correct for dis-cretisation effects. We make the further systematic im-provement here, introduced in [24], to include coefficientsof O ( v ) operators that have been matched to continuumQCD through O ( α s ). We call this improved-NRQCD(iNRQCD). This action has already been used to suc-cessfully determine bottomonium S , P and D wave masssplittings [24, 25], precise hyperfine splittings [26, 27], B meson decay constants [28], Υ and Υ (cid:48) leptonic widths[29], B , D meson mass splittings [27] and hindered M1radiative decays between bottomonium states [30].The iNRQCD Hamiltonian evolution equations can be The spin-independent O ( v ) terms are subleading effects for the¯ b ¯ bbb energies relevant to this study. TABLE III. Details of the gauge ensembles used in this study. β is the gauge coupling. a (fm) is the lattice spacing [24, 31], am q are the sea quark masses, N s × N T gives the spatial andtemporal extent of the lattices in lattice units and n cfg is thenumber of configurations used for each ensemble. We use16 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 am s am c N s × N T n cfg .
00 0 . . . .
635 24 ×
64 10522 6 .
00 0 . . . .
628 48 ×
64 10003 6 .
30 0 . . .
037 0 .
440 32 ×
96 10084 6 .
72 0 . . .
024 0 .
286 48 ×
144 400 written as G ( x , t + 1) = e − aH G ( x , t ) G ( x , t src ) = δ x , x (11)with e − aH = (cid:18) − aδH | t +1 (cid:19) (cid:18) − aH | t +1 n (cid:19) n U † t ( x ) × (cid:18) − aH | t n (cid:19) n (cid:18) − aδH | t (cid:19) (12) aH = − ∆ (2) am b ,aδH = aδH v + aδH v ; aδH v = − c (∆ (2) ) am b ) + c i am b ) (cid:16) ∇ · ˜E − ˜E · ∇ (cid:17) − c am b ) σ · (cid:16) ˜ ∇ × ˜E − ˜E × ˜ ∇ (cid:17) − c am b σ · ˜B + c ∆ (4) am b − c (∆ (2) ) n ( am b ) aδH v = − c am b ) (cid:110) ∆ (2) , σ · ˜B (cid:111) − c am b ) (cid:110) ∆ (2) , σ · (cid:16) ˜ ∇ × ˜E − ˜E × ˜ ∇ (cid:17)(cid:111) − c i am b ) σ · ˜E × ˜E . (13)The parameter n is used to prevent instabilities at largemomentum from the kinetic energy operator. A valueof n = 4 is chosen for all am b values. We evaluate thepropagator using local sources Eq. (11). Here, am b isthe bare b quark mass, ∇ is the symmetric lattice deriva-tive, with ˜ ∇ the improved version, and ∆ (2) , ∆ (4) arethe lattice discretisations of Σ i D i , Σ i D i respectively. ˜E , TABLE IV. Parameters used for the valence quarks. am b is the bare b -quark mass in lattice units, u L is the tadpoleparameter and the c i are coefficients of terms in the NRQCDHamiltonian (see Eq. 13). Details of their calculation can befound in [24, 32]. c , c , c and c are included at tree-level.Set am b u L c , c c c c .
73 0 . .
31 1 .
02 1 .
19 1 .
162 2 .
66 0 . .
31 1 .
02 1 .
19 1 .
163 1 .
95 0 . .
21 1 .
29 1 .
18 1 .
124 1 .
22 0 . .
15 1 .
00 1 .
12 1 . ˜B are the improved chromoelectric and chromomagneticfields, details of which can be found in [24]. Each ofthese fields, as well as the covariant derivatives, must betadpole-improved. We take the mean trace of the gluonfield in Landau gauge, u L = (cid:104) Tr U µ ( x ) (cid:105) , as the tadpoleparameter, calculated in [24, 28]. The matching coeffi-cients c , c , c , c , c in the above Hamiltonian havebeen computed perturbatively to one-loop [24, 32]. c was found to be close to the tree-level value of one non-perturbatively [24] and we set it, as well as the rest ofthe matching coefficients, to their tree-level values. Thequark 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 termsin the action). The above Hamiltonian neglects the four-fermion operators which appear at O ( α s v ), as well asother operators which appear at higher order in the non-relativistic expansion. Simple power-counting estimates[24, 26] would lead us to expect contributions of order afew percent (or a few MeV) at most to binding energiesfrom these terms. In the case where a tetraquark boundstate is observed then we can estimate the systematiceffect from neglecting these contributions.The parameters used in this study are summarised inTable IV. There, c , c and c are the correct values for an O ( v ) iNRQCD action [24]. For Set 4, all parameters arethose for the O ( v ) action. The small changes to thesecoefficients in going to an O ( v ) iNRQCD action (whichare similar in magnitude to the two-loop corrections) arenot appreciable for the purpose of this work: whetheror not a tetraquark candidate exists below the lowestbottomonium-pair threshold. All other parameters listedin Table IV are taken from [26, 30].Within iNRQCD the single-particle energy-eigenstatescan be decomposed in the standard non-relativistic ex-pansion as E ( P ) = M S + | P | M K + . . . (14)with M S , M K the static and kinetic masses respectively[24]. Because the quark mass term is removed fromthe iNRQCD Hamiltonian, the static mass is unphysi-cal, differing by a constant shift from the physical (ki-netic) mass. This means that only static mass differ-ences determined fully non-perturbatively can be com-pared to experimental results. However, this constantshift can be calculated in lattice perturbation theory ifrequired [33], or by using an additional experimental in-put. The kinetic mass does not suffer this problem asit acquires the quark mass contributions from the quarkkinetic 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 K − ( x | y ) = (cid:32) G ψ ( x | y ) 00 − G χ ( x | y ) (cid:33) (15)where G ψ ( x | y ) is the two-spinor component quark prop-agator and G χ ( x | y ) is the two-spinor component anti-quark propagator. Taken together, we can now computethe b -quark propagator via the iNRQCD evolution equa-tions on the gluon ensembles listed in Table III. The lastpiece needed to calculate the two-point correlators, andhence energies, are the discretised finite-volume versionsof the interpolating operators. C. Discrete Finite-Volume Operators
Together, the isotropic discretisation and the peri-odic 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. IIhave well-defined J P C quantum numbers associated with SO (3), we need to construct operators which transformwithin the irreps of the O h symmetry group (relevantfor lattice calculation). This can be achieved by themethod of subduced representations, where an operatorwith a specific J P C can be taken to a specific lattice ir-rep Λ P C by using the subduction coefficients found inAppendix A of [18]. At rest, each of our J P C = 0 ++ and J P C = 1 + − operators trivially subduce into a sin-gle lattice irrep labelled by A ++1 and T + − respectively.However, the J P C = 2 ++ case is slightly more compli-cated and subduces into two lattice irreps labelled by T ++2 and E ++ (which are three- and two-dimensional).We construct both the T /E operators as √ O [2] T /E = O J =2 ,m =2 ± O J =2 ,m = − , which are correctly subducedfrom the J = 2 ++ operators defined in Sec. II. In prin-ciple, each lattice irrep allows mixing between different J states, e.g., the A irrep contains not only the J = 0states but also the J = 4 [34]. However in practice since The conserved quantum numbers of a symmetry group are de-termined using the little group, which for SO (3) and O h dependon the momenta type [35]. In this study we focus on states atrest. TABLE V. The ¯ b ¯ bbb interpolating operators used in thisstudy. Operators in each column are subduced from theinfinite-volume continuum quantum numbers J PC given inthe first row. The superscript on each operator denotes thelattice irrep of that operator and the subscript denotes thebuilding blocks of the operator, as explained in the text. Wegenerate each operator with three different spatial configu-rations as shown in Eq. (4): where the building blocks areseparated by a distance r x = 0 , x -direction. 0 ++ + − ++ source sink source/sink source/sink O A ( η b ,η b ) O A ( η b ,η b ) O T (Υ ,η b ) O T (Υ , Υ) O A ( η b ,η b ) O A (Υ , Υ) O T ( D ¯3 c ,A c ) O E (Υ , Υ) O A (Υ , Υ) O A ( η b ,η b ) O T ( D ¯3 c ,A c ) O A (Υ , Υ) O A (Υ , Υ) O E ( D ¯3 c ,A c ) O A ( D ¯3 c ,A c ) O A ( D ¯3 c ,A c ) O A ( D c ,A ¯6 c ) O A ( D c ,A ¯6 c ) we are only looking for the ground state of the ¯ b ¯ bbb cor-relators we are not sensitive to these mixing effects.A complete list of ¯ b ¯ bbb interpolating operators usedto produce the correlator data herein is given in TableV. In fact, this is an over-constrained set due to theFierz identities (shown in Table II) which relate the two-meson and diquark-antidiquark correlators. We includethis over constrained system and ensure that we repro-duce the Fierz relations to numerical precision, perform-ing a non-trivial check on our data. Additionally, wealso reproduce the relations between the 8 c × c colourcombination and the others [36] on a subset of the data.It has been found [37] that separating each hadronwithin the two-hadron interpolating operator by a spe-cific distance r can significantly aid in the extractionof the (ground) state energy. In this direction, we usethree different spatial configurations of the ¯ b ¯ bbb correla-tors where the individual building blocks are separatedby a distance of r x = 0 , x -direction .Finally, the subduced lattice interpolating operatorsare defined in terms of the Dirac fields as in Eq. (1), andthe correlators are defined analogously, as in Eq. (10).We can then use the decomposition of K − ( x | y ) givenin Eq. (15) with suitable boundary conditions to writethe correlator in terms of the iNRQCD quark propagator G ψ ( x | y ). Due to our use of iNRQCD, there are no back-ward propagating valence anti-quarks in our calculation. The η b and Υ energies used to determine the non-interacting 2 η b and 2Υ thresholds (needed to determine if a state exists belowthem) are found from locally smeared meson correlators only. Consequently, the appreciable finite-temporal effects seenin relativistic two-meson lattice correlators [38] do notarise in our calculation, simplifying the analysis. Withthis methodology, it is now possible to compute the low-est energy levels associated with the ¯ b ¯ bbb system. IV. THE LOW-LYING ENERGY EIGENSTATESOF THE ++ , + − AND ++ ¯ b ¯ bbb SYSTEM
In order to determine if there is an energy eigenstatebelow the 2 η b threshold, we first need to find the non-interacting thresholds on each ensemble listed in TableIII. This is achieved by computing the bottomonium η b and Υ two-point correlators as described above, thenfitting them to the functional form given in Eq. (6) toextract the single particle energies. As in the range ofstudies listed in Sec. III B, amongst others, we utilise thewell-established Bayesian fitting methodology [39] in thiswork and refer the reader to [24] for technical details.Although we fit the correlator data in order to extractfit parameters, in the following we display effective massplots so the reader can visualise the data. The singleparticle effective mass is constructed as aE eff J PC = log (cid:32) C J PC i,j ( t ) C J PC i,j ( t + 1) (cid:33) (16)= aE J PC + Z i Z j, ∗ Z i Z j, ∗ e − ( E − E ) t (1 − e − ( E − E ) ) + . . . (17) t →∞ −−−→ aE J PC . (18)As can be inferred from Eq. (17), the greater the massgap E − E or the larger the ground state overlap Z then the quicker aE eff J PC converges to a plateau, whichgives aE J PC . The η b and Υ effective masses are shown inFigure 3, where the returned fitted ground-state energyfrom the correlator fits is also shown overlaid in black.The large difference in energy between the η (cid:48) b (Υ (cid:48) ) andthe ground state η b (Υ) means that the effective massplots fall rapidly to a plateau given by the ground stateenergy, indicating that the fit to the correlator data willextract the ground state energy precisely.Also evident from the effective mass plot is the con-stant signal-to-noise in the η b data, as might be ex-pected from straightforward application of the Parisi-Lepage arguments [40, 41] for noise growth in a sys-tem where all the quarks are the same and no 0 ++ bound tetraquark exists. The argument specifies thatthe noise of the two-point correlator should behave like These arise in lattice two-meson calculations when a relativisticformulation of valence quark is used due to one of the mesonspropagating forward in time while the other propagates aroundthe temporal boundary backwards in time [38]. t / a . . . . . . a E e ff ( t ) t ( fm ) a . N s ϒη b E fit ϒ E fit η b FIG. 3. The effective mass plot for the η b and Υ on the“superfine” ensemble (Set 4 listed in Table III). The effectivemass plots on the other ensembles are qualitatively identical.(colour online) exp { ( E J PC − E GS / t } where E J PC is the lowest en-ergy eigenstate of the bottomonium operator O J PC con-structed to have the quantum numbers J P C and E GS isthe lowest energy eigenstate of the mean squared correla-tor which controls the noise. Thus, it would be surprisingif a tetraquark candidate did exist below the 2 η b thresh-old from the lattice perspective alone as then E GS < E η b and the noise of the η b data would grow exponentially.However, the lattice calculation still needs to be per-formed for a conclusive statement to be made about theexistence of this tetraquark candidate since the Parisi-Lepage arguments do not allow for raw crossed Wickcontractions that would contribute to either the full two-meson or tetraquark correlator.Lattice correlators are affected by both the discrete na-ture of the space-time lattice and separately by its finitevolume. Each has a separate but calculable effect on theextracted lattice energies. Corrections in energies dueto discretisation effects are proportional to a k , where k depends on the level of improvement. Here systematicdiscretisation errors are reduced to α s a by the improve-ments made, as for those studies listed in Sec. III B, andwe expect this to be small enough to have little impact.We can assess this from our results with different valuesof the lattice spacing.Finite-volume effects for single-particle energies (aris-ing from the lightest particle in the sea propagatingaround the spatial boundaries) are known to behave likeexp( − aM π N s ) [42] and are not appreciable for the en-sembles used here which have aM π N s ≈ N s on Set 1 and 2 differ by a factor of two, giving a basictest of volume-dependence. However, finite-volume inter-actions can shift a two-particle energy by an amount thatdepends on the infinite-volume scattering matrix. Fur-ther, these shifts are non-trivial to parameterise (see forexample [44–50]). As the specific purpose of this study isto search for a hypothetical tetraquark bound state be-
20 40 60 80 100 120 140 t / a a E e ff = l og ( C ( t ) / C ( t + )) a . N s Mock Data With Free Mesonshhhh M lat ϒ M lat η b E ++ FIG. 4. Assuming a tetraquark exists 100 MeV below the2 η b threshold, different normally distributed couplings in ¯ b ¯ bbb correlator mock-data produces different effective mass curveson the “superfine” ensemble (Set 4 listed in Table III) as de-scribed in the text. (colour online) low the lowest threshold, we will not attempt to quantifythese finite-volume interactions.Eq. (7) describes the non-relativistic two-particle con-tribution to the correlator after t = 1 fm (as shown inAppendix A). In this case, because of the additional1 /t time-dependence, the effective mass formula forthese contributions differs from their single-particle coun-terparts. Removing the leading order time-dependenceyields an effective mass defined as aE eff ,tJ PC = log (cid:32) t C J PC i,j ( t )( t + 1) C J PC i,j ( t + 1) (cid:33) . (19)For the 0 ++ , 1 + − and 2 ++ operators that are con-structed, if a stable tetraquark exists below the 2 η b , η b +Υor 2Υ thresholds then it would show up as the groundstate of each correlator and hence also in the effectivemasses. Otherwise, each threshold would be the low-est energy eigenstate. Higher energy states will also ap-pear in each correlator. For example, the 2Υ and η b η (cid:48) b in the 0 ++ case, the Υ η (cid:48) b in the 1 + − and the ΥΥ (cid:48) and χ b χ b (1 P ) in the 2 ++ . Of these, if no tetraquark state ispresent, only the 2Υ in the 0 ++ might be noticeable whilestudying the ground state, as it is O (100) MeV abovethe 2 η b . All other excited states have similar energy dif-ferences to those appearing in the η b / Υ effective massesshown in Figure 3, which rapidly falls to the ground state.It is helpful at this stage to generate mock correlatordata to illustrate how we might expect the ¯ b ¯ bbb correlatorresults to behave in the presence or absence of a low-lyingtetraquark state (neglecting two-particle finite-volume ef-fects). Using the extracted lattice η b and Υ masses, wecan compute the non-interacting 2 η b and 2Υ thresholdson our ensembles. Further, for a fixed value of Z X , wecan also compute their leading order two-particle contri-bution to the correlator from Eq. (7). Additionally wecan infer the values of the non-interacting η b η (cid:48) b and ΥΥ (cid:48) masses on our ensembles using the experimental PDG[4] values as input in order to include their two-particlecontributions in the mock correlator data also. Further,in the 0 ++ channel, if we assume a tetraquark bound-state exists 100 MeV below the 2 η b threshold , for afixed value of the tetraquarks non-perturbative overlap, Z b , this hypothetical state’s contribution to the corre-lator is given by Eq. (6). Then, given the effective massformula defined in Eq. (16), for each different choice ofthe non-perturbative coefficients we can generate a sep-arate effective mass curve. In practice, we choose dif-ferent values of the coefficients from a normal distribu-tion with zero mean and unit variance. Figure 4 showssuch a plot for the “superfine” ensemble (Set 4 in Ta-ble III) where the solid blue curves represent differentvalues of the normally distributed coefficients. As thenon-perturbative coefficients show up through the ratio Z X /Z b in the effective mass formula, analogously toEq. (17), once the energy difference E − E is set theeffective mass is only sensitive to the relative size of theof the tetraquark overlap to that of the lowest threshold(once the contribution of excited states has become neg-ligible). With this knowledge, the lower red dotted curvegives mock-data in the situation where there is only a newstate present in the correlator data ( Z X = 0), the middlered dashed curve indicates the case when the tetraquarkand two-meson states have the same value of coefficient( Z X = Z b ), while the upper dot-dashed red curve is themock-data in the case where there is no new state in thedata ( Z b = 0). The blue curves below the middle redone have a larger overlap onto the new state while thoseabove have an increasingly vanishing one. With our over-constrained colour-spin basis of S -wave operators at leastone operator should have an appreciable overlap onto atetraquark state below threshold (if it exists) and, as il-lustrated by the mock-data, give a easy/clean signal sim-ilar to the middle red dashed curve where the effectivemass drops below the 2 η b threshold. Even though Figure4 is mock-data, the upper red curve with no tetraquarkstate present looks very like the real data shown in Fig-ures 5 and 6.One important point to note is the additional factor of1 /t appearing in the two-particle contribution in Eq. (7)relative to the single-particle one in Eq. (6). This factorsuppresses the two-particle contribution relative to thesingle-particle state, e.g., t = 100 gives a suppressionof the two-particle state of (0 . . . This is one rea-son why the middle red dashed curve has a particularlyrapid fall to the new state which lies only 100 MeV belowthe 2 η b (compared to Figure 5d where the 2Υ is O (100)MeV above the 2 η b ). Overall this effect would producean enhancement of the stable tetraquark state if it ex-ists. Further, while we used the “superfine” ensemble for The smallest binding of a below threshold ¯ b ¯ bbb tetraquark fromthe phenomenological studies (which are shown in Figure 17) was108 MeV [9]. t / a . . . . . a E e ff η b → η b ( t ) t ( fm ) a . N s E eff , t ++ ( t ) : r x = E eff0 ++ ( t ) : r x = E eff , t ++ ( t ) : r x = E eff0 ++ ( t ) : r x = E eff , t ++ ( t ) : r x = E eff0 ++ ( t ) : r x = M lat ϒ M lat η b E fit0 ++ (a) 2 η b → η b t / a . . . . . a E e ff η b → ϒ ( t ) t ( fm ) a . N s E eff , t ++ ( t ) : r x = E eff0 ++ ( t ) : r x = E eff , t ++ ( t ) : r x = E eff0 ++ ( t ) : r x = E eff , t ++ ( t ) : r x = E eff0 ++ ( t ) : r x = M lat ϒ M lat η b E fit0 ++ (b) 2 η b → t / a . . . . . a E e ff ϒ → η b ( t ) t ( fm ) a . N s E eff , t ++ ( t ) : r x = E eff0 ++ ( t ) : r x = E eff , t ++ ( t ) : r x = E eff0 ++ ( t ) : r x = E eff , t ++ ( t ) : r x = E eff0 ++ ( t ) : r x = M lat ϒ M lat η b E fit0 ++ (c) 2Υ → η b t / a . . . . . a E e ff ϒ → ϒ ( t ) t ( fm ) a . N s E eff , t ++ ( t ) : r x = E eff0 ++ ( t ) : r x = E eff , t ++ ( t ) : r x = E eff0 ++ ( t ) : r x = E eff , t ++ ( t ) : r x = E eff0 ++ ( t ) : r x = M lat ϒ M lat η b E fit0 ++ (d) 2Υ → FIG. 5. The ¯ b ¯ bbb effective masses for the 0 ++ ( M , M ) → ( M , M ) correlators where M , M are the η b or Υ. E eff and E eff ,t are the single- and two-particle effective masses defined in Eq. (18) and Eq. (19) respectively. The mesons are separated by adistance r x in the x -direction when constructing the two-meson interpolating operator as given in Eq. (4). Gray points are notused when fitting the data (c.f., App. A). (color online)TABLE VI. The ground state static masses extracted from the lattice ¯ b ¯ bbb correlator data as described in the text.Set aM η b aM Υ aM ++ aM + − aM ++ T aM ++ E . . . . . . . . . . . . . . . . . . . . . . . . illustrative purposes, the other ensembles with larger lat-tice spacings have similar features but with less temporalresolution due to the larger numerical value of the latticespacing.If one is searching for a bound state below thresholdthen both single- and two-particle contributions appear in the correlator and one may use the single particle ef-fective mass formula given in Eq. (18) in order to high-light the bound ground state mass (as done in the mockdata above). However if no bound state exists in thedata, which we do not know a priori , then not remov-ing the 1 /t dependence from the two-particle contri-1butions has an important effect: an additional factor of1 . /t ) is introduced into the effective mass for-mula in Eq. (18), giving a contamination that vanishesslowly as t → ∞ . This would produce a confusing pic-ture of what is actually contributing. The two-particleeffective mass formula Eq. (19) removes this contribution.In the results to be reported now, we will overlay bothsingle- and two-particle effective masses on the same plotfor the reader’s convenience.We generate the ¯ b ¯ bbb correlator data for the operatorsgiven in Table V using the ensembles listed in Table IIIand fit this data simultaneously with the bottomoniummeson data so to include correlations between data sets.All the ¯ b ¯ bbb data within a specific irrep and those whichare unrelated by a Fierz relation are fit using Eq. (7)for the two-particle contributions and Eq. (6) for a hy-pothetical tetraquark state below threshold. The meanof the prior energy of the 2 η b , η b + Υ and 2Υ thresholdsare roughly estimated based on the effective masses andthen given a suitably wide prior width of 100 MeV whilea 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 thresh-old, no energy eigenstate is found below this thresholdand variations of the tetraquark prior energy are insignif-icant. Similar behaviour is seen with the other quantumnumbers.Again, while we fit the correlator data in order to ex-tract particle energies, so that the reader can visualisethis data we display effective mass plots on the differ-ent ensembles. The “superfine” ensemble (Set 4) 0 ++ two-meson effective masses are shown in Figures 5, whilethe 0 ++ diquark-antidiquark are given in Figure 6, the“physical coarse” (Set 2) 1 + − two-meson and diquark-antidiquark in Figure 7 and the “fine” (Set 3) 2 ++ two-meson and diquark-antidiquark correlators subduced intothe T lattice irrep are shown in Figure 8. Each plot hasthe fitted ground state energy overlaid in black for com-parison. The 2 ++ subduced into the E irrep is similar tothe T case. Further, the behaviour of the lattice data onall ensembles is qualitatively similar to those shown. Theextracted ground state energies in each channel are givenin Table VI and a comparison of the energies is shown ingraphical form in Figure 9.It should also be noted that the numerical value ofthe effective mass (shift) plateau in two-hadron correla-tors has been shown, in certain cases, to be sensitive tothe choice of interpolating operators [51, 52]. There, theauthors found that when the noise growth in the cor-relator data restricts the study of effective masses to amaximum propagation time of approximately 2 fm, fakeplateaus can appear. These fake plateaus can be a conse-quence of different choices of source and sink (smeared) Simultaneously fitting data sets related by a Fierz identity wouldmean the correlation matrix would have a zero eigenvalue andthus not be invertable for use in a least-squares minimisation. t / a . . . . . a E e ff D × ¯ → × ¯ ( t ) t ( fm ) a . N s E eff , t ++ ( t ) : r x = E eff0 ++ ( t ) : r x = E eff , t ++ ( t ) : r x = E eff0 ++ ( t ) : r x = E eff , t ++ ( t ) : r x = E eff0 ++ ( t ) : r x = M lat ϒ M lat η b E fit0 ++ t / a . . . . . a E e ff D × ¯ → × ¯ ( t ) t ( fm ) a . N s E eff , t ++ ( t ) : r x = E eff0 ++ ( t ) : r x = E eff , t ++ ( t ) : r x = E eff0 ++ ( t ) : r x = E eff , t ++ ( t ) : r x = E eff0 ++ ( t ) : r x = M lat ϒ M lat η b E fit0 ++ FIG. 6. The ¯ b ¯ bbb effective masses for the 0 ++ diquark-antidiquark ¯3 × ¯3 and 6 × ¯6 correlators. E eff and E eff ,t are thesingle- and two-particle effective masses defined in Eq. (18)and Eq. (19) respectively. The diquarks are separated by adistance r x in the x -direction when constructing the diquark-antidiquark interpolating operator as given in Eq. (4). (coloronline) operators: this can cause the a negative sign in the Z term in the effective mass formula Eq. (17) and a dipbelow threshold can appear for a short time range whichcan be misinterpreted as a bound state. Wall-sourceswere shown to be particularly prone to this behaviour ofproducing a “false dip” and obtaining an appreciably dif-ferent “plateau” than a Gaussian source. Here we onlyuse local quark sources. In addition, the elastic scatteringstates can also have a dependence on the choice of opera-tor, which can cause a slow decay to the ground state andmimic a slowly varying effective mass that can be mis-taken for a plateau over a short time range. As noted inthese studies, a necessary check for a real effective massplateau when using different source and sink operators isthe convergence of all data to a single plateau at timeslarger than approximately 2 fm. As we separate the op-erators by r x = 0 , t > t / a . . . . . . . . a E e ff ϒ η b → ϒ η b ( t ) t ( fm ) a . N s E eff , t + − ( t ) : r x = E eff1 + − ( t ) : r x = E eff , t + − ( t ) : r x = E eff1 + − ( t ) : r x = E eff , t + − ( t ) : r x = E eff1 + − ( t ) : r x = M lat ϒ + M lat η b E fit1 + − t / a . . . . . . . . a E e ff D × ¯ → × ¯ ( t ) t ( fm ) a . N s E eff , t + − ( t ) : r x = E eff1 + − ( t ) : r x = E eff , t + − ( t ) : r x = E eff1 + − ( t ) : r x = M lat ϒ + M lat η b E fit1 + − FIG. 7. The ¯ b ¯ bbb effective masses for the 1 + − Υ η b anddiquark-antidiquark ¯3 × ¯3 correlators. E eff and E eff ,t are thesingle- and two-particle effective masses defined in Eq. (18)and Eq. (19) respectively. (color online) A few notable features of the ¯ b ¯ bbb effective mass plotsare evident. First and foremost, no value of the effec-tive mass is observed below the lowest non-interactingbottomonium-pair threshold in any channel, in line withwhat one would expect if no stable tetraquark candi-date existed below threshold. Indeed, the ¯ b ¯ bbb effectivemass plots are strikingly similar to the upper dot-dashedcurve in the mock data in Figure 4 where no boundtetraquark state is present. Additionally, the 2 η b → η b effective mass shown in Figure 5a plateaus very earlydue to the larger overlap onto the 2 η b threshold, whilethe 0 ++ →
2Υ effective mass shown in Figure 5dfalls more slowly to the 2 η b threshold due to the largeroverlap onto the nearby 2Υ threshold. The cross correla-tors 2 η b →
2Υ show how the different operators convergeto a single plateau at a time greater than t ≈ ++ correlator data is a linearcombination of the 2 η b and 2Υ two-meson data, relatedby the Fierz identities given in Table II, and the effec-tive masses shown in Figure 6 reflects this. It is em-pirically observed that separating the diquark from the t / a . . . . . . . . a E e ff ϒϒ → ϒϒ ( t ) t ( fm ) a . N s E eff , t ++ ( t ) : r x = E eff2 ++ ( t ) : r x = E eff , t ++ ( t ) : r x = E eff2 ++ ( t ) : r x = E eff , t ++ ( t ) : r x = E eff2 ++ ( t ) : r x = M lat ϒ M lat η b E fit2 ++ t / a . . . . . . . . a E e ff D × ¯ → × ¯ ( t ) t ( fm ) a . N s E eff , t ++ ( t ) : r x = E eff2 ++ ( t ) : r x = E eff , t ++ ( t ) : r x = E eff2 ++ ( t ) : r x = M lat ϒ M lat η b E fit2 ++ FIG. 8. The ¯ b ¯ bbb effective masses for the 2 ++ ΥΥ anddiquark-antidiquark ¯3 × ¯3 correlators subduced into the T irrep. E eff and E eff ,t are the single- and two-particle effectivemasses defined in Eq. (18) and Eq. (19) respectively. (coloronline) anti-diquark by too large a distance r x results in largernoise due to the separation of colour sources. The 1 + − Υ η b and diquark-antidiquark effective masses are shownin Figure 7, where the noise starts to increase after t ≈ η b threshold.Based on this, one would also expect the signal-to-noiseto be worse for the 2 ++ data, which is also evident fromthe correlator data subduced into the T irrep shown inFig. 8. As Set 2 has physical m l /m s corresponding toa pion mass of O (131) MeV, while the other ensembleshave nonphysical m l /m s corresponding to pion massesof O (300) MeV [53], no sensitivity to light sea-quarks isobserved. This would be expected from the smallness ofthe Van-der-Waals potential generated by the two-pionexchange between two 1 S bottomonium mesons [54]. Ascan be seen in Figures 8 and 9, the 2 ++ ground stateobtained from the lattice is slightly higher than that ofthe non-interacting threshold. However, this is the statewhich has the largest signal-to-noise, restricting the datato shorter time regions and it is possible that we are sen-3 − − −
20 0 20 40 60 E GS − E th (MeV) a . N s a . N s a . N s a . N s ++ E ++ T + − ++ FIG. 9. A summary of the ¯ b ¯ bbb ground state energies withthe lowest non-interacting bottomonium-pair threshold sub-tracted, across the different lattice ensembles listed in TableIII. Statistical error only. Note, as shown in Table III, fewerconfigurations were used on the a = 0 .
06 fm ensemble thanon the others. (color online) sitive to the same aforementioned issue of a slowly vary-ing fake plateau. Alternatively, this positive shift in thetwo-particle energy could potentially indicate appreciableinfinite-volume continuum scattering arising from finite-volume interactions [44], but quantifying these phase-shifts is outside the remit of this study. Regardless, theseeffects do not indicate that a bound tetraquark state ex-ists in this channel.For illustration purposes, we also show the effectivemasses of the individual Wick contractions contributingto the 2 η b → η b correlator in Figure 10. As is evident, ineach individual Wick contraction the effective mass dropsbelow the 2 η b threshold but then rises slowly to thresh-old. However, importantly, when all Wick contractionsare added together to yield the full correlator (shownFigure 5a) the effective mass falls rapidly to thresholdfrom above. This behaviour will be discussed further inSec. VI.After analysing all our data, as hinted by the effectivemass plots, there is no indication of a bound tetraquarkstate below the non-interacting thresholds on any ensem-ble, as shown in Fig. 9. We see no evidence of any changein the ground-state energies (with respect to the thresh-olds) as we vary lattice spacing or sea quark masses.Searching for a new tetraquark candidate has at least atwo-dimensional parameter space: the hypothetical statewould have an energy and also an overlap onto a spe-cific operator. Using the lattice data presented here, we t / a . . . . . a E e ff η b → η b ( t ) t ( fm ) a . N s E eff , t ++ ( t ) : r x = E eff0 ++ ( t ) : r x = E eff , t ++ ( t ) : r x = E eff0 ++ ( t ) : r x = E eff , t ++ ( t ) : r x = E eff0 ++ ( t ) : r x = M lat ϒ M lat η b t / a . . . . . a E e ff η b → η b ( t ) t ( fm ) a . N s E eff , t ++ ( t ) : r x = E eff0 ++ ( t ) : r x = E eff , t ++ ( t ) : r x = E eff0 ++ ( t ) : r x = E eff , t ++ ( t ) : r x = E eff0 ++ ( t ) : r x = M lat ϒ M lat η b FIG. 10. The individual Wick contraction effective masses ofthe 2 η b → η b correlators. The upper figure is the Direct1and the lower is the Xchange2 contraction (each shown di-agrammatically in Figs. 1a and 1b). E eff and E eff ,t are thesingle- and two-particle effective masses defined in Eq. (18)and Eq. (19) respectively. (color online) can determine a relationship between these parameters.Assuming that a tetraquark does exist below the low-est bottomonium-pair threshold in our data, at a certaintime t ∗ the correlator can be modeled with a two-stateansatz. Specifically for the 0 ++ channel, given that thetetraquark has an energy E b and an overlap Z b ontoa particular operator, at a large enough time the onlyother appreciable contribution will come from the higher2 η b threshold which has an overlap Z η b with the sameoperator. In this case the correlator is given by C ( t ∗ ) = Z b e − aE b t ∗ + Z η b e − aE ηb t ∗ (cid:18) aM η b πt ∗ (cid:19) . (20)Using this ansatz in the effective mass formula Eq. (16)and rewriting the equation in terms of the non-4 − − − −
100 0 E b − E η b (MeV) − − − − − − − − Z b / Z η b a . N s xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx σ − σ σ − σ FIG. 11. The excluded region for the ratio of tetraquark/2 η b overlaps, Z b /Z η b , onto the O ( η b ,η b ) operator, assuming atetraquark with mass, E b , lying below the 2 η b threshold, E η b . The red hashed region is excluded at 5 σ by the dataas described in the text. The 1 σ − σ and 3 σ − σ exclusionbands are also shown for reference. (color online) perturbative overlaps yields the constraint Z b Z η b = − (cid:16) t ∗ t ∗ +1 (cid:17) exp ( aE eff ( t ∗ ) − aE η b )exp( aE eff ( t ∗ ) − aE b ) − × e − a ∆ Et ∗ (cid:18) aM η b πt ∗ (cid:19) (21)with a ∆ E = aE η b − aE b >
0. As can be seen, if thetetraquark is not observed by a time t ∗ then the over-lap onto this new state must be (at least) exponentiallysuppressed with the binding of the tetraquark state, e.g,with − ∆ E .This point illustrates that if a tetraquark did exist with E b < E η b then it is possible that it was not observed inour data because Z b ≈ t ∗ where the two-stateansatz is valid, aE eff ( t ∗ ) from correlator data constructedwith a specific operator, as well as aE η b = 2 aM η b . Forthe local O ( η b ,η b ) operator on the a = 0 .
06 fm ensemble,by examining Figure 5a, a choice of t ∗ = 143 ensures thatthe two-state ansatz is valid (given the long plateau atthe 2 η b threshold). Here, aE eff ( t ∗ = 143) = 0 . aM η b , givenin Table VI, is found from the η b -meson data. Then, us-ing this data in the constraint, for a certain choice of E b ,a numerical value of the ratio of overlaps is found suchthat it is consistent with zero within its small 1 σ statis-tical error. Any value of the ratio of overlaps larger thanthis 1 σ error is inconsistent, at this level of confidence,with our data observing a tetraquark at this value of E b . We use this model to estimate how small the hypo-thetical tetraquark overlap would need to be so that thetetraquark was not observed within our statistical preci-sion. A 1 σ , 3 σ and 5 σ exclusion plot of the parameterspace is given in Figure 11. As the input data into theconstraint has a long propagation time past t ∗ > aE eff ( t ∗ ) which doesnot fall below the threshold, a significant amount of pa-rameter space is excluded. It should be understood thatthis figure is only valid for a particular operator in acertain channel. The given 0 ++ channel in Figure 11 ex-cludes the largest amount of parameter space as it is themost statistically precise. Also Figure 11 is constructedfrom data on the “superfine” ensemble alone, where dis-cretisation effects are smallest and would not change thequantitative behaviour significantly.To conclude this section, we find no evidence of a stabletetraquark candidate below the non-interacting thresh-olds by studying a full S -wave colour-spin basis of QCDoperators. In the next section we will perform an ex-ploratory and complementary study of an alternative ap-proach so to ensure the robustness of our conclusions. V. NRQCD WITH A HARMONIC OSCILLATORPOTENTIAL
A stable tetraquark state in the 0 ++ , 1 + − and 2 ++ channels would overlap with the full basis of S -wavecolour-spin operators utilised above. In this section, wego one additional step by exploring an alternative ap-proach in order to further investigate the possibility of atetraquark state. Adding a central confining potential tothe quark interactions can produce a more deeply boundtetraquark relative to the threshold as the strength ofthis interaction is increased (as we will see below). Fur-thermore, an appropriate choice of additional interactioncan reduce the fiducial volume of the lattice and thusthin the allowed discrete momenta states of the two me-son degrees of freedom. Adding an external attractivescalar central-potential to the QCD interactions yieldsthese desired effects.The harmonic oscillator potential is a particularly suit-able choice of scalar interaction between quarks. For aparticle of mass m at position x away from the centre x the potential is just κr / ≡ κ | x − x | /
2. Defin-ing ω = (cid:112) ( κ/m ), the ground state energy and wave-function are E = 3 ω/ ψ ( r ) = C exp ( − mωr / . Thus we expect that for Defining r = | x − x | and R cm = (( x + x ) / − x ), we can . . . . . . . ω (GeV) . . . . . . r r m s (f m ) η b η b : Set 12 η b : Set 22 η b : Set 32 η b : Set 4 FIG. 12. The root mean square distance r rms of the η b andthe 2 η b as a function of the harmonic oscillator strength cal-culated from a potential model with the different ensembleparameters listed in Table III as discussed in the text. (coloronline) small values of ω the ground state for 2 η b mesons is ap-proximately 2 M η b ( ω ) + 3 ω . The 3 ω term comes fromtwo colour-singlet η b mesons in the harmonic oscillatorpotential and the mass of the η b is shifted slightly fromthe value at ω = 0 because of the additional harmonicoscillator interaction combined with the QCD interac-tions that bind the two heavy quarks into a η b . How-ever, for a compact tetraquark state the mass would be M (¯ b ¯ bbb )( ω ) + 3 ω/
2, as there is only one colour-singletstate in the central harmonic oscillator potential (3 ω/ η b and much less than the effective latticevolume) then its mass will also receive only a modestpositive correction due to ω .Hence if there were a tetraquark state near thresholdthen this additional interaction could drive it further be-low threshold, giving a much cleaner and distinct signalfor the tetraquark candidate in our calculation. As thepotential model framework describing the η b has beenlargely successful, we can use this framework as a generalguide for the exploratory non-perturbative lattice calcu-lation when including the harmonic oscillator potential.Of course, we are mainly interested in QCD (and not theharmonic oscillator) and, as such, if we do find a stabletetraquark state when including the harmonic oscillatorpotential then we must take the ω → separate κ [( x − x ) + ( x − x ) ] / κ/ r + (2 κ ) R cm ] / For sufficiently small ω the shift in η b is directly related to therms radius of the η b state since from perturbation theory it isgiven by (cid:104) η b ( ω = 0) | kr / | η b ( ω = 0) (cid:105) . The harmonic oscillator potential is defined as δH HO = m b ω | x − x | (22)where the quarks are pulled towards the fixed point x with a strength ω . We choose x to be the same positionas the source. Intuitively, as the quarks/ η b ’s start topropagate further away from the source, the harmonicoscillator 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 thequarks/ η b ’s are confined into by the addition of the har-monic oscillator potential. The root mean square dis-tance, r rms , gives an indication of this. We determine the r rms of the η b based on solutions of the Schrodinger equa-tion using a Cornell potential . The Matthieu equationcan describe the behaviour of two free η b ’s in a harmonicoscillator potential on a periodic box, and the solutionsof which can yield r rms for the 2 η b state. The resultsfrom such a calculation are plotted in Figure 12.Based on this, in order to confine the quarks suffi-ciently so that the two η b ’s overlap, and also to studythe dependence on ω , values of ω = (75 , , , ω = (75 , , , κ/ am b )( aω ) / . , . , . , . . , . , . , . e − a ˜ H = (cid:18) − aδH HO l (cid:19) l e − aH (cid:18) − aδH HO l (cid:19) l (23)where e − aH is the purely NRQCD evolution equationdefined in Eq. (12). This implementation was chosen sothat the evolution equation is still time-reversal symmet-ric. 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 cal-culations on Set 1 and 3 respectively. Following thesedetails, we are now able to present results from the non-perturbative lattice calculations. A. Numerical Results
All correlator data from Set 1 and 3 discussed in Sec-tion IV was generated again with the inclusion of theharmonic oscillator potential at the four different ω val-ues given above. The harmonic oscillator alters both the In a periodic box of length L the harmonic oscillator is alsoperiodic. The potential form is − α s r + rb with α s = 0 . b = 2 .
34 GeVand reduced mass 2 .
59 GeV . . . M e V a . N s (’hhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhh’,) E fit η b E eff , ξη b ( t ) E eff η b ( t ) . . . M e V . . . M e V t / a . . . ω = M e V FIG. 13. The effective mass plot for the η b when includingthe harmonic oscillator potential on Set 3. E eff is given byEq. (16) while E eff ,ξ removes the leading 1+ e − ωt dependencefrom the correlator (24) to enable a better comparison withthe data when no harmonic oscillator potential is included.(colour online) single- and two-particle contributions to the correlator sothat they become (as derived in Appendix B) dependenton ω as C J PC i,j ( t, ω ) = (cid:88) n Z in Z j, ∗ n (1 + e − ωt ) e − ( M ( ω ) n + ω ) t (24) C J PC i,j ( t, ω ) = (cid:88) X Z iX Z j, ∗ X (cid:18) ωµ r π − − e − ωt (cid:19) × e − ( M S ( ω )+ M S ( ω )+3 ω ) t + · · · . (25)First, Figure 13 shows the effective masses, as defined inEq. (16), of the η b when including the harmonic oscilla-tor. Also overlaid are the effective masses when removingthe 1 + e − ωt dependence to enable a better comparisonwith the data when no harmonic oscillator is included.As can be seen, the dip in the harmonic oscillator ef-fective masses is from this additional time dependence.Physically, this can be understood to be due to the b -quarks travelling non-relativistically and so it takes timefor the harmonic oscillator to have an effect.The η b correlator data when including the harmonicoscillator potential is fit to the functional form given byEq. (24) in order to extract the lowest energy eigenstate M ( ω ) η b +3 ω/ η b mass with no harmonic oscillator included ( M ( ω = 0) η b )and then plot the energy differences against ω , as shownin Figure 14. Good qualitative agreement between the . . . . . . . . ω (GeV) . . . . . . . ∆ M ( ω ) η b + ω / ( G e V ) Model a ≈ .
12 fm a ≈ .
09 fm
FIG. 14. The lattice η b energy M ( ω ) η b +3 ω/ M ( ω = 0) η b subtractedcompared to the model predictions as discussed in the text.(color online) lattice results and the potential model predictions is ob-served.For the ¯ b ¯ bbb system, we show the 0 ++ effective masseson Set 3 in Figure 15. It is evident that the 0 ++ andthe η b data contains more noise when a harmonic oscil-lator potential is included. While fitting the data to theform in Eq. (25) can be performed, it is not necessary asthe purpose of this exploratory work is to determine ifa stable tetraquark exists when ω (cid:54) = 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 Direct1and Xchange2 Wick contractions of the 2 η b → η b corre-lator in Figure 16. As before, the effective masses of theindividual Wick contractions drop below the 2 η b thresh-old, even though importantly, when added together toyield the full correlator shown in Figure 15a the effectivemass is always above threshold. As this was also seenin the pure NRQCD data shown in Sec. IV, it may be aproblematic feature of models that utilise a phenomeno-logically motivated four-body potential for this system.To conclude this section, despite adding an auxiliarypotential into the QCD interactions that should push anear threshold tetraquark candidate increasingly lowerwe find no indication of any state below the 2 η b threshold.The conclusions of this section then agree with those ofSec. IV. VI. DISCUSSION AND CONCLUSIONS
In this work we have studied the low-lying spectrumof the ¯ b ¯ bbb system using the first-principles lattice non-relativistic QCD methodology in order to search for astable tetraquark state below the lowest non-interactingbottomonium-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Υ.7In Section III we describe our numerical methodology.Four gluon ensembles were employed with lattice spac-ings ranging from a = 0 . − .
12 fm, and one ensem-ble which has physical light-quark masses. All ensembleshave u , d , s and c quarks in the sea.In Sec. IV we presented the majority of the resultsin this work. Here, we determined the lowest energyeigenstate of the ¯ b ¯ bbb system with the quantum num-bers 0 ++ , 1 + − and 2 ++ using an over-constrained S -wave colour/spin basis (arising from Fierz relations be-tween the diquark-antidiquark and two-meson systems asshown in Table II). We did not observe any state belowthe lowest non-interacting bottomonium-pair thresholdin any channel, as can be seen in Figures 5, 6, 7 and 8,and a summary of our results from this section is givenin Figure 9.In Sec. V, to ensure the robustness of our conclu-sions, we performed an exploratory calculation of a novelmethod which added an auxiliary scalar potential into theQCD interactions with the objective of pushing a nearthreshold tetraquark increasingly lower than the thresh-old. This would give a more distinct and cleaner signalfor its presence in our calculation. The harmonic oscil-lator was found to be a suitable central scalar potential.For the η b -meson with this potential, we first verifiedagreement between the non-perturbative lattice calcu-lations and a potential model (as shown in Figure 14)and then used this potential model as a general guideto choose multiple appropriate values of the potentialstrength. Despite studying the ¯ b ¯ bbb system with this ad-ditional scalar potential on the lattice, no indication ofthe QCD tetraquark was observed.This work is the only first-principles study of the low-lying ¯ b ¯ bbb spectrum in the literature. However, there areothers which utilise different methodologies. For exam-ple, [8] predicts the tetraquark mass by solving the two-particle Schrodinger equation with a phenomenologicallymotivated non-confining potential between the point-likediquark and anti-diquark, and finds a 0 ++ , + − and 2 ++ tetraquark to be bound by 44, 51 and 5 MeV respec-tively . The authors of [14] used a diquark model includ-ing 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-mean-square distance between the diquark and anti-diquarkinside the tetraquark in this model is similar in mag-nitude to the distance between the quarks inside eachdiquark [9]. Consequently, such an approach is inter-nally inconsistent. More recently, [10] used a Hamilto-nian including only spin-spin interactions mediated by aone-gluon exchange. Here, all other effects such as the A model in which the diquarks are taken to be fundamentalparticles cannot determine the two-meson threshold from theSchrodinger equation because the diquarks cannot recombineinto mesons. Thus, the experimental meson masses [4] are usedto determine the lowest threshold. chromoelectric interactions, colour confinement and the b -quark mass need to be set separately. The authorsset these additional contributions in two ways: by es-timating 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 statecould be either below the 2 η b threshold or lie inbetweenthe 2 η b and 2Υ thresholds (where in both cases thethresholds were determined using the experimental me-son masses). In [55] the author also uses a model in-cluding only chromomagnetic interactions and finds anunbound tetraquark. Within the QCD sum-rules frame-work, [11] finds a tetraquark candidate approximately300 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 thethresholds. Indeed, in the limit of very heavy quarkswhere the force proceeds through one-gluon exchangecontaining only colour-Coulomb contributions (safely ne-glecting spin and long distance effects), the authors of[7] used a variational methodology to determine thata bound tetraquark exists for a QQ ¯ Q (cid:48) ¯ Q (cid:48) system when m Q /m Q (cid:48) < ∼ .
15 (where both m Q and m Q (cid:48) are heavy rel-ative to Λ QCD ). However, if m Q /m Q (cid:48) is varied and thetetraquark becomes unbound, then as the free two-mesoneigenstate becomes the ground state of the system thisnumerical methodology has an increasingly slow conver-gence to a solution [7] (being numerically ill-posed) dueto a redundant degree of freedom in the minimisationprocedure. Indeed [7] indicates that for m Q /m Q (cid:48) = 1the solution is unstable, a hint that no bound tetraquarkexists for all identical quarks in the very heavy masslimit. The authors of [14] assume that the b -quark issufficiently heavy so to use the one-gluon exchange withonly colour-Coulomb contributions, and by also neglect-ing the mixing between different colour-components ofthe 2 × η b mass to de-termine the threshold). In an orthogonal direction tothe above work, the authors of [15] only include a linearstring contribution in the one-gluon exchange (neglectingspin effects and the appreciable short-distance Coulombcontributions), and without mixing between the differentcolour components of the potential matrix, find a boundtetraquark when m Q /m Q (cid:48) = 1. However, in subsequentwork [56, 57], by modelling the aforementioned mixingthey concluded that no bound tetraquark exists. Per-haps the most sophisticated non first-principles method-ology used to study the four-body ¯ b ¯ bbb tetraquark is thediffusion Monte Carlo method utilised by [9]. Here, onedetermines the ground state of a phenomenologically mo-tivated Hamiltonian by solving the Schrodinger equa-tion and examining the stability of (cid:80) n e − ( E n − E ) t Ψ n ( x )to determine E , where E n (Ψ n ) is the n -th energy-eigenstate (eigenfunction). The authors include both thecolour-Coulomb and linear contributions in the gluonexchange but neglect the mixing between the different8 . . . M e V a . N s (’hhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhh’,)2 M lat η b M lat η b E eff , ‡ η b ( t ) E eff η b ( t ) . . . M e V . . . M e V t / a . . . ω = M e V (a) 2 η b → η b . . . M e V a . N s (’hhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhh’,)2 M lat η b M lat η b E eff , ‡ η b ( t ) E eff η b ( t ) . . . M e V . . . M e V t / a . . . ω = M e V (b) 2Υ → . . . M e V a . N s (’hhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhh’,)2 M lat η b M lat η b E eff , ‡ η b ( t ) E eff η b ( t ) . . . M e V . . . M e V t / a . . . ω = M e V (c) 3 c × ¯3 c → c × ¯3 c . . . M e V a . N s (’hhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhh’,)2 M lat η b M lat η b E eff , ‡ η b ( t ) E eff η b ( t ) . . . M e V . . . M e V t / a . . . ω = M e V (d) 6 c × ¯6 c → c × ¯6 c FIG. 15. The ¯ b ¯ bbb effective masses for the 0 ++ correlators when including the harmonic oscillator potential on Set 3. E eff isgiven by Eq. (16), while E eff , ‡ removes the leading 1 − e − ωt dependence from the correlator (25) to enable a better comparisonwith the data when no harmonic oscillator potential is included. (color online) colour components in the potential matrix, and find a sta-ble tetraquark candidate 108 MeV below the 2 η b thresh-old (determined from the experimental meson mass).Consequently, there is no study in the literature whichis not from first-principles that includes all the appre-ciable effects relevant for the bb ¯ b ¯ b system: treating thebottom-quarks as fundamental particles, including bothshort and long distance effects in the gluon exchange andincluding the mixing between the different colour compo-nents in the 2 × η b → η b correlator are shown in Figure 10. As is evident there,in each individual Wick contraction the effective massdrops below the 2 η b threshold but rises slowly to thresh-old even though when all Wick contractions are addedtogether 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 withthe additional scalar potential, shown in Figure 16, pos-sibly indicating that this may be a problematic featureof models that utilise a phenomenologically motivatedfour-body potential: a subset of the interactions show9 . . . M e V a . N s (’hhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhh’,)2 M lat η b M lat η b E eff , ‡ η b ( t ) E eff η b ( t ) . . . M e V . . . M e V t / a . . . ω = M e V . . . M e V a . N s (’hhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhh’,)2 M lat η b M lat η b E eff , ‡ η b ( t ) E eff η b ( t ) . . . M e V . . . M e V t / a . . . ω = M e V FIG. 16. The effective masses for the individual 2 η b → η b Wick contraction correlator data when including a harmonicoscillator potential on Set 3. The upper figure is the Direct1and the lower is the Xchange2 contraction (each shown dia-grammatically in Figs. 1a and 1b). (color online) behaviour that may be misinterpreted as a bound statebelow threshold, while when all interactions are includedno bound state is seen. Particularly, the slow rise tothreshold from below could make the diffusion MonteCarlo method practically difficult due to the slowly vary-ing stability condition combined with the fact that a longevolution time (greater than 8 fm) is necessary.In conclusion: we find no evidence of a ¯ b ¯ bbb tetraquark with a mass below the lowest non-interactingbottomonium-pair thresholds in the 0 ++ , 1 + − or 2 ++ channels. We give a constraint in Eq. (21) that futurephenomenological models must satisfy if such QCD statesare postulated. For the 0 ++ channel, we use this con-straint to estimate how small the non-perturbative over-lap of the hypothetical tetraquark (onto a particular op-erator) would need to be, relative to the 2 η b , so thatit was not observed within our statistical precision. A − − −
200 0 200 E b − E η b (MeV) This Calc.Lattice QCDPheno.Models CLSZ17BLO16AFGSZ17KNR17W17WLCZ16 † BLN11AFGSZ17VVR13 ξ FIG. 17. A comparison of our result for the ¯ b ¯ bbb ground stateenergy in the 0 ++ channel (stat. error only) and predictionsfrom phenomenological models. The hatched region indicatesthe exclusion of a bound tetraquark with an energy E b sub-ject to the value of its non-perturbative overlap as given inFigure 11. In this comparison, we take our ground state en-ergy obtained on the “superfine” ensemble (Set 4 in Table III)as a representative because it has the smallest discretisaioneffects and the statistical error encompasses the results on theother ensembles as shown in Figure 9. The y -axis labels re-sults from different phenomenological models [7–13] by lastinitial of authors and year of publication. An error is plottedif given in the reference. The two results from WLCZ16 † dif-fer by how the mass scale was set. The result VVR13 ξ findsno bound tetraquark and we indicate this by placing theirresult on threshold. (color online) σ , 3 σ and 5 σ exclusion plot of the parameter space isshown in Figure 11, and discussed in Sec. IV. As we havepropagation times longer than 8 fm and statistically pre-cise data, we can exclude all but the most finely-tunedparameter space. Our lattice results then rule out thephenomenological models discussed above that predict atetraquark below the lowest bottomonium-pair thresh-olds which have a value of non-perturbative overlap thatis excluded by Fig. 11. A comparison of these resultswith ours is shown in Figure 17.Further studies of possible heavy tetraquark channelsthat include orbital angular momentum either betweenthe mesons in the tetraquark or between the quarks in themeson could be performed with the methodology used0here . Similarly one could also study whether stable¯ c ¯ ccc , ¯ b ¯ cbc or ¯ b ¯ bcc tetraquarks exist or not. Additionally,two-hadron systems receive a finite-volume energy shiftwhich depends on the infinite-volume scattering ampli-tude which is non-trivial to parameterise. Here we donot calculate these finite-volume energy shifts. Doing soin a more extended study would allow statements to bemade about the existence of resonant tetraquark statesabove the lowest thresholds, that likely do exist in nature.Quantifying these shifts would be an exciting avenue forfuture work.Finally, recent work based on heavy-quark symmetry[6] and phenomenological arguments [58] indicates that a J P = 1 + ¯ b ¯ bud tetraquark will be stable in QCD. In fact,by extracting a potential from the lattice in the staticheavy-quark limit and solving the Schrodinger equation[59–61] also finds binding in this channel. Initial latticecalculations hint that such a state exists [62] but calcu-lations are difficult because of a signal-to-noise problemfor heavy-light states [63]. Lattice QCD calculations inthis direction are essential for a conclusive first-principlesstatement to be made and to give further motivation fora targeted experimental search for these tetraquark con-figurations of nature. The unaveraged correlator data, that have been anal-ysed to produce the results in this work, are publicly avail-able in a SQLite database from any of the authors uponrequest . ACKNOWLEDGMENTS
We would like to thank William Bardeen and ZhenLiu for the many insightful discussions on tetraquarks,as well as Gavin Cheung who in addition gave guidanceon the tetraquark implementation. We are also grate-ful to the MILC collaboration for the use of their gaugeconfigurations. This manuscript has been authored byFermi Research Alliance, LLC under Contract No. DE-AC02-07CH11359 with the U. S. Department of Energy,Office of Science, Office of High Energy Physics. TheUnited States Government retains and the publisher, byaccepting the article for publication, acknowledges thatthe United States Government retains a non-exclusive,paid-up, irrevocable, world-wide license to publish or re-produce the published form of this manuscript, or allowothers to do so, for United States Government purposes. It should be noted that although we focused on S -wave com-binations of quarks, the channels we study also exclude certaincombinations of orbital angular momentum from producing abound tetraquark. For example, the 0 ++ overlaps with 2Υ in anorbital angular momentum D -wave configuration. If this stateproduced a low-lying bound tetraquark it would also show up inour calculation. The unaveraged correlator data is too large to be hosted for freeon a remote server.
The results described here were obtained using the Dar-win Supercomputer of the University of Cambridge HighPerformance Computing Service as part of the DiRACfacility jointly funded by STFC, the Large Facilities Cap-ital Fund of BIS and the Universities of Cambridge andGlasgow.
Appendix A: Two-Point Correlator Fit Functions
Here we derive the non-relativistic two-particle con-tribution to the correlator on our ensembles. To begin,the correlator is given in Eq. (5). For clarity, the i, j subscripts are dropped. The completeness relation for atwo-hadron system is [64] I = (cid:88) X (cid:90) d P tot (2 π ) d k (2 π ) E ( X ) | X P tot , k ) (cid:105)(cid:104) X P tot , k ) | (A1)where | X P tot , k ) (cid:105) = | M ( k ) M ( P tot − k ) (cid:105) is a two-hadronstate (with quantum numbers suppressed) and to avoidsuperfluous notation, we will also set P tot = 0. A keydifference from the one-hadron system is the internalrelative momentum, k , which contributes an additionalthree-integral. Substituting the completeness relationEq. (A1) into the correlator Eq. (5) and performing themomentum conserving integrals yields C ( t ) = (cid:88) X (cid:90) d k (2 π ) Z X ( k ) e − E ( X ) t (A2)where Z X ( k ) is a non-perturbative coefficient.However, on a discrete finite-volume the above inte-gral over elastic states is replaced by a finite sum with k m ∈ ( − N s / , . . . , , N s /
2) in units of (2 π/aN s ). Inturn, for P tot = 0, Eq. (A2) becomes a sum over back-to-back hadronic states which have values of the discretemomenta that are equal in magnitude but opposite in di-rection. One can expand the two-particle energy using anon-relativistic dispersion relation, appropriate since weare using NRQCD, as E ( X ) = (cid:113) M + | k | + (cid:113) M + | k | (A3) ≈ M S + M S + | k | µ r (A4)where we have defined the static, kinetic and reducedmasses by M S , M K and µ r = M K M K / ( M K + M K )respectively. In a finite-volume there will be an addi-tional contribution to Eq. (A4) dependent on the infinite-volume scattering phase shift, which will be discussedfurther below. Eq. (A4) also illustrates the density ofback-to-back states on our ensembles. As an exam-ple, examining the a = 0 .
09 fm ensemble, and taking M η b = 9 . | k | / µ r ≈
20 MeV or 0 . | k | / µ r (with the methodology used in [38]).Practically, this would be overly computationally expen-sive and instead, the fact that the states with k (cid:54) = arerelated by the dispersion relation (and are not indepen-dent as the sum would assume) should be included.This can be achieved by first expanding the non-perturbative coefficient Z X ( k ) as a polynomial in | k | /µ r , as dictated by rotational symmetry and by en-suring the Taylor coefficients have the same dimension,then keeping all terms needed to a certain precision. Af-ter this the correlator can be written as C ( t ) = (cid:88) X e − ( M S + M S ) t (cid:88) k (cid:110) ∞ (cid:88) i =0 Z lX | k | l µ lr (cid:111) e − | k | µr t (A5)= (cid:88) X e − ( M S + M S ) t (cid:90) πa − πa d k (2 π ) (cid:110) ∞ (cid:88) i =0 Z lX | k | l µ lr (cid:111) e − | k | µr t (A6)where going from the first to the second line we havereplaced the finite sum by an integral. Taking the limitsof the integral to ±∞ and performing the integrals over k analytically yields the fit function given in Eq. (7). Onceit is shown that it is possible to replace the finite sumby the indefinite integral within our statistical precisionthen it is valid to use the above fit function with ourdata. − − | k | (GeV) . . . . . . . . | k | e xp (cid:16) − | k | µ r t (cid:17) a . N s hhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhh ± π / at = t = t = t = t = t = t = t = FIG. 18. The integrands of the moments given in Eq. (A8)at multiple times. The crosses represent the discrete finite-volume momentum contributions on the coarse (Set 1) en-semble as discussed in the text. Due to the Gaussian timedependence, the integrand peak moves towards the origin forlarger times. (color online)
To do so, using spherical coordinates in Eq. (A6), wedefine the quantities that we need to compare as I ( l ) ( t ) = 1 µ lr (cid:90) ∞−∞ d | k || k | l +2 e − | k | µr t (A7) D ( l ) ( t ) = 1 µ lr (cid:88) | k | | k | l +2 e − | k | µr t . (A8)The integrands of both are shown diagrammatically inFigure 18, where it is observed that due to the Gaus-sian time-dependence the peaks of the integrand movetowards the origin with larger t . As such, one objectiveis to choose a large enough ˆ t such that a sufficient major-ity of the integrand is within the maximum momentum π/a . We can replace the discrete finite-volume fit func-tion with it’s infinite-volume counterpart if the relativedifference between them is less than our statistical errors.Specifically if (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:80) l max l =0 Z , ( l ) I ( l ) ( t ) − (cid:80) ∞ l =0 Z , ( l ) D ( l ) ( t ) (cid:80) l max l =0 Z , ( l ) I ( l ) ( t ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (A9) ≤ (cid:12)(cid:12)(cid:12)(cid:80) l max l =0 Z , ( l ) I ( l ) ( t ) − (cid:80) ∞ l =0 Z , ( l ) D ( l ) ( t ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) Z , (0) I (0) (cid:12)(cid:12) (A10) ≤ l max (cid:88) l =0 Z , ( l ) Z , (0) (cid:12)(cid:12) I ( l ) ( t ) − D ( l ) ( t ) (cid:12)(cid:12) I (0) + ∞ (cid:88) l = l max +1 Z , ( l ) Z , (0) D ( l ) ( t ) I (0) ≤ l max (cid:88) l =0 (cid:12)(cid:12) I ( l ) ( t ) − D ( l ) ( t ) (cid:12)(cid:12) I (0) + ∞ (cid:88) l = l max +1 D ( l ) ( t ) I (0) (A11) ≤ δC ( t ) C ( t ) (A12)where l max is the maximum number of moments to be in-cluded in the fit function, the inequality in the second lineholds as the moment integrands are positive (shown dia-grammatically in Figure 18), in the third line the Cauchyinequality has been used, and in the fourth line it is as-sumed that the leading moment gives the dominant con-tribution ( Z , ( l ) ≤ Z , (0) ). Studying Eq. (A11) insteadof Eq. (A9) is a conservative option.Each part of the first term in Eq. (A11) representshow similar I ( l ) ( t ) and D ( l ) ( t ) need to be in order to beconsidered equivalent within statistical precision. This isshown in Figure 19. For a particular l max , the secondterm represents when the higher moments look like noisewithin statistical precision, also shown in Figure 19. EachFigure was generated with the coarse ensemble parame-ters (listed as Set 1 in Table III) as this ensemble hasthe largest lattice spacing (and hence smallest π/a value − the upper limit on the integral of interest) and alsothe smallest N s (the number of discrete momenta usedin the finite-volume sum). As such, the other ensembleswill give better approximations and studying Set 1 is con-servative. Overlaid on each plot is the smallest relativestatistical error from the data on any ensemble. Due tothe constant signal-to-noise ratio, the number of config-urations and the size of the lattice spacing, the smallest2 − − − − − − − I ( l ) ( t ) − D ( l ) ( t ) I ( ) ( t ) a . N s Stat.Err. t ∗ t (fm) − − − − − − − D ( l ) ( t ) / I ( ) ( t ) hhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhh l = l = l = l = l = FIG. 19. The difference between the discrete finite-volumeand infinite-volume continuum moments (upper) and whichmoments can be neglected compared to our statistical preci-sion (lower) as discussed in the text. (color online) statistical error was the 2 η b correlator on the fine ensem-ble. Only examining situations below this curve is themost conservative option for all data generated. As canbe observed in Figure 19, the discrete finite-volume sumsare well represented by the indefinite integrals. Addition-ally, in order to neglect the higher moments within ourstatistical precision, a choice of ˆ t = 1 fm and l max = 2 issufficient.Expanding the finite-volume two-particle energy non-relativistically in Eq. (A4) neglected a possible finite-volume energy shift which depends on the infinite-volumescattering amplitude. In the small scattering-lengthlimit, the energy shift is known to be volume suppressed[40]. The two-particle systems under study are in thislimit as the η b and Υ are compact due to the heavy-quark mass, with a size of 0 . − . Appendix B: Two-Point Correlator Fit FunctionsWith A Harmonic Oscillator
In the non-relativistic limit the free propagator (toleading order) is∆( x , t ) = (cid:90) d p (2 π ) exp ( i p · x ) exp (cid:18) − (cid:26) m + p m (cid:27) t (cid:19) . (B1)The free two-meson propagator with P tot = 0, bothstarting at common origin x = ( ,
0) and ending at time t is given by ˜∆( t ) = (cid:90) d x ∆ ( x , t )∆ ( x , t ) . (B2)Using Eq. (B1) in Eq. (B2) produces the large-time be-haviour of the free two-meson propagator as˜∆( t ) = (cid:16) µ r πt (cid:17) / e − ( m + m ) t . (B3)This agrees with the leading behaviour derived in App. A.Next, for the harmonic oscillator case, the one dimen-sional Hamiltonian is ∂ψ∂t = Eψ = − m ∂ ψ∂x + κ x ψ. (B4)Solutions of this system can be related to the solutionsof the Mehler differential equation via ∂φ∂ ˜ t = ∂ φ∂ρ − ρ φ (B5)with the identifications ω = (cid:112) κ/m , t = 2˜ t/ω and r =( κm ) − / ρ . The Greens function (propagator) for theMehler differential equation is given by∆( ρ , ρ , ˜ t ) = 1 (cid:112) π sinh(2˜ t ) exp (cid:16) − coth (2˜ t ) ( ρ + ρ )2+ csch (2˜ t ) ρ ρ (cid:17) . (B6)To normalize this propagator one can first com-pare the large t behaviour of this solution with theknown behaviour of the harmonic oscillator propagator,lim t →∞ G ( t ) = | Ψ( ) | e − ωt , where the wavefunction atthe origin is given by Ψ( ) = ( mω/π ) / . As such, theharmonic oscillator solution is∆( x , , t ) = √ mω (cid:112) π sinh( ωt ) exp (cid:18) − mω x ωt ) (cid:19) . (B7)The three-dimensional solution can then be obtained us-ing the separability of each spatial direction, so thatthe zero spatial-momentum single-particle correlator ina harmonic oscillator potential is given by (cid:90) d x ∆( x , , t ) = (cid:18) ωt ) (cid:19) / e − mt . (B8)3Finally, the equal mass two-particle propagator startingat a common origin x = ( ,
0) and ending at a time t with P tot = , in the presence of an external har-monic oscillator potential, is found from using Eq. (B7) in Eq. (B2), to give˜∆( t ) = (cid:18) mω π ωt ) (cid:19) / e − mt . (B9)By comparing (B9) to (B3), and noting that m = 2 µ r , wesee that the free two-meson harmonic oscillator propaga-tor reduces to the non harmonic oscillator case as ω → [1] R. J. Jaffe, Phys. Rev. D , 267 (1977).[2] J. J. Dudek, R. G. Edwards, and D. J. Wilson(Hadron Spectrum Collaboration), Phys. Rev. D93 ,094506 (2016), arXiv:1602.05122 [hep-ph].[3] R. A. Briceno, J. J. Dudek, R. G. Edwards,and D. J. Wilson (Hadron Spectrum Collaboration),arXiv:1708.06667 [hep-lat].[4] C. Patrignani et al. (Particle Data Group), Chin. Phys.
C40 , 100001 (2016).[5] N. Brambilla, S. Eidelman, B. Heltsley, R. Vogt, G. Bod-win, et al. (Quarkonium Working Group), Eur.Phys.J.
C71 , 1534 (2011), arXiv:1010.5827 [hep-ph].[6] E. J. Eichten and C. Quigg, (2017), arXiv:1707.09575[hep-ph].[7] A. Czarnecki, B. Leng, and M. B. Voloshin, (2017),arXiv:1708.04594 [hep-ph].[8] A. V. Berezhnoy, A. V. Luchinsky, and A. A. Novoselov,Phys. Rev.
D86 , 034004 (2012), arXiv:1111.1867 [hep-ph].[9] Y. Bai, S. Lu, and J. Osborne, (2016), arXiv:1612.00012[hep-ph].[10] J. Wu, Y.-R. Liu, K. Chen, X. Liu, and S.-L. Zhu,(2016), arXiv:1605.01134 [hep-ph].[11] W. Chen, H.-X. Chen, X. Liu, T. G. Steele, and S.-L.Zhu, (2016), arXiv:1605.01647 [hep-ph].[12] Z.-G. Wang, Eur. Phys. J.
C77 , 432 (2017),arXiv:1701.04285 [hep-ph].[13] M. Karliner, S. Nussinov, and J. L. Rosner, Phys. Rev.
D95 , 034011 (2017), arXiv:1611.00348 [hep-ph].[14] M. N. Anwar, J. Ferretti, F.-K. Guo, E. Santopinto, andB.-S. Zou, (2017), arXiv:1710.02540 [hep-ph].[15] J. Vijande, A. Valcarce, and J. M. Richard, Phys. Rev.
D76 , 114013 (2007), arXiv:0707.3996 [hep-ph].[16] E. Eichten and Z. Liu, arXiv:1709.09605 [hep-ph].[17] R. Vega-Morales, (2017), arXiv:1710.02738 [hep-ph].[18] J. J. Dudek, R. G. Edwards, M. J. Peardon, D. G.Richards, and C. E. Thomas (Hadron Spectrum Collabo-ration), Phys. Rev. D , 034508 (2010), arXiv:1004.4930[hep-ph].[19] G. T. Bodwin, E. Braaten, and G. P. Lepage,Phys. Rev. D51 , 1125 (1995), [Erratum: Phys.Rev.D55,5853(1997)], arXiv:hep-ph/9407339 [hep-ph].[20] A. Bazavov et al. (MILC Collaboration), Phys. Rev.
D82 , 074501 (2010), arXiv:1004.0342 [hep-lat].[21] A. Hart, G. M. von Hippel, and R. R. Horgan(HPQCD Collaboration), Phys. Rev. D , 074008(2009), arXiv:0812.0503 [hep-lat].[22] E. Follana et al. (HPQCD Collaboration, UKQCD Col-laboration), Phys. Rev. D , 054502 (2007), arXiv:hep-lat/0610092 [hep-lat]. [23] G. P. Lepage, L. Magnea, C. Nakhleh, U. Magnea, andK. Hornbostel, Phys. Rev. D , 4052 (1992), arXiv:hep-lat/9205007 [hep-lat].[24] R. J. Dowdall et al. (HPQCD Collaboration), Phys. Rev.D , 054509 (2012), arXiv:1110.6887 [hep-lat].[25] J. O. Daldrop, C. T. H. Davies, and R. J. Dowdall(HPQCD Collaboration), Phys. Rev. Lett. , 102003(2012), arXiv:1112.2590 [hep-lat].[26] R. J. Dowdall, C. T. H. Davies, T. Hammant, R. R.Horgan, and C. Hughes (HPQCD Collaboration),Phys. Rev. D89 , 031502 (2014), [Erratum: Phys.Rev.D92,039904(2015)], arXiv:1309.5797 [hep-lat].[27] R. J. Dowdall, C. T. H. Davies, T. C. Hammant, andR. R. Horgan (HPQCD Collaboration), Phys. Rev. D ,094510 (2012), arXiv:1207.5149 [hep-lat].[28] R. J. Dowdall, C. T. H. Davies, R. R. Horgan, C. J. Mon-ahan, and J. Shigemitsu (HPQCD Collaboration), Phys.Rev. Lett. , 222003 (2013), arXiv:1302.2644 [hep-lat].[29] B. Colquhoun, R. J. Dowdall, C. T. H. Davies, K. Horn-bostel, and G. P. Lepage (HPQCD Collaboration), Phys.Rev. D , 074514 (2015).[30] C. Hughes, R. J. Dowdall, C. T. H. Davies, R. R. Horgan,G. von Hippel, and M. Wingate (HPQCD Collabora-tion), Phys. Rev. D92 , 094501 (2015), arXiv:1508.01694[hep-lat].[31] B. Chakraborty, C. T. H. Davies, B. Galloway, P. Knecht,J. Koponen, G. C. Donald, R. J. Dowdall, G. P. Lep-age, and C. McNeile (HPQCD Collaboration), Phys.Rev.
D91 , 054508 (2015), arXiv:1408.4169 [hep-lat].[32] T. C. Hammant, A. G. Hart, G. M. von Hippel, R. R.Horgan, and C. J. Monahan (HPQCD Collaboration),Phys. Rev. D , 014505 (2013), arXiv:1303.3234 [hep-lat].[33] C. J. Morningstar, Phys. Rev. D50 , 5902 (1994),arXiv:hep-lat/9406002 [hep-lat].[34] D. C. Moore and G. T. Fleming, Phys. Rev. D , 014504(2006), arXiv:hep-lat/0507018 [hep-lat].[35] C. E. Thomas, R. G. Edwards, and J. J. Dudek(Hadron Spectrum Collaboration), Phys. Rev. D ,014507 (2012), arXiv:1107.1930 [hep-lat].[36] J. Vijande and A. Valcarce, Symmetry , 155 (2009),arXiv:0912.3605 [hep-ph].[37] E. Berkowitz, T. Kurth, A. Nicholson, B. Joo, E. Ri-naldi, M. Strother, P. M. Vranas, and A. Walker-Loud(CalLat Collaboration), Phys. Lett. B765 , 285 (2017),arXiv:1508.00886 [hep-lat].[38] J. J. Dudek, R. G. Edwards, and C. E. Thomas(Hadron Spectrum Collaboration), Phys. Rev.
D86 ,034031 (2012), arXiv:1203.6041 [hep-ph]. [39] G. Lepage, B. Clark, C. Davies, K. Hornbostel,P. Mackenzie, et al. , Nucl. Phys. Proc. Suppl. , 12(2002), arXiv:hep-lat/0110175 [hep-lat].[40] G. Parisi, Physics Reports , 203 (1984).[41] G. P. Lepage, in Boulder ASI 1989:97-120 (1989) pp.97–120.[42] M. L¨uscher, Communications in Mathematical Physics , 177.[43] R. J. Dowdall, C. T. H. Davies, G. P. Lepage, and C. Mc-Neile (HPQCD Collaboration), Phys. Rev.
D88 , 074504(2013), arXiv:1303.1670 [hep-lat].[44] M. L¨uscher, Nuclear Physics B , 531 (1991).[45] R. A. Briceno and Z. Davoudi, Phys. Rev.
D88 , 094507(2013), arXiv:1204.1110 [hep-lat].[46] R. A. Briceno, Phys. Rev.
D89 , 074507 (2014),arXiv:1401.3312 [hep-lat].[47] P. Guo, J. Dudek, R. Edwards, and A. P. Szczepaniak,Phys. Rev.
D88 , 014501 (2013), arXiv:1211.0929 [hep-lat].[48] M. Gockeler, R. Horsley, M. Lage, U. G. Meissner,P. E. L. Rakow, A. Rusetsky, G. Schierholz, and J. M.Zanotti, Phys. Rev.
D86 , 094513 (2012), arXiv:1206.4141[hep-lat].[49] C. h. Kim, C. T. Sachrajda, and S. R. Sharpe, Nucl.Phys.
B727 , 218 (2005), arXiv:hep-lat/0507006 [hep-lat].[50] S. He, X. Feng, and C. Liu, JHEP , 011 (2005),arXiv:hep-lat/0504019 [hep-lat].[51] T. Iritani et al. (HAL QCD Collaboration), JHEP ,101 (2016), arXiv:1607.06371 [hep-lat].[52] T. Yamazaki, K.-I. Ishikawa, Y. Kuramashi, andA. Ukawa (PACS Collaboration), PoS LATTICE2016 , 108, arXiv:1702.00541 [hep-lat].[53] A. Bazavov et al. (MILC Collaboration), Phys. Rev.
D87 , 054505 (2013), arXiv:1212.4768 [hep-lat].[54] N. Brambilla, G. Krein, J. Tarrs Castell, and A. Vairo,Phys. Rev.
D93 , 054002 (2016), arXiv:1510.05895 [hep-ph].[55] B. Silvestre-Brac, Phys. Rev. D , 2179 (1992).[56] J. Vijande, A. Valcarce, and J.-M. Richard, Phys. Rev. D87 , 034040 (2013), arXiv:1301.6212 [hep-ph].[57] J.-M. Richard, A. Valcarce, and J. Vijande, Phys. Rev.
D95 , 054019 (2017), arXiv:1703.00783 [hep-ph].[58] M. Karliner and J. L. Rosner, (2017), arXiv:1707.07666[hep-ph].[59] P. Bicudo and M. Wagner (European Twisted MassCollaboration), Phys. Rev.
D87 , 114511 (2013),arXiv:1209.6274 [hep-ph].[60] P. Bicudo, K. Cichy, A. Peters, B. Wagenbach,and M. Wagner, Phys. Rev.
D92 , 014507 (2015),arXiv:1505.00613 [hep-lat].[61] P. Bicudo, M. Cardoso, A. Peters, M. Pflaumer, andM. Wagner, (2017), arXiv:1704.02383 [hep-lat].[62] A. Francis, R. J. Hudspith, R. Lewis, and K. Maltman,Phys. Rev. Lett. , 142001 (2017), arXiv:1607.05214[hep-lat].[63] C. T. H. Davies, C. McNeile, E. Follana, G. P. Lep-age, H. Na, and J. Shigemitsu (HPQCD Collaboration),Phys. Rev.
D82 , 114504 (2010), arXiv:1008.4018 [hep-lat].[64] M. Srednicki,