Search for ¯ b ¯ bus and ¯ b ¯ cud tetraquark bound states using lattice QCD

We use lattice QCD to investigate the existence of strong-interaction-stable antiheavy-antiheavy-light-light tetraquarks. We study the ¯ b ¯ bus system with quantum numbers J P = 1 + as well as the ¯ b ¯ cud systems with quantum numbers I ( J P ) = 0(0 + ) and I ( J P ) = 0(1 + ). We carry out computations on ﬁve gauge-link ensembles with 2 + 1 ﬂavors of domain-wall fermions, including one at the physical pion mass. The bottom quarks are implemented using lattice nonrelativistic QCD, and the charm quarks using an anisotropic clover action. In addition to local diquark-antidiquark and local meson-meson interpolating operators, we include nonlocal meson-meson operators at the sink, which facilitates the reliable determination of the low-lying energy levels. We ﬁnd clear evidence for the existence of a strong-interaction-stable ¯ b ¯ bus tetraquark with binding energy ( − 86 ± 22 ± 10) MeV and mass (10609 ± 22 ± 10) MeV. For the ¯ b ¯ cud systems we do not ﬁnd any indication for the existence of bound states, but cannot rule out their existence either.


I. INTRODUCTION
Hadrons with integer spin, in particular those corresponding to low-lying states in the respective spectra, are typically ordinary mesons composed of a single valence quark and a single valence antiquark.They might, however, also contain two valence quarks and two valence antiquarks.Such so-called tetraquarks 1 were discovered only recently, primarily in the heavy-quark sector [1][2][3][4][5][6][7][8].Of particular importance is the recent discovery of an anticharm-anticharmlight-light tetraquark T cc by the LHCb collaboration with isospin I = 0 and mass slightly below the lowest two-meson threshold corresponding to DD * [9,10].Such antiheavy-antiheavy-light-light systems Q Qqq are manifestly flavorexotic and are simpler to investigate theoretically than their QQqq counterparts, because the lowest relevant decay threshold consists of a pair of heavy-light mesons, typically with similar mass, and not the significantly lighter scattering states containing a light meson and ordinary quarkonium (or even the annihilation products of the quarkonium).Moreover, strong-interaction-stable Q Qqq tetraquarks are expected to exist for sufficiently large heavy quark masses m Q [11][12][13].In this limit, the two heavy antiquarks form a color-triplet with size of order (α s m Q ) −1 and binding energy of order α 2 s m Q due to the attractive Coulomb potential at small Q Q separations.Q Qqq tetraquarks are then quite similar to heavy-light-light baryons Qqq, just like heavy-heavy-light baryons Q Qq are related to heavy-light mesons Qq [14][15][16][17].Thus, the question is whether the physical heavy quark mass m c or m b is sufficiently large for Q Qqq bound states to exist below the corresponding lowest Qq-Qq two-meson thresholds.
Note that Q Qqq tetraquarks with heavy b quarks have not yet been observed experimentally.However, possible search strategies are discussed in Refs.[57][58][59].As mentioned above, the closely related T cc tetraquark with quark content ccud was recently discovered by the LHCb collaboration [9,10].A first lattice-QCD study of this system at a heavier-than-physical pion mass can be found in Ref. [60].
In this work we focus on the bb us system with quantum numbers J P = 1 + and the bcud systems with quantum numbers I(J P ) = 0(0 + ) and I(J P ) = 0(1 + ).We employ the same lattice QCD setup as in our previous study of the bb ud tetraquark with quantum numbers I(J P ) = 0(1 + ) [41], i.e. we use NRQCD to discretize b quarks and domain-wall light quarks.The charm quarks, which were not part of our previous study, are implemented using an anisotropic clover action with three parameters tuned nonperturbatively to eliminate heavy-quark discretization errors.In the construction of the two-point correlation functions, we consider not only local interpolating operators (in which the four quarks are jointly projected to zero momentum, i.e.where each quark is centered around the same point in space), but also non-local interpolating operators (in which each of the two quark-antiquark pairs forming a color-singlet is projected to zero momentum individually).It has been shown in previous studies of other four-quark systems that including both types of interpolating operators is required to reliably determine ground state energies in exotic channels [41,61,62].In this way we expand on the works of Refs.[39,40,[48][49][50], where non-local interpolating operators were not considered.This article is organized in the following way.In Section II we briefly summarize our lattice setup.In Section III we discuss the interpolating operators for the three systems we investigate and the corresponding correlation functions.In Section IV we give the lattice results for the single heavy-light meson energies.Section V is the main section, where we present our numerical results for the antiheavy-antiheavy-light-light four-quark systems.We explore the importance of each of our interpolating operators, extract finite-volume energy levels for all ensembles, and formulate conclusions concerning the existence of antiheavy-antiheavy-light-light tetraquarks at the physical u and d quark mass and in infinite spatial volume.We summarize the main points of our work in Section VI and give a brief outlook.Note that results obtained at an early stage of this project were presented at recent conferences [63,64].

A. Gauge link configurations, light quark and bottom quark propagators
The computations presented in this work were carried out on five ensembles of gauge link configurations generated by the RBC and UKQCD collaborations [65,66] using the Iwasaki gauge action [67] and N f = 2 + 1 flavors of domain-wall fermions [68][69][70][71].The ensembles differ in the lattice spacing, the lattice size and the pion mass and are summarized in Table I.Further details can be found in our previous lattice QCD study of a bb ud tetraquark [41], where we used exactly the same ensembles.: sea-quark mass of flavor q; am (val) q : valence-quark mass of flavor q; mπ: pion mass.We use all-mode-averaging [72,73] with 32 or 64 sloppy (sl) and 1 or 2 exact (ex) samples per configuration, leading to the total numbers of samples given in the last column of the table.
We use point-to-all propagators with Gaussian-smeared sources (cf.Sec.III A 4).We employ the all-mode averaging technique [72,73] with 32 or 64 sloppy (sl) and 1 or 2 exact (ex) samples per configuration, where the sloppy correlationfunction samples differ from the exact samples in that they use light and strange propagators computed with a reduced solver iteration count.The light-quark propagators are identical to those used in Ref. [41].The valence strange-quark masses are close to the physical value [66].For the bottom quarks we use lattice NRQCD [74,75]; also here the setup is the same as in Ref. [41].

B. Charm quark propagators
For the charm quarks we use an anisotropic clover action, following the approach developed in Refs.[76][77][78][79][80][81][82], which allows the removal of discretization errors of order |ap|, (am) n , and |ap|(am) n for all non-negative integers n.Specifically, our action is of the same form as in Ref. [82], and we tuned the mass am c (denoted as am 0 in Ref. [82]), anisotropy parameter ζ, and clover coefficient c P nonperturbatively such that the D s meson rest mass, kinetic mass, and hyperfine splitting extracted from two-point functions on each ensemble match the experimental values [83].These observables calculated on each ensemble are found to agree with experiment within 0.4%, 1.0%, and 1.4% (or better) precision, respectively.The values of the action parameters are given in Table II.[82], where mc is denoted by m0.

III. INTERPOLATING OPERATORS AND CORRELATION FUNCTIONS
A. Four-quark systems The main goal of this work is to compute low-lying energy levels of antiheavy-antiheavy-light-light four-quark systems with quark content bb us and bcud and to explore whether the ground-state energies are below the lowest corresponding meson-meson thresholds.A ground-state energy sufficiently far below threshold (compared to the expected size of finite-volume effects) would indicate a four-quark bound state, i.e. the existence of a strong-interactionstable tetraquark.In the bb us case we consider exclusively the J P = 1 + channel, which is the only channel where one can expect sufficiently strong attractive forces to generate a bound state (see the symmetry arguments given in Section III.B of Ref. [37]).In the bcud case we focus on I = 0, again because of the related stronger attraction of the four quarks [13,35].There are two promising I = 0 channels, because the heavy antiquark pair bc can be either flavor symmetric or flavor antisymmetric.The symmetric I(J P ) = 0(1 + ) channel is conceptually similar to the J P = 1 + channel for bb us (and also for bb ud, as investigated in detail within the same setup in our previous work [41]), while the antisymmetric I(J P ) = 0(0 + ) channel is different.
To be able to resolve possibly existing four-quark bound states as well as meson-meson scattering states, we employ both local interpolating operators and non-local ("scattering") interpolating operators.Local operators are constructed from products of four quark fields at the same point in space, followed by projection of the product to total momentum zero.Scattering operators, on the other hand, resemble two heavy-light mesons with independent spatial locations and individual projection of each meson to momentum zero.Local interpolating operators can be categorized further into meson-meson and diquark-antidiquark operators.The local meson-meson operators (as well as the scattering operators) resemble pairs of mesons with overall quantum numbers identical to those of the fourquark system of interest.For each local meson-meson operator we also consider a corresponding scattering operator, which differs only in the momentum projection.The importance of diquark-antidiquark pairs was pointed out in Refs.[84][85][86].Following Jaffe's notation of "good" and "bad" diquarks [84], our diquark-antidiquark operators are designed in such a way that the light diquark (us or ud) is a "good" diquark.If possible, we choose for the heavy diquark also a "good" configuration (in the case of bcud with I(J P ) = 0(0 + )), otherwise we use a "bad" heavy diquark (for bb us and for bcud with I(J P ) = 0(1 + )).
As we demonstrated in our previous work [41], scattering operators play an important role in extracting lowlying energy levels, because they generate sizable overlaps to energy eigenstates close to two-meson thresholds.In particular, if a four-quark bound state exists, scattering operators can eliminate contamination in the fit result for the corresponding energy level caused by nearby scattering states.In contrast to the bb ud system discussed in Ref. [41], which has an SU(2) isospin symmetry, there is no such symmetry for the light us quarks for the bb us system.The consequence is that there are not only two, but three relevant meson-meson thresholds, which are rather close, within around 50 MeV.They correspond to BB * s , B * B s (which is around 5 MeV above BB * s ) and B * B * s (which is around 50 MeV above BB * s ).The corresponding four local interpolating operators (three meson-meson operators and one diquark-antidiquark operator) are and the three scattering operators are Above, a, b are color indices, j, k, l are spatial indices, and C = γ 0 γ 2 is the charge conjugation matrix.We note that the operators O 3 , O 4 and O 7 are antisymmetric in the light quark flavors.The operators O 1 and O 2 as well as the operators O 5 and O 6 can be linearly combined in such a way that there is one symmetric and one antisymmetric light flavor combination.
2. Interpolating operators for bcud with I(J P ) = 0(0 + ) The lowest meson-meson thresholds in this channel are BD and B * D * .Their energy difference is, however, sizable, approximately 190 MeV.Thus, we expect that resolving energy levels close to the B * D * threshold is not of central importance when studying this channel and exploring the possible existence of a four-quark bound state below the BD threshold.Consequently, we only consider a single meson-meson structure of BD type.The corresponding two local operators are and the only scattering operator is The quantum number I = 0 implies the antisymmetric light flavor combination ud − du (as in our previous study [41] of the bb ud system).The heavy quark flavors bc are also in an antisymmetric combination, allowing J = 0, which is not possible for heavy quark flavors bb .
3. Interpolating operators for bcud with I(J P ) = 0(1 + ) For total angular momentum J = 1 the lowest meson-meson thresholds are B * D, BD * and B * D * .We follow a similar strategy as in the previous subsection and do not consider a B * D * meson-meson structure.The other two thresholds are separated by approximately 100 MeV.Thus, we use the three local operators and the two scattering operators

Quark propagators and correlation functions
As in our previous work [41] we apply standard smearing techniques to improve the overlap generated by the interpolating operators to the low-lying energy eigenstates.All quark-fields in Eq. (1) to Eq. ( 15) are Gaussiansmeared, where ∆ is the nearest-neighbor gauge-covariant spatial Laplacian.For the Gaussian smearing of the up, down, and strange quarks we use APE-smeared spatial gauge links [87] 2 , while for the charm quarks we use stout-smeared spatial gauge links [89].The reason for using different types of smearing is that we reuse quark propagators computed previously for other projects.No link smearing is used in the bottom quarks.All smearing parameters are listed in Table III.For each of the three systems discussed in Section III A 1 to Section III A 3 we computed temporal correlation matrices ( . . .denotes the expectation value of the lattice QCD path integral, and j, k now label different operator structures), from which we determine the low-lying energy eigenvalues and obtain information about the quark composition of the corresponding eigenstates, as discussed in detail in Section V.
All computations are based on point-to-all propagators with sources smeared as discussed above.For the light quarks we used the same propagators as in our previous work [41], where further technical details are discussed.As a consequence, we are restricted to correlation functions with a local interpolating operators at the source, for which one can use translational invariance to replace the spatial sum by a simple multiplication with the spatial volume.At the sink, however, both local and non-local interpolating operators are used.Thus, our correlation matrices are non-square matrices of sizes 7 × 4, 3 × 2 and 5 × 3, respectively, for the systems discussed in Section III A 1 to Section III A 3. It is straightforward to show that all three correlation matrices are real-valued and that the square sub-matrices are symmetric.We verified that our numerical results are consistent with these properties and exploited them to increase statistical precision.Similarly, we used the time reversal symmetry to relate C jk (t) and C jk (−t), which reduces statistical uncertainties even further.

B. B, Bs and D mesons
In Section V we will compare the resulting ground state energies of the bb us and bcud four-quark systems discussed above to the respective lowest meson-meson thresholds.To this end, we also computed the energies of the pseudoscalar and vector B, B s , and D mesons using exactly the same setup.The corresponding interpolating operators are

IV. ENERGIES OF PSEUDOSCALAR AND VECTOR B, Bs AND D MESONS
We determined the ground-state energies of pseudoscalar and vector B, B s and D mesons via uncorrelated χ 2minimizing fits of constants to the corresponding effective-energy functions at sufficiently large temporal separations combined with a jackknife analysis.As usual, these effective energies are defined as where C(t) is a temporal correlation function of one of the interpolating operators (18) to (23).The results for all six mesons for each of the five ensembles are listed in Table IV.As a cross-check we also determined these meson energies by correlated exponential fitting as in our previous work [41] and found consistent results.To exemplify the quality of our numerical data, we show in Fig. 1 effective-energy plots for ensemble C005 together with the corresponding plateau fits.Note that the energies of the B, B * , B s and B * s mesons listed in Table IV do not correspond to the full meson masses, as e.g.measured in experiment.The reason is the use of NRQCD, resulting in negative energy shifts proportional to n b , the number of b quarks present in the corresponding states.At tree level, this shift amounts to −n b m b , where m b is the b-quark mass.Since we exclusively consider energy differences between four-quark states and meson-meson thresholds with the same n b , these energy shifts cancel and there is no need to determine them.

V. RESULTS ON ANTIHEAVY-ANTIHEAVY-LIGHT-LIGHT FOUR-QUARK SYSTEMS
The correlation matrix (17) with interpolating operators from Section III A 1, Section III A 2 or Section III A 3 can be written as a sum over the energy eigenstates |n of the respective flavor and J P sector, with real valued and |Ω denoting the vacuum.To extract the energy levels E n and overlap factors Z n j from the numerical lattice-QCD results for C jk (t), we carry out correlated χ 2 -minimizing multi-exponential fits of a truncated version of the right hand side of Eq. ( 25), in a suitably chosen range t min ≤ t ≤ t max .For further technical details concerning this multi-exponential fitting we refer to Section V A of our previous work [41].To check for and to exclude systematic errors as well as to minimize statistical errors, we also consider submatrices of the correlation matrices defined in Section III and vary the temporal fit range.
A. bb us with J P = 1 + 1. Reduction of the size of the correlation matrix from 7 × 4 to 6 × 3 In a preparatory step we replace the local interpolating operators (1) to ( 4) by linear combinations of these operators, The coefficients vn j were determined by solving generalized eigenvalue problems where C jk (t) is the lattice-QCD result for the 4 × 4 correlation matrix containing the local operators O 1 , O 2 , O 3 and O 4 .We normalized the eigenvector components such that j |v n j (t)| 2 = 1 and show them for ensemble C01 in Fig. 2, where one can see that the eigenvector components v n j (t) are fairly independent of t, in particular for larger values of t.Thus we defined the coefficients in Eq. ( 28) as vn j = v n k (t/a = 8), where t/a = 8 was selected because the v n k (t/a = 8) have rather small statistical uncertainties and are already consistent with the plateaus formed at larger values of t (for ensemble C01 the coefficients vn j are collected in Table V; for the other four ensembles they are quite similar).With this definition, operator O j , when applied to the vacuum, should create a trial state with large overlap to energy eigenstate |j − 1 .Thus, this new set of operators offers the possibility to discard some of them (e.g.O 4 or even O 4 and O 3 ) to keep the corresponding correlation matrix small, while retaining at the same time the overlap to the low-lying energy eigenstates of interest.This is beneficial for the precision of the numerical analyses discussed below.Since we are mainly interested in the energy level of the ground state, O 1 is of particular importance.In practice, it turned out that using in addition also O 2 and O 3 is favorable with respect to a precise determination of energy levels.O 4 , however, does not seem to be advantageous in our context and is therefore discarded.Altogether our analysis is based on the three local interpolating operators O 1 , O 2 , O 3 and the three non-local interpolating operators defined in Eqs. ( 5) to (7).Thus, in the following we will study a 6 × 3 correlation matrix and its submatrices.

Energy levels
To reliably determine the lowest energy levels, in particular that of the ground state, we carried out multi-exponential fits as discussed at the beginning of this section.We considered various submatrices, numbers of exponentials N , and fit ranges t min ≤ t ≤ t max .The corresponding results with correlated χ 2 /d.o.f.< 2 are summarized for ensemble C01 in Fig. 3, while those for the other ensembles are collected in Appendix A. The boxes at the bottom of Fig. 3 indicate, for each fit, which interpolating operators were included.A filled/empty box represents an operator that was included/excluded.From bottom to top, the boxes represent O 1 , O 2 , ..., O 6 .Local operators are colored in black, scattering operators in red.The fit results for E 0 and E 1 are shown as blue and green points with error bars, where the energy of the lowest threshold, E B + E B * s , is subtracted (this threshold is represented by the horizontal dashed line).Above the plot, further details are provided for each fit: the number of exponentials, the temporal fit range, and the resulting correlated χ 2 /d.o.f. .
2. Squared normalized eigenvector components |v n j | 2 as functions of t for ensemble C01, obtained by solving generalized eigenvalue problems as defined in Eq. ( 29).The corresponding 4 × 4 correlation matrix contains the four local interpolating operators (1) to (4).The dashed horizontal lines represent the squares of the coefficients vn j , where vn The first seven columns from the left represent fits in which only local interpolating operators were considered.Each of the three local operators seems to be associated with a specific energy, O 1 with ≈ 0 MeV, O 2 with ≈ 130 MeV and O 3 with ≈ 200 MeV.This is not surprising, given that these operators were constructed in such a way that the corresponding 3 × 3 correlation matrix is approximately diagonal in the region of t separations that enter the multi-exponential fits.Clearly, O 1 is of particular importance for a precise determination of the energy of the ground state.Thus, O 1 was included in all further fits, where in addition to local operators also scattering operators were used.
It is crucial to note that for all fits that include at least O 1 and one of the scattering operators O 4 to O 6 , the fit result for E 0 is around 100 MeV below the BB * s threshold.This is a first clear indication that the ground state in the bb us and J P = 1 + sector is a strong-interaction-stable tetraquark.One can also see that the fit result for E 1 is in many cases close to 0 MeV, which is consistent with the expectation that the first excitation is a meson-meson scattering state close to the BB * s threshold.We note that the results for the other four ensembles are comparable, i.e.E 0 is around 100 MeV below the BB * s threshold, and E 1 is around 0 MeV for several fits (see Appendix A).As one can see from Fig. 3, the results for E 0 from fits including at least O 1 and one of the scattering operators O 4 to O 6 (represented by the filled blue data points) agree within the statistical uncertainties.Thus, these fit results seem to be suited to estimate the ground state energy and its uncertainty.We computed such an estimate by a weighted average of these fit results, assuming 100% correlation, using a standard method also employed by the FLAG collaboration [90] (see Appendix B for a brief summary).The estimated ground-state energies are also plotted in the corresponding figures, e.g. for ensemble C01 in Fig. 3 (the blue horizontal line and the light blue error band).Concerning the energy of the first excitation, Fig. 3 suggests that it is somewhere around the BB * s threshold.We refrain from estimating this energy in a quantitative way by computing a weighted average of selected fit results for E 1 .The reason is that it is hard to decide whether E 1 obtained by a particular fit indeed corresponds to the energy of the first excitation.There are several states that could be close to the BB * s threshold, e.g. a BB * s or a B * B s scattering state.Additionally, there might also be a B * B * s scattering state in that energy region because of the finite spatial volume and the attractive interaction of the two mesons [37].The low-lying excitations could correspond to superpositions of these structures and are expected to have similar energies.Thus, a fit result for E 1 close to the BB * s threshold could, for example, reflect the energy of the first or the second excitation or a mix of both.In principle, one could try to disentangle these excitations by studying the resulting overlap factors Z n j for each fit in detail.Since we only need the ground-state energy for our final analysis in Section V D, we discuss the overlap factors just for a single fit with N = 3 exponentials to the full 6 × 3 correlation matrix (see the following subsection).

Overlap factors
A trial state O j † |Ω can be expanded according to which shows that the overlap factors Z n j contain information about the composition and quark arrangement of the energy eigenstates |n .For example, an overlap factor |Z In Fig. 4 we show normalized overlap factors obtained via a multi-exponential fit with N = 3 in the range 16 ≤ t/a ≤ 24 to the full 6 × 3 correlation matrix of ensemble F004.Corresponding results for the other ensembles are qualitatively identical.We start with an extensive discussion of the overlap factors Z 0 j associated with the ground state |0 and then briefly comment on the overlap factors Z n j with n > 0 related to the excitations.operator O 2 operator O 6 FIG. 4. Normalized overlap factors Zn j for the bb us system obtained via a multi-exponential fit with N = 3 in the range 16 ≤ t/a ≤ 24 to the full 6 × 3 correlation matrix of ensemble F004.The index of the operator above each plot is identical to the index j, while the labels of the energy eigenstates below each plot correspond to the index n.O 3 ) of the same order of magnitude.Such a meson-meson composition is expected from existing static-light lattice QCD results [37] on the strong-interaction-stable bb ud tetraquark with I(J P ) = 0(1 + ), a closely related four-quark system (same quantum numbers J P , and because of isospin I = 0 antisymmetric in the light flavors), where it was found that it is a roughly even mixture of BB * and B * B * .The bb us system also has a sizable diquark-antidiquark component (operator O 4 ), albeit somewhat smaller than the aforementioned meson-meson components.This, too, is expected and is consistent with recent static-light lattice-QCD results on the bb ud tetraquark, where the meson-meson to diquark-antidiquark ratio was estimated to be around 60%/40% [86].

