Importance of meson-meson and of diquark-antidiquark creation operators for a \bar{b} \bar{b} u d tetraquark
IImportance of meson-meson and of diquark-antidiquark creation operators for a ¯ b ¯ bud tetraquark (1) Pedro Bicudo, ∗ (2) Antje Peters, † (3) Sebastian Velten, ‡ and (3) , (4) Marc Wagner § (1) CFTP, Dep. Física, Instituto Superior Técnico,Universidade de Lisboa, Av. Rovisco Pais, 1049-001 Lisboa, Portugal (2)
Westfälische Wilhelms-Universität Münster, Institut für Medizinische Psychologie und Systemneurowissenschaften,Von-Esmarch-Straße 52, D-48149 Münster, Germany (3)
Johann Wolfgang Goethe-Universität Frankfurt am Main, Institut für Theoretische Physik,Max-von-Laue-Straße 1, D-60438 Frankfurt am Main, Germany and (4)
Helmholtz Research Academy Hesse for FAIR, Campus Riedberg,Max-von-Laue-Straße 12, D-60438 Frankfurt am Main, Germany
In recent years, the existence of a hadronically stable ¯ b ¯ bud tetraquark with quantum numbers I ( J P ) = 0(1 + ) was confirmed by first principles lattice QCD computations. In this work we uselattice QCD to compare two frequently discussed competing structures for this tetraquark by consid-ering meson-meson as well as diquark-antidiquark creation operators. We use the static-light approx-imation, where the two ¯ b quarks are assumed to be infinitely heavy with frozen positions, while thelight u and d quarks are fully relativistic. By minimizing effective energies and by solving generalizedeigenvalue problems we determine the importance of the meson-meson and the diquark-antidiquarkcreation operators with respect to the ground state. It turns out, that the diquark-antidiquarkstructure dominates for ¯ b ¯ b separations r < ∼ . fm, whereas it becomes increasingly more irrelevantfor larger separations, where the I ( J P ) = 0(1 + ) tetraquark is mostly a meson-meson state. We alsoestimate the meson-meson to diquark-antidiquark ratio of this tetraquark and find around / . PACS numbers: 12.38.Gc, 13.75.Lb, 14.40.Rt, 14.65.Fy. ∗ [email protected] † [email protected] ‡ [email protected] § [email protected] a r X i v : . [ h e p - l a t ] J a n I. INTRODUCTION
A long standing problem in particle physics is to understand exotic hadrons, i.e. hadrons which have a structuremore complicated than a quark-antiquark pair or a triplet of quarks [1]. Such exotic hadrons turned out to be notonly difficult to observe experimentally, but also technical to address in quark models [2]. In the last couple ofyears, several exotic hadrons, the majority pertaining to the class of tetraquarks with at least two heavy quarks, wereclearly confirmed by the BELLE, BES-III and LHCb experimental collaborations. The observed exotic hadrons areresonances high in the spectrum. Studying them theoretically from first principles with lattice QCD requires theapplication and development of specific techniques from lattice hadron spectrosopy and scattering theory (see e.g.Refs. [3, 4]).There are two classes of doubly-heavy tetraquarks. Tetraquarks with one heavy quark and one heavy antiquark ¯ QQ ¯ qq including the Z c and Z b states are easier to detect experimentally. Their observation at Belle [5–7], Cleo-C[8], BESIII [9–13] and LHCb [14] collaborations turned tetraquarks into a main highlight of particle physics in recentyears. But since they are resonances with more than one decay channel, we study in this paper tetraquarks withtwo heavy antiquarks ¯ Q ¯ Qqq (or equivalently two heavy quarks, i.e. QQ ¯ q ¯ q ), which are theoretically simpler, becausethey are either hadronically stable or can only decay into a pair of heavy-light mesons. Moreover, with the recentobservation of hadronic systems with two heavy quarks [15, 16] at LHCb we expect this second class of tetraquarksto be observed in the near future. Their discovery potential is discussed in Refs. [17, 18].These ¯ Q ¯ Qqq tetraquarks are expected to form bound states, when the antiquarks are sufficiently heavy [19–29].Recently this was confirmed with lattice QCD computations. One of the approaches uses the Born-Oppenheimerapproximation [30, 31], where the problem is split into two steps. The first step is to compute the potentials of twostatic antiquarks in the presence of two light quarks using state-of-the-art lattice QCD techniques (see e.g. Refs. [32–37]). Then, in the second step, the heavy quark dynamics is studied using a quantum mechanical Hamiltonian withthe previously computed lattice QCD potentials. Using this approach, a ¯ b ¯ bud tetraquark bound state with quantumnumbers I ( J P ) = 0(1 + ) was predicted [36–40]. Shortly afterwards, this was confirmed by several full lattice QCDcomputations using four quarks of finite mass [41–45].Our present goal is to explore the structure of this ¯ b ¯ bud tetraquark with quantum numbers I ( J P ) = 0(1 + ) usingfirst principles lattice QCD computations, and to answer a hotly debated theoretical question [46–70]: Is it a diquark-antidiquark system (denoted in the following as Dd ) or rather a meson-meson system (denoted in the following as BB )?The lattice QCD result for the static potential relevant for the ¯ b ¯ bud tetraquark with quantum numbers I ( J P ) =0(1 + ) can be parameterized by a screened Coulomb potential (see left plot of Figure 1 and Refs. [32, 33, 35–37]).This is consistent with the following. (i) At ¯ b ¯ b separations larger than the typical meson radius each of the two heavyantiquarks forms a bound state with one of the light quark. Thus we have a system composed of two separated B mesons with interactions well-known to decay exponentially [71] with Yukawa-like potentials [72]. (ii) At small ¯ b ¯ b separations the heavy antiquarks interact directly via gluon exchange and are immersed in a light quark cloud. Thusat small distances the potential is a Coulomb potential as expected from the asymptotic freedom of perturbative QCD(see e.g. Ref. [73] and references therein).To compare quantitatively the importance of a Dd structure and of a BB structure for given ¯ b ¯ b separation, we utilizea set of lattice QCD creation operators, both of Dd type and of BB type with static ¯ b quarks. A similar approach waspreviously used to explore systems of four quarks of finite mass. For example, studies of four-quark systems including aheavy ¯ cc pair were to some extent inconclusive. Neither clear evidence for the existence of a tetraquark was found, nora signal improvement was observed for diquark-antidiquark operators [53, 74]. Also the tetraquark candidates a (980) and the D ∗ s (2317) were investigated in that way, employing sets of different creation operators, including quark-antiquark, meson-meson and/or diquark-antidiquark type [75–79]. The focus of these studies was more to distinguishbetween a quark-antiquark and a four-quark structure, and since computations of tetraquark correlation functions,where all quarks have a finite mass, are extremely challenging [80], no definite conclusion concerning meson-meson ordiquark-antidiquark dominance was reached. On the other hand, studies of potentials and the corresponding gluonfield distributions were performed with four static quarks, which are close to a system of four bottom quarks ¯ b ¯ bbb [81–83]. In this case it was possible to clearly distinguish between a diquark-antidiquark structure and a meson-mesonstructure for the ground state depending on the geometric positions of the static sources.This paper is structured as follows. In section II we review important technical steps from our previous work [38]and discuss in detail the lattice QCD creation operators of Dd type and of BB type. We detail our lattice QCD setupin section III. In Section IV we present our numerical results concerning the relative importance of a Dd structureand of a BB structure at given ¯ b ¯ b separation. At the end of this section we use these results to crudely estimate thepercentage of Dd and of BB in the ¯ b ¯ bud tetraquark with quantum numbers I ( J P ) = 0(1 + ) . We conclude in section V. II. POTENTIALS OF TWO STATIC ANTIQUARKS IN THE PRESENCE OF TWO LIGHT QUARKSAND CORRESPONDING CREATION OPERATORS
In preceding papers [33, 35, 37, 39] we have computed potentials V qq,j z , P , P x ( r ) of two static antiquarks ¯ Q ¯ Q atseparation r in the presence of two light quarks qq using lattice QCD. The computations have been carried out formany different sectors characterized by the following quantum numbers: light flavor qq with q ∈ { u, d, s, c } , totalangular momentum of the light quarks and gluons j z with respect to the ¯ Q ¯ Q separation axis, parity P and reflectionalong an axis perpendicular to the ¯ Q ¯ Q separation axis P x . There are both attractive and repulsive sectors. Mostpromising with respect to the existence of ¯ Q ¯ Qqq tetraquark bound states or resonances are attractive potentials withlight quarks q ∈ { u, d } , since they are rather wide and deep.Using the Born-Oppenheimer approximation, which amounts to solving the Schrödinger equation for the radialcoordinate of the two heavy quarks ¯ Q ¯ Q = ¯ b ¯ b with the computed potentials V qq,j z , P , P x ( r ) , (cid:18) m b (cid:18) − d dr + L ( L + 1) r (cid:19) + V qq,j z , P , P x ( r ) − m sl (cid:19) R ( r ) = ER ( r ) , (1)one can explore the existence of hadronically stable tetraquarks. m b is the b quark mass, L is the relative orbitalangular momentum of ¯ b ¯ b pair and m sl is the mass of the lightest static-light meson (computed within the same latticeQCD setup as V qq,j z , P , P x ( r ) ; see e.g. Refs. [84, 85]). There is one particular potential V ( r ) = V ud − du, , − , + ( r ) (shownin the left plot of Figure 1), which has quantum numbers ( I, j z , P , P x ) = (0 , , − , +) , leading for L = 0 to a stable ¯ b ¯ bud tetraquark with quantum numbers I ( J P ) = 0(1 + ) and binding energy − E = 38(18) MeV [86]. The probability densityof the ¯ b ¯ b separation π | R ( r ) | (shown in the right plot of Figure 1) indicates that one typically finds separations inthe range . fm < ∼ r < ∼ . fm. None of the other potentials is sufficiently wide or deep to host a bound state [39]. Figure 1. (left)
Lattice QCD results for the potential V ( r ) = V ud − du, , − , + ( r ) together with the parameterization − ( α/r ) e − ( r/d ) p with α = 0 . , d = 0 . fm and p = 2 . . (right) Probability density of the ¯ b ¯ b separation π | R ( r ) | . (The results shown inthe two plots are taken from Ref. [38].) As discussed in the introduction, the main goal of this work is to investigate the structure of the predicted ¯ b ¯ bud tetraquark. In particular, we explore, whether the tetraquark is more similar to a meson-meson state BB or to adiquark-antidiquark state Dd , two scenarios frequently discussed in the literature and at conferences also for othertetraquark candidates (see the discussion in section I). To this end, we refine the existing lattice QCD computationof the potential V ( r ) by using two types of creation operators.The first type of creation operators, which we used already in our previous computations [33, 35, 37, 39], excitestwo B mesons at separation r , O BB, Γ = 2 N BB ( C Γ) AB ( C ˜Γ) CD (cid:16) ¯ Q aC ( r ) ψ ( f ) aA ( r ) (cid:17)(cid:16) ¯ Q bD ( r ) ψ ( f (cid:48) ) bB ( r ) (cid:17) (2)with r = | r − r | , color indices a, b , spin indices A, B, C, D and ψ ( f ) ψ ( f (cid:48) ) = ud − du . N BB is a normalization, which willbe discussed later. There are two independent choices for the light spin matrix Γ consistent with ( j z , P , P x ) = (0 , − , +) . Γ = (1 + γ ) γ predominantly excites two negative parity ground state mesons B ( ∗ ) B ( ∗ ) , while Γ = (1 − γ ) γ mostlygenerates two positive parity excited mesons B ∗ , B ∗ , , as one can see e.g. by applying a Fierz transformation to O BB, Γ (see also Ref. [37]). Since static spins have no effect on energy levels, the heavy spin matrix is irrelevant and can bechosen arbitrarily, ˜Γ ∈ { (1 − γ ) γ , (1 − γ ) γ j } .The second type of creation operators, which we use here for the first time, resembles a diquark-antidiquark pairwith heavy quarks separated by r and connected by a gluonic string, O Dd, Γ = − N Dd (cid:15) abc (cid:16) Ψ ( f ) bA ( z )( C Γ) AB Ψ ( f (cid:48) ) cB ( z ) (cid:17) (cid:15) ade (cid:16) ¯ Q fC ( r ) U fd ( r ; z )( C ˜Γ) CD ¯ Q gD ( r ) U ge ( r ; z ) (cid:17) . (3)Again N Dd is a normalization and the allowed light and heavy spin matrices are the same as for the operator O BB, Γ ,i.e. Γ ∈ { (1 − γ ) γ , (1+ γ ) γ } and ˜Γ ∈ { (1 − γ ) γ , (1 − γ ) γ j } . Since the heavy spins are not part of the Hamiltonian,it is important to use the same ˜Γ for the operators O BB, Γ and O Dd, Γ , whenever they are part of the same correlationmatrix. For definiteness, we choose ˜Γ = (1 − γ ) γ . r and r are always separated along one of the lattice axes and z = ( r + r ) / . For odd r/a ( a denotes the lattice spacing), z does not coincide with one of the lattice sites. In thesecases we take the average of the operator (3) with z = ( r + r ) / a ˆ r / and with z = ( r + r ) / − a ˆ r / , where ˆ r = ( r − r ) / | r − r | .We use smeared light quark and gluon fields, which implies that the operator ψ ( f ) aA ( r ) does not generate a point-likeexcitation at r , but rather a cloud-like excitation of spherical shape with diameter ≈ . fm. Similarly, the gauge linksconnecting the two static antiquarks in the operator O Dd, Γ generate a flux tube with a certain thickness. For detailssee section III, where our lattice setup is discussed.Note that the creation operators O BB, Γ and O Dd, Γ do not generate orthogonal states, when applied to the vacuum | Ω (cid:105) . For r = 0 , the operators are even identical, when properly normalized, i.e. O BB, Γ = O Dd, Γ , if N BB = N Dd . Thiscan easily be shown by using the identity (cid:15) abc (cid:15) ade = δ bd δ ce − δ be δ cd . For increasing r , however, they become more andmore different, as we will show numerically in section IV A. One obvious reason for that is that O Dd, Γ creates a fluxtube of length r , whereas O BB, Γ does not. In section IV B, section IV C and section IV D we will explore, whether theground state of the ( I, j z , P , P x ) = (0 , , − , +) sector at a given separation r is more similar to a BB state O Dd, Γ | Ω (cid:105) or to a DD state O Dd, Γ | Ω (cid:105) . Since the energy of this ground state is the potential V ( r ) used in the Born-Oppenheimerprediction of the ¯ b ¯ bud tetraquark with quantum numbers I ( J P ) = 0(1 + ) , our investigation will provide informationon the structure of that tetraquark. We will discuss that in section IV E. III. LATTICE QCD SETUP AND TECHNIQUESA. Lattice actions
The light quark action used in this work is the Wilson twisted mass action [87, 88], S F [ χ, ¯ χ, U ] = (cid:88) x ¯ χ ( x ) (cid:16) D W ( m ) + iµγ τ (cid:17) χ ( x ) , (4)where D W is the standard Wilson Dirac operator, D W ( m ) = 12 (cid:16) γ µ ( ∇ µ + ∇ ∗ µ ) − ∇ ∗ µ ∇ µ (cid:17) + m , (5)with the forward and backward covariant derivatives ∇ µ and ∇ ∗ µ . χ = ( χ ( u ) , χ ( d ) ) is the light quark doublet in theso-called twisted basis, which is related to quark fields in the usual “physical basis” via the twist rotation ψ = e iγ τ ω/ χ, (6)where ω is the twist angle.The gluon action used in this work is the tree-level Symanzik improved action [89], S G [ U ] = β (cid:88) x (cid:18) b (cid:88) µ,ν =1 Re Tr (cid:16) − P × µν ( x ) (cid:17) + b (cid:88) µ (cid:54) = ν ReTr (cid:16) − P × µν ( x ) (cid:17)(cid:19) , (7)where b = − / , b = 1 − b , β = 6 /g , g is the bare coupling and P × µν and P × µν are the plaquette and a × Wilson loop, respectively.To achieve automatic O ( a ) improvement, the hopping parameter κ = (8 + 2 am ) − is tuned to its critical value, atwhich the PCAC quark mass vanishes [87, 90–93]. This corresponds to maximal twist, i.e. to tuning ω in Eq. (6) to π/ . B. Ensembles of gauge link configurations
We have performed computations on two ensembles of gauge link configurations generated by the European TwistedMass Collaboration (ETMC) [94–96]. These ensembles have different lattice spacings, a ≈ . fm and a ≈ . fm(set via the pion mass and pion decay constant and chiral perturbation theory [96]), which allows to study a finerresolution in the static antiquark-antiquark separation r and to check for discretization errors. The pion masses m PS and spatial and temporal extents L and T are, however, very similar. The parameters and details of the ensemblesare collected in Table I. ensemble name β a in fm ( L/a ) × T /a κ µ m PS in MeV .
90 0 . ×
48 0 . .
004 340(13) 108
C30.32 .
05 0 . ×
64 0 . .
003 325(10) 98
Table I. Ensembles of gauge link configurations.
C. Smearing techniques
As already mentioned in section II, we have used several smearing techniques.The spatial gauge links appearing in the Dd creation operator (3) are APE smeared [97] with parameters N APE = 30 and α APE = 0 . (see Ref. [84], Eq. (23)). The quark fields appearing in the BB and Dd creation operators (2) and (3)are Gaussian smeared [98] with parameters N Gauss = 50 and κ Gauss = 0 . (see Ref. [84], Eq. (25)) and APE-smearedspatial gauge links. The intention of both APE and Gaussian smearing is to increase the ground state overlapsgenerated by the creation operators.In addition to that we also use HYP2-smeared gauge links in temporal direction [99–101] with parameters α = α = 1 . and α = 0 . . This reduces the self energy of the static quarks and, thus, improves the signal-to-noiseratio of the computed correlation functions. IV. NUMERICAL RESULTS
We consider the following four creation operators, O j , j ∈ (cid:110) [ BB, (1 + γ ) γ ] , [ BB, γ ] , [ Dd, (1 + γ ) γ ] , [ Dd, γ ] (cid:111) , (8)and define the corresponding trial states as | Φ j (cid:105) = O j | Ω (cid:105) . (9) O BB, (1+ γ ) γ predominantly excites two negative parity ground state mesons. Thus, | Φ BB, (1+ γ ) γ (cid:105) is expected tohave the largest overlap to the ground state of the ( I, j z , P , P x ) = (0 , , − , +) sector at large r . At small r , however,the diquark-antidiquark operators might be advantageous. We use | Φ Dd,γ (cid:105) , i.e. a light diquark with just γ , astypically discussed in the literature. Since O BB, Γ ∝ O Dd, Γ for r = 0 (see section II), the diquark-antidiquark trialstate | Φ Dd, (1+ γ ) γ (cid:105) might even be a better candidate for having a large ground state overlap. For completeness wealso include O BB,γ , the mesonic molecule counterpart of O Dd,γ . O BB,γ excites a linear combination of two negativeparity ground state mesons and two significantly heavier positive parity excited mesons [37].With these operators we computed the × correlation matrix, C jk ( t ) = (cid:68) O † j ( t ) O k ( t ) (cid:69) = (cid:104) Ω |O † j ( t ) O k ( t ) | Ω (cid:105) = (cid:104) Φ j ( t ) | Φ k ( t ) (cid:105) (10)with (cid:104) . . . (cid:105) denoting the path integral expectation value and t/a = ( t − t ) /a ≥ . For the second equality we assumedthat in the spectral decomposition propagation over temporal separation T − t is suppressed for all states except forthe vacuum. To cross-check our computations, we checked the numerical results with respect to the symmetries γ hermiticity, parity, time reversal, charge conjugation and cubic rotations around the axis of separation (for details seeRef. [37]). In a second step we averaged elements of the correlation matrices related by these symmetries, to reducestatistical errors. A. Squared overlaps of the normalized BB and Dd trial states In this subsection we study α jk ( t ) = | C jk ( t ) | C jj ( t ) C kk ( t ) . (11)For t → this quantity is the squared normalized overlap of trial state | Φ j (cid:105) = O j | Ω (cid:105) and trial state | Φ k (cid:105) = O k | Ω (cid:105) ,i.e. α jk = lim t → α jk ( t ) = |(cid:104) Φ j | Φ k (cid:105)| . Clearly, ≤ α jk ≤ . For an arbitrary state | Φ (cid:105) and an orthonormal basis | Φ j (cid:105) , j = 1 , , , . . . ∞ (cid:88) k =1 α k = 1 . (12)Thus for two trial states | Φ j (cid:105) and | Φ k (cid:105) , α jk can be interpreted as a measure of their orthogonality, where α jk ≈ indicates almost orthogonal and α jk ≈ almost parallel states. For large t all α jk approach , because the groundstate dominates in that limit.In the left plot of Figure 2 we show α jk for j = BB, (1 + γ ) γ and k = Dd, (1 + γ ) γ as function of t for severalfixed r [102] for ensemble B40.24. The corresponding trial states become more similar for smaller r , as indicatedby larger values of α jk close to . This is not surprising, because for r = 0 the two operators O BB, Γ and O Dd, Γ are identical, if normalized according to N BB = N Dd , as discussed in section II. For larger heavy quark separations r , the two trial states clearly differ. O BB, (1+ γ ) γ generates a pair of spatially separated B mesons, which interactonly by weak residual hadronic forces. In contrast to that, the antiquarks ¯ Q ¯ Q forming the heavy diquark in theoperator O Dd, (1+ γ ) γ are connected by a gluonic flux tube of length r . The data points for ≤ t/a ≤ can be fittedconsistently with degree-2 polynomials, which are also shown in Figure 2. These fits represent crude extrapolationsof α jk to t = 0 , which are the squared overlaps of the corresponding normalized BB and Dd trial states. Figure 2. α jk as a function of t for several fixed r for ensemble B40.24. (left) j = BB, (1 + γ ) γ , k = Dd, (1 + γ ) γ . (right) j = BB, γ , k = Dd, γ . For t → , α jk is the squared overlap of the corresponding normalized trial states. The curves,which are fits with degree-2 polynomials, represent crude extrapolations to t = 0 . In the right plot of Figure 2 we show the corresponding results for j = BB, γ and k = Dd, γ . They are quitesimilar to those for j = BB, (1 + γ ) γ and k = Dd, (1 + γ ) γ and can be interpreted in the same way. The mainpoint of these two plots is to demonstrate that BB and Dd trial states, even though not orthogonal, are not linearlydependent either. In the following subsections we explore, whether the ground state of the ( I, j z , P , P x ) = (0 , , − , +) sector for given r is more similar to a BB trial state or to a Dd trial state.Analog plots of α jk for ensemble C30.32 are very similar to those of Figure 2 and, thus, not shown. B. Effective energies corresponding to diagonal elements of the correlation matrix at small and largetemporal separations
Now we consider effective energies corresponding to diagonal elements of the correlation matrix, i.e. V eff j ( r, t ) = − a log (cid:18) C jj ( t ) C jj ( t − a ) (cid:19) (no sum over j ) . (13)As discussed in section II, all four operators probe the ( I, j z , P , P x ) = (0 , , − , +) sector. Thus, for fixed r and atsufficiently large t , all four V eff j ( r, t ) should approach the same constant, which is the ground state energy V ( r ) . Forseparations r < ∼ . fm, our numerical results confirm that expectation (see e.g. the left plot of Figure 3, where V eff j ( r, t ) is shown for r ≈ . fm as a function of t ). For larger r , however, the two Dd operators seem to generate only littleoverlap to the ground state. The consequence is that the corresponding effective energies V eff Dd, (1+ γ ) γ and V eff Dd,γ do not convincingly converge to the plateau at V ( r ) in the t region, where we carried out computations, and haverather large statistical errors for larger t separations (see e.g. the right plot of Figure 3, where V eff j ( r, t ) is shown for r ≈ . fm as a function of t ). In contrast to that, the two BB operators still lead to clear effective energy plateaus. Figure 3. Effective energies V eff j corresponding to diagonal elements of the correlation matrix for fixed r as functions of t forensemble B40.24. (left) r = 2 a ≈ . fm. (right) r = 10 a ≈ . fm. The dotted gray line in both plots represents the BB threshold at m B , where m B is the mass of the lightest static light meson taken from our previous work [85]. A first indicator concerning the structure of the ground state are effective energies at small temporal separations,i.e. V eff j ( r, t = 2 a ) . Since there is little suppression of excited states by the Euclidean time evolution, a small valueof V eff j ( r, t = 2 a ) close to the ground state energy V ( r ) implies an operator, which predominantly excites the groundstate. A larger value of V eff j ( r, t = 2 a ) , on the other hand, is a sign that the corresponding operator creates a trialstate less similar to the ground state.We start with a comparison of the two operators O BB, (1+ γ ) γ and O BB,γ for ensemble B40.24 by showing thedifference of their effective energies, V eff BB, (1+ γ ) γ ( r, t = 2 a ) − V eff BB,γ ( r, t = 2 a ) , in the upper left plot of Figure 4.This difference is clearly negative for all separations r , which is not surprising. O BB, (1+ γ ) γ creates predominantly apair of ground state static-light mesons, while O BB,γ creates roughly a / superposition of a pair of negativeparity ground state mesons and a pair of significantly heavier positive parity static-light mesons, as discussed aboveand as can be shown e.g. by a Fierz transformation (for details see Ref. [37]). Results from an analog comparison ofthe two diquark-antidiquark operators are very similar (see upper right plot of Figure 4). This is interesting, becauseit indicates that a light diquark with spin structure given by (1 + γ ) γ is energetically preferred over a light diquarkwith spin structure given just by γ .Most interesting, of course, is the comparison of a meson-meson and a diquark-antidiquark operator, specifically of O BB, (1+ γ ) γ and O Dd, (1+ γ ) γ , which we have just identified as being superior to O BB,γ and O Dd,γ , respectively.We show the difference of their effective masses, V eff BB, (1+ γ ) γ ( r, t = 2 a ) − V eff Dd, (1+ γ ) γ ( r, t = 2 a ) , in the lower plotof Figure 4. For r < ∼ . a ≈ . fm the difference is positive, indicating that for small separations the diquark-antidiquark operator generates a trial state more similar to the ground state. For larger separations, r > ∼ . fm the Figure 4. V eff j ( r, t = 2 a ) − V eff k ( r, t = 2 a ) , i.e. differences of effective energies corresponding to diagonal elements of thecorrelation matrix for t/a = 2 , as functions of r for ensemble B40.24. (top left) j = BB, (1 + γ ) γ , k = BB, γ . (topright) j = Dd, (1 + γ ) γ , k = Dd, γ . (bottom) j = BB, (1 + γ ) γ , k = Dd, (1 + γ ) γ . difference becomes negative and strongly points towards a meson-meson structure. While the latter is expected (atlarge r the flux tube present in a heavy diquark-antidiquark is energetically disfavored), the former is a first hinttowards a diquark-antidiquark dominance at smaller r . We continue to investigate this in more detail in the followingsubsections.Analog plots of differences of effective energies for ensemble C30.32 are very similar to those of Figure 4 and, thus,not shown. C. Optimizing trial states by minimizing effective energies
Now we consider the 2-dimensional space spanned by the states | Φ BB, (1+ γ ) γ (cid:105) and | Φ Dd, (1+ γ ) γ (cid:105) . Any trial statefrom that space can be written as | Φ b,d (cid:105) = b | Φ BB, (1+ γ ) γ (cid:105) + d | Φ Dd, (1+ γ ) γ (cid:105) (14)with coefficients b, d ∈ C . To identify a trial state as similar to the ground state as possible, i.e. with large overlap tothe ground state and little overlap to excitations, we minimize the corresponding effective energy V eff b,d ( r, t ) = − a log (cid:18) C [ b,d ][ b,d ] ( t ) C [ b,d ][ b,d ] ( t − a ) (cid:19) (15)with respect to b and d . Since, V eff b,d is independent of the norm and the phase of | Φ b,d (cid:105) , we can fix b = 1 and minimize V eff b,d for given r and t with respect to the real and the imaginary part of d . Such a 2-dimensional minimization isnumerically straightforward. In particular, no further lattice QCD computations are needed, because C [ b,d ][ b,d ] can beexpressed in terms of the correlation matrix introduced in Eq. (10), C [ b,d ][ b,d ] ( t ) = (cid:32) bd (cid:33) † j C jk ( t ) (cid:32) bd (cid:33) k . (16)In the following we consider w BB = | b | | b | + | d | , w Dd = | d | | b | + | d | (17)with b = 1 and d minimizing V eff b,d . w BB and w Dd are the normalized absolute squares of the coefficients of theoptimized trial states appearing in Eq. (14). These quantities exhibit only a weak dependence on t . For ≤ t/a ≤ and ensemble B40.24 ( ≤ t/a ≤ and ensemble C30.32) they are consistent with a constant. For t/a ≥ ( t/a ≥ ),statistical fluctuations and errors become large and the signal is quickly lost in noise. The latter is not surprising,because w BB and w Dd are subtle quantities depending on the amount of excited states in BB and Dd correlationfunctions, which are exponentially suppressed in t . In Figure 5 we show example plots of w BB and w Dd as functionsof t for selected separations r/a = 2 , r/a = 5 and r/a = 8 for ensemble B. Figure 5. w BB and w Dd = 1 − w BB , the normalized absolute squares of the coefficients of the optimized trial states forseveral fixed r as functions of t for ensemble B40.24. The horizontal red lines indicate the fit results ¯ w BB and ¯ w Dd and thecorresponding statistical errors. We determine each plateau value by a χ minimizing fit of a constant in the range ≤ t/a ≤ ( ≤ t/a ≤ ). Theresulting numbers, ¯ w BB ( r ) and ¯ w Dd ( r ) = 1 − ¯ w BB ( r ) , can be interpreted as the relative weight of a meson-mesonand a diquark-antidiquark structure at ¯ b ¯ b separation r in the ground state, which corresponds to the potential V ( r ) of two static antiquarks and is, thus, closely related to the ¯ b ¯ bud tetraquark with quantum numbers I ( J P ) = 0(1 + ) . InFigure 6 we plot ¯ w BB and ¯ w Dd as functions of r . One can clearly see that there is a diquark-antidiquark dominancefor ¯ b ¯ b separations r < ∼ . fm. For r > ∼ . fm, the meson-meson structure is more prominent and for r > ∼ . fm, thediquark-antidiquark contribution is negligible, i.e. the system is exclusively composed of two B mesons.We note that the quantities w BB and w Dd as well as the fitted ¯ w BB and ¯ w Dd depend on the normalization ofthe operators O BB, (1+ γ ) γ and O Dd, (1+ γ ) γ or, equivalently, on the normalization and the corresponding states | Φ BB, (1+ γ ) γ (cid:105) and | Φ Dd, (1+ γ ) γ (cid:105) . To allow a meaningful interpretation in terms of the relative weight of a BB anda Dd structure, the norms of these two states has to be similar. In this work we use N BB = N Dd , which implies O BB, Γ = O Dd, Γ for r = 0 and, thus, | Φ BB, Γ (cid:105) = | Φ Dd, Γ (cid:105) for r = 0 (see section II). We expect that N BB = N Dd alsoresults in similar norms for r > . A common alternative, which we have used in previous lattice QCD projects, is tonormalize the operators O j such that C jj ( t = a ) = 1 (no sum over j ; see e.g. Refs. [78, 79, 103]). As a cross-checkwe also explored this normalization in our current work and found almost identical results to those obtained with N BB = N Dd .0 Figure 6. ¯ w BB and ¯ w Dd = 1 − ¯ w BB , the fitted normalized absolute squares of the coefficients of the optimized trial state, asfunctions of r for both ensembles. D. Eigenvector components obtained by solving a generalized eigenvalue problem
To further explore the structure of the ground state in the ( I, j z , P , P x ) = (0 , , − , +) sector, we now use N × N correlation matrices C jk as defined in Eq. (10) and solve the generalized eigenvalue problem (GEVP) C jk ( t ) v ( n ) k ( t ) = λ ( n ) ( t ) C jk ( t ) v ( n ) k ( t ) , n = 0 , . . . , N − (18)for t /a ≥ and t/a > t /a (for a detailed discussion of the GEVP in lattice field theory see e.g. Ref. [104]). Effectiveenergies for the lowest N energy eigenstates are then given by V eff , ( n ) ( r, t ) = − a log (cid:18) λ ( n ) ( t ) λ ( n ) ( t − a ) (cid:19) . (19)These are generalizations of the effective energy defined in Eq. (13), since V eff , (0) ( r, t ) = V eff j ( r, t ) for N = 1 . For N > , V eff , (0) ( r, t ) approaches the same constant V ( r ) for large t , but plateaus can typically be identified atsomewhat smaller t , because of an elimination of excitations (see also the second next paragraph and appendix A,where a minimization of effective energies is related to the GEVP).The eigenvector components v ( n ) j ( t ) , which we always normalize according to (cid:80) j | v ( n ) j ( t ) | = 1 , contain informationabout the relative importance of the creation operators included in the correlation matrix and, thus, hints about thestructure of the corresponding energy eigenstates. For large t and t , | n (cid:105) ≈ (cid:88) j v ( n ) j ( t ) | Φ j (cid:105) , (20)where the ≈ sign denotes an approximate expansion of the energy eigenstate | n (cid:105) in terms of the trial states | Φ j (cid:105) . Forsuch values of t and t , the squared eigenvector components as functions of t form plateaus and we determine thecorresponding asymptotic values of | v (0) j ( t ) | by χ minimizing fits of constants. The results of these fits are denotedby | ¯ v (0) j | . The squared eigenvector components | v ( n ) j ( t ) | as well as the fitted | ¯ v (0) j | depend on the normalization ofthe creation operators, as it is the case for w BB , w Dd , ¯ w BB and ¯ w Dd (see the discussion at the end of section IV C).As before, we use N BB = N Dd , which amounts to having trial states | Φ j (cid:105) with similar norm.It is interesting to note that for a × correlation matrix with trial states | Φ (cid:105) = | Φ BB, (1+ γ ) γ (cid:105) and | Φ (cid:105) = | Φ Dd, (1+ γ ) γ (cid:105) and t /a = t/a − , the eigenvector components v (0) j are proportional to the coefficients b and d minimizing V eff b,d ( r, t ) defined in Eq. (15). Moreover, ( | v (0) BB, (1+ γ ) γ | , | v (0) Dd, (1+ γ ) γ | ) = ( w BB , w Dd ) and, thus, ( | ¯ v (0) BB, (1+ γ ) γ | , | ¯ v (0) Dd, (1+ γ ) γ | ) = ( ¯ w BB , ¯ w Dd ) , as we show in appendix A. In other words, the results on the structureof the ¯ b ¯ bud ground state obtained in section IV C by determining a trial state, which minimizes an effective energy,are identical to results from a specific corresponding GEVP. This allows to understand and interpret the GEVP1eigenvector components from another perspective. Compared to the trial state optimization from section IV C, theGEVP, however, offers further possibilities, for example to choose t independent of t (i.e. not as t /a = t/a − ) orto study the full correlation matrix with N = 4 .As discussed in the previous paragraph, solving the GEVP with a × correlation matrix including the operators O BB, (1+ γ ) γ and O Dd, (1+ γ ) γ and using t /a = t/a − yields exactly the same results as shown in Figure 6 (one justhas to replace labels according to ¯ w BB → | ¯ v (0) BB, (1+ γ ) γ | and ¯ w Dd → | ¯ v (0) Dd, (1+ γ ) γ | ). However, many lattice QCDpapers using the GEVP fix t to a rather small value independent of t . For small values of t , statistical errors aresomewhat reduced, which in turn allows to consider larger values of t . Thus we also computed | ¯ v (0) BB, (1+ γ ) γ | and | ¯ v (0) Dd, (1+ γ ) γ | for t /a = 1 . As already observed in the previous subsection for the related quantities w BB and w Dd ,the resulting squared eigenvector components | v (0) BB, (1+ γ ) γ | and | v (0) Dd, (1+ γ ) γ | exhibit only a mild t dependence andthe majority is consistent with a constant for large t . Results for selected ¯ b ¯ b separations r are shown in Figure 7.For ensemble C30.32 we use the fit range ≤ t/a ≤ to determine | ¯ v (0) BB, (1+ γ ) γ | and | ¯ v (0) Dd, (1+ γ ) γ | . For ensembleB40.24 the analysis and interpretation of the data is less obvious. While for t/a ≤ plateaus are indicated, for t/a ≥ there is a weak, almost linearly increasing deviation from these plateaus for some separations r (see e.g. the plot inthe center of the upper line of Figure 7). Since errors are also increasing at larger t , we interpret this as statisticalfluctuation and use the fit range ≤ t/a ≤ . Figure 7. The squared eigenvector components | v (0) BB, (1+ γ ) γ | and | v (0) Dd, (1+ γ ) γ | = 1 − | v (0) BB, (1+ γ ) γ | for several fixed r asfunctions of t . (upper line) Ensemble B40.24. (lower line)
Ensemble C30.32. The horizontal red lines indicate the fit results | ¯ v (0) BB, (1+ γ ) γ | and | ¯ v (0) Dd, (1+ γ ) γ | and the corresponding statistical errors. In the left plot of Figure 8 we show | ¯ v (0) BB, (1+ γ ) γ | and | ¯ v (0) Dd, (1+ γ ) γ | for both ensembles. These curves for t /a = 1 are quite similar to those corresponding to t /a = t/a − (and shown in Figure 6). Again, one can see a clear dominanceof the diquark-antidiquark operator for separations r < ∼ . fm. In the range . fm ≤ r ≤ . fm there is a rapidchange towards a meson-meson structure. For r > ∼ . fm there is almost no diquark-antidiquark contribution anymoreand the ¯ b ¯ bud four-quark system seems to be composed exclusively of two B mesons. This plot confirms our resultsfrom section IV B obtained by minimizing effective energies.It should be noted that the quantities ¯ w BB and ¯ w Dd as well as the fitted squared eigenvector components2 Figure 8. Fitted squared eigenvector components | ¯ v (0) j | as functions of r . (left) × correlation matrix including the operators O BB, (1+ γ ) γ and O Dd, (1+ γ ) γ for both ensembles. (right) × correlation matrix including all operators (see Eq. (8)) forensemble C30.32. | ¯ v (0) BB, (1+ γ ) γ | and | ¯ v (0) Dd, (1+ γ ) γ | depend on the creation operators used, e.g. the details of the quark and gaugefield smearing or their normalization. Besides discretization errors, this might be part of the reason for the slightdifferences between the results obtained for ensemble B40.24 and ensemble C30.32 shown in Figure 6 and the leftplot of Figure 8. Since these differences are quite small and since we found almost identical results for differentnormalization prescriptions, we consider ¯ w BB and ¯ w Dd as well as | ¯ v (0) j | as reliable indicators characterizing the quarkand gluon structure of the ground state of the ( I, j z , P , P x ) = (0 , , − , +) sector.Finally we performed the same GEVP analysis using the full × correlation matrix including all operators definedin Eq. (8). We determined | ¯ v (0) j | using the same fit ranges as for the previous × analyses (results for ensembleC30.32 are shown in the right plot of Figure 8). Again, there is a diquark-antidiquark dominance for r < ∼ . fm, thistime represented by the eigenvector components corresponding to two operators O Dd, (1+ γ ) γ and O Dd,γ , while thereis meson-meson dominance for r > ∼ . fm, reflected by the eigenvector components corresponding to two operators O BB, (1+ γ ) γ and O BB,γ . When adding the two diquark eigenvector components as well as the two meson-mesoneigenvector components, i.e. when considering | ¯ v (0) Dd, (1+ γ ) γ | + | ¯ v (0) Dd,γ | and | ¯ v (0) BB, (1+ γ ) γ | + | ¯ v (0) BB,γ | , we find withinerrors the same curves as obtained above by the × analysis. This is reassuring, because the results concerningthe meson-meson percentage and the diquark-antidiquark percentage do not change, even though we increased thedimension of the basis of trial states used to approximate the ground state from to . E. Meson-meson and diquark-antidiquark percentages of the ¯ b ¯ bud tetraquark with I ( J P ) = 0(1 + ) During the lattice QCD computation of the potential V ( r ) the heavy antiquarks ¯ b ¯ b are considered as static, i.e.their positions are fixed and only the light quarks ud and the gluons are dynamical degrees of freedom. The dynamicsof the heavy quarks can, however, be studied in a second step, by inserting the potential V ( r ) into the Schrödingerequation (1). This two step approach is widely known as the Born-Oppenheimer approximation (see e.g. Ref. [31]for a detailed discussion in the context of exotic mesons). In Ref. [38] we solved the Schrödinger equation (1) using m b = m B = 5279 MeV [105] and found a single bound state with binding energy − E = 38(18) MeV indicatingthe existence of a hadronically stable ¯ b ¯ bud tetraquark with quantum numbers I ( J P ) = 0(1 + ) . Moreover, the wavefunction R ( r ) /r gives the probability density of the ¯ b ¯ b separation, p r ( r ) = 4 π | R ( r ) | , which is shown in the right plotof Figure 1.The quantities ¯ w BB and ¯ w Dd (see Figure 6) as well as | ¯ v (0) BB, (1+ γ ) γ | and | ¯ v (0) Dd, (1+ γ ) γ | (see left plot of Figure 8)can be interpreted as meson-meson and diquark-antidiquark percentages for fixed r . We define p BB,w ( r ) = ¯ w BB , p Dd,w ( r ) = ¯ w Dd (21)and p BB,v ( r ) = | ¯ v (0) BB, (1+ γ ) γ | , p Dd,v ( r ) = | ¯ v (0) Dd, (1+ γ ) γ | , (22)3where p BB,w ( r ) ≈ p BB,v ( r ) and p Dd,w ( r ) ≈ p Dd,v ( r ) . These percentages together with the probabilty density p r ( r ) can be used to crudely estimate the total meson-meson and diquark-antidiquark percentages of the ¯ b ¯ bud tetraquark,which we denote by % BB j and % Dd j . To this end, we use the parameterization p BB,j ( r ) = 12 (cid:16) tanh ( α j ( r − r ,j ) + 1 (cid:17) , p Dd,j ( r ) = 1 − p BB,j ( r ) , (23)where the parameters r ,j and α j are determined by χ minimizing fits to the results of both ensembles as shown inFigure 6 and Figure 8. Then we compute % BB and % Dd via % BB j = (cid:90) dr p r ( r ) p BB,j ( r ) , % Dd j = (cid:90) dr p r ( r ) p Dd,j ( r ) = 1 − % BB j . (24)We find % BB w = 0 . , % Dd w = 0 . and % BB v = 0 . , % Dd v = 0 . , i.e. almost the same result, when using ¯ w BB and ¯ w Dd and when using | ¯ v (0) BB, (1+ γ ) γ | and | ¯ v (0) Dd, (1+ γ ) γ | . Since the eigenvector components are slightly operatordependent, as discussed in section IV D, these numbers should only be considered as crude estimates with a systematicerror of maybe ± . Still it seems to be clear that the ¯ b ¯ bud tetraquark with quantum numbers I ( J P ) = 0(1 + ) isneither strongly meson-meson dominated nor strongly diquark-antidiquark dominated, but rather an approximatelyequal linear combination of both structures. Finally we note that these results for % BB j and for % Dd j are fullyconsistent with squared eigenvector components obtained during a recent lattice QCD study [44] of the same tetraquarkusing four quarks of finite mass. There, a meson-meson component of . and a diquark-antidiquark componentof . was found [106]. V. CONCLUSIONS AND OUTLOOK
In this work we used lattice QCD to study a recurrent question on the nature of tetraquarks in the context of a ¯ b ¯ bud tetraquark with quantum numbers I ( J P ) = 0(1 + ) : Are they more similar to meson-meson systems or ratherto diquark-antidiquark pairs? Moreover we addressed the Dirac structure of the light quarks, comparing γ with (1 + γ ) γ , which are both consistent with I ( J P ) = 0(1 + ) .We implemented four different lattice QCD creation operators, two of meson-meson type and two of diquark-antidiquark type. We solved the GEVP both for × and × correlation matrices, to determine quantitatively,which of the implemented structures is preponderant in the ground state at ¯ b ¯ b separation r . Moreover, we optimizedtrial states by minimizing effective energies and proved that this is equivalent to solving a specific correspondingGEVP.Notice that the question we are addressing is quite subtle. We first showed that the BB and Dd trial states are notorthogonal. Nevertheless, since they are not linearly dependent either, we were able to determine, which one is moresimilar to the ground state. In what concerns light spin we found the (1 + γ ) γ is the dominant Dirac structure at allseparations r , both for BB and for Dd . This is what we expected, since the less favorable γ structure generates notonly negative parity mesons, but also excited positive parity mesons. As for color we showed that at small separations r < ∼ . fm the diquark-antidiquark structure dominates, whereas at larger separations the meson-meson trial state isclearly more similar to the ground state of the ¯ b ¯ bud system. For r > ∼ . fm, the percentage of BB is already largerthan and approaches for even larger r .Using these results as well as the wave function of the ¯ b ¯ b separation already obtained in Ref. [38], we estimated themeson-meson to diquark-antidiquark ratio of the ¯ b ¯ bud tetraquark with quantum numbers I ( J P ) = 0(1 + ) and foundaround / . Thus BB and Dd components seem to be present in comparable parts. Appendix A: Equivalence of optimizing a trial state by minimizing an effective energy and of solving a GEVP
We start with the N × N GEVP defined in Eq. (18), choose t = t − a and use the simplified notation v ( n ) ≡ ( v ( n )1 ( t ) , v ( n )2 ( t ) , . . . , v ( n ) N ( t )) and λ ( n ) ≡ λ ( n ) ( t ) , C ( t ) v ( n ) = λ ( n ) C ( t − a ) v ( n ) , n = 0 , . . . , N − . (A1)The correlation matrix C defined in Eq. (10) is hermitean, i.e. C † = C , the eigenvalues λ ( n ) are real and we assumethem to be positive and non-degenerate, i.e. λ (0) > λ (1) > . . . > λ ( N − > .We can rewrite v ( m ) † C ( t ) v ( n ) in two ways by using Eq. (A1) and its hermitean conjugate, v ( m ) † C ( t ) v ( n ) = v ( m ) † C ( t − a ) v ( n ) λ ( n ) (A2) v ( m ) † C ( t ) v ( n ) = v ( m ) † C ( t − a ) v ( n ) λ ( m ) . (A3)4Since λ ( m ) (cid:54) = λ ( n ) for m (cid:54) = n , we conclude v ( m ) † C ( t − a ) v ( n ) = 0 for m (cid:54) = n . In the same way one can show v ( m ) † C ( t ) v ( n ) = 0 for m (cid:54) = n . Moreover, it is convenient to normalize the eigenvectors according to v ( n ) † C ( t − a ) v ( n ) = 1 .Now we consider an arbitrary trial state | Ψ (cid:105) from the N -dimensional space spanned by the states | Φ j (cid:105) included inthe correlation matrix C , | Ψ (cid:105) = (cid:88) j a j | Φ j (cid:105) . (A4)The corresponding correlation function is C Ψ ( t ) = a † C ( t ) a . (A5)Since a (as well as any other complex N -component vector) can be expanded in terms of the eigenvectors accordingto a = (cid:88) n µ ( n ) v ( n ) (A6)with µ ( n ) ∈ C , we can write for temporal separation tC Ψ ( t ) = (cid:88) m,n µ ( m ) ∗ v ( m ) † C ( t ) v ( n ) µ ( n ) = (cid:88) n | µ ( n ) | v ( n ) † C ( t ) v ( n ) = (cid:88) n | µ ( n ) | λ ( n ) v ( n ) † C ( t − a ) v ( n ) == (cid:88) n | µ ( n ) | λ ( n ) (A7)and analogously for temporal separation t − aC Ψ ( t − a ) = (cid:88) n | µ ( n ) | v ( n ) † C ( t − a ) v ( n ) = (cid:88) n | µ ( n ) | . (A8)Now we consider the effective energy corresponding to the correlation function C Ψ , E eff (cid:126)µ ( t ) = − a log (cid:18) C Ψ ( t ) C Ψ ( t − a ) (cid:19) = − a log (cid:18) (cid:88) m | µ ( m ) | (cid:80) n | µ ( n ) | λ ( n ) (cid:19) . (A9)Since ≤ | µ ( m ) | / (cid:80) n | µ ( n ) | ≤ and (cid:80) m ( | µ ( m ) | / (cid:80) n | µ ( n ) | ) = 1 , the argument of the logarithm in Eq. (A9)is a weighted sum of the eigenvalues, i.e. can assume values between the maximal eigenvalue λ (0) and the minimaleigenvalue λ ( N − . Minimizing E eff (cid:126)µ ( t ) with respect to (cid:126)µ is equivalent to maximizing the argument of the logarithm,which corresponds to arbitrary µ (0) and µ (1) = µ (2) = . . . = µ ( N − = 0 . Thus, E eff (cid:126)µ ( t ) is minimized for a ∝ v (0) ,which implies | Ψ (cid:105) ∝ (cid:80) j v (0) j | Φ j (cid:105) .This is, what we wanted to show: The coefficients of the trial state (A4) with minimal effective energy at temporalseparation t are identical to the components of the eigenvector v (0) from the GEVP with t = t − a . ACKNOWLEDGMENTS
We acknowledge useful discussions with A. Ali, J. Kämper, M. Pflaumer and G. Schierholz.M.W. acknowledges support by the Heisenberg Programme of the Deutsche Forschungsgemeinschaft (DFG, GermanResearch Foundation) – project number 399217702.Calculations on the GOETHE-HLR and on the on the FUCHS-CSC high-performance computers of the FrankfurtUniversity were conducted for this research. We would like to thank HPC-Hessen, funded by the State Ministry ofHigher Education, Research and the Arts, for programming advice. [1] R. L. Jaffe, Phys. Rev.
D15 , 267 (1977).[2] P. Bicudo and M. Cardoso, Phys. Rev.
D94 , 094032 (2016), arXiv:1509.04943 [hep-ph].[3] S. Prelovsek, C. Lang, and D. Mohler, in
Mini-Workshop Bled 2011: Understanding Hadronic Spectra (2011) pp. 73–81,arXiv:1110.4520 [hep-ph]. [4] C. Morningstar, J. Bulava, B. Fahy, J. Fallica, A. Hanlon, B. Hörz, K. Juge, and C. H. Wong, Acta Phys. Polon. Supp. , 421 (2016).[5] A. Bondar et al. (Belle), Phys. Rev. Lett. , 122001 (2012), arXiv:1110.2251 [hep-ex].[6] Z. Q. Liu et al. (Belle), Phys. Rev. Lett. , 252002 (2013), arXiv:1304.0121 [hep-ex].[7] K. Chilikin et al. (Belle), Phys. Rev. D90 , 112009 (2014), arXiv:1408.6457 [hep-ex].[8] T. Xiao, S. Dobbs, A. Tomaradze, and K. K. Seth, Phys. Lett.
B727 , 366 (2013), arXiv:1304.3036 [hep-ex].[9] M. Ablikim et al. (BESIII), Phys. Rev. Lett. , 252001 (2013), arXiv:1303.5949 [hep-ex].[10] M. Ablikim et al. (BESIII), Phys. Rev. Lett. , 132001 (2014), arXiv:1308.2760 [hep-ex].[11] M. Ablikim et al. (BESIII), Phys. Rev. Lett. , 242001 (2013), arXiv:1309.1896 [hep-ex].[12] M. Ablikim et al. (BESIII), Phys. Rev. Lett. , 022001 (2014), arXiv:1310.1163 [hep-ex].[13] M. Ablikim et al. (BESIII), Phys. Rev. Lett. , 212002 (2014), arXiv:1409.6577 [hep-ex].[14] R. Aaij et al. (LHCb), Phys. Rev. Lett. , 222002 (2014), arXiv:1404.1903 [hep-ex].[15] R. Maciuła, V. A. Saleev, A. V. Shipilova, and A. Szczurek, Phys. Lett. B , 458 (2016), arXiv:1601.06981 [hep-ph].[16] R. Aaij et al. (LHCb), Phys. Rev. Lett. , 112001 (2017), arXiv:1707.01621 [hep-ex].[17] A. Ali, A. Y. Parkhomenko, Q. Qin, and W. Wang, Phys. Lett. B , 412 (2018), arXiv:1805.02535 [hep-ph].[18] A. Ali, Q. Qin, and W. Wang, Phys. Lett. B , 605 (2018), arXiv:1806.09288 [hep-ph].[19] J. P. Ader, J. M. Richard, and P. Taxil, Phys. Rev.
D25 , 2370 (1982).[20] J. l. Ballot and J. M. Richard, Phys. Lett.
B123 , 449 (1983).[21] L. Heller and J. A. Tjon, Phys. Rev.
D35 , 969 (1987).[22] J. Carlson, L. Heller, and J. A. Tjon, Phys. Rev.
D37 , 744 (1988).[23] H. J. Lipkin, Phys. Lett.
B172 , 242 (1986).[24] D. M. Brink and F. Stancu, Phys. Rev.
D57 , 6778 (1998).[25] B. A. Gelman and S. Nussinov, Phys. Lett.
B551 , 296 (2003), arXiv:hep-ph/0209095 [hep-ph].[26] J. Vijande, F. Fernandez, A. Valcarce, and B. Silvestre-Brac, Eur. Phys. J.
A19 , 383 (2004), arXiv:hep-ph/0310007[hep-ph].[27] D. Janc and M. Rosina, Few Body Syst. , 175 (2004), arXiv:hep-ph/0405208 [hep-ph].[28] T. D. Cohen and P. M. Hohler, Phys. Rev. D74 , 094003 (2006), arXiv:hep-ph/0606084 [hep-ph].[29] J. Vijande, A. Valcarce, and J. M. Richard, Phys. Rev.
D76 , 114013 (2007), arXiv:0707.3996 [hep-ph].[30] M. Born and R. Oppenheimer, Annalen der Physik , 457 (1927).[31] E. Braaten, C. Langmack, and D. H. Smith, Phys. Rev. D , 014044 (2014), arXiv:1402.0438 [hep-ph].[32] W. Detmold, K. Orginos, and M. J. Savage, Phys. Rev. D76 , 114503 (2007), arXiv:hep-lat/0703009 [hep-lat].[33] M. Wagner (ETM), PoS
LATTICE2010 , 162 (2010), arXiv:1008.1538 [hep-lat].[34] G. Bali and M. Hetzenegger (QCDSF), PoS
LATTICE2010 , 142 (2010), arXiv:1011.0571 [hep-lat].[35] M. Wagner (ETM), Acta Phys. Polon. Supp. , 747 (2011), arXiv:1103.5147 [hep-lat].[36] Z. S. Brown and K. Orginos, Phys. Rev. D86 , 114506 (2012), arXiv:1210.1953 [hep-lat].[37] P. Bicudo, K. Cichy, A. Peters, and M. Wagner, Phys. Rev.
D93 , 034501 (2016), arXiv:1510.03441 [hep-lat].[38] P. Bicudo and M. Wagner, Phys. Rev.
D87 , 114511 (2013), arXiv:1209.6274 [hep-ph].[39] P. Bicudo, K. Cichy, A. Peters, B. Wagenbach, and M. Wagner, Phys. Rev.
D92 , 014507 (2015), arXiv:1505.00613[hep-lat].[40] P. Bicudo, J. Scheunert, and M. Wagner, Phys. Rev.
D95 , 034502 (2017), arXiv:1612.02758 [hep-lat].[41] A. Francis, R. J. Hudspith, R. Lewis, and K. Maltman, Phys. Rev. Lett. , 142001 (2017), arXiv:1607.05214 [hep-lat].[42] A. Francis, R. J. Hudspith, R. Lewis, and K. Maltman, Phys. Rev. D , 054505 (2019), arXiv:1810.10550 [hep-lat].[43] P. Junnarkar, N. Mathur, and M. Padmanath, Phys. Rev. D , 034507 (2019), arXiv:1810.12285 [hep-lat].[44] L. Leskovec, S. Meinel, M. Pflaumer, and M. Wagner, Phys. Rev. D , 014503 (2019), arXiv:1904.04197 [hep-lat].[45] R. Hudspith, B. Colquhoun, A. Francis, R. Lewis, and K. Maltman, Phys. Rev. D , 114506 (2020), arXiv:2006.14294[hep-lat].[46] P. Bicudo, Nucl. Phys. A , 537 (2005), arXiv:hep-ph/0401106.[47] A. K. Rai, J. Pandya, and P. Vinodkummar, Nucl. Phys. A , 406 (2007), arXiv:hep-ph/0612244.[48] B. Chakrabarti, A. Bhattachary, and S. Mani, Phys. Scripta , 025103 (2009), arXiv:0806.4036 [hep-ph].[49] S. Godfrey and S. L. Olsen, Ann. Rev. Nucl. Part. Sci. , 51 (2008), arXiv:0801.3867 [hep-ph].[50] S. H. Lee, K. Morita, and M. Nielsen, Phys. Rev. D , 076001 (2008), arXiv:0808.3168 [hep-ph].[51] N. Brambilla et al. , Eur. Phys. J. C , 1534 (2011), arXiv:1010.5827 [hep-ph].[52] R. Faccini, A. Pilloni, and A. D. Polosa, Mod. Phys. Lett. A , 1230025 (2012), arXiv:1209.0107 [hep-ph].[53] A. Esposito, M. Papinutto, A. Pilloni, A. Polosa, and N. Tantalo, Phys. Rev. D , 054029 (2013), arXiv:1307.2873[hep-ph].[54] F.-K. Guo, L. Liu, U.-G. Meissner, and P. Wang, Phys. Rev. D , 074506 (2013), arXiv:1308.2545 [hep-lat].[55] W. Xie, L. Mo, P. Wang, and S. R. Cotanch, Phys. Lett. B , 148 (2013), arXiv:1302.5737 [hep-ph].[56] A. Esposito, A. Guerrieri, and A. Pilloni, Phys. Lett. B , 194 (2015), arXiv:1409.3551 [hep-ph].[57] S. R. Ignjatović and V. B. Jovanović, Facta Univ. Ser. Phys. Chem. Tech. , 151 (2014).[58] G. Eichmann, C. S. Fischer, and W. Heupel, Phys. Lett. B , 282 (2016), arXiv:1508.07178 [hep-ph].[59] S. J. Brodsky and R. F. Lebed, Phys. Rev. D , 114025 (2015), arXiv:1505.00803 [hep-ph].[60] L. Maiani, A. Polosa, and V. Riquer, JHEP , 160 (2016), arXiv:1605.04839 [hep-ph].[61] A. Esposito, A. Pilloni, and A. Polosa, Phys. Lett. B , 292 (2016), arXiv:1603.07667 [hep-ph].[62] G. Rupp and E. van Beveren, Chin. Phys. C , 053104 (2017), arXiv:1611.00793 [hep-ph]. [63] A. Ali, J. S. Lange, and S. Stone, Prog. Part. Nucl. Phys. , 123 (2017), arXiv:1706.00610 [hep-ph].[64] J. Sonnenschein and D. Weissman, Eur. Phys. J. C , 326 (2019), arXiv:1812.01619 [hep-ph].[65] J.-M. Richard, A. Valcarce, and J. Vijande, Annals Phys. , 168009 (2020), arXiv:1910.08295 [nucl-th].[66] N. Brambilla, S. Eidelman, C. Hanhart, A. Nefediev, C.-P. Shen, C. E. Thomas, A. Vairo, and C.-Z. Yuan, Phys. Rept. , 1 (2020), arXiv:1907.07583 [hep-ex].[67] L. Maiani, A. D. Polosa, and V. Riquer, Phys. Rev. D , 074002 (2019), arXiv:1908.03244 [hep-ph].[68] Y. Liu, M. A. Nowak, and I. Zahed, Phys. Rev. D , 126023 (2019), arXiv:1904.05189 [hep-ph].[69] H.-X. Chen, W. Chen, R.-R. Dong, and N. Su, Chin. Phys. Lett. , 101201 (2020), arXiv:2008.07516 [hep-ph].[70] G. Yang, J. Ping, and J. Segovia, Phys. Rev. D , 054023 (2020), arXiv:2007.05190 [hep-ph].[71] M. Gavela, A. Le Yaouanc, L. Oliver, O. Pene, J. Raynal, and S. Sood, Phys. Lett. B , 431 (1979).[72] H. Yukawa, Proc. Phys. Math. Soc. Jap. , 48 (1935).[73] F. Karbstein, A. Peters, and M. Wagner, JHEP , 114 (2014), arXiv:1407.7503 [hep-ph].[74] S. Prelovsek, C. Lang, L. Leskovec, and D. Mohler, Phys. Rev. D , 014504 (2015), arXiv:1405.7623 [hep-lat].[75] D. Mohler, C. Lang, L. Leskovec, S. Prelovsek, and R. Woloshyn, Phys. Rev. Lett. , 222001 (2013), arXiv:1308.3175[hep-lat].[76] C. Lang, L. Leskovec, D. Mohler, S. Prelovsek, and R. Woloshyn, Phys. Rev. D , 034510 (2014), arXiv:1403.8103[hep-lat].[77] G. S. Bali, S. Collins, A. Cox, and A. Schäfer, Phys. Rev. D , 074501 (2017), arXiv:1706.01247 [hep-lat].[78] C. Alexandrou, J. Berlin, M. Dalla Brida, J. Finkenrath, T. Leontiou, and M. Wagner, Phys. Rev. D , 034506 (2018),arXiv:1711.09815 [hep-lat].[79] C. Alexandrou, J. Berlin, J. Finkenrath, T. Leontiou, and M. Wagner, Phys. Rev. D , 034502 (2020), arXiv:1911.08435[hep-lat].[80] A. Abdel-Rehim, C. Alexandrou, J. Berlin, M. Dalla Brida, J. Finkenrath, and M. Wagner, Comput. Phys. Commun. , 97 (2017), arXiv:1701.07228 [hep-lat].[81] N. Cardoso, M. Cardoso, and P. Bicudo, Phys. Rev. D , 054508 (2011), arXiv:1107.1355 [hep-lat].[82] M. Cardoso, N. Cardoso, and P. Bicudo, Phys. Rev. D , 014503 (2012), arXiv:1204.5131 [hep-lat].[83] P. Bicudo, M. Cardoso, O. Oliveira, and P. Silva, Phys. Rev. D , 074508 (2017), arXiv:1702.07789 [hep-lat].[84] K. Jansen, C. Michael, A. Shindler, and M. Wagner (ETM), JHEP , 058 (2008), arXiv:0810.1843 [hep-lat].[85] C. Michael, A. Shindler, and M. Wagner (ETM), JHEP , 009 (2010), arXiv:1004.4235 [hep-lat].[86] In our most recent and refined computation, which includes an extrapolation to physically light u/d quark masses and acoupled channel Schrödinger equation with a 2-component wave function, we find for the tetraquark mass
10 545 +38 − MeV.This corresponds to binding energy − E = 59 +30 − MeV with respect to the BB ∗ threshold (see Ref. [40] for details).[87] R. Frezzotti, P. A. Grassi, S. Sint, and P. Weisz (Alpha), JHEP , 058 (2001), arXiv:hep-lat/0101001.[88] R. Frezzotti and G. Rossi, JHEP , 007 (2004), arXiv:hep-lat/0306014.[89] P. Weisz, Nucl. Phys. B , 1 (1983).[90] F. Farchioni, C. Urbach, R. Frezzotti, K. Jansen, I. Montvay, G. Rossi, E. Scholz, A. Shindler, N. Ukita, and I. Wetzorke,Nucl. Phys. B Proc. Suppl. , 240 (2005), arXiv:hep-lat/0409098.[91] F. Farchioni, K. Jansen, I. Montvay, E. Scholz, L. Scorzato, A. Shindler, N. Ukita, C. Urbach, and I. Wetzorke, Eur.Phys. J. C , 73 (2005), arXiv:hep-lat/0410031.[92] R. Frezzotti, G. Martinelli, M. Papinutto, and G. Rossi, JHEP , 038 (2006), arXiv:hep-lat/0503034.[93] K. Jansen, M. Papinutto, A. Shindler, C. Urbach, and I. Wetzorke (XLF), JHEP , 071 (2005), arXiv:hep-lat/0507010.[94] P. Boucaud et al. (ETM), Phys. Lett. B , 304 (2007), arXiv:hep-lat/0701012.[95] P. Boucaud et al. (ETM), Comput. Phys. Commun. , 695 (2008), arXiv:0803.0224 [hep-lat].[96] R. Baron et al. (ETM), JHEP , 097 (2010), arXiv:0911.5061 [hep-lat].[97] M. Albanese et al. (APE), Phys. Lett. B , 163 (1987).[98] S. Gusken, Nucl. Phys. B Proc. Suppl. , 361 (1990).[99] A. Hasenfratz and F. Knechtli, Phys. Rev. D , 034504 (2001), arXiv:hep-lat/0103029.[100] M. Della Morte, S. Durr, J. Heitger, H. Molke, J. Rolf, A. Shindler, and R. Sommer (ALPHA), Phys. Lett. B , 93(2004), [Erratum: Phys.Lett.B 612, 313–314 (2005)], arXiv:hep-lat/0307021.[101] M. Della Morte, A. Shindler, and R. Sommer, JHEP , 051 (2005), arXiv:hep-lat/0506008.[102] Static-static correlation functions computed with lattice QCD exhibit strong discretization errors for separations r/a < ∼ ,as is e.g. well known in the context of computations of Wilson loops and the ordinary static potential. Thus, throughoutthis work we only show and discuss results for r/a ≥ .[103] M. Kalinowski and M. Wagner, Phys. Rev. D , 094508 (2015), arXiv:1509.02396 [hep-lat].[104] B. Blossier, M. Della Morte, G. von Hippel, T. Mendes, and R. Sommer, JHEP , 094 (2009), arXiv:0902.1265 [hep-lat].[105] P. Zyla et al. (Particle Data Group), PTEP2020