Searching for beauty-fully bound tetraquarks using lattice nonrelativistic QCD Searching for beauty-fully bound tetraquarks using lattice nonrelativistic QCD Ciaran Hughes,1,* Estia Eichten,1,† and Christine T. H. Davies2,‡ 1Fermi National Accelerator Laboratory, Batavia, Illinois 60510, USA 2SUPA, School of Physics and Astronomy, University of Glasgow, Glasgow G12 8QQ, United Kingdom (Received 16 October 2017; published 14 March 2018) Motivated by multiple phenomenological considerations, we perform the first search for the existence of a b̄b̄bb tetraquark bound state with a mass below the lowest noninteracting bottomonium-pair threshold using the first-principles lattice nonrelativistic QCD methodology. We use a full S-wave color/spin basis for the b̄b̄bb operators in the three 0þþ, 1þ− and 2þþ channels. We employ four gluon field ensembles at multiple lattice spacing values ranging from a ¼ 0.06–0.12 fm, all of which include u, d, s and c quarks in the sea, and one ensemble which has physical light-quark masses. Additionally, we perform novel exploratory work with the objective of highlighting any signal of a near threshold tetraquark, if it existed, by adding an auxiliary potential into the QCD interactions. With our results we find no evidence of a QCD bound tetraquark below the lowest noninteracting thresholds in the channels studied. DOI: 10.1103/PhysRevD.97.054505 I. INTRODUCTION Tetraquarks were first considered theoretically decades ago in the context of light-quark physics in order to explain, amongst other experimental features, the a0ð980Þ and f0ð980Þ broad resonances [1].1 More recently, there has been exciting experimental evidence indicating the poten- tial existence of tetraquark candidates amongst the so- called XYZ states—states whose behavior differs from predictions of the heavy quark-antiquark potential model. The observed XYZ states apparently contain two heavy quarks, ðcc̄Þ or ðbb̄), and two light quarks [4]. The dynamics of these systems involves both the short-distance and long-distance behavior of QCD and hence theoretical predictions are difficult. Consequently, many competing phenomenological models currently exist for these states [5]. Lattice QCD studies of the observed XYZ states are also difficult because these states are high up in the spectrum as well as being in the threshold region for strong decays into two heavy flavor mesons. While there are theoretical arguments that some tetraquark states with doubly heavy flavor (e.g., bbūd̄, bbūs̄ and bbd̄s̄) should be bound and stable against all strong decays [6], no general arguments exist for tetraquarks with heavy quark-antiquark content such as QQ̄0qq̄0 states. A tetraquark system composed of four heavy quarks is a much cleaner system to study theoretically as long-distance effects from light quarks are expected not to be appreciable, as opposed to systems which are a mixture of heavy and light quarks. In the limit of very heavy quarks perturbative QCD single-gluon exchange dominates [7] and so the dynamics are relatively simple. This makes these systems particularly useful to study in order to shed light on the aforementioned XYZ states. In fact, there is a multitude of phenomenological models (with a quark mass ranging from the bottom to the very heavy limit) which predict the existence of a Q̄Q̄QQ bound tetraquark [8–15]. However, these are not calculations from first principles and have an unquantifiable systematic error associated with the choice of four-body potential. In reality, the heaviest possible tetraquark system in nature would be a b̄b̄bb tetraquark. For this, nonperturbative QCD cannot be ignored, making a first-principles lattice QCD study essential. If such a bound b̄b̄bb tetraquark did exist, how it would be observed at the LHC has already been addressed [16,17]. Given these pressing theoretical motivations, in this work we perform the first lattice QCD study of the b̄b̄bb system. The sole objective of this exploratory work is to determine if the dynamics of QCD generates enough binding force between the b̄b̄bb to produce a tetraquark state with a mass below the lowest noninteracting botto- monium-pair threshold, ensuring that it is stable against simple strong decays. Searching for such a bound b̄b̄bb *chughes@fnal.gov †eichten@fnal.gov ‡christine.davies@glasgow.ac.uk Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article’s title, journal citation, and DOI. Funded by SCOAP3. 1Recent lattice studies of the scattering amplitude pole do indicate that these states are in fact resonances as opposed to cusp effects etc. [2,3]. PHYSICAL REVIEW D 97, 054505 (2018) 2470-0010=2018=97(5)=054505(22) 054505-1 Published by the American Physical Society https://crossmark.crossref.org/dialog/?doi=10.1103/PhysRevD.97.054505&domain=pdf&date_stamp=2018-03-14 https://doi.org/10.1103/PhysRevD.97.054505 https://doi.org/10.1103/PhysRevD.97.054505 https://doi.org/10.1103/PhysRevD.97.054505 https://doi.org/10.1103/PhysRevD.97.054505 https://creativecommons.org/licenses/by/4.0/ https://creativecommons.org/licenses/by/4.0/ tetraquark candidate is particularly well suited to the first- principles lattice QCD methodology because this state, if it existed, would be the ground state in the b̄b̄bb system. This means that it should be relatively easy to identify. Further, b̄b̄bb annihilation effects are strongly suppressed by the heavy-quark mass, as in the bottomonium system, and so can be ignored. This paper is organized as follows: in Sec. II the inter- polating operators used in this study are discussed, in Sec. III the computational methodology is given, while the majority of the results are presented in Sec. IV. In Sec. V we explore a novel method of adding an auxiliary potential into QCD with the objective of highlighting a possible tetraquark signal. We then discuss our conclusions in Sec. VI. II. INFINITE-VOLUME CONTINUUM EIGENSTATES, OPERATORS AND TWO-POINT CORRELATORS The QCD Fock space contains all color-singlet single- particle states such as the conventional mesons jηbðkÞi, jΥðkÞi etc. and, if a b̄b̄bb bound state also exists, a tetraquark state jT4bðkÞi. In addition, there are also the two-particle states which can be labeled by appropriate quantum numbers as jPtot; JPC; jkrelj; JP1C11 ; JP2C22 ; Lreli where Ptot (J PC) is the total (angular) momentum of the two-particle system, with JPiCii being the quantum numbers of the individual particles and krel (Lrel) the relative (orbital angular) momentum between the two particles. The sole motivation of this work is to search for a possible b̄b̄bb tetraquark candidate within QCD that couples to a bottomonium pair and which lies below the lowest threshold. The bottomonium mesons we study can be classified as JPC ¼ 2Sþ1LJ. As any orbital angular momenta are expected to raise the internal kinetic energy of the state (and hence its rest mass) we focus on two-body S-wave systems (L ¼ 0) with no orbital angular momen- tum between them (Lrel ¼ 0). With the quantum numbers of the Υ=ηb being JPC ¼ 1−−=0−þ, the S-wave 2ηb and 2Υ can have (through the addition of angular momenta) a quantum number of 0þþ, while the Υηb has 1þ− and the 2Υ can also be in a 2þþ configuration. We now want to construct a full basis of S- wave color/spin interpolating operators that has overlap with these quantum numbers. To do so, we start by forming all possible color combinations that the 2b̄ and 2b can be in. These are specified in Table I. We can construct meson interpolating operators as O 1ð8Þ M ðt; xÞ ¼ G 1ð8Þ efg b̄fΓMbgðt; xÞ; ð1Þ where ΓM ¼ iγ5; γk projects onto the quantum numbers of the ηb and Υ respectively, and G 1ð8Þ efg is the color projection 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̄ Ĉ f ΓDbgðt; xÞ; ð2Þ O 3ð6̄Þ A ðt; xÞ ¼ G 3ð6̄Þ efg b̄fΓAb Ĉ g ðt; xÞ; ð3Þ where ðbĈÞα ¼ Cαβb̄β is the charge-conjugated field with C ¼ −iγ0γ2.2 As the two quarks have the same flavor, the Pauli-exclusion principle applies and the wave function has to be completely antisymmetric. With our choice to focus on S-wave combinations of particles, the spatial wave function must be symmetric. As the color (triplet) sextet has a (anti) symmetric color wave function, this forces the spin-wave function to be in a (triplet) singlet with (Γ ¼ γk) Γ ¼ iγ5. With these building blocks, we can form four classes of b̄b̄bb color singlets by contracting the color factors G in any irreducible representation (irrep) with its conjugate color factors, i.e., 1c × 1c, 8c × 8c, 3 3c × 3̄c and 6c × 6̄c. These SUð3Þ invariant color contractions are given in Table I. After doing this, we need to project the operators onto a specific angular momentum JP by using the standard SOð3Þ Clebsch-Gordan coefficients (using a spherical basis of spin matrices [18]) as OJ;mðP;QÞðt; xÞ ¼ X m1;m2 hJ; mjJ1; m1; J2; m2i × OJ1;m1P ðt; x þ rÞOJ2;m2Q ðt; xÞ ð4Þ with ðP; QÞ describing the blocks this configuration is built from, i.e., ðηb; ηbÞ, ðΥ; ΥÞ, ðD; AÞ, etc. We also allow the possibility of the two blocks being separated by a distance r. For the r ¼ 0 case, the operators project onto a definite total angular momentum J. For the r ≠ 0 case, one can Taylor expand OJ1;m1P ðt; x þ rÞ around r ¼ 0 to notice that the operator projects onto a superposition of quantum TABLE I. The color representations of the different quark combinations. Note that, as described in the text, once the color representation of the (anti-) diquark is chosen, the Pauli-exclusion principle enforces certain spin combinations in S-wave. Also given are the SUð3Þ color contractions needed for the b̄b̄bb operators. b b̄ b̄b b̄b̄ bb Color irrep 3c 3̄c 1c, 8c 3c, 6̄c 3̄c; 6c G1efgG 1 ef0g0 δfgδf0g0 G8efgG 8 ef0g0 2δfg0δf0g − 2δfgδf0g0=3 G3efgG 3 ef0g0 ðδff0δgg0 − δfg0δgf0Þ=2 G6efgG 6 ef0g0 ðδff0δgg0 þ δfg0δgf0Þ=2 2γ0 ¼ diagð1; −1Þ in the convention used by NRQCD. 3Since the 8c irrep of SUð3Þ is the adjoint, it is similar to its conjugate 8̄c. HUGHES, EICHTEN, and DAVIES PHYS. REV. D 97, 054505 (2018) 054505-2 numbers. Consequently, it is possible to utilize this to further search for the lowest ground state of the four quark system. When dealing with the diquark components, to project onto a definite value of charge conjugation in Eq. (4) one can form the linear combination O1;m1D O 1;m2 A � O1;m2D O1;m1A . In fact, not all of these color combinations are indepen- dent. Fierz relations constrain the number of independent color-spin operators that are possible. For the local oper- ators in S-wave, the relations between the two-meson and diquark-antidiquark states are given in Table II. The simplest quantity that can be calculated on the lattice in order to extract particle masses is the Euclidean two- point correlator. This is defined as CJ PC i;j ðt; Ptot ¼ 0Þ ¼ Z d3xhOJ;mii ðt; xÞO J;mj j ð0; 0Þ†i ð5Þ where we choose to project to zero spatial momentum and i, j label potentially different operators at the source and sink with the same JPC, e.g., i ¼ ðηb; ηbÞ, j ¼ ðΥ; ΥÞ. The single-particle contributions to the correlator are deter- mined by inserting a complete set of single-particle states in the Hilbert-space formalism into Eq. (5) to yield CJ PC i;j ðt; Ptot ¼ 0Þ ¼ X n ZinZ j;� n e−Ent ð6Þ with Zin ¼ h0jOJ;mii jni being the nonperturbative overlap of the operator to the eigenstate jni and Enjni ¼ Hjni the energy eigenvalue. Note that all states jni with the same quantum numbers contribute to this correlator, e.g., for the bottomonium 0−þ pseudoscalar correlator the jηbi as well as all radial excitations contribute. The two-particle contributions to the correlator are slightly more complicated. In this case, as derived in Appendix A, the nonrelativistic two-particle states give a contribution to the correlator that is CJ PC i;j ðt;Ptot ¼ 0Þ ¼ � μr 2πt �3 2 X X2 e−ðM S 1 þMS 2 Þt × � Z0X2 þZ2X2 3 ðtμrÞ þZ4X2 15 ðtμrÞ2 þ��� � ð7Þ where the sum is over all distinct two-particle states X2 with quantum numbers JPC and Ptot ¼ 0, MSi (MKi ) is the static (kinetic) mass of the particles, as defined in Eq. (14), μr ¼ MK1 M K 2 =ðMK1 þ MK2 Þ is the reduced mass and Z2lX2 are nonperturbative coefficients. Energies of states can be extracted using the above functional form once the correlator has been computed. Examining Eq. (5) in the path-integral formalism we can perform the connected Wick contractions4 so that the correlator can be written as an integral over the gluon fields with the integrand consisting of products of b-quark propagators. For each two-meson-type operator, e.g., O1cηbO 1c ηb, as all quarks have the same flavor there are four connected Wick contractions. These are shown diagram- matically in Fig. 1. The first Wick contraction for the two-meson correlator, called Direct1 and shown in Fig. 1(a), has the expression GRAdeG R Ad0e0G R BghG R Bg0h0 X x Tr½ΓM1K−1ðt; x; 0; zÞe0g × Γ†M1K −1ðt; x; 0; zÞ†hd0�Tr½ΓM2K−1ðt; x0; 0; z0Þeg0 × Γ†M2K −1ðt; x0; 0; z0Þ†h0d� ð8Þ TABLE II. Fierz relations in the b̄b̄bb system relating the two- meson and the diquark-antidiquark bilinears. JPC Diquark antidiquark Two meson 0þþ 3̄c × 3c − 1 2 j0; ΥΥi þ ffiffi 3 p 2 j0; ηbηbi 0þþ 6c × 6̄c ffiffi 3 p 2 j0; ΥΥi þ 1 2 j0; ηbηbi 1þ− 3̄c × 3c 1ffiffi 2 p ðj1; Υηbi þ j1; ηbΥiÞ 2þþ 3̄c × 3c j2; ΥΥi (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 flavor. The grey region represents a color 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 antiquark is exchanged between the meson pair, (c) Direct3 where each meson propagates to the other, and (d) Xchange4 where a quark is exchanged between the meson pair. 4Annihilation diagrams are suppressed by powers of the heavy-quark mass [19] and are expected to be negligible. SEARCHING FOR BEAUTY-FULLY BOUND TETRAQUARKS … PHYS. REV. D 97, 054505 (2018) 054505-3 while the second, called Xchange2, is given by − GRAdeG R Ad0e0G R BghG R Bg0h0 X x Tr½ΓM1K−1ðt; x; 0; zÞe0g × Γ†M1K −1ðt; x0; 0; zÞ†hdΓM2K−1ðt; x0; 0; z0Þeg0 × Γ†M2K −1ðt; x; 0; z0Þ†h0d0� ð9Þ where x0 ¼ x þ r. The other diagrams, Direct3 and Xchange4, have similar expressions. For the diquark- antidiquark-type operators, e.g., O6cD O 6̄c A , there are also four Wick contractions which can be combined into one expression as CJ PCðt; Ptot ¼ 0Þ ¼ ½1 � sgnðCΓDÞT � sgnðCΓAÞT þ sgnðCΓDÞTsgnðCΓAÞT�GRAdeGRAd0e0GRBghGRBg0h0 × X x Tr½CΓDK−1ðt; x; 0; zÞeg0Γ†DCK−1ðt; x; 0; zÞTdh0� × Tr½ΓACK−1ðt; x0; 0; z0Þ�e0gCΓ†AK−1ðt; x0; 0; z0Þ†hd0� ð10Þ where the sign function is defined by sgnðXÞT ¼ �1 if XT ¼ �X. The � in Eq. (10) corresponds to the 3 or 6 color representation. It is this prefactor with the signs which enforces the Pauli-exclusion principle: the sum cancels for spin combinations that do not make the wave function overall antisymmetric. The spin-triplet/singlet configura- tions we consider here obey ðCγkÞT ¼ þðCγkÞ and ðCγ5ÞT ¼ −ðCγ5Þ. Diagrammatically the four connected Wick contractions contributing to the diquark correlator are shown in Fig. 2. To calculate the two-point correlators described above within the first-principles Feynman path-integral approach to QCD needs the methodology of lattice QCD. We now discuss our lattice QCD approach. III. LATTICE QCD METHODOLOGY A. Second generation Nf = 2 + 1 + 1 gluon ensembles Our lattice calculation uses gauge field configurations generated by the MILC collaboration [20]. For the gauge fields, theyused the tadpole-improvedLüscher-Weisz gauge action correct to Oðαsa2Þ [21] and include 2 þ 1 þ 1 flavors in the sea, the up and down quarks (treated as two degenerate light quarks with mass ml), the strange quark, and the charm quark. The sea quarks are included using the highly improved staggered quark formulation [22]. Four ensembles are chosen in this study. As one might expect roughly double the discretization errors in a 2b̄b system relative to the b̄b one, to ensure that the heavy-quark potential is accurately represented (where short-distance details may be important for a compact tetraquark candidate) we utilize three ensembles that span relatively fine lattice spacings ranging from a ¼ 0.06–0.12 fm. Details of the ensembles are given in Table III. Due to the computational expense, most of the ensembles use heavier ml than in the real world. However, to test ml dependence, we use one ensemble (set 2 in Table III) that has physical aml=ams. Additionally, the ensembles have been fixed to Coulomb gauge to allow nongauge invariant nonlocal operators to be used [as constructed in Eq. (4)]. B. b-quarks using iNRQCD A nonrelativistic effective field theory is appropriate for physical systems where the relative velocity of the (a) (b) (c) (d) FIG. 2. There are four connected Wick contractions for the diquark-antidiquark-type correlator when the quarks have the same flavor. The blue shaded region represents a diquark, the red shaded region the antidiquark, a blue line a quark and the red line an antiquark. The uncrossing of the lines in Fig. 2(b) to produce Fig. 2(a) gives a �, as discussed in the text, which enforces the Pauli-exclusion principle. TABLE III. Details of the gauge ensembles used in this study. β is the gauge coupling. a (fm) is the lattice spacing [23,24], amq are the sea quark masses, Ns × NT gives the spatial and temporal extent of the lattices in lattice units and ncfg is the number of configurations used for each ensemble. We use 16 time sources on each configuration to increase statistics. Ensembles 1 and 2 are referred to as “coarse,” 3 as “fine,” and 4 as “superfine”. Set β a (fm) aml ams amc Ns × NT ncfg 1 6.00 0.1219(9) 0.0102 0.0509 0.635 24 × 64 1052 2 6.00 0.1189(9) 0.00184 0.0507 0.628 48 × 64 1000 3 6.30 0.0884(6) 0.0074 0.037 0.440 32 × 96 1008 4 6.72 0.0592(3) 0.0048 0.024 0.286 48 × 144 400 HUGHES, EICHTEN, and DAVIES PHYS. REV. D 97, 054505 (2018) 054505-4 constituent particles inside the bound state is much smaller than 1 (in Planck units). When applied to the strong force, this framework is called nonrelativistic QCD (NRQCD). It is well known that b-quarks are very nonrelativistic inside their low-lying bound states, where v2 ≈ 0.1 [25]. For lattice NRQCD, the continuum NRQCD action is discretized onto the lattice [25] with operators included to a predetermined level in v2. Here we use an action accurate through Oðv4Þ with additional spin-dependent terms at Oðv6Þ.5 Operators are also added to correct for discretiza- tion effects. We make the further systematic improvement here, introduced in [23], to include coefficients of Oðv4Þ operators that have been matched to continuum QCD through OðαsÞ. We call this improved NRQCD (iNRQCD). This action has already been used to success- fully determine bottomonium S, P and D wave mass splittings [23,26], precise hyperfine splittings [27,28], B meson decay constants [29], Υ and Υ0 leptonic widths [30], B, D meson mass splittings [28] and hindered M1 radiative decays between bottomonium states [31]. The iNRQCD Hamiltonian evolution equations can be written as Gðx; t þ 1Þ ¼ e−aHGðx; tÞ Gðx; tsrcÞ ¼ δx;x0 ð11Þ with e−aH ¼ � 1 − aδHjtþ1 2 �� 1 − aH0jtþ1 2n � n U†t ðxÞ × � 1 − aH0jt 2n � n � 1 − aδHjt 2 � ð12Þ aH0 ¼ − Δð2Þ 2amb ; aδH ¼ aδHv4 þ aδHv6; aδHv4 ¼ −c1 ðΔð2ÞÞ2 8ðambÞ3 þ c2 i 8ðambÞ2 ð∇ · ~E − ~E · ∇Þ − c3 1 8ðambÞ2 σ · ð ~∇ × ~E − ~E × ~∇Þ − c4 1 2amb σ · ~B þ c5 Δð4Þ 24amb − c6 ðΔð2ÞÞ2 16nðambÞ2 aδHv6 ¼ −c7 1 8ðambÞ3 fΔð2Þ; σ · ~Bg − c8 3 64ðambÞ4 fΔð2Þ; σ · ð ~∇ × ~E − ~E × ~∇Þg − c9 i 8ðambÞ3 σ · ~E × ~E: ð13Þ The parameter n is used to prevent instabilities at large momentum from the kinetic energy operator. A value of n ¼ 4 is chosen for all amb values. We evaluate the propagator using local sources Eq. (11). Here, amb is the bare b quark mass, ∇ is the symmetric lattice derivative, with ~∇ being the improved version, and Δð2Þ, Δð4Þ are the lattice discretizations of ΣiD2i , ΣiD 4 i respectively. ~E, ~B are the improved chromoelectric and chromomagnetic fields, details of which can be found in [23]. Each of these fields, as well as the covariant derivatives, must be tadpole improved. We take the mean trace of the gluon field in Landau gauge, u0L ¼ h13 TrUμðxÞi, as the tadpole param- eter, calculated in [23,29]. The matching coefficients c1, c2, c4, c5, c6 in the above Hamiltonian have been computed perturbatively to one loop [23,32]. c3 was found to be close to the tree-level value of one nonperturbatively [23] and we set it, as well as the rest of the matching coefficients, to their tree-level values. The quark mass amb is tuned fully nonperturbatively in iNRQCD [27] using the spin-averaged kinetic mass of the Υ and ηb (which is less sensitive to spin- dependent terms in the action). The above Hamiltonian neglects the four-fermion operators which appear at Oðα2sv3Þ, as well as other operators which appear at higher order in the nonrelativistic expansion. Simple power- counting estimates [23,27] would lead us to expect con- tributions of order a few percent (or a few MeV) at most to binding energies from these terms. In the case where a tetraquark bound state is observed then we can estimate the systematic effect from neglecting these contributions. The parameters used in this study are summarized in Table IV. There, c1, c5 and c6 are the correct values for an Oðv4Þ iNRQCD action [23]. For set 4, all parameters are those for the Oðv4Þ action. The small changes to these coefficients in going to an Oðv6Þ iNRQCD action (which are similar in magnitude to the two-loop corrections) are not appreciable for the purpose of this work: whether or not a tetraquark candidate exists below the lowest bottomo- nium-pair threshold. All other parameters listed in Table IV are taken from [27,31]. Within iNRQCD the single-particle energy eigenstates can be decomposed in the standard nonrelativistic expan- sion as TABLE IV. Parameters used for the valence quarks. amb is the bare b-quark mass in lattice units, u0L is the tadpole parameter and the ci are coefficients of terms in the NRQCD Hamiltonian [see Eq. (13))]. Details of their calculation can be found in [23,32]. c3, c7, c8 and c9 are included at tree level. Set amb u0L c1, c6 c2 c4 c5 1 2.73 0.8346 1.31 1.02 1.19 1.16 2 2.66 0.8350 1.31 1.02 1.19 1.16 3 1.95 0.8525 1.21 1.29 1.18 1.12 4 1.22 0.8709 1.15 1.00 1.12 1.10 5The spin-independent Oðv6Þ terms are subleading effects for the b̄ b̄ bb energies relevant to this study. SEARCHING FOR BEAUTY-FULLY BOUND TETRAQUARKS … PHYS. REV. D 97, 054505 (2018) 054505-5 EðPÞ ¼ MS þ jPj 2 2MK þ � � � ð14Þ with MS, MK being the static and kinetic masses respec- tively [23]. Because the quark mass term is removed from the iNRQCD Hamiltonian, the static mass is unphysical, differing by a constant shift from the physical (kinetic) mass. This means that only static mass differences deter- mined fully nonperturbatively can be compared to exper- imental results. However, this constant shift can be calculated in lattice perturbation theory if required [33], or by using an additional experimental input. The kinetic mass does not suffer this problem as it acquires the quark mass contributions from the quark kinetic terms [23]. Also, within iNRQCD the Dirac field Ψ can be written in terms of the quark ψ and antiquark χ as Ψ ¼ ðψ; χÞT. The propa- gator is then K−1ðxjyÞ ¼ � GψðxjyÞ 0 0 −GχðxjyÞ � ð15Þ where GψðxjyÞ is the two-spinor component quark propa- gator and GχðxjyÞ is the two-spinor component antiquark propagator. Taken together, we can now compute the b- quark propagator via the iNRQCD evolution equations on the gluon ensembles listed in Table III. The last piece needed to calculate the two-point correlators, and hence energies, is the discretized finite-volume versions of the interpolating operators. C. Discrete finite-volume operators Together, the isotropic discretization and the periodic finite volume break the infinite-volume continuum SOð3Þ symmetry of NRQCD to the octahedral group, Oh [34]. Thus, while the operators constructed in Sec. II have well- defined JPC quantum numbers associated with SOð3Þ, we need to construct operators which transform within the irreps of the Oh symmetry group (relevant for lattice calculation). This can be achieved by the method of subduced representations, where an operator with a specific JPC can be taken to a specific lattice irrep6 ΛPC by using the subduction coefficients found in Appendix A of [18]. At rest, each of our JPC ¼ 0þþ and JPC ¼ 1þ− operators trivially subduce into a single lattice irrep labeled by Aþþ1 and Tþ−1 respectively. However, the J PC ¼ 2þþ case is slightly more complicated and subduces into two lattice irreps labeled by Tþþ2 and E þþ (which are three and two dimensional). We construct both the T2=E operators asffiffiffi 2 p O ½2� T2=E ¼ OJ¼2;m¼2 � OJ¼2;m¼−2, which are correctly subduced from the J ¼ 2þþ operators defined in Sec. II. In principle, each lattice irrep allows mixing between different J states, e.g., the A1 irrep contains not only the J ¼ 0 states but also the J ¼ 4 [34]. However in practice since we are only looking for the ground state of the b̄b̄bb correlators we are not sensitive to these mixing effects. A complete list of b̄b̄bb interpolating operators used to produce the correlator data herein is given in Table V. In fact, this is an overconstrained set due to the Fierz identities (shown in Table II) which relate the two-meson and diquark-antidiquark correlators. We include this overcon- strained system and ensure that we reproduce the Fierz relations to numerical precision, performing a nontrivial check on our data. Additionally, we also reproduce the relations between the 8c × 8c color combination and the others [36] on a subset of the data. It has been found [37] that separating each hadron within the two-hadron interpolating operator by a specific distance r can significantly aid in the extraction of the (ground) state energy. In this direction, we use three different spatial configurations of the b̄b̄bb correlators where the individual building blocks are separated by a distance of rx ¼ 0, 1 or 2 lattice units in the x-direction.7 Finally, the subduced lattice interpolating operators are defined in terms of the Dirac fields as in Eq. (1), and the correlators are defined analogously, as in Eq. (10). We can then use the decomposition of K−1ðxjyÞ given in Eq. (15) with suitable boundary conditions to write the correlator in terms of the iNRQCD quark propagator GψðxjyÞ. Due to TABLE V. The b̄b̄bb interpolating operators used in this study. Operators in each column are subduced from the infinite-volume continuum quantum numbers JPC given in the first row. The superscript on each operator denotes the lattice irrep of that operator and the subscript denotes the building blocks of the operator, as explained in the text. We generate each operator with three different spatial configurations as shown in Eq. (4): where the building blocks are separated by a distance rx ¼ 0, 1 or 2 lattice units in the x-direction. 0þþ 1þ− 2þþ Source Sink Source/sink Source/sink OA1ðηb;ηbÞ O A1 ðηb;ηbÞ O T1 ðΥ;ηbÞ O T2 ðΥ;ΥÞ OA1ðηb;ηbÞ O A1 ðΥ;ΥÞ O T1 ðD3̄c ;A3c Þ OEðΥ;ΥÞ OA1ðΥ;ΥÞ O A1 ðηb;ηbÞ O T2 ðD3̄c ;A3c Þ OA1ðΥ;ΥÞ O A1 ðΥ;ΥÞ O E ðD3̄c ;A3c Þ OA1ðD3̄c ;A3c Þ OA1ðD3̄c ;A3c Þ OA1ðD6c ;A6̄c Þ OA1ðD6c ;A6̄c Þ 6The conserved quantum numbers of a symmetry group are determined using the little group, which for SOð3Þ and Oh depend on the momenta type [35]. In this study we focus on states at rest. 7The ηb and Υ energies used to determine the noninteracting 2ηb and 2Υ thresholds (needed to determine if a state exists below them) are found from locally smeared meson correlators only. HUGHES, EICHTEN, and DAVIES PHYS. REV. D 97, 054505 (2018) 054505-6 our use of iNRQCD, there are no backward propagating valence antiquarks in our calculation. Consequently, the appreciable finite-temporal effects seen in relativistic two- meson lattice correlators8 [38] do not arise in our calcu- lation, simplifying the analysis. With this methodology, it is now possible to compute the lowest energy levels asso- ciated with the b̄b̄bb system. IV. THE LOW-LYING ENERGY EIGENSTATES OF THE 0+ + , 1+ − AND 2+ + b̄b̄bb SYSTEM In order to determineif there is an energy eigenstate below the 2ηb threshold, we first need to find the noninteracting thresholds on each ensemble listed in Table III. This is achieved by computing the bottomonium ηb and Υ two- point correlators as described above, then fitting them to the functional form given in Eq. (6) to extract the single-particle energies. As in the range of studies listed in Sec. III B, amongst others, we utilize the well-established Bayesian fitting methodology [39] in this work and refer the reader to [23] for technical details. Although we fit the correlator data in order to extract fit parameters, in the following we display effective mass plots so the reader can visualize the data. The single-particle effective mass is constructed as aEeff JPC ¼ log � CJ PC i;j ðtÞ CJ PC i;j ðt þ 1Þ � ð16Þ ¼aEJPC þ Zi1Z j;� 1 Zi0Z j;� 0 e−ðE1−E0Þtð1−e−ðE1−E0ÞÞþ��� ð17Þ ⟶ t→∞ aEJPC: ð18Þ As can be inferred from Eq. (17), the greater the mass gap E1 − E0 or the larger the ground state overlap Z0, the quicker aEeff JPC converges to a plateau, which gives aEJPC. The ηb and Υ effective masses are shown in Fig. 3, where the returned fitted ground-state energy from the correlator fits is also shown overlaid in black. The large difference in energy betweentheη0b (Υ 0) andthegroundstateηb (Υ)meansthatthe effective mass plots fall rapidly to a plateau given by the ground-state energy, indicating that the fit to the correlator data extracts the ground-state energy precisely. Also evident from the effective mass plot is the constant signal to noise in the ηb data, as might be expected from straightforward application of the Parisi-Lepage arguments [40,41] for noise growth in a system where all the quarks are the same and no 0þþ bound tetraquark exists. The argument specifies that the noise of the two-point correlator should behave like expfðEJPC − EGS=2Þtg where EJPC is the lowest energy eigenstate of the bottomonium operator OJPC con- structed to have the quantum numbers JPC and EGS is the lowest energy eigenstate of the mean squared correlator which controls the noise. Thus, it would be surprising if a tetraquark candidate did exist below the 2ηb threshold from thelatticeperspectivealoneasthenEGS < 2Eηb andthenoise of the ηb data would grow exponentially. However, the lattice calculation still needs to be performed for a conclusive statement to be made about the existence of this tetraquark candidate since the Parisi-Lepage arguments do not allow for raw crossed Wick contractions that wouldcontribute toeither the full two-meson or tetraquark correlator. Lattice correlators are affected by both the discrete nature of the space-time lattice and separately by its finite volume. Each has a separate but calculable effect on the extracted lattice energies. Corrections in energies due to discretization effects are proportional to ak, where k depends on the level of improvement. Here systematic discretization errors are reduced to α2sa 2 by the improvements made, as for those studies listed in Sec. III B, and we expect this to be small enough to have little impact. We can assess this from our results with different values of the lattice spacing. Finite-volume effects for single-particle energies (arising from the lightest particle in the sea propagating around the spatial boundaries) are known to behave like expð−aMπNsÞ [42] and are not appreciable for the ensembles used here which have aMπNs ≈ 4 [43]. In fact, Ns on set 1 and 2 differ by a factor of 2, giving a basic test of volume dependence. However, finite-volume interactions can shift a two-particle energy by an amount that depends on the infinite-volume scattering matrix. Further, these shifts are nontrivial to parameterise (see for example [44–50]). As the specific purpose of this study is to search for a hypothetical tetraquark bound state below the lowest threshold, we do not attempt to quantify these finite-volume interactions. Equation (7) describes the nonrelativistic two-particle contribution to the correlator after t ¼ 1 fm (as shown in FIG. 3. The effective mass plot for the ηb and Υ on the superfine ensemble (set 4 listed in Table III). The effective mass plots on the other ensembles are qualitatively identical. 8These arise in lattice two-meson calculations when a rela- tivistic formulation of valence quark is used due to one of the mesons propagating forward in time while the other propagates around the temporal boundary backwards in time [38]. SEARCHING FOR BEAUTY-FULLY BOUND TETRAQUARKS … PHYS. REV. D 97, 054505 (2018) 054505-7 Appendix A). In this case, because of the additional 1=t 3 2 time dependence, the effective mass formula for these contributions differs from their single-particle counterparts. Removing the leading order time dependence yields an effective mass defined as aEeff;t JPC ¼ log � t 3 2CJ PC i;j ðtÞ ðt þ 1Þ32CJPCi;j ðt þ 1Þ � : ð19Þ For the 0þþ, 1þ− and 2þþ operators that are constructed, if a stable tetraquark exists below the 2ηb, ηb þ Υ or 2Υ thresholds then it would show up as the ground state of each correlator and hence also in the effective masses. Otherwise, each threshold would be the lowest energy eigenstate. Higher energy states also appear in each correlator. For example, the 2Υ and ηbη0b in the 0 þþ case, the Υη0b in the 1 þ− and the ΥΥ0 and χb0χb2ð1PÞ in the 2þþ. Of these, if no tetraquark state is present, only the 2Υ in the 0þþ might be noticeable while studying the ground state, as it is Oð100Þ MeV above the 2ηb. All other excited states have similar energy differences to those appearing in the ηb=Υ effective masses shown in Fig. 3, which rapidly falls to the ground state. It is helpful at this stage to generate mock correlator data to illustrate how we might expect the b̄b̄bb correlator results to behave in the presence or absence of a low-lying tetraquark state (neglecting two-particle finite-volume effects). Using the extracted lattice ηb and Υ masses, we can compute the noninteracting 2ηb and 2Υ thresholds on our ensembles. Further, for a fixed value of Z0X2, we can also compute their leading order two-particle contribution to the correlator from Eq. (7). Additionally we can infer the values of the noninteracting ηbη 0 b and ΥΥ 0 masses on our ensembles using the experimental PDG [4] values as input in order to include their two-particle contributions in the mock correlator data also. Further, in the 0þþ channel, if we assume a tetraquark bound state exists 100 MeV below the 2ηb threshold, 9 for a fixed value of the tetraquarks non- perturbative overlap, Z4b, this hypothetical state’s contri- bution to the correlator is given by Eq. (6). Then, given the effective mass formula defined in Eq. (16), for each different choice of the nonperturbative coefficients we can generate a separate effective mass curve. In practice, we choose different values of the coefficients from a normal distribution with zero mean and unit variance. Figure 4 shows such a plot for the superfine ensemble (set 4 in Table III) where the solid blue curves represent different values of the normally distributed coefficients. As the nonperturbative coefficients show up through the ratio Z0X2=Z4b in the effective mass formula, analogously to Eq. (17), once the energy difference E1 − E0 is set the effective mass is only sensitive to the relative size of the tetraquark overlap to that of the lowest threshold (once the contribution of excited states has become negligible). With this knowledge, the lower red dotted curve gives mock data in the situation where there is only a new state present in the correlator data (Z0X2 ¼ 0), the middle red dashed curve indicates the case when the tetraquark and two-meson states have the same value of coefficient (Z0X2 ¼ Z4b), while the upper dot-dashed red curve is the mock data in the case where there is no new state in the data (Z4b ¼ 0). The blue curves below the middle red one have a larger overlap onto the new state while those above have an increasingly vanishing one. With our overconstrained color-spin basis of S-wave operators at least one operator should have an appreciable overlap onto a tetraquark state below threshold (if it exists) and, as illustrated by the mock data, give an easy/clean signal similar to the middle red dashed curve where the effective mass drops below the 2ηb threshold. Even though Fig. 4 is mock data, the upper red curve with no tetraquark state present looks very like the real data shown in Figs. 5 and 6. One important point to note is the additional factor of 1=t 3 2 appearing in the two-particle contribution in Eq. (7) relative to the single-particle one in Eq. (6). This factor suppresses the two-particle contribution relative to the single-particle state; e.g., t ¼ 100 gives a suppression of the two-particle state of ð0.01Þ1.5. This is one reason why the middle red dashed curve has a particularly rapid fall to the new state which lies only 100 MeV below the 2ηb [compared to Fig. 5(d) where the 2Υ is Oð100Þ MeV above the 2ηb]. Overall this effect would produce an enhancement of the stable tetraquark state if it exists. Further, while we used the superfine ensemble for illustrative purposes, the other ensembles with larger lattice spacings have similar features but with less temporal resolution due to the larger numerical value of the lattice spacing. FIG. 4. Assuming a tetraquark exists 100 MeV below the 2ηb threshold, different normally distributed couplings in b̄b̄bb correlator mock data produce different effective mass curves on the superfine ensemble (set 4 listed in Table III) as described in the text. 9The smallest binding of a below threshold b̄b̄bb tetraquark from the phenomenological studies (which are shown in Fig. 17) was 108 MeV [9]. HUGHES, EICHTEN, and DAVIES PHYS. REV. D 97, 054505 (2018) 054505-8 If one is searching for a bound state below threshold then both single- and two-particle contributions appear in the correlator and one may use the single-particle effective mass formula given in Eq. (18) in order to highlight the bound ground-state mass (as done in the mock data above). However if no bound state exists in the data, which we do not know a priori, then not removing the 1=t 3 2 dependence from the two-particle contributions has an important effect: an additional factor of 1.5 logð1 þ 1=tÞ is introduced into the effective mass formula in Eq. (18), giving a contami- nation that vanishes slowly as t → ∞. This would produce a confusing picture of what is actually contributing. The two-particle effective mass formula Eq. (19) removes this contribution. In the results to be reported now, we overlay both single- and two-particle effective masses on the same plot for the reader’s convenience. We generate the b̄b̄bb correlator data for the operators given in Table V using the ensembles listed in Table III and fit these data simultaneously with the bottomonium meson data to include correlations between data sets. All the b̄b̄bb data within a specific irrep and those which are unrelated by a Fierz relation10 are fit using Eq. (7) for the two-particle contributions and Eq. (6) for a hypothetical tetraquark state below threshold. The mean of the prior energy of the 2ηb, ηb þ Υ and 2Υ thresholds is roughly estimated based on the effective masses and then given a suitably wide prior width of 100 MeV while a tetraquark state prior energy is taken to be 250(100) MeV below each threshold. As can be seen in Fig. 5, since the data plateau to the noninteracting 2ηb threshold, no energy eigenstate is found below this thresh- old and variations of the tetraquark prior energy are insignificant. Similar behavior is seen with the other quantum numbers. (a) (c) (d) (b) FIG. 5. The b̄b̄bb effective masses for the 0þþ ðM1; M1Þ → ðM2; M2Þ correlators where M1, M2 are the ηb or Υ. Eeff and Eeff;t are the single- and two-particle effective masses defined in Eqs. (18) and (19) respectively. The mesons are separated by a distance rx in the x-direction when constructing the two-meson interpolating operator as given in Eq. (4). Gray points are not used when fitting the data (cf., Appendix A). 10Simultaneously fitting data sets related by a Fierz identity would mean the correlation matrix would have a zero eigenvalue and thus not be invertable for use in a least-squares minimization. SEARCHING FOR BEAUTY-FULLY BOUND TETRAQUARKS … PHYS. REV. D 97, 054505 (2018) 054505-9 Again, while we fit the correlator data in order to extract particle energies, so that the reader can visualize these data, we display effective mass plots on the different ensembles. The superfine ensemble (set 4) 0þþ two-meson effective masses are shown in Figs. 5, while the 0þþ diquark- antidiquark effective masses are given in Fig. 6, the physical coarse (set 2) 1þ− two-meson and diquark-anti- diquark effective masses in Fig. 7 and the fine (set 3) 2þþ two-meson and diquark-antidiquark effective masses sub- duced into the T2 lattice irrep are shown in Fig. 8. Each plot has the fitted ground-state energy overlaid in black for comparison. The 2þþ subduced into the E irrep is similar to the T2 case. Further, the behavior of the lattice data on all ensembles is qualitatively similar to those shown. The extracted ground-state energies in each channel are given in Table VI and a comparison of the energies is shown in graphical form in Fig. 9. It should also be noted that the numerical value of the effective mass (shift) plateau in two-hadron correlators has been shown, in certain cases, to be sensitive to the choice of interpolating operators [51,52]. There, the authors found that when the noise growth in the correlator data restricts the study of effective masses to a maximum propagation time of approximately 2 fm, fake plateaus can appear. These fake plateaus can be a consequence of different choices of source and sink (smeared) operators: this can cause a negative sign in the Z1 term in the effective mass formula Eq. (17) and a dip below threshold can appear for a short time range which can be misinterpreted as a bound state. Wall sources were shown to be particularly prone to this behavior of producing a “false dip” and obtaining an appreciably different plateau than a Gaussian source. Here we only use local quark sources. In addition, the elastic scattering states can also have a depend- ence on the choice of operator, which can cause a slow decay to the ground state andmimic a slowly varying effective mass that can be mistaken for a plateau over a short time range. As noted in these studies, a necessary check for a real effective mass plateau when using different source and sink operators istheconvergenceofalldatatoasingleplateauattimeslarger than approximately 2 fm. As we separate the operators by rx ¼ 0, 1 and 2 lattice units (as described in Sec. III C) and propagate to t > 8 fm, this is a consistency check we satisfy. A few notable features of the b̄b̄bb effective mass plots are evident. First and foremost, no value of the FIG. 6. The b̄b̄bb effective masses for the 0þþ diquark- antidiquark 3̄ × 3̄ and 6 × 6̄ correlators. Eeff and Eeff;t are the single- and two-particle effective masses defined in Eqs. (18) and (19) respectively. The diquarks are separated by a distance rx in the x-direction when constructing the diquark-antidiquark inter- polating operator as given in Eq. (4). FIG. 7. The b̄b̄bb effective masses for the 1þ− Υηb and diquark-antidiquark 3̄ × 3̄ correlators. Eeff and Eeff;t are the single- and two-particle effective masses defined in Eqs. (18) and (19) respectively. HUGHES, EICHTEN, and DAVIES PHYS. REV. D 97, 054505 (2018) 054505-10 effective mass is observed below the lowest noninteracting bottomonium-pair threshold in any channel, in line with what one would expect if no stable tetraquark candidate existed below threshold. Indeed, the b̄b̄bb effective mass plots are strikingly similar to the upper dot-dashed curve in the mock data in Fig. 4 where no bound tetraquark state is present. Additionally, the 2ηb → 2ηb effective mass shown in Fig. 5(a) plateaus very early due to the larger overlap onto the 2ηb threshold, while the 0 þþ 2Υ → 2Υ effective mass shown in Fig. 5(d) falls more slowly to the 2ηb threshold due to the larger overlap onto the nearby 2Υ threshold. The cross correlators 2ηb → 2Υ show how the different operators converge to a single plateau at a time greater than t ≈ 4 fm, a necessity for a true plateau as discussed above. The local diquark-antidiquark 0þþ corre- lator data are a linear combination of the 2ηb and 2Υ two- meson data, related by the Fierz identities given in Table II, and the effective masses shown in Fig. 6 reflect this. It is empirically observed that separating the diquark from the antidiquark by too large a distance rx results in larger noise due to the separation of color sources. The 1þ− Υηb and diquark-antidiquark effective masses are shown in Fig. 7, where the noise starts to increase after t ≈ 7 fm due to the Parisi-Lepage argument mentioned above with the noise being set by at least the 2ηb threshold. Based on this, one would also expect the signal to noise to be worse for the 2þþ data, which is also evident from the correlator data subduced into the T2 irrep shown in Fig. 8. As set 2 has physical ml=ms corresponding to a pion mass of Oð131Þ MeV, while the other ensembles have nonphysical ml=ms corresponding to pion masses of Oð300Þ MeV [53], no sensitivity to light sea quarks is observed. This would be expected from the smallness of the Van-der-Waals potential generated by the two-pion exchange between two 1S bottomonium mesons [54]. As can be seen in Figs. 8 and 9, the 2þþ ground state obtained from the lattice is slightly higher than that of the noninteracting threshold. However, this is the state which has the largest signal to noise, restricting the data to shorter time regions and it is possible that we are sensitive to the same aforementioned issue of a slowly varying fake plateau. Alternatively, this positive shift in the two-particle energy could potentially indicate appreciable infinite-volume continuum scattering arising from finite-volume interactions [44], but quantify- ing these phase shifts is outside the remit of this study. Regardless, these effects do not indicate that a bound tetraquark state exists in this channel. For illustration purposes, we also show the effective masses of the individual Wick contractions contributing to the 2ηb → 2ηb correlator in Fig. 10. As is evident, in each individual Wick contraction the effective mass drops below the 2ηb threshold but then rises slowly to threshold. However, importantly, when all Wick contractions are added together to yield the full correlator [shown Fig. 5(a)] the effective mass falls rapidly to threshold from above. This behavior is discussed further in Sec. VI. After analyzing all our data, as hinted by the effective mass plots, there is no indication of a bound tetraquark state below the noninteracting thresholds on any ensemble, as FIG. 8. The b̄b̄bb effective masses for the 2þþ ΥΥ and diquark- antidiquark 3̄ × 3̄ correlators subduced into the T2 irrep. E eff and Eeff;t are the single- and two-particle effective masses defined in Eqs. (18) and (19) respectively. TABLE VI. The ground-state static masses extracted from the lattice b̄b̄bb correlator data as described in the text. Set aMηb aMΥ aM0þþ aM1þ− aM2þþT2 aM2þþE 1 0.25548(13) 0.29180(22) 0.5121(7) 0.5500(12) 0.5840(22) 0.5863(29) 2 0.25741(19) 0.29365(39) 0.5162(12) 0.5534(21) 0.5921(20) 0.5911(28) 3 0.26570(8) 0.29288(13) 0.5321(16) 0.5594(11) 0.5899(21) 0.5888(18) 4 0.43288(8) 0.45209(12) 0.8658(10) 0.8865(23) 0.9083(36) 0.9079(35) SEARCHING FOR BEAUTY-FULLY BOUND TETRAQUARKS … PHYS. REV. D 97, 054505 (2018) 054505-11 shown in Fig. 9. We see no evidence of any change in the ground-state energies (with respect to the thresholds) as we vary lattice spacing or sea quark masses. Searching for a new tetraquark candidate has at least a two-dimensional parameter space: the hypothetical state would have an energy and also an overlap onto a specific operator. Using the lattice data presented here, we can determine a relationship between these parameters. Assuming that a tetraquark does exist below the lowest bottomonium-pair threshold in our data, at a certain time t� the correlator can be modeled with a two-state ansatz. Specifically for the 0þþ channel, given that the tetraquark has an energy E4b and an overlap Z4b onto a particular operator, at a large enough time the only other appreciable contribution will come from the higher 2ηb threshold which has an overlap Z2ηb with the same operator. In this case the correlator is given by Cðt�Þ ¼ Z24be−aE4bt � þ Z22ηbe −aE2ηb t � � aMηb 4πt� �3 2 : ð20Þ Using this ansatz in the effective mass formula Eq. (16) and rewriting the equation in terms of the nonperturbative overlaps yields the constraint Z24b Z22ηb ¼ � 1 − ð t�t�þ1Þ 3 2 expðaEeffðt�Þ − aE2ηbÞ expðaEeffðt�Þ − aE4bÞ − 1 � × e−aΔEt � � aMηb 4πt� �3 2 ð21Þ with aΔE ¼ aE2ηb − aE4b > 0. As can be seen, if the tetraquark is not observed by a time t� then the overlap onto this new state must be (at least) exponentially suppressed with the binding of the tetraquark state, e.g., with −ΔE. This point illustrates that if a tetraquark did exist with E4b < E2ηb then it is possible that it was not observed in our data because Z4b ≈ 0 within statistical precision. In this scenario, we can use the constraint Eq. (21) to estimate an upper bound on the magnitude of the overlaps given that no clear evidence of the tetraquark is observed within our statistical precision. The needed inputs for the constraint include the value of t� where the two-state ansatz is valid, aEeffðt�Þ from correlator data constructed with a specific operator, as well as aE2ηb ¼ 2aMηb. For the local Oðηb;ηbÞ operator on the a ¼ 0.06 fm ensemble, by examining Fig. 5(a), a choice of t� ¼ 143 ensures that the two-state ansatz is valid (given the long plateau at the 2ηb threshold). Here, aEeffðt� ¼ 143Þ ¼ 0.87634ð61Þ can also be precisely obtained. The value of aMηb, given in Table VI, is found from the ηb-meson data. Then, using this data in the constraint, for a certain choice of E4b, a numerical value FIG. 9. A summary of the b̄b̄bb ground-state energies with the lowest noninteracting bottomonium-pair threshold subtracted, across the different lattice ensembles listed in Table III, statistical error only. Note, as shown in Table III, fewer configurations were used on the a ¼ 0.06 fm ensemble than on the others. FIG. 10. The individual Wick contraction effective masses of the 2ηb → 2ηb correlators. The upper figure is the Direct1 and the lower is the Xchange2 contraction [each shown diagrammatically inFigs. 1(a)and 1(b)]. Eeff and Eeff;t arethesingle-and two-particle effective masses defined in Eqs. (18) and (19) respectively. HUGHES, EICHTEN, and DAVIES PHYS. REV. D 97, 054505 (2018) 054505-12 of the ratio of overlaps is found such that it is consistent with 0 within its small 1σ statistical error. Any value of the ratio of overlaps larger than this 1σ error is inconsistent, at this level of confidence, with our data observing a tetra- quark at this value of E4b. We use this model to estimate how small the hypo- thetical tetraquark overlap would need to be so that the tetraquark was not observed within our statistical pre- cision. A 1σ, 3σ and 5σ exclusion plot of the parameter space is given in Fig. 11. As the input data into the constraint has a long propagation time past t� > 8 fm and a statistically precise value of aEeffðt�Þ which does not fall below the threshold, a significant amount of param- eter space is excluded. It should be understood that this figure is only valid for a particular operator in a certain channel. The given 0þþ channel in Fig. 11 excludes the largest amount of parameter space as it is the most statistically precise. Also Fig. 11 is constructed from data on the superfine ensemble alone, where discretiza- tion effects are smallest and would not change the quantitative behavior significantly. To conclude this section, we find no evidence of a stable tetraquark candidate below the noninteracting thresholds by studying a full S-wave color-spin basis of QCD operators. In the next section we perform an exploratory and complementary study of an alternative approach to ensure the robustness of our conclusions. V. NRQCD WITH A HARMONIC OSCILLATOR POTENTIAL A stable tetraquark state in the 0þþ, 1þ− and 2þþ channels would overlap with the full basis of S-wave color-spin operators utilized above. In this section, we go one additional step by exploring an alternative approach in order to further investigate the possibility of a tetraquark state. Adding a central confining potential to the quark interactions can produce a more deeply bound tetraquark relative to the threshold as the strength of this interaction is increased (as we see below). Furthermore, an appropriate choice of additional interaction can reduce the fiducial volume of the lattice and thus thin the allowed discrete momenta states of the two-meson degrees of freedom. Adding an external attractive scalar central potential to the QCD interactions yields these desired effects. The harmonic oscillator potential is a particularly suitable choice of scalar interaction between quarks. For a particle of mass m at position x away from the center x0 the potential is just κr2=2 ≡ κjx − x0j2=2. Defining ω ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffi ðκ=mÞ p , the ground-state energy and wave function are E0 ¼ 3ω=2 and ψðrÞ ¼ C exp ð−mωr2=2Þ. Additionally, the separabil- ity of the combined QCD and harmonic oscillator potential into total and relative coordinates ensures that solutions of multiquark systems can be split into two parts with the total coordinate piece analytically solvable. This follows from the nature of the harmonic oscillator potential.11 Thus we expect that for small values of ω the ground state for 2ηb mesons is approximately 2MηbðωÞ þ 3ω. The 3ω term comes from two color-singlet ηb mesons in the harmonic oscillator potential and the mass of the ηb is shifted slightly from the value at ω ¼ 0 because of the additional harmonic oscillator inter- action combined with the QCD interactions that bind the two heavy quarks into a ηb. 12 However, for a compact tetraquark state the mass would be Mðb̄b̄bbÞðωÞ þ 3ω=2, as there is only one color-singlet state in the central harmonic oscillator potential (3ω=2) and if the tetraquark state is also a compact state (on the scale of ηb and much less than the effective lattice volume), then its mass will also receive only a modest positive correction due to ω. Hence if there were a tetraquark state near threshold then this additional interaction could drive it further below threshold, giving a much cleaner and distinct signal for the tetraquark candidate in our calculation. As the potential model framework describing the ηb has been largely successful, we can use this framework as a general guide for the exploratory nonperturbative lattice calculation when including the harmonic oscillator potential. Of course, we are mainly interested in QCD (and not the harmonic oscillator) and, as such, if we do find a stable tetraquark state when including the harmonic oscillator potential then we must take the ω → 0 limit. Thus, the objective of this section is to determine if a stable tetraquark state exists when the quarks are exposed to an auxiliary potential (which could push the tetraquark increasingly lower than the threshold) and if it does, whether it will survive the QCD limit. FIG. 11. The excluded region for the ratio of tetraquark/2ηb overlaps, Z4b=Z2ηb, onto the Oðηb;ηbÞ operator, assuming a tetraquark with mass, E4b, lying below the 2ηb threshold, E2ηb. The red hashed region is excluded at 5σ by the data as described in the text. The 1σ − 3σ and 3σ − 5σ exclusion bands are also shown for reference. 11Defining r ¼ jx1 − x2j and Rcm ¼ ððx1 þ x2Þ=2 − x0Þ, we can separate κ½ðx1 − x0Þ2 þ ðx2 − x0Þ2�=2 into relative and center-of-mass coordinates ½ðκ=2Þr2 þ ð2κÞR2cm�=2. 12For sufficiently small ω the shift in ηb is directly related to the rms radius of the ηb state since from perturbation theory it is given by hηbðω ¼ 0Þjkr2=2jηbðω ¼ 0Þi. SEARCHING FOR BEAUTY-FULLY BOUND TETRAQUARKS … PHYS. REV. D 97, 054505 (2018) 054505-13 The harmonic oscillator potential is defined as13 δHHO ¼ mbω 2 2 jx − x0j2 ð22Þ where the quarks are pulled towards the fixed point x0 with a strength ω. We choose x0 to be the same position as the source. Intuitively, as the quarks/ηb’s start to propagate further away from the source, the harmonic oscillator potential pulls them closer together. In turn, this restricts the quarks/ηb’s to be in a certain volume. First, it is necessary to determine which volumes the quarks/ηb’s are confined into by the addition of the harmonic oscillator potential. The root mean square dis- tance, rrms, gives an indication of this. We determine the rrms of the ηb based on solutions of the Schrodinger equation using a Cornell potential.14 The Matthieu equation can describe the behavior of two free ηb’s in a harmonic oscillator potential on a periodic box, and the solutions of which can yield rrms for the 2ηb state. The results from such a calculation are plotted in Fig. 12. Based on this, in order to confine the quarks sufficiently so that the two ηb’s overlap, and also to study thedependence on ω, values of ω ¼ ð75; 150; 300; 350Þ MeV and ω ¼ ð75; 150; 350; 500Þ MeV are chosen for the lattice ensembles called set 1 and set 3 in Table III. In lattice units, the simulated values of κ=2 ¼ ðambÞðaωÞ2=2 are (0.0029,0.0117,0.0469,0.638) and (0.0011,0.0044,0.0240, 0.0489) respectively. The harmonic oscillator is implemented through a minor modification of the NRQCD evolution equations via e−a ~H ¼ � 1 − aδHHO 2l � l e−aH � 1 − aδHHO 2l � l ð23Þ where e−aH is the purely NRQCD evolution equation defined in Eq. (12). This implementation was chosen so that the evolution equation is still time-reversal symmetric. Here, l is a stability parameter akin to n in Eq. (12) which is used to prevent possible numerical instabilities [25]. Values of l ¼ 13 and 10 were chosen for the calculations on set 1 and 3 respectively. Following these details, we are now able to present results from the nonperturbative lattice calculations. A. Numerical results All correlator data from set 1 and 3 discussed in Sec. IV were generated again with the inclusion of the harmonic oscillator potential at the four different ω values given above. The harmonic oscillator alters both the single- and two-particle contributions to the correlator so that they become (as derived in Appendix B) dependent on ω as CJ PC i;j ðt; ωÞ ¼ X n ZinZ j;� n ð1 þ e−2ωtÞ32 e−ðMðωÞnþ 3 2 ωÞt ð24Þ CJ PC i;j ðt; ωÞ ¼ X X2 ZiX2Z j;� X2 � 2ωμrπ −1 1 − e−4ωt �3 2 e−ðM S 1 ðωÞþMS 2 ðωÞþ3ωÞt þ � � � : ð25Þ First, Fig. 13 shows the effective masses, as defined in Eq. (16), of the ηb when including the harmonic oscillator. Also overlaid are the effective masses when removing the 1 þ e−2ωt dependence to enable a better comparison with the data when no harmonic oscillator is included. As can be seen, the dip in the harmonic oscillator effective masses is from this additional time dependence. Physically, this can be understood to be due to the b-quarks travelling non- relativistically and so it takes time for the harmonic oscillator to have an effect. FIG. 12. The root mean square distance rrms of the ηb and the 2ηb as a function of the harmonic oscillator strength calculated from a potential model with the different ensemble parameters listed in Table III as discussed in the text. FIG. 13. The effective mass plot for the ηb when including the harmonic oscillator potential on set 3. Eeff is given by Eq. (16) while Eeff;ξ removes the leading 1 þ e−2ωt dependence from the correlator (24) to enable a better comparison with the data when no harmonic oscillator potential is included. 13In a periodic box of length L the harmonic oscillator is also periodic. 14The potential form is − 4αs 3r þ rb2 with αs ¼ 0.36, b ¼ 2.34 GeV and reduced mass 2.59 GeV. HUGHES, EICHTEN, and DAVIES PHYS. REV. D 97, 054505 (2018) 054505-14 The ηb correlator data when including the harmonic oscillator potential are fit to the functional form given by Eq. (24) in order to extract the lowest energy eigenstate MðωÞηb þ 3ω=2 from the asymptotic behavior. We show the fitted result overlaid on the effective mass plot in Fig. 13. As before, the long plateau indicates that the ground state will be extracted accurately. To compare to the potential model predictions, we subtract the ηb mass with no harmonic oscillator included [Mðω ¼ 0Þηb] and then plot the energy differences against ω, as shown in Fig. 14. Good qualitative agreement between the lattice results and the potential model predictions is observed. Forthe b̄b̄bb system,weshow the 0þþ effective masseson set 3 in Fig. 15. It is evident that the 0þþ and the ηb data contain more noise when a harmonic oscillator potential is included. While fitting the data to the form in Eq. (25) can be performed, it is not necessary as the purpose of this FIG. 14. The lattice ηb energy MðωÞηb þ 3ω=2 when including the harmonic oscillator potential with Mðω ¼ 0Þηb subtracted compared to the model predictions as discussed in the text. (a) (b) (c) (d) FIG. 15. The b̄b̄bb effective masses for the 0þþ correlators when including the harmonic oscillator potential on set 3. Eeff is given by Eq. (16), while Eeff;‡ removes the leading 1 − e−4ωt dependence from the correlator (25) to enable a better comparison with the data when no harmonic oscillator potential is included. SEARCHING FOR BEAUTY-FULLY BOUND TETRAQUARKS … PHYS. REV. D 97, 054505 (2018) 054505-15 exploratory work is to determine if a stable tetraquark exists when ω ≠ 0. As can be seen, there is no fall below the 2ηb threshold for any valueofω. Similarbehavior isseen withthe data on set 1. We show the effective masses for the individual Direct1 and Xchange2 Wick contractions of the 2ηb → 2ηb corre- lator in Fig. 16. As before, the effective masses of the individual Wick contractions drop below the 2ηb threshold, even though importantly, when added together to yield the full correlator shown in Fig. 15(a) the effective mass is always above threshold. As this was also seen in the pure NRQCD data shown in Sec. IV, it may be a problematic feature of models that utilize a phenomenologically moti- vated four-body potential for this system. To conclude this section, despite adding an auxiliary potential into the QCD interactions that should push a near threshold tetraquark candidate increasingly lower we find no indication of any state below the 2ηb threshold. The conclusions of this section then agree with those of Sec. IV. VI. DISCUSSION AND CONCLUSIONS In this work we have studied the low-lying spectrum of the b̄b̄bb system using the first-principles lattice nonrelativistic QCD methodology in order to search for a stable tetraquark state below the lowest noninteracting bottomonium-pair threshold in three different channels: the 0þþ which couples to the 2ηb and 2Υ, the 1þ− which couples to Υηb and the 2þþ which couples to 2Υ. In Sec. III we describe our numerical methodology. Four gluon ensembles were employed with lattice spacings ranging from a ¼ 0.06–0.12 fm, and one ensemble which has physical light-quark masses. All ensem- bles have u, d, s and c quarks in the sea. In Sec. IV we presented the majority of the results in this work. Here, we determined the lowest energy eigenstate of the b̄b̄bb system with the quantum numbers 0þþ, 1þ− and 2þþ using an overconstrained S-wave color/spin basis (arising fromFierzrelations betweenthediquark-antidiquark and two-meson systems as shown in Table II). We did not observe any state below the lowest noninteracting bottomo- nium-pair threshold in any channel, as can be seen in Figs. 5–8, and a summary of our results from this section is given in Fig. 9. In Sec. V, to ensure the robustness of our conclusions, we performed an exploratory calculation of a novel method which added an auxiliary scalar potential into the QCD interactions with the objective of pushing a near threshold tetraquark increasingly lower than the threshold. This would give a more distinct and cleaner signal for its presence in our calculation. The harmonic oscillator was found to be a suitable central scalar potential. For the ηb-meson with this potential, we first verified agreement between the non- perturbative lattice calculations and a potential model (as shown in Fig. 14) and then used this potential model as a general guide to choose multiple appropriate values of the potential strength. Despite studying the b̄b̄bb system with this additional scalar potential on the lattice, no indication of the QCD tetraquark was observed. This work is the only first-principles study of the low- lying b̄b̄bb spectrum in the literature. However, there are others which utilize different methodologies. For example, Ref. [8] predicts the tetraquark mass by solving the two- particle Schrodinger equation with a phenomenologically motivated nonconfining potential between the pointlike diquark and antidiquark, and finds a 0þþ; 1þ− and 2þþ tetraquark to be bound by 44, 51 and 5 MeV respectively.15 The authors of [14] used a diquark model including a confining linear potential, but neglected spin effects, and found a 0þþ tetraquark to be bound by 48 MeV. However, it has also been found that the root mean square distance FIG. 16. The effective masses for the individual 2ηb → 2ηb Wick contraction correlator data when including a harmonic oscillator potential on set 3. The upper figure is the Direct1 and the lower is the Xchange2 contraction [each shown diagram- matically in Figs. 1(a) and 1(b)]. 15A model in which the diquarks are taken to be fundamental particles cannot determine the two-meson threshold from the Schrodinger equation because the diquarks cannot recombine into mesons. Thus, the experimental meson masses [4] are used to determine the lowest threshold. HUGHES, EICHTEN, and DAVIES PHYS. REV. D 97, 054505 (2018) 054505-16 between the diquark and antidiquark inside the tetraquark in this model is similar in magnitude to the distance between the quarks inside each diquark [9]. Consequently, such an approach is internally inconsistent. More recently, [10] used a Hamiltonian including only spin-spin interactions medi- ated by a one-gluon exchange. Here, all other effects such as the chromoelectric interactions, color confinement and b- quark mass need to be set separately. The authors set these additional contributions in two ways: by estimating the effects using an effective heavy quark [10] or by using the experimental meson masses as input. In this way, the authors find that the tetraquark state could be either below the 2ηb threshold or lie in between the 2ηb and 2Υ thresholds (where in both cases the thresholds were determined using the experimental meson masses). In [55] the author also uses a model including only chromomagnetic interactions and finds an unbound tetraquark. Within the QCD sum-rules framework, [11] finds a tetraquark candidate approximately 300 MeV below the experimental 2ηb threshold while [12] finds a tetraquark lying in between the experimental 2ηb and 2Υ thresholds. Using phenomenological arguments, [13] also finds a tetraquark candidate lying in between the thresholds. Indeed, in the limit of very heavy quarks where the force proceeds through one-gluon exchange containing only color-Coulomb contributions (safely neglecting spin and long-distance effects), the authors of [7] used a variational methodology to determine that a bound tetra- quark exists for a QQQ̄0Q̄0 system when mQ=mQ0 ≲ 0.15 (where both mQ and mQ0 are heavy relative to ΛQCD). However, if mQ=mQ0 is varied and the tetraquark becomes unbound, then as the free two-meson eigenstate becomes the ground state of the system this numerical methodology has an increasingly slow convergence to a solution [7] (being numerically ill posed) due to a redundant degree of freedom in the minimization procedure. Indeed [7] indicates that for mQ=mQ0 ¼ 1 the solution is unstable, a hint that no bound tetraquark exists for all identical quarks in the very heavy mass limit. The authors of [14] assume that the b-quark is sufficiently heavy to use the one-gluon exchange with only color-Coulomb contributions, and by also neglecting the mixing between different color components of the 2 × 2 potential matrix, find a tetraquark bound by 78(20) MeV (by using the experimental ηb mass to determine the threshold). In an orthogonal direction to the above work, the authors of [15] only include a linear string contribution in the one- gluon exchange (neglecting spin effects and the appreciable short-distance Coulomb contributions), and without mixing between the different color components of the potential matrix, find a bound tetraquark when mQ=mQ0 ¼ 1. However, in subsequent work [56,57], by modelling the aforementioned mixing they concluded that no bound tetraquark exists. Perhaps the most sophisticated non first-principles methodology used to study the four-body b̄b̄bb tetraquark is the diffusion Monte Carlo method utilized by [9]. Here, one determines the ground state of a phenomenologically motivated Hamiltonian by solving the Schrodinger equation and examining the stability ofP ne −ðEn−E0ÞtΨnðxÞ to determine E0, where En (Ψn) is the n-th energy-eigenstate (eigenfunction). The authors include boththe color-Coulomband linear contributionsin thegluon exchange but neglect the mixing between the different color components in the potential matrix, and find a stable tetraquark candidate 108 MeV below the 2ηb threshold (determined from the experimental meson mass). Consequently, there is no study in the literature which is not from first principles that includes all the appreciable effects relevant for the bbb̄b̄ system: treating the bottom quarks as fundamental particles, including both short- and long-distance effects in the gluon exchange and including the mixing between the different color components in the 2 × 2 potential matrix. It should be emphasized however that these studies, unlike ours, are not from first principles and thus have an unquantifiable systematic error associated with the choice of four-body potential. To emphasize this further, thinking of each Wick contraction (shown diagrammatically in Fig. 1) as a different potential contributing to the QCD dynamics, then only studying a subset of these interactions can change the energies of states. This can lead to the misidentification of a new state below threshold. For example, the effective masses of the individual Wick contractions contributing to the 2ηb → 2ηb correlator are shown in Fig. 10. As is evident there, in each individual Wick contraction the effective mass drops below the 2ηb threshold but rises slowly to threshold even though when all Wick contractions are added together to yield the full correlator [shown Fig. 5(a)] the effective mass falls rapidly to threshold from above. This behavior is even more pronounced in the data with the additional scalar potential, shown in Fig. 16, possibly indicating that this may be a problematic feature of models that utilize a phenomeno- logically motivated four-body potential: a subset of the interactions show behavior that may be misinterpreted as a bound state below threshold, while when all interactions are included no bound state is seen. Particularly, the slow rise to threshold from below could make the diffusion Monte Carlo method practically difficult due to the slowly varying stability condition combined with the fact that a long evolution time (greater than 8 fm) is necessary. In conclusion, we find no evidence of a b̄b̄bb tetraquark with a mass below the lowest noninteracting bottomonium- pair thresholds in the 0þþ, 1þ− or 2þþ channels. We give a constraint in Eq. (21) that future phenomenological models must satisfy if such QCD states are postulated. For the 0þþ channel, we use this constraint to estimate how small the nonperturbative overlap of the hypothetical tetraquark (onto a particular operator) would need to be, relative to the 2ηb, so that it was not observed within our statistical precision. A 1σ, 3σ and 5σ exclusion plot of the parameter space is shown in Fig. 11, and discussed in Sec. IV. As we SEARCHING FOR BEAUTY-FULLY BOUND TETRAQUARKS … PHYS. REV. D 97, 054505 (2018) 054505-17 have propagation times longer than 8 fm and statistically precise data, we can exclude all but the most finely tuned parameter space. Our lattice results then rule out the phenomenological models discussed above that predict a tetraquark below the lowest bottomonium-pair thresholds which have a value of nonperturbative overlap that is excluded by Fig. 11. A comparison of these results with ours is shown in Fig. 17. Further studies of possible heavy tetraquark channels that includeorbital angular momentum either between themesons in the tetraquark or between the quarks in the meson could be performed with the methodology used here.16 Similarly one could also study whether stable c̄c̄cc, b̄c̄bc or b̄b̄cc tetraquarks exist or not. Additionally, two-hadron systems receive a finite-volume energy shift which depends on the infinite-volume scattering amplitude which is nontrivial to parametrize. Here we do not calculate these finite-volume energy shifts. Doing so in a more extended study would allow statements to be made about the existence of resonant tetraquark states above the lowest thresholds, that likely do exist in nature. Quantifying these shifts would be an exciting avenue for future work. Finally, recent work based on heavy-quark symmetry [6] and phenomenological arguments [58] indicates that a JP ¼ 1þ b̄b̄ud tetraquark is stable in QCD. In fact, by extracting a potential from the lattice in the static heavy-quark limit and solving the Schrodinger equation [59–61] also find binding in this channel. Initial lattice calculations hint that such a state exists [62], but calcu- lations are difficult because of a signal-to-noise problem for heavy-light states [63]. Lattice QCD calculations in this direction are essential for a conclusive first-principles statement to be made and to give further motivation for a targeted experimental search for these tetraquark configurations of nature. ACKNOWLEDGMENTS We would like to thank William Bardeen and Zhen Liu for the many insightful discussions on tetraquarks, as well as Gavin Cheung who in addition gave guidance on the tetraquark implementation. We are also grateful to the MILC collaboration for the use of their gauge configura- tions. This manuscript has been authored by Fermi Research Alliance, LLC under Contract No. DE-AC02- 07CH11359 with the U.S. Department of Energy, Office of Science, Office of High Energy Physics. The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a nonexclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this manuscript, or allow others to do so, for United States Government purposes. The results described here were obtained using the Darwin Supercomputer of the University of Cambridge High Performance Computing Service as part of the DiRAC facility jointly funded by STFC, the Large Facilities Capital Fund of United Kingdom Department of Business, Innovation and Skills (BIS), and the Universities of Cambridge and Glasgow. APPENDIX A: TWO-POINT CORRELATOR FIT FUNCTIONS Here we derive the nonrelativistic two-particle contri- bution 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 a two-hadron system is [64] −600 −400 −200 0 200 E4b − E2ηb (MeV) This Calc. Lattice QCD Pheno. Models CLSZ17 BLO16 AFGSZ17 KNR17 W17 WLCZ16† BLN11 AFGSZ17 VVR13ξ FIG. 17. A comparison of our result for the b̄b̄bb ground-state energyinthe0þþ channel(statisticalerroronly)andpredictionsfrom phenomenological models. The hatched region indicates the ex- clusion of a bound tetraquark with an energy E4b subject to the value of its nonperturbative overlap as given in Fig. 11. In this comparison, we take our ground-state energy obtained on the superfine ensemble (set 4 in Table III) as a representative because it has the smallest discretization effects and the statistical error encompasses the results on the other ensembles as shown in Fig. 9. The y-axis labels results from different phenomenological models [7–13] by last initial of authors and year of publication. An error is plotted if given in the reference. The two results from WLCZ16† differ by how the mass scale was set. The result VVR13ξ finds no bound tetraquark and we indicate this by placing their result on threshold. 16It should be noted that although we focused on S-wave combinations of quarks, the channels we study also exclude certain combinationsoforbitalangularmomentumfromproducingabound tetraquark. For example, the 0þþ overlaps with 2Υ in an orbital angular momentum D-wave configuration. If this state produced a low-lying bound tetraquark itwould alsoshowupinourcalculation. HUGHES, EICHTEN, and DAVIES PHYS. REV. D 97, 054505 (2018) 054505-18 I ¼ X X2 Z d3Ptot ð2πÞ3 d3k ð2πÞ3 1 2EðX2Þ jX 2 ðPtot;kÞihX 2 ðPtot;kÞj ðA1Þ where jX2ðPtot;kÞi ¼ jM1ðkÞM2ðPtot − kÞi is a two-hadron state (with quantum numbers suppressed) and to avoid superfluous notation, we also set Ptot ¼ 0. A key difference from the one-hadron system is the internal relative momen- tum, k, which contributes an additional three-integral. Substituting the completeness relation Eq. (A1) into the correlator Eq. (5) and performing the momentum conserv- ing integrals yields CðtÞ ¼ X X2 Z d3k ð2πÞ3 ZX2ðkÞ 2e−EðX 2Þt ðA2Þ where ZX2ðkÞ is a nonperturbative coefficient. However, on a discrete finite volume the above integral over elastic states is replaced by a finite sum with km ∈ ð−Ns=2 þ 1; …; ; Ns=2Þ in units of ð2π=aNsÞ. In turn, for Ptot ¼ 0, Eq. (A2) becomes a sum over back-to-back hadronic states which have values of the discrete momenta that are equal in magnitude but opposite in direction. One can expand the two-particle energy using a nonrelativistic dispersion relation, appropriate since we are using NRQCD, as EðX2Þ ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi M21 þ jkj2 q þ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi M22 þ jkj2 q ðA3Þ ≈ MS1 þ MS2 þ jkj2 2μr ðA4Þ wherewe have defined the static, kinetic and reduced masses by MS, MK and μr ¼ MK1 MK2 =ðMK1 þ MK2 Þ respectively. In a finite volume there is an additional contribution to Eq. (A4) dependent on the infinite-volume scattering phase shift, which is discussed further below. Equation (A4) also illustrates the density of back-to-back states on our ensem- bles. As an example, examining the a ¼ 0.09 fm ensemble, and taking Mηb ¼ 9.399ð2Þ GeV from the PDG [4], the smallest allowed jkj2=2μr ≈ 20 MeV or 0.0092 in lattice units with all other back-to-back states separated by multi- ples of this value. Consequently, due to the bottomonium mass being large compared to the smallest allowed momen- tum, adjacent back-to-back states are sufficiently close in energy that fitting the momentum states as a discrete sum would require a vast set of correlators projected onto each separate jkj2=2μr (with the methodology used in [38]). Practically, this would be overly computationally expensive and instead, the fact that the states with k ≠ 0 are related by the dispersion relation (and are not independent as the sum would assume) should be included. This can be achieved by first expanding the nonpertur- bative coefficient ZX2ðkÞ as a polynomial in jkj2=μ2r, as dictated by rotational symmetry and by ensuring the Taylor coefficients have the same dimension, then keeping all terms needed to a certain precision. After this the correlator can be written as CðtÞ ¼ X X2 e−ðM S 1 þMS 2 Þt X k �X∞ i¼0 Z2l X2 jkj2l μ2lr � e− jkj2 2μr t ðA5Þ ¼ X X2 e−ðM S 1 þMS 2 Þt Z π a −πa d3k ð2πÞ3 �X∞ i¼0 Z2l X2 jkj2l μ2lr � e− jkj2 2μr t ðA6Þ where going from the first to the second line we have replaced the finite sum by an integral. Taking the limits of the integral to �∞ and performing the integrals over k analytically yields the fit function given in Eq. (7). Once it is shown that it is possible to replace the finite sum by the indefinite integral within our statistical precision then it is valid to use the above fit function with our data. To do so, using spherical coordinates in Eq. (A6), we define the quantities that we need to compare as IðlÞðtÞ ¼ 1 μ2lr Z ∞ −∞ djkjjkj2lþ2e− jkj2 2μr t; ðA7Þ DðlÞðtÞ ¼ 1 μ2lr X jkj jkj2lþ2e− jkj2 2μr t: ðA8Þ The integrands of both are shown diagrammatically in Fig. 18, where it is observed that due to the Gaussian time dependence the peaks of the integrand move towards the origin withlarger t. As such,one objective istochoosea large enough t̂ such that a sufficient majority of the integrand is within the maximum momentum π=a. We can replace the discrete finite-volume fit function with its infinite-volume 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) ensemble as discussed in the text. Due to the Gaussian time dependence, the integrand peak moves towards the origin for larger times. SEARCHING FOR BEAUTY-FULLY BOUND TETRAQUARKS … PHYS. REV. D 97, 054505 (2018) 054505-19 counterpart if the relative difference between them is less than our statistical errors. Specifically if ���� Plmax l¼0 Z 2;ðlÞIðlÞðtÞ − P∞l¼0 Z2;ðlÞDðlÞðtÞPlmax l¼0 Z 2;ðlÞIðlÞðtÞ ���� ðA9Þ ≤ j Plmaxl¼0 Z2;ðlÞIðlÞðtÞ − P∞l¼0 Z2;ðlÞDðlÞðtÞj jZ2;ð0ÞIð0Þj ðA10Þ ≤ Xlmax l¼0 Z2;ðlÞ Z2;ð0Þ jIðlÞðtÞ − DðlÞðtÞj Ið0Þ þ X∞ l¼lmaxþ1 Z2;ðlÞ Z2;ð0Þ DðlÞðtÞ Ið0Þ ≤ Xlmax l¼0 jIðlÞðtÞ − DðlÞðtÞj Ið0Þ þ X∞ l¼lmaxþ1 DðlÞðtÞ Ið0Þ ðA11Þ ≤ δCðtÞ CðtÞ ðA12Þ where lmax is the maximum number of moments to be included in the fit function, the inequality in the second line holds as the moment integrands are positive (shown dia- grammatically in Fig. 18), in the third line the Cauchy inequality has been used, and in the fourth line it is assumed that the leading moment gives the dominant contribution (Z2;ðlÞ ≤ Z2;ð0Þ). Studying Eq. (A11) instead of Eq. (A9) is a conservative option. Each part of the first term in Eq. (A11) represents how similar IðlÞðtÞ and DðlÞðtÞ need to be in order to be considered equivalent within statistical precision. This is shown in Fig. 19. For a particular lmax, the second term represents when the higher moments look like noise within statistical precision, also shown in Fig. 19. Each figure was generated with the coarse ensemble parameters (listed as set 1 in Table III) as this ensemble has the largest lattice spacing (and hence smallest π=a value—the upper limit on the integral of interest) and also the smallest Ns (the number of discrete momenta used in the finite-volume sum). As such, the other ensembles give better approx- imations and studying set 1 is conservative. Overlaid on each plot is the smallest relative statistical error from the data on any ensemble. Due to the constant signal-to-noise ratio, the number of configurations and the size of the lattice spacing, the smallest statistical error was the 2ηb correlator on the fine ensemble. Only examining situations below this curve is the most conservative option for all data generated. As can be observed in Fig. 19, the discrete finite- volume sums are well represented by the indefinite inte- grals. Additionally, in order to neglect the higher moments within our statistical precision, a choice of t̂ ¼ 1 fm and lmax ¼ 2 is sufficient. Expanding the finite-volume two-particle energy non- relativistically in Eq. (A4) neglected a possible finite-volume energy shift which depends on the infinite-volume scattering amplitude. In the small scattering-length limit, the energy shift is known to be volume suppressed [40]. The two- particle systems under study are in this limit as the ηb and Υ are compact due to the heavy-quark mass, with a size of 0.2– 0.3 fm. As such, the low-momentum energy shifts are not expected to be appreciable given the large volume ensembles we employ. Energy shifts in higher momentum states from Eq. (A4) are exponentially suppressed due to the Gaussian integral in Eq. (A6). Consequently, these too are not appreciable and no large influence of finite-volume energy shifts is seen (cf., the effective mass figures in Sec. IV). Quantifying these finite-volume scattering shifts is outside the remit of this study. Regardless, the scattering shifts would be positive and push the finite-volume two-particle energy higher and not contribute to a misidentification of a bound tetraquark below the noninteracting threshold. APPENDIX B: TWO-POINT CORRELATOR FIT FUNCTIONS WITH A HARMONIC OSCILLATOR In the nonrelativistic limit the free propagator (to leading order) is Δðx; tÞ ¼ Z d3p ð2πÞ3 exp ðip · xÞ exp � − � m þ p 2 2m � t � : ðB1Þ The free two-meson propagator with Ptot ¼ 0, both starting at common origin x0 ¼ ð0; 0Þ and ending at time t, is given by ~ΔðtÞ ¼ Z d3xΔ1ðx; tÞΔ2ðx; tÞ: ðB2Þ Using Eq. (B1) in Eq. (B2) produces the large-time behavior of the free two-meson propagator as FIG. 19. The difference between the discrete finite-volume and infinite-volume continuum moments (upper) and which moments can be neglected compared to our statistical precision (lower) as discussed in the text. HUGHES, EICHTEN, and DAVIES PHYS. REV. D 97, 054505 (2018) 054505-20 ~ΔðtÞ ¼ � μr 2πt � 3=2 e−ðm1þm2Þt: ðB3Þ ThisagreeswiththeleadingbehaviorderivedinAppendixA. Next, for the harmonic oscillator case, the one-dimensional Hamiltonian is ∂ψ ∂t ¼ Eψ ¼ − 1 2m ∂2ψ ∂x2 þ κ 2 x2ψ: ðB4Þ Solutions of this system can be related to the solutions of the Mehler differential equation via ∂ϕ ∂~t ¼ ∂2ϕ ∂ρ2 − ρ 2ϕ ðB5Þ with the identifications ω ¼ ffiffiffiffiffiffiffiffiffi κ=m p , t ¼ 2~t=ω and r ¼ ðκmÞ−1=4ρ. The Greens function (propagator) for the Mehler differential equation is given by Δðρ1; ρ2; ~tÞ ¼ 1ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2π sinhð2~tÞ p exp � − coth ð2~tÞ ðρ 2 1 þ ρ22Þ 2 þ cschð2~tÞρ1ρ2 � : ðB6Þ To normalize this propagator one can first compare the large t behavior of this solution with the known behavior of the harmonic oscillator propagator, limt→∞GðtÞ ¼ jΨð0Þj2e−12ωt, where the wave function at the origin is given by Ψð0Þ ¼ ðmω=πÞ1=4. As such, the harmonic oscillator solution is Δðx;0;tÞ ¼ ffiffiffiffiffiffiffi mω p ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2πsinhðωtÞ p exp � − mωx2 2 cothðωtÞ � : ðB7Þ The three-dimensional solution can then be obtained using the separability of each spatial direction, so that the zero spatial-momentum single-particle correlator in a harmonic oscillator potential is given by Z d3xΔðx; 0; tÞ ¼ � 1 coshðωtÞ � 3=2 e−mt: ðB8Þ Finally, the equal mass two-particle propagator starting at a common origin x0 ¼ ð0; 0Þ and ending at a time t with Ptot ¼ 0, in the presence of an external harmonic oscillator potential, is found from using Eq. (B7) in Eq. (B2), to give ~ΔðtÞ ¼ � mω 2π 1 sinhð2ωtÞ � 3=2 e−2mt: ðB9Þ By comparing (B9) to (B3), and noting that m ¼ 2μr, we see that the free two-meson harmonic oscillator propagator reduces to the nonharmonic oscillator case as ω → 0. [1] R. J. Jaffe, Phys. Rev. D 15, 267 (1977). [2] J. J. Dudek, R. G. Edwards, and D. J. Wilson (Hadron Spectrum Collaboration), Phys. Rev. D 93, 094506 (2016). [3] R. A. Briceno, J. J. Dudek, R. G. Edwards, and D. J. Wilson (Hadron Spectrum Collaboration), arXiv:1708.06667. [4] C. Patrignani et al. (Particle Data Group Collaboration), Chin. Phys. C 40, 100001 (2016). [5] N. Brambilla, S. Eidelman, B. Heltsley, R. Vogt, G. Bodwin et al. (Quarkonium Working Group Collaboration), Eur. Phys. J. C 71, 1534 (2011). [6] E. J. Eichten and C. Quigg, Phys. Rev. Lett. 119, 202002 (2017). [7] A. Czarnecki, B. Leng, and M. B. Voloshin, arXiv: 1708.04594. [8] A. V. Berezhnoy, A. V. Luchinsky, and A. A. Novoselov, Phys. Rev. D 86, 034004 (2012). [9] Y. Bai, S. Lu, and J. Osborne, arXiv:1612.00012. [10] J. Wu, Y.-R. Liu, K. Chen, X. Liu, and S.-L. Zhu, arXiv: 1605.01134. [11] W. Chen, H.-X. Chen, X. Liu, T. G. Steele, and S.-L. Zhu, Phys. Lett. B 773, 247 (2017). [12] Z.-G. Wang, Eur. Phys. J. C 77, 432 (2017). [13] M. Karliner, S. Nussinov, and J. L. Rosner, Phys. Rev. D 95, 034011 (2017). [14] M. N. Anwar, J. Ferretti, F.-K. Guo, E. Santopinto, and B.-S. Zou, arXiv:1710.02540. [15] J. Vijande, A. Valcarce, and J. M. Richard, Phys. Rev. D 76, 114013 (2007). [16] E. Eichten and Z. Liu, arXiv:1709.09605. [17] R. Vega-Morales, arXiv:1710.02738. [18] J. J. Dudek, R. G. Edwards, M. J. Peardon, D. G. Richards, and C. E. Thomas (Hadron Spectrum Collaboration), Phys. Rev. D 82, 034508 (2010). [19] G. T. Bodwin, E. Braaten, and G. P. Lepage, Phys. Rev. D 51, 1125 (1995); 55, 5853(E) (1997). [20] A. Bazavov et al. (MILC Collaboration), Phys. Rev. D 82, 074501 (2010). [21] A. Hart, G. M. von Hippel, and R. R. Horgan (HPQCD Collaboration), Phys. Rev. D 79, 074008 (2009). [22] E. Follana, Q. Mason, C. Davies, K. Hornbostel, G. P. Lepage, J. Shigemitsu, H. Trottier, and K. Wong (HPQCD Collaboration and UKQCD Collaboration), Phys. Rev. D 75, 054502 (2007). [23] R. J. Dowdall et al. (HPQCD Collaboration), Phys. Rev. D 85, 054509 (2012). SEARCHING FOR BEAUTY-FULLY BOUND TETRAQUARKS … PHYS. REV. D 97, 054505 (2018) 054505-21 https://doi.org/10.1103/PhysRevD.15.267 https://doi.org/10.1103/PhysRevD.93.094506 https://doi.org/10.1103/PhysRevD.93.094506 http://arXiv.org/abs/1708.06667 https://doi.org/10.1088/1674-1137/40/10/100001 https://doi.org/10.1140/epjc/s10052-010-1534-9 https://doi.org/10.1140/epjc/s10052-010-1534-9 https://doi.org/10.1103/PhysRevLett.119.202002 https://doi.org/10.1103/PhysRevLett.119.202002 http://arXiv.org/abs/1708.04594 http://arXiv.org/abs/1708.04594 https://doi.org/10.1103/PhysRevD.86.034004 http://arXiv.org/abs/1612.00012 http://arXiv.org/abs/1605.01134 http://arXiv.org/abs/1605.01134 https://doi.org/10.1016/j.physletb.2017.08.034 https://doi.org/10.1140/epjc/s10052-017-4997-0 https://doi.org/10.1103/PhysRevD.95.034011 https://doi.org/10.1103/PhysRevD.95.034011 http://arXiv.org/abs/1710.02540 https://doi.org/10.1103/PhysRevD.76.114013 https://doi.org/10.1103/PhysRevD.76.114013 http://arXiv.org/abs/1709.09605 http://arXiv.org/abs/1710.02738 https://doi.org/10.1103/PhysRevD.82.034508 https://doi.org/10.1103/PhysRevD.82.034508 https://doi.org/10.1103/PhysRevD.51.1125 https://doi.org/10.1103/PhysRevD.51.1125 https://doi.org/10.1103/PhysRevD.55.5853 https://doi.org/10.1103/PhysRevD.82.074501 https://doi.org/10.1103/PhysRevD.82.074501 https://doi.org/10.1103/PhysRevD.79.074008 https://doi.org/10.1103/PhysRevD.75.054502 https://doi.org/10.1103/PhysRevD.75.054502 https://doi.org/10.1103/PhysRevD.85.054509 https://doi.org/10.1103/PhysRevD.85.054509 [24] B. Chakraborty, C. T. H. Davies, B. Galloway, P. Knecht, J. Koponen, G. C. Donald, R. J. Dowdall, G. P. Lepage, and C. McNeile (HPQCD Collaboration), Phys. Rev. D 91, 054508 (2015). [25] G. P. Lepage, L. Magnea, C. Nakhleh, U. Magnea, and K. Hornbostel, Phys. Rev. D 46, 4052 (1992). [26] J. O. Daldrop, C. T. H. Davies, and R. J. Dowdall (HPQCD Collaboration), Phys. Rev. Lett. 108, 102003 (2012). [27] R. J. Dowdall, C. T. H. Davies, T. Hammant, R. R. Horgan, and C. Hughes (HPQCD Collaboration), Phys. Rev. D 89, 031502 (2014); 92, 039904(E) (2015). [28] R. J. Dowdall, C. T. H. Davies, T. C. Hammant, and R. R. Horgan (HPQCD Collaboration), Phys. Rev. D 86, 094510 (2012). [29] R. J. Dowdall, C. T. H. Davies, R. R. Horgan, C. J. Monahan, and J. Shigemitsu (HPQCD Collaboration), Phys. Rev. Lett. 110, 222003 (2013). [30] B. Colquhoun, R. J. Dowdall, C. T. H. Davies, K. Hornbostel, and G. P. Lepage (HPQCD Collaboration), Phys. Rev. D 91, 074514 (2015). [31] C. Hughes, R. J. Dowdall, C. T. H. Davies, R. R. Horgan, G. von Hippel, and M. Wingate (HPQCD Collaboration), Phys. Rev. D 92, 094501 (2015). [32] T. C. Hammant, A. G. Hart, G. M. von Hippel, R. R. Horgan, and C. J. Monahan (HPQCD Collaboration), Phys. Rev. D 88, 014505 (2013). [33] C. J. Morningstar, Phys. Rev. D 50, 5902 (1994). [34] D. C. Moore and G. T. Fleming, Phys. Rev. D 73, 014504 (2006). [35] C. E. Thomas, R. G. Edwards, and J. J. Dudek (Hadron Spectrum Collaboration), Phys. Rev. D 85, 014507 (2012). [36] J. Vijande and A. Valcarce, Symmetry 1, 155 (2009). [37] E. Berkowitz, T. Kurth, A. Nicholson, B. Joo, E. Rinaldi, M. Strother, P. M. Vranas, and A. Walker-Loud (CalLat Collaboration), Phys. Lett. B 765, 285 (2017). [38] J. J. Dudek, R. G. Edwards, and C. E. Thomas (Hadron Spectrum Collaboration), Phys. Rev. D 86, 034031 (2012). [39] G. P. Lepage, B. Clark, C. T. H. Davies, K. Hornbostel, P. B. Mackenzie, C. Morningstar, and H. Trottier, Nucl. Phys. B, Proc. Suppl. 106, 12 (2002). [40] G. Parisi, Phys. Rep. 103, 203 (1984). [41] G. P. Lepage, Report No. CLNS-89-971. [42] M. Lüscher, Commun. Math. Phys. 104, 177 (1986). [43] R. J. Dowdall, C. T. H. Davies, G. P. Lepage, and C. McNeile (HPQCD Collaboration), Phys. Rev. D 88, 074504 (2013). [44] M. Lüscher, Nucl. Phys. B354, 531 (1991). [45] R. A. Briceno and Z. Davoudi, Phys. Rev. D 88, 094507 (2013). [46] R. A. Briceno, Phys. Rev. D 89, 074507 (2014). [47] P. Guo, J. Dudek, R. Edwards, and A. P. Szczepaniak, Phys. Rev. D 88, 014501 (2013). [48] M. Gockeler, R. Horsley, M. Lage, U. G. Meissner, P. E. L. Rakow, A. Rusetsky, G. Schierholz, and J. M. Zanotti, Phys. Rev. D 86, 094513 (2012). [49] C. h. Kim, C. T. Sachrajda, and S. R. Sharpe, Nucl. Phys. B727, 218 (2005). [50] S. He, X. Feng, and C. Liu, J. High Energy Phys. 07 (2005) 011. [51] T. Iritani et al. (HAL QCD Collaboration), J. High Energy Phys. 10 (2016) 101. [52] T. Yamazaki, K.-I. Ishikawa, Y. Kuramashi, and A. Ukawa (PACS Collaboration), Proc. Sci., LATTICE2016 (2017) 108. [53] A. Bazavov et al. (MILC Collaboration), Phys. Rev. D 87, 054505 (2013). [54] N. Brambilla, G. Krein, J. T. Castell, and A. Vairo, Phys. Rev. D 93, 054002 (2016). [55] B. Silvestre-Brac, Phys. Rev. D 46, 2179 (1992). [56] J. Vijande, A. Valcarce, and J.-M. Richard, Phys. Rev. D 87, 034040 (2013). [57] J.-M. Richard, A. Valcarce, and J. Vijande, Phys. Rev. D 95, 054019 (2017). [58] M. Karliner and J. L. Rosner, Phys. Rev. Lett. 119, 202001 (2017). [59] P. Bicudo and M. Wagner (European Twisted Mass Collaboration), Phys. Rev. D 87, 114511 (2013). [60] P. Bicudo, K. Cichy, A. Peters, B. Wagenbach, and M. Wagner, Phys. Rev. D 92, 014507 (2015). [61] P. Bicudo, M. Cardoso, A. Peters, M. Pflaumer, and M. Wagner, Phys. Rev. D 96, 054510 (2017). [62] A. Francis, R. J. Hudspith, R. Lewis, and K. Maltman, Phys. Rev. Lett. 118, 142001 (2017). [63] C. T. H. Davies, C. McNeile, E. Follana, G. P. Lepage, H. Na, and J. Shigemitsu (HPQCD Collaboration), Phys. Rev. D 82, 114504 (2010). [64] M. Srednicki, Quantum Field Theory, 1st ed. (Cambridge University Press, Cambridge, 2007). HUGHES, EICHTEN, and DAVIES PHYS. REV. D 97, 054505 (2018) 054505-22 https://doi.org/10.1103/PhysRevD.91.054508 https://doi.org/10.1103/PhysRevD.91.054508 https://doi.org/10.1103/PhysRevD.46.4052 https://doi.org/10.1103/PhysRevLett.108.102003 https://doi.org/10.1103/PhysRevD.89.031502 https://doi.org/10.1103/PhysRevD.89.031502 https://doi.org/10.1103/PhysRevD.92.039904 https://doi.org/10.1103/PhysRevD.86.094510 https://doi.org/10.1103/PhysRevD.86.094510 https://doi.org/10.1103/PhysRevLett.110.222003 https://doi.org/10.1103/PhysRevLett.110.222003 https://doi.org/10.1103/PhysRevD.91.074514 https://doi.org/10.1103/PhysRevD.91.074514 https://doi.org/10.1103/PhysRevD.92.094501 https://doi.org/10.1103/PhysRevD.92.094501 https://doi.org/10.1103/PhysRevD.88.014505 https://doi.org/10.1103/PhysRevD.88.014505 https://doi.org/10.1103/PhysRevD.50.5902 https://doi.org/10.1103/PhysRevD.73.014504 https://doi.org/10.1103/PhysRevD.73.014504 https://doi.org/10.1103/PhysRevD.85.014507 https://doi.org/10.1103/PhysRevD.85.014507 https://doi.org/10.3390/sym1020155 https://doi.org/10.1016/j.physletb.2016.12.024 https://doi.org/10.1103/PhysRevD.86.034031 https://doi.org/10.1103/PhysRevD.86.034031 https://doi.org/10.1016/S0920-5632(01)01638-3 https://doi.org/10.1016/S0920-5632(01)01638-3 https://doi.org/10.1016/0370-1573(84)90081-4 https://doi.org/10.1007/BF01211589 https://doi.org/10.1103/PhysRevD.88.074504 https://doi.org/10.1103/PhysRevD.88.074504 https://doi.org/10.1016/0550-3213(91)90366-6 https://doi.org/10.1103/PhysRevD.88.094507 https://doi.org/10.1103/PhysRevD.88.094507 https://doi.org/10.1103/PhysRevD.89.074507 https://doi.org/10.1103/PhysRevD.88.014501 https://doi.org/10.1103/PhysRevD.88.014501 https://doi.org/10.1103/PhysRevD.86.094513 https://doi.org/10.1103/PhysRevD.86.094513 https://doi.org/10.1016/j.nuclphysb.2005.08.029 https://doi.org/10.1016/j.nuclphysb.2005.08.029 https://doi.org/10.1088/1126-6708/2005/07/011 https://doi.org/10.1088/1126-6708/2005/07/011 https://doi.org/10.1007/JHEP10(2016)101 https://doi.org/10.1007/JHEP10(2016)101 https://doi.org/10.1103/PhysRevD.87.054505 https://doi.org/10.1103/PhysRevD.87.054505 https://doi.org/10.1103/PhysRevD.93.054002 https://doi.org/10.1103/PhysRevD.93.054002 https://doi.org/10.1103/PhysRevD.46.2179 https://doi.org/10.1103/PhysRevD.87.034040 https://doi.org/10.1103/PhysRevD.87.034040 https://doi.org/10.1103/PhysRevD.95.054019 https://doi.org/10.1103/PhysRevD.95.054019 https://doi.org/10.1103/PhysRevLett.119.202001 https://doi.org/10.1103/PhysRevLett.119.202001 https://doi.org/10.1103/PhysRevD.87.114511 https://doi.org/10.1103/PhysRevD.92.014507 https://doi.org/10.1103/PhysRevD.96.054510 https://doi.org/10.1103/PhysRevLett.118.142001 https://doi.org/10.1103/PhysRevLett.118.142001 https://doi.org/10.1103/PhysRevD.82.114504 https://doi.org/10.1103/PhysRevD.82.114504