Importance of meson-meson and of diquark-antidiquark creation operators for a $\bar{b} \bar{b} u d$ tetraquark

In recent years, the existence of a hadronically stable $\bar{b} \bar{b} u d$ tetraquark with quantum numbers $I(J^P) = 0(1^+)$ was confirmed by first principles lattice QCD computations. In this work we use lattice QCD to compare two frequently discussed competing structures for this tetraquark by considering meson-meson as well as diquark-antidiquark creation operators. We use the static-light approximation, where the two $\bar{b}$ quarks are assumed to be infinitely heavy with frozen positions, while the light $u$ and $d$ quarks are fully relativistic. By minimizing effective energies and by solving generalized eigenvalue problems we determine the importance of the meson-meson and the diquark-antidiquark creation operators with respect to the ground state. It turns out, that the diquark-antidiquark structure dominates for $\bar{b} \bar{b}$ separations $r<0.25 \, \text{fm}$, whereas it becomes increasingly more irrelevant for larger separations, where the $I(J^P) = 0(1^+)$ tetraquark is mostly a meson-meson state. We also estimate the meson-meson to diquark-antidiquark ratio of this tetraquark and find around $60\% / 40\%$.


I. INTRODUCTION
A long standing problem in particle physics is to understand exotic hadrons, i.e. hadrons which have a structure more complicated than a quark-antiquark pair or a triplet of quarks [1]. Such exotic hadrons turned out to be not only difficult to observe experimentally, but also technical to address in quark models [2]. In the last couple of years, several exotic hadrons, the majority pertaining to the class of tetraquarks with at least two heavy quarks, were clearly confirmed by the BELLE, BES-III and LHCb experimental collaborations. The observed exotic hadrons are resonances high in the spectrum. Studying them theoretically from first principles with lattice QCD requires the application 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 QQqq including the Z c and Z b states are easier to detect experimentally. Their observation at Belle [5][6][7], Cleo-C [8], BESIII [9][10][11][12][13] and LHCb [14] collaborations turned tetraquarks into a main highlight of particle physics in recent years. But since they are resonances with more than one decay channel, we study in this paper tetraquarks with two heavy antiquarksQQqq (or equivalently two heavy quarks, i.e. QQqq), which are theoretically simpler, because they are either hadronically stable or can only decay into a pair of heavy-light mesons. Moreover, with the recent observation of hadronic systems with two heavy quarks [15,16] at LHCb we expect this second class of tetraquarks to be observed in the near future. Their discovery potential is discussed in Refs. [17,18].
Recently this was confirmed with lattice QCD computations. One of the approaches uses the Born-Oppenheimer approximation [30,31], where the problem is split into two steps. The first step is to compute the potentials of two static antiquarks in the presence of two light quarks using state-of-the-art lattice QCD techniques (see e.g. Refs. [32][33][34][35][36][37]). Then, in the second step, the heavy quark dynamics is studied using a quantum mechanical Hamiltonian with the previously computed lattice QCD potentials. Using this approach, abbud tetraquark bound state with quantum numbers I(J P ) = 0(1 + ) was predicted [36][37][38][39][40]. Shortly afterwards, this was confirmed by several full lattice QCD computations using four quarks of finite mass [41][42][43][44][45].
Our present goal is to explore the structure of thisbbud tetraquark with quantum numbers I(J P ) = 0(1 + ) using first principles lattice QCD computations, and to answer a hotly debated theoretical question : Is it a diquarkantidiquark 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 thebbud 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][36][37]). This is consistent with the following. (i) Atbb separations larger than the typical meson radius each of the two heavy antiquarks 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 smallbb separations the heavy antiquarks interact directly via gluon exchange and are immersed in a light quark cloud. Thus at 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 givenbb separation, we utilize a set of lattice QCD creation operators, both of Dd type and of BB type with staticb quarks. A similar approach was previously used to explore systems of four quarks of finite mass. For example, studies of four-quark systems including a heavycc pair were to some extent inconclusive. Neither clear evidence for the existence of a tetraquark was found, nor a signal improvement was observed for diquark-antidiquark operators [53,74]. Also the tetraquark candidates a 0 (980) and the D * s0 (2317) were investigated in that way, employing sets of different creation operators, including quarkantiquark, meson-meson and/or diquark-antidiquark type [75][76][77][78][79]. The focus of these studies was more to distinguish between 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 or diquark-antidiquark dominance was reached. On the other hand, studies of potentials and the corresponding gluon field distributions were performed with four static quarks, which are close to a system of four bottom quarksbbbb [81][82][83]. In this case it was possible to clearly distinguish between a diquark-antidiquark structure and a meson-meson structure 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 setup in section III. In Section IV we present our numerical results concerning the relative importance of a Dd structure and of a BB structure at givenbb separation. At the end of this section we use these results to crudely estimate the percentage of Dd and of BB in thebbud 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 QUARKS AND CORRESPONDING CREATION OPERATORS
In preceding papers [33,35,37,39] we have computed potentials V qq,jz,P,Px (r) of two static antiquarksQQ at separation r in the presence of two light quarks qq using lattice QCD. The computations have been carried out for many different sectors characterized by the following quantum numbers: light flavor qq with q ∈ {u, d, s, c}, total angular momentum of the light quarks and gluons j z with respect to theQQ separation axis, parity P and reflection along an axis perpendicular to theQQ separation axis P x . There are both attractive and repulsive sectors. Most promising with respect to the existence ofQQqq tetraquark bound states or resonances are attractive potentials with light 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 radial coordinate of the two heavy quarksQQ =bb with the computed potentials V qq,jz,P,Px (r), one can explore the existence of hadronically stable tetraquarks. m b is the b quark mass, L is the relative orbital angular momentum ofbb pair and m sl is the mass of the lightest static-light meson (computed within the same lattice QCD setup as V qq,jz,P,Px (r); see e.g. Refs. [84,85]). There is one particular potential V (r) = V ud−du,0,−,+ (r) (shown in the left plot of Figure 1), which has quantum numbers (I, j z , P, P x ) = (0, 0, −, +), leading for L = 0 to a stablebbud tetraquark with quantum numbers I(J P ) = 0(1 + ) and binding energy −E = 38 (18) MeV [86]. The probability density of thebb separation 4π|R(r)| 2 (shown in the right plot of Figure 1) indicates that one typically finds separations in the range 0.1 fm < ∼ r < ∼ 0.6 fm. None of the other potentials is sufficiently wide or deep to host a bound state [39]. As discussed in the introduction, the main goal of this work is to investigate the structure of the predictedbbud tetraquark. In particular, we explore, whether the tetraquark is more similar to a meson-meson state BB or to a diquark-antidiquark state Dd, two scenarios frequently discussed in the literature and at conferences also for other tetraquark candidates (see the discussion in section I). To this end, we refine the existing lattice QCD computation of 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], excites two B mesons at separation r, with r = |r 2 −r 1 |, color indices a, b, spin indices A, B, C, D and ψ (f ) ψ (f ) = ud−du. N BB is a normalization, which will be discussed later. There are two independent choices for the light spin matrix Γ consistent with (j z , P, P x ) = (0, −, +). Γ = (1 + γ 0 )γ 5 predominantly excites two negative parity ground state mesons B ( * ) B ( * ) , while Γ = (1 − γ 0 )γ 5 mostly generates two positive parity excited mesons B * 0,1 B * 0,1 , 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 be chosen arbitrarily, The second type of creation operators, which we use here for the first time, resembles a diquark-antidiquark pair with heavy quarks separated by r and connected by a gluonic string, Again N Dd is a normalization and the allowed light and heavy spin matrices are the same as for the operator 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 correlation matrix. For definiteness, we chooseΓ = (1 − γ 0 )γ 3 . r 1 and r 2 are always separated along one of the lattice axes and z = (r 1 + r 2 )/2. For odd r/a (a denotes the lattice spacing), z does not coincide with one of the lattice sites. In these cases we take the average of the operator (3) with z = (r 1 + r 2 )/2 + ar/2 and with z = (r 1 + r 2 )/2 − ar/2, wherê r = (r 2 − r 1 )/|r 2 − r 1 |. We use smeared light quark and gluon fields, which implies that the operator ψ (f )a A (r) does not generate a point-like excitation at r, but rather a cloud-like excitation of spherical shape with diameter ≈ 0.5 fm. Similarly, the gauge links connecting the two static antiquarks in the operator O Dd,Γ generate a flux tube with a certain thickness. For details see 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 |Ω . For r = 0, the operators are even identical, when properly normalized, i.e. O BB,Γ = O Dd,Γ , if N BB = N Dd . This can easily be shown by using the identity abc ade = δ bd δ ce − δ be δ cd . For increasing r, however, they become more and more different, as we will show numerically in section IV A. One obvious reason for that is that O Dd,Γ creates a flux tube of length r, whereas O BB,Γ does not. In section IV B, section IV C and section IV D we will explore, whether the ground state of the (I, j z , P, P x ) = (0, 0, −, +) sector at a given separation r is more similar to a BB state O Dd,Γ |Ω or to a DD state O Dd,Γ |Ω . Since the energy of this ground state is the potential V (r) used in the Born-Oppenheimer prediction of thebbud tetraquark with quantum numbers I(J P ) = 0(1 + ), our investigation will provide information on the structure of that tetraquark. We will discuss that in section IV E.