The result |Z
The ) that is symmetric in the light flavors us, i.e. the analog of an I = 1 operator for light flavors ud.This confirms that the bb us ground state is antisymmetric in the light flavors and indicates that it is the counterpart of the bb ud tetraquark with I(J P ) = 0(1 + ).While the operator O 3 is flavor antisymmetric, it was constructed via the GEVP in a way to generate almost no overlap with the ground state and with the lowest flavor-symmetric excitation |2 .Thus it is not surprising that Z0  [37]) and the finite spatial volume could lead to an energy level of the first excitation close to the BB * s threshold, as indicated by Fig. 3.
Finally, the overlap factors Z 2 j represent almost exclusively symmetric light flavor combinations.This indicates that also for the scattering states in our finite spatial lattice volume, SU(3) flavor symmetry is approximately preserved.Thus, the second excitation seems to be the analog of the ground state in the bb ud four-quark sector with I = 1, where no strong-interaction-stable four-quark state was found in a static-light lattice-QCD study [35,37].
B. bcud with I(J P ) = 0(0 + ) As discussed in Section III A 2, we consider three interpolating operators for this system: two local operators and one scattering operator.Thus, the corresponding correlation matrix has size 3 × 2. Since this is a rather small matrix, there is no need to further reduce the number of operators in a preparatory step, as done for the bb us system.
To determine the energy of the ground state we proceed as in Section III A 2 and carry out multi-exponential fits.Again we consider various submatrices, numbers of exponentials N , and fit ranges t min ≤ t ≤ t max .The corresponding results with correlated χ 2 /d.o.f.< 2 are summarized for ensemble C01 in Fig. 5, while those for the other ensembles are collected in Appendix A. Like for the bb us system, we again find significantly lower values for E 0 once the scattering operator O 3 (see Eq. ( 10)) is included, compared to fits in which only local interpolating operators are used.Averaging over the fits that include the scattering operator leads to an estimate for the ground-state energy, which is slightly above, but within its uncertainty compatible with, the BD threshold.We find similar results for the other four ensembles (see Appendix A).This suggest that there is no strong-interaction-stable four-quark state in this channel.The lowest energy eigenstate rather seems to be a BD scattering state.
In Fig. 6 we show the normalized overlap factors Zn j obtained via a multi-exponential fit with N = 3 in the range 6 ≤ t/a ≤ 10 to the full 3 × 2 correlation matrix of ensemble F004.Corresponding results for the other ensembles are qualitatively identical.It is obvious that the BD scattering trial state O 3 † |Ω has large overlap to the ground state and almost negligible overlap to the first and second excitation.This supports our above conclusion that the ground state is a meson-meson scattering state.operator O 3 FIG.6. Normalized overlap factors Zn j for the bcud system with I(J P ) = 0(0 + ) obtained via a multi-exponential fit with N = 3 in the range 6 ≤ t/a ≤ 10 to the full 3 × 2 correlation matrix of ensemble F004.The index of the operator above each plot is identical to the index j, while the labels of the energy eigenstates below each plot correspond to the index n.
C. bcud with I(J P ) = 0(1 + ) According to Section III A 3 we consider five interpolating operators here: three local operators and two scattering operators.Thus, the corresponding correlation matrix has size 5 × 3. We do not reduce the number of operators in a preparatory step as done for the bb us system.
To determine the energies of the ground state and of the first excitation, we again carry out multi-exponential fits and consider various submatrices, numbers of exponentials N , and fit ranges t min ≤ t ≤ t max .The corresponding results with correlated χ 2 /d.o.f.< 2 are summarized for ensemble C01 in Fig. 3, while those for the other ensembles are collected in Appendix A. As for the previously investigated four-quark systems, we find significantly lower values for E 0 and E 1 as soon as the scattering operators O 4 and O 5 (see Eq. ( 14) and Eq. ( 15)) are included.In particular, operator O 4 , which has a B * D-like meson-meson structure, favors small values for E 0 close to the B * D threshold.Since the local operator O 1 is also of B * D type, we estimate the ground state energy by averaging over the fits that include both O 1 and O 4 .The result is slightly above, but within its uncertainty compatible with, the B * D threshold.As before, we do not estimate the energy of the first excitation quantitatively by computing a weighted average of selected fit results for E 1 .We note, however, that this energy level seems to be close to the BD * threshold, which is around 100 MeV above the B * D threshold.We found similar results for the other four ensembles (see Appendix A).In summary, this suggests that there is no strong-interaction-stable four-quark state in this channel.The lowest energy eigenstate rather seems to be a B * D scattering state.
In Fig. 8  operator O 5 FIG. 8. Normalized overlap factors Zn j for the bcud system with I(J P ) = 0(1 + ) obtained via a multi-exponential fit with N = 3 in the range 14 ≤ t/a ≤ 20 to the full 5 × 3 correlation matrix of ensemble F004.The index of the operator above each plot is identical to the index j, while the labels of the energy eigenstates below each plot correspond to the index n.

D. Final results for the bb us and bcud ground-state energies
We list the final results for the ground-state energies relative to the lowest meson-meson thresholds for the three investigated four-quark systems and for all five ensembles in Table VI.These energies correspond to the horizontal blue lines and light blue error bands in Fig. 3, Fig. 5, Fig. 7, and Fig. 11 to Fig. 22.In Fig. 9, we plot these results as a function of m 2 π .
1. bb us with J P = 1 + For the bb us system we found ground-state energies around 70 MeV to 100 MeV below the BB * s threshold.These are the energies in a finite periodic spatial volume of linear extent N s a ≈ 2.7 fm for ensembles C005, C01, F004 and F006 and N s a ≈ 5.3 fm for ensemble C00078.To extrapolate to infinite volume, we could, in principle, proceed as in our previous work [41] on the bb ud tetraquark with I(J P ) = 0(1 + ) and use Lüscher's finite volume method [91,92].For the bb us system this is, however, technically more complicated, because one has to take into account at least two scattering channels, BB * s and B * B s , which have almost the same threshold energy.Moreover, the energy levels of the corresponding excitations are difficult to determine, as discussed in Section V A 2. However, since the finite-volume ground-state energies are significantly below these thresholds, we expect only mild finite-volume corrections, much smaller than our current statistical errors.This expectation is supported by our infinite-volume extrapolations of bb ud results in Ref. [41], where the finite-volume ground-state energies turned out to be essentially identical to their infinite-volume counterparts.Thus, we do not carry out an infinite-volume extrapolation in this work, but postpone such an analysis until we have improved lattice data available, in particular correlation functions with scattering operators at both the source and the sink.
Our five ensembles differ in the light-quark mass, corresponding to pion masses in the range 139 MeV < ∼ m π < ∼ 431 MeV, which allows us to perform an extrapolation of ∆E s to the physical point (note that one of our ensembles, C00078, has a light quark mass that is almost physical).Since the observed dependence on the light-quark mass is mild (in fact, consistent with no dependence), a fit that is linear in m u;d and hence quadratic in m 2 π is sufficient.We performed a χ 2 -minimizing fit using the ansatz where ∆E 0 (m π,phys ) and c are fit parameters and m π,phys = 135 MeV.The resulting values for these parameters are with χ 2 /d.o.f.= 0.81, indicating consistency of the lattice data with our linear ansatz.The data points and the fit are shown in the upper plot of Fig. 9.
There are also systematic errors due to the finite lattice spacing and the NRQCD action.We expect these errors to be of the same order as for the related bb ud system with I(J P ) = 0(1 + ).We have discussed these errors in detail in Section VII our previous work [41] and estimated them to be not larger than 10 MeV.Thus, our final results for the bb us tetraquark binding energy and mass are ∆E 0 (m π,phys ) = (−86 ± 22 ± 10) MeV , mbb us tetraquark (m π,phys ) = (10609 ± 22 ± 10) MeV, (35) where mbb us tetraquark is obtained by adding the experimental results of the B and B * s masses [83] to ∆E 0 .
2. bcud with I(J P ) = 0(0 + ) and I(J P ) = 0(1 + ) For both bcud systems, the finite-volume ground-state energies are compatible with the corresponding lowest mesonmeson thresholds.Thus, there is no indication that strong-interaction-stable tetraquarks exist in these channels.However, because of the statistical uncertainties of order 20 MeV . . .50 MeV (see Table VI), we cannot exclude the existence of a shallow bound state with binding energy of only a few MeV below the respective threshold.
Since we are not in a position to quantify finite-volume corrections, which might be sizable in particular for states close to the threshold, we also refrain from extrapolating our lattice results to physical pion mass.To summarize our finite-volume results in a graphical way, we nevertheless plot them in Fig. 9 in the same style as their bb us counterparts together with the relevant meson-meson thresholds.