A. Lattice actions
The light quark action used in this work is the Wilson twisted mass action [87,88], where D W is the standard Wilson Dirac operator, with the forward and backward covariant derivatives ∇ µ and ∇ * µ . χ = (χ (u) , χ (d) ) is the light quark doublet in the so-called twisted basis, which is related to quark fields in the usual "physical basis" via the twist rotation where ω is the twist angle. The gluon action used in this work is the tree-level Symanzik improved action [89], where b 1 = −1/12, b 0 = 1 − 8b 1 , β = 6/g 2 0 , g 0 is the bare coupling and P 1×1 µν and P 1×2 µν are the plaquette and a 1 × 2 Wilson loop, respectively.

B. Ensembles of gauge link configurations
We have performed computations on two ensembles of gauge link configurations generated by the European Twisted Mass Collaboration (ETMC) [94][95][96]. These ensembles have different lattice spacings, a ≈ 0.079 fm and a ≈ 0.063 fm (set via the pion mass and pion decay constant and chiral perturbation theory [96]), which allows to study a finer resolution 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 ensembles are collected in Table I.  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.5 (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.5 (see Ref. [84], Eq. (25)) and APE-smeared spatial gauge links. The intention of both APE and Gaussian smearing is to increase the ground state overlaps generated by the creation operators.
In addition to that we also use HYP2-smeared gauge links in temporal direction [99][100][101] with parameters α 1 = α 2 = 1.0 and α 3 = 0.5. This reduces the self energy of the static quarks and, thus, improves the signal-to-noise ratio of the computed correlation functions.

IV. NUMERICAL RESULTS
We consider the following four creation operators, and define the corresponding trial states as O BB,(1+γ0)γ5 predominantly excites two negative parity ground state mesons. Thus, |Φ BB,(1+γ0)γ5 is expected to have the largest overlap to the ground state of the (I, j z , P, P x ) = (0, 0, −, +) sector at large r. At small r, however, the diquark-antidiquark operators might be advantageous. We use |Φ Dd,γ5 , i.e. a light diquark with just γ 5 , as typically discussed in the literature. Since O BB,Γ ∝ O Dd,Γ for r = 0 (see section II), the diquark-antidiquark trial state |Φ Dd,(1+γ0)γ5 might even be a better candidate for having a large ground state overlap. For completeness we also include O BB,γ5 , the mesonic molecule counterpart of O Dd,γ5 . O BB,γ5 excites a linear combination of two negative parity ground state mesons and two significantly heavier positive parity excited mesons [37]. With these operators we computed the 4 × 4 correlation matrix, with . . . denoting the path integral expectation value and t/a = (t 2 − t 1 )/a ≥ 1. For the second equality we assumed that in the spectral decomposition propagation over temporal separation T − t is suppressed for all states except for the vacuum. To cross-check our computations, we checked the numerical results with respect to the symmetries γ 5 hermiticity, parity, time reversal, charge conjugation and cubic rotations around the axis of separation (for details see Ref. [37]). In a second step we averaged elements of the correlation matrices related by these symmetries, to reduce statistical errors.
A. Squared overlaps of the normalized BB and Dd trial states In this subsection we study For t → 0 this quantity is the squared normalized overlap of trial state Clearly, 0 ≤ α 0 jk ≤ 1. Note that for an arbitrary state |Ψ and an orthonormal basis |k , k = 1, 2, 3, . . .
Thus for two trial states |Φ j and |Φ k , α 0 jk can be interpreted as a measure of their orthogonality, where α 0 jk ≈ 0 indicates almost orthogonal and α 0 jk ≈ 1 almost parallel states. For large t all α jk approach 1, because the ground state dominates in that limit.
In the left plot of Figure 2 we show α jk for j = BB, (1 + γ 0 )γ 5 and k = Dd, (1 + γ 0 )γ 5 as function of t for several fixed r [102] for ensemble B40. 24. The corresponding trial states become more similar for smaller r, as indicated by larger values of α jk close to 1. 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+γ0)γ5 generates a pair of spatially separated B mesons, which interact only by weak residual hadronic forces. In contrast to that, the antiquarksQQ forming the heavy diquark in the operator O Dd,(1+γ0)γ5 are connected by a gluonic flux tube of length r. The data points for 1 ≤ t/a ≤ 5 can be fitted consistently with degree-2 polynomials, which are also shown in Figure 2. These fits represent crude extrapolations of α jk to t = 0, which are the squared overlaps of the corresponding normalized BB and Dd trial states. (right) j = BB, γ5, k = Dd, γ5. For t → 0, α 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, γ 5 and k = Dd, γ 5 . They are quite similar to those for j = BB, (1 + γ 0 )γ 5 and k = Dd, (1 + γ 0 )γ 5 and can be interpreted in the same way. The main point of these two plots is to demonstrate that BB and Dd trial states, even though not orthogonal, are not linearly dependent either. In the following subsections we explore, whether the ground state of the (I, j z , P, P x ) = (0, 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 large temporal separations Now we consider effective energies corresponding to diagonal elements of the correlation matrix, i.e.
(no sum over j).
As discussed in section II, all four operators probe the (I, j z , P, P x ) = (0, 0, −, +) sector. Thus, for fixed r and at sufficiently large t, all four V eff j (r, t) should approach the same constant, which is the ground state energy V (r). For separations r < ∼ 0.3 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 ≈ 0.16 fm as a function of t). For larger r, however, the two Dd operators seem to generate only little overlap to the ground state. The consequence is that the corresponding effective energies V eff Dd,(1+γ0)γ5 and V eff Dd,γ5 do not convincingly converge to the plateau at V (r) in the t region, where we carried out computations, and have rather 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 ≈ 0.79 fm as a function of t). In contrast to that, the two BB operators still lead to clear effective energy plateaus. Thus, if one is just interested to determine V (r), it is sufficient to implement BB operators. This is, what we did in previous work, e.g. in Ref. [38] (see also the left plot of Figure 1).  Figure 4. This difference is clearly negative for all separations r, which is not surprising. O BB,(1+γ0)γ5 creates predominantly a pair of ground state static-light mesons, while O BB,γ5 creates roughly a 50%/50% superposition of a pair of negative parity ground state mesons and a pair of significantly heavier positive parity static-light mesons, as discussed above and as can be shown e.g. by a Fierz transformation (for details see Ref. [37]). Results from an analog comparison of the two diquark-antidiquark operators are very similar (see upper right plot of Figure 4). This is interesting, because it indicates that a light diquark with spin structure given by (1 + γ 0 )γ 5 is energetically preferred over a light diquark with spin structure given just by γ 5 .
Most interesting, of course, is the comparison of a meson-meson and a diquark-antidiquark operator, specifically of O BB,(1+γ0)γ5 and O Dd,(1+γ0)γ5 , which we have just identified as being superior to O BB,γ5 and O Dd,γ5 , respectively. We show the difference of their effective masses, V eff BB,(1+γ0)γ5 (r, t = 2a) − V eff Dd,(1+γ0)γ5 (r, t = 2a), in the lower plot of Figure 4. For r < ∼ 3.15 a ≈ 0.25 fm the difference is positive, indicating that for small separations the diquarkantidiquark operator generates a trial state more similar to the ground state. For larger separations, r > ∼ 0.25 fm the difference becomes negative and strongly points towards a meson-meson structure. While the latter is expected (at large r the flux tube present in a heavy diquark-antidiquark is energetically disfavored), the former is a first hint towards a diquark-antidiquark dominance at smaller r. We continue to investigate this in more detail in the following subsections.
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+γ0)γ5 and |Φ Dd,(1+γ0)γ5 . Any trial state from that space can be written as with coefficients b, d ∈ C. To identify a trial state as similar to the ground state as possible, i.e. with large overlap to the ground state and little overlap to excitations, we minimize the corresponding effective energy with respect to b and d. Since, V eff b,d is independent of the norm and the phase of |Φ b,d , 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 is numerically straightforward. In particular, no further lattice QCD computations are needed, because C [b,d] [b,d] can be expressed in terms of the correlation matrix introduced in Eq. (10), In the following we consider with b = 1 and d minimizing V eff b,d . w BB and w Dd are the normalized absolute squares of the coefficients of the optimized trial states appearing in Eq. (15). These quantities exhibit only a weak dependence on t. For 3 ≤ t/a ≤ 5 and ensemble B40.24 (4 ≤ t/a ≤ 6 and ensemble C30.32) they are consistent with a constant. For t/a ≥ 6 (t/a ≥ 7), 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 correlation functions, which are exponentially suppressed in t. In Figure 5 we show example plots of w BB and w Dd as functions of t for selected separations r/a = 2, r/a = 5 and r/a = 8 for ensemble B. We determine each plateau value by a χ 2 minimizing fit of a constant in the range 3 ≤ t/a ≤ 5 (4 ≤ t/a ≤ 6). The resulting numbers,w BB (r) andw Dd (r) = 1 −w BB (r), can be interpreted as the relative weight of a meson-meson and a diquark-antidiquark structure atbb separation r in the ground state, which corresponds to the potential V (r) of two static antiquarks and is, thus, closely related to thebbud tetraquark with quantum numbers I(J P ) = 0(1 + ). In Figure 6 we plotw BB andw Dd as functions of r. One can clearly see that there is a diquark-antidiquark dominance forbb separations r < ∼ 0.20 fm. For r > ∼ 0.30 fm, the meson-meson structure is more prominent and for r > ∼ 0.50 fm, the diquark-antidiquark contribution is negligible, i.e. the system is exclusively composed of two B mesons. It is interesting to note that the separation r ≈ 0.3 fm, where the meson-meson structure starts to dominate, is of the same order as the size of a B meson. A precise comparison, however, seems to be difficult, because the size of a B meson is model dependent and not precisely known (see e.g. Refs. [103,104]).
We note that the quantities w BB and w Dd as well as the fittedw BB andw Dd depend on the normalization of the operators O BB,(1+γ0)γ5 and O Dd,(1+γ0)γ5 or, equivalently, on the normalization and the corresponding states |Φ BB,(1+γ0)γ5 and |Φ Dd,(1+γ0)γ5 . To allow a meaningful interpretation in terms of the relative weight of a BB and a 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,Γ = |Φ Dd,Γ for r = 0 (see section II). We expect that N BB = N Dd also results in similar norms for r > 0. A common alternative, which we have used in previous lattice QCD projects, is to normalize the operators O j such that C jj (t = a) = 1 (no sum over j; see e.g. Refs. [78,79,105]). As a cross-check we also explored this normalization in our current work and found almost identical results to those obtained with N BB = N Dd .

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, 0, −, +) sector, we now use N × N correlation matrices C jk as defined in Eq. (10) and solve the generalized eigenvalue problem (GEVP) for t 0 /a ≥ 1 and t/a > t 0 /a (for detailed discussions of the GEVP in lattice field theory see e.g. Refs. [106][107][108][109][110][111][112][113]). Effective energies for the lowest N energy eigenstates are then given by . (20) These are generalizations of the effective energy defined in Eq. (14), since V eff,(0) (r, t) = V eff j (r, t) for N = 1. For N > 1, V eff,(0) (r, t) approaches the same constant V (r) for large t, but plateaus can typically be identified at somewhat 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 where the ≈ sign denotes an approximate expansion of the energy eigenstate |n in terms of the trial states |Φ j . For such values of t and t 0 , the squared eigenvector components as functions of t form plateaus and we determine the corresponding asymptotic values of |v  j | 2 depend on the normalization of the creation operators, as it is the case for w BB , w Dd ,w BB andw 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 with similar norm.
It is interesting to note that for a 2 × 2 correlation matrix with trial states |Φ 1 = |Φ BB,(1+γ0)γ5 and |Φ 2 = |Φ Dd,(1+γ0)γ5 and t 0 /a = t/a − 1, the eigenvector components v Dd,(1+γ0)γ5 | 2 ) = (w BB , w Dd ) and, thus, (|v Dd,(1+γ0)γ5 | 2 ) = (w BB ,w Dd ), as we show in appendix A. In other words, the results on the structure of thebbud 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 GEVP eigenvector components from another perspective. Compared to the trial state optimization from section IV C, the GEVP, however, offers further possibilities, for example to choose t 0 independent of t (i.e. not as t 0 /a = t/a − 1) or to study the full correlation matrix with N = 4.
As discussed in the previous paragraph, solving the GEVP with a 2 × 2 correlation matrix including the operators O BB,(1+γ0)γ5 and O Dd,(1+γ0)γ5 and using t 0 /a = t/a − 1 yields exactly the same results as shown in Figure 6 (one just has to replace labels according tow BB → |v Dd,(1+γ0)γ5 | 2 ). However, many lattice QCD papers using the GEVP fix t 0 to a rather small value independent of t. For small values of t 0 , statistical errors are somewhat reduced, which in turn allows to consider larger values of t. Thus we also computed |v Dd,(1+γ0)γ5 | 2 for t 0 /a = 1. As already observed in the previous subsection for the related quantities w BB and w Dd , the resulting squared eigenvector components |v Dd,(1+γ0)γ5 | 2 . For ensemble B40.24 the analysis and interpretation of the data is less obvious. While for t/a ≤ 5 plateaus are indicated, for t/a ≥ 6 there is a weak, almost linearly increasing deviation from these plateaus for some separations r (see e.g. the plot in the center of the upper line of Figure 7). Note, however, that errors are also increasing at larger t. For example, the data points at 8 ≤ t/a ≤ 9 are just around 2σ away from those at 4 ≤ t/a ≤ 5. Also the monotonic almost linear behavior is not necessarily an indication of a systematic deviation from a constant, because all data points were computed on the same gauge link configurations and neighboring points in t are, thus, correlated. Consequently, we interpret the observed deviations as statistical fluctuations and use the fit range 4 ≤ t/a ≤ 5.
In the left plot of Figure 8 we show |v Dd,(1+γ0)γ5 | 2 for both ensembles. These curves for t 0 /a = 1 are quite similar to those corresponding to t 0 /a = t/a−1 (and shown in Figure 6). Again, one can see a clear dominance of the diquark-antidiquark operator for separations r < ∼ 0.20 fm. In the range 0.20 fm ≤ r ≤ 0.30 fm there is a rapid change towards a meson-meson structure. For r > ∼ 0.50 fm there is almost no diquark-antidiquark contribution anymore and thebbud four-quark system seems to be composed exclusively of two B mesons. This plot confirms our results from section IV B obtained by minimizing effective energies.
It should be noted that the quantitiesw BB andw Dd as well as the fitted squared eigenvector components |v (0) BB,(1+γ0)γ5 | 2 and |v (0) Dd,(1+γ0)γ5 | 2 depend on the creation operators used, e.g. the details of the quark and gauge field smearing or their normalization. Besides discretization errors, this might be part of the reason for the slight differences between the results obtained for ensemble B40.24 and ensemble C30.32 shown in Figure 6 and the left plot of Figure 8. Since these differences are quite small and since we found almost identical results for different normalization prescriptions, we considerw BB andw Dd as well as |v (0) j | 2 as reliable indicators characterizing the quark and gluon structure of the ground state of the (I, j z , P, P x ) = (0, 0, −, +) sector.
Finally we performed the same GEVP analysis using the full 4 × 4 correlation matrix including all operators defined in Eq. (8). We determined |v (0) j | 2 using the same fit ranges as for the previous 2 × 2 analyses (results for ensemble C30.32 are shown in the right plot of Figure 8). Again, there is a diquark-antidiquark dominance for r < ∼ 0.20 fm, this time represented by the eigenvector components corresponding to two operators O Dd,(1+γ0)γ5 and O Dd,γ5 , while there is meson-meson dominance for r > ∼ 0.30 fm, reflected by the eigenvector components corresponding to two operators O BB,(1+γ0)γ5 and O BB,γ5 . When adding the two diquark eigenvector components as well as the two meson-meson eigenvector components, i.e. when considering |v BB,γ5 | 2 , we find within errors the same curves as obtained above by the 2 × 2 analysis. This is reassuring, because the results concerning the meson-meson percentage and the diquark-antidiquark percentage do not change, even though we increased the dimension of the basis of trial states used to approximate the ground state from 2 to 4.
We note that the GEVP can also be used to study excited states. In our case, i.e. for quantum numbers (I, j z , P, P x ) = (0, 0, −, +), the first and second excitation correspond to the repulsive potential of a negative and a positive parity B meson and to the attractive potential of two positive parity B mesons. We computed these potentials in a pevious work [37] (see in particular Figure 4 in Ref. [37], "singlet A"), where operators O BB,Γ with a larger set of matrices Γ were used. To study the structure of these excitations would require also a comparable set of operators During the lattice QCD computation of the potential V (r) the heavy antiquarksbb 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 dynamics of the heavy quarks can, however, be studied in a second step, by inserting the potential V (r) into the Schrödinger equation (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 [114] and found a single bound state with binding energy −E = 38 (18) MeV indicating the existence of a hadronically stablebbud tetraquark with quantum numbers I(J P ) = 0(1 + ). Moreover, the wave function R(r)/r gives the probability density of thebb separation, p r (r) = 4π|R(r)| 2 , which is shown in the right plot of Figure 1.
The quantitiesw BB andw Dd (see Figure 6) as well as |v Dd,(1+γ0)γ5 | 2 (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 (22) and where 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 thebbud tetraquark, which we denote by %BB j and %Dd j . To this end, we use a parameterization, which is a simple mathematical function able to consistently describe our lattice QCD results, where the parameters r 0,j and α j are determined by χ 2 minimizing fits to the results of both ensembles as shown in Figure 6 and Figure 8. Then we compute %BB and %Dd via We find %BB w = 0.58, %Dd w = 0.42 and %BB v = 0.60, %Dd v = 0.40, i.e. almost the same result, when usingw BB andw Dd and when using |v Dd,(1+γ0)γ5 | 2 . Note that we have lattice QCD results for p BB,j (r) and p Dd,j (r) only for separations r > ∼ 0.1 fm. Thus, it is unclear, whether the parameterization (24) is a valid description also for r < ∼ 0.1 fm. The coresponding systematic error is, however, quite small, because the probability to find thebb pair at separation r < ∼ 0.1 fm is also rather small (see the plot of p r (r) = 4π|R(r)| 2 in Figure 1). To quote a crude and very conservative upper bound for this systematic error, we solved again the integral in Eq. (25) replacing in the interval 0 ≤ r ≤ 0.1 fm the almost vanishing p BB,j (r) = (tanh(α j (r − r 0,j ) + 1)/2 ≈ 0 by p BB,j (r) = 1. The results are quite similar, %BB w = 0.63, %Dd w = 0.37 and %BB v = 0.65, %Dd v = 0.35, indicating that the corresponding systematic error is well below 0.05. Moreover, the eigenvector components are slightly operator dependent, as discussed in section IV D. Thus, the percentages %BB j and %Dd j should only be considered as crude estimates. As total systematic error, reflecting both the parameterization and the operator dependence, we estimate ≈ 0.10. Still it seems to be clear that thebbud tetraquark with quantum numbers I(J P ) = 0(1 + ) is neither strongly meson-meson dominated nor strongly diquark-antidiquark dominated, but rather an approximately equal linear combination of both structures. Finally we note that these results for %BB j and for %Dd j are fully consistent with squared eigenvector components obtained during a recent lattice QCD study [44] of the same tetraquark using four quarks of finite mass. There, a meson-meson component of 0.65(4) and a diquark-antidiquark component of 0.35(4) was found [115].

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 ā bbud tetraquark with quantum numbers I(J P ) = 0(1 + ): Are they more similar to meson-meson systems or rather to diquark-antidiquark pairs? Moreover we addressed the Dirac structure of the light quarks, comparing γ 5 with (1 + γ 0 )γ 5 , 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 diquarkantidiquark type. We solved the GEVP both for 2 × 2 and 4 × 4 correlation matrices, to determine quantitatively, which of the implemented structures is preponderant in the ground state atbb separation r. Moreover, we optimized trial states by minimizing effective energies and proved that this is equivalent to solving a specific corresponding GEVP.
Notice that the question we are addressing is quite subtle. We first showed that the BB and Dd trial states are not orthogonal. Nevertheless, since they are not linearly dependent either, we were able to determine, which one is more similar to the ground state. In what concerns light spin we found the (1 + γ 0 )γ 5 is the dominant Dirac structure at all separations r, both for BB and for Dd. This is what we expected, since the less favorable γ 5 structure generates not only negative parity mesons, but also excited positive parity mesons. As for color we showed that at small separations r < ∼ 0.25 fm the diquark-antidiquark structure dominates, whereas at larger separations the meson-meson trial state is clearly more similar to the ground state of thebbud system. For r > ∼ 0.5 fm, the percentage of BB is already larger than 95% and approaches 100% for even larger r.
Using these results as well as the wave function of thebb separation already obtained in Ref. [38], we estimated the meson-meson to diquark-antidiquark ratio of thebbud tetraquark with quantum numbers I(J P ) = 0(1 + ) and found around 60%/40%. 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. (19), choose t 0 = t − a and use the simplified notation v (n) ≡ (v The correlation matrix C defined in Eq. (10) is hermitean, i.e. C † = C, the eigenvalues λ (n) are real and we assume them to be positive and non-degenerate, i.e. λ (0) > λ (1) > . . . > λ (N −1) > 0. We can rewrite v (m) † C(t)v (n) in two ways by using Eq. (A1) and its hermitean conjugate, Since λ (m) = λ (n) for m = n, we conclude v (m) † C(t − a)v (n) = 0 for m = n. In the same way one can show v (m) † C(t)v (n) = 0 for m = 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 |Ψ from the N -dimensional space spanned by the states |Φ j included in the correlation matrix C, |Ψ = j a j |Φ j . (A4) The corresponding correlation function is Since a (as well as any other complex N -component vector) can be expanded in terms of the eigenvectors according to with µ (n) ∈ C, we can write for temporal separation t C Ψ (t) = Since 0 ≤ |µ (m) | 2 / n |µ (n) | 2 ≤ 1 and m (|µ (m) | 2 / n |µ (n) | 2 ) = 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 minimal eigenvalue λ (N −1) . Minimizing E eff µ (t) with respect to µ is equivalent to maximizing the argument of the logarithm, which corresponds to arbitrary µ (0) and µ (1) = µ (2) = . . . = µ (N −1) = 0. Thus, E eff µ (t) is minimized for a ∝ v (0) , which implies |Ψ ∝ j v (0) j |Φ j . This is, what we wanted to show: The coefficients of the trial state (A4) with minimal effective energy at temporal separation t are identical to the components of the eigenvector v (0) from the GEVP with t 0 = t − a.