VI. CONCLUSIONS AND OUTLOOK
We investigated a bb us and two bcud four-quark systems using lattice QCD with dynamical domain-wall u, d, and s quarks.The charm quarks were implemented using an anisotropic clover action with parameters tuned to remove heavy-quark discretization errors, while the b quarks were discretized within the framework of NRQCD.Our work improves upon existing similar studies [39,40,[48][49][50] by including also non-local (scattering) interpolating operators.In the bb us sector with quantum numbers J P = 1 + , we find clear evidence for a strong-interaction-stable tetraquark.The binding energy with respect to the BB * s threshold is (−86±22±10) MeV, which is consistent with previous lattice QCD results from Refs.[39,40].In Fig. 10 we summarize and compare these lattice QCD results with results obtained using different approaches, e.g.quark models, phenomenological considerations, or sum rules [13,18,24,26,28,30,[43][44][45][46][47].As discussed in the introduction, there are strong discrepancies, even on a qualitative level, between these nonlattice results.Thus, it is important to have multiple independent first-principles lattice-QCD computations, and the agreement of the lattice results from different groups, as shown with the blue and black data points in Fig. 10, increase the confidence in these results.
For the bcud systems with quantum numbers I(J P ) = 0(0 + ) and I(J P ) = 0(1 + ) the situation is less clear.We find finite-volume ground-state energies that are compatible with the lowest thresholds corresponding to BD and B * D, respectively.To decide whether there is a shallow bound state, more precise data and infinite-volume extrapolations will be needed.Results from previous lattice QCD studies [48][49][50] are mostly consistent with our results, but are also inconclusive.It is interesting to note that Ref. [50] reports a ground-state energy for I(J P ) = 0(1 + ) below the B * D threshold for a fine lattice spacing a ≈ 0.06 fm, but not for the coarse lattice spacing a ≈ 0.12 fm.The authors of Ref. [50] conclude that taking the continuum limit might be essential for the bcud system.We do not observe such a trend (see Fig. 9), but it should be kept in mind that the types of lattice actions used here differ from Ref. [50], except for the bottom quarks.As discussed in the introduction, also non-lattice studies do not clarify the possible existence of a strong-interaction-stable bcud tetraquark, since they exhibit strong discrepancies (Refs.[26,27,43,[51][52][53][54][55][56] predict the existence of a stable tetraquark, while Refs.[13,24,30,44,46] claim the opposite).
correlation matrices (rather than just the sinks as done here).We expect that this will allow us to determine the lowlying energy levels, in particular those associated with scattering states, more reliably and more precisely.We could then carry out infinite-volume extrapolations for the bcud systems using Lüscher's method [91] and possibly clarify the existence or non-existence of a strong-interaction-stable bcud tetraquark.Another interesting direction could be to explore heavy-heavy-light-light four-quark systems with other quantum numbers for which stable tetraquarks are not expected, but for which resonances could exist.A clear candidate is the bb ud system with I(J P ) = 0(1 − ), where such a resonance around 15 MeV above the BB threshold was predicted using static-static-light-light potentials computed with lattice QCD and the Born-Oppenheimer approximation [38].
Appendix A: Summary plots of multi-exponential fits to determine energy levels for ensembles C00078, C005, F004 and F006 In this appendix we show the results of multi-exponential fits to determine E 0 and E 1 for the ensembles C00078, C005, F004, and F006: • bb us with J P = 1 + : Fig. 11 to Fig. 14.
The style of these figures is identical to Fig. 3, Fig. 5 and Fig. 7, respectively, where the same quantities are shown for ensemble C01, and which are discussed in detail in Section V.        FIG. 17. Fit results for E0 for the bcud system with I(J P ) = 0(0 + ) relative to the BD threshold for ensemble F004.

1 .
Interpolating operators for bb us with J P = 1 +

FIG. 1 .
FIG. 1. Effective energies for pseudoscalar and vector B, Bs and D mesons for ensemble C005.The horizontal lines represent the corresponding plateau fits in the range t/a = 7 . . .20.

66 FIG. 3 .
FIG.3.Fit results for E0 (blue) and E1 (green) for the bb us system relative to the BB * s threshold for ensemble C01.

FIG. 5 .
FIG.5.Fit results for E0 for the bcud system with I(J P ) = 0(0 + ) relative to the BD threshold for ensemble C01.
we show the normalized overlap factors Zn j obtained via a multi-exponential fit with N = 3 in the range 14 ≤ t/a ≤ 20 to the full 5 × 3 correlation matrix of ensemble F004.Corresponding results for the other ensembles are qualitatively identical.One can see that the B * D scattering trial state O 4 † |Ω almost exclusively overlaps with the ground state, i.e. .e. the BD * scattering trial state O 5 † |Ω almost exclusively overlaps with the first excitation.This supports our interpretation of the ground state and the first excitation as B * D and BD * scattering states.
FIG.9.Ground-state energy as function of the squared pion mass for the bb us system (top), the bcud system with I(J P ) = 0(0 + ) (bottom left) and the bcud system with I(J P ) = 0(1 + ) (bottom right).For bb us we also show the fit and linear extrapolation to the physical point at m π,phys = 135 MeV [see Eq.(33) and Eq.(34)].Horizontal dashed lines indicate the lowest corresponding thresholds: the BB * s threshold for bb us, the BD threshold for bcud with I(J P ) = 0(0 + ), and the B * D threshold for bcud with I(J P ) = 0(1 + ).

10 FIG. 11 .
FIG. 11.Fit results for E0 (blue) and E1 for the bb us system relative to the BB * s threshold for ensemble C00078.

85 FIG. 13 .
FIG.13.Fit results for E0 (blue) and E1 (green) for the bb us system relative to the BB * s threshold for ensemble F004.

FIG. 16 .
FIG.16.Fit results for E0 for the bcud system with I(J P ) = 0(0 + ) relative to the BD threshold for ensemble C005.

TABLE I .
[65,66]ink ensembles[65,66]and light quark propagators used in this work.Ns, Nt: number of lattice sites in spatial and temporal directions; a: lattice spacing; am

TABLE II .
Parameters used in the anisotropic clover action for the charm quarks.The form of the heavy-quark action is given in Ref.
2 1 | suggests that the trial state O 1 † |Ω has a large ground-state overlap, i.e., is rather similar to the ground state.Recall that O 1 is a weighted sum of four local operators O 1 to O 4 (Eq.(28) with coefficients v0 j as listed in Table V for ensemble C01).Since v0 1 ≈ −v 0 2 , there is a local BB * s and B * B s component (operators O 1 and O 2 ) that is antisymmetric in the light flavors us.There is also a local antisymmetric B * B * s component (operator overlap factors Z n 2 and Z n 3 clearly show that the trial states O 2 † |Ω and O 3 † |Ω are essentially orthogonal to the ground state |0 .According to Table V, the operator O 2 is a local combination of BB * s and B * B s (operators O 1 and O 2 |Ω both have overlaps to the ground state |0 , but also sizable overlaps to the first and second excitations.Thus, one should not infer that the ground state is quite similar to a scattering state.Since the scattering operators O 4 and O 5 contain all terms present in the local operators O 1 and O 2 , the non-vanishing overlaps Z 0 4 and Z 0 5 rather support our conclusions above, namely that the bb us ground state is a four-quark bound state with a large local flavor antisymmetric BB * s and B * B s component.As already discussed above in the context of energy levels and the fit parameter E 1 , one should be cautious in formulating conclusions concerning the excitations based on our multi-exponential fit results.Still, it seems noteworthy to mention that the trial states O 3 † |Ω and O 6 † |Ω have large overlap with the first excitation |1 and only little overlap with |0 and |2 .Since O 6 is a B * B * s scattering operator and the dominant component of O 3 is a local B * B * s structure (see Table V), this might be a hint that the first excitation is of B * B * s type or at least contains a significant B * B * s component.Even though the B * B * s threshold is around 50 MeV above the BB * s threshold, the expected attraction of a B * meson and a B * s meson (see Ref.