Lattice QCD investigation of a doubly-bottom $\bar{b} \bar{b} u d$ tetraquark with quantum numbers $I(J^P) = 0(1^+)$

We use lattice QCD to investigate the spectrum of the $\bar{b} \bar{b} u d$ four-quark system with quantum numbers $I(J^P) = 0(1^+)$. We use five different gauge-link ensembles with $2+1$ flavors of domain-wall fermions, including one at the physical pion mass, and treat the heavy $\bar{b}$ quark within the framework of lattice nonrelativistic QCD. Our work improves upon previous similar computations by considering in addition to local four-quark interpolators also nonlocal two-meson interpolators and by performing a L\"uscher analysis to extrapolate our results to infinite volume. We obtain a binding energy of $(-128 \pm 24 \pm 10) \, \textrm{MeV}$, corresponding to the mass $(10476 \pm 24 \pm 10) \, \textrm{MeV}$, which confirms the existence of a $\bar{b} \bar{b} u d$ tetraquark that is stable with respect to the strong and electromagnetic interactions.

exotic channels. In this way we expand on the works of Refs. [25][26][27], where nonlocal interpolating fields were not considered. Having both local and nonlocal interpolating fields allows us to determine the ground-state and firstexcited-state energy in the I(J P ) = 0(1 + ) channel and perform a Lüscher analysis of BB * scattering.
The paper is structured as follows. In Sec. II we summarize our lattice setup, including the computation of quark propagators. In Sec. III we discuss the interpolating operators and the corresponding correlation functions. The extraction of the energy levels on the lattice is discussed in Secs. IV and V. In Sec. VI we present the scattering analysis, and in Sec. VII we perform a fit of the pion-mass dependence of the binding energy and estimate systematic uncertainties. Our conclusions are given in Sec. VIII.

A. Gauge-link configurations and light-quark propagators
We performed the computations presented here using gauge-link configurations generated by the RBC and UKQCD collaborations [31,32] with 2+1 flavors of domain-wall fermions [33][34][35][36] and the Iwasaki gauge action [37]. We use the five ensembles listed in Table I, which differ in the lattice spacing a ≈ 0.083 fm . . . 0.114 fm, the lattice size (spatial extent ≈ 2.65 fm . . . 5.48 fm) and the pion mass m π ≈ 139 MeV . . . 431 MeV. Ensemble C00078 uses the Möbius domain-wall action [36] with length of the fifth dimension N 5 = 24, while the other ensembles use the Shamir action [35] with N 5 = 16. The lattice spacings listed in the Table were determined in Ref. [32]. : bare strange sea quark mass; mπ: pion mass. We use all-mode-averaging [38,39] with 32 or 64 sloppy (sl) and 1 or 2 exact (ex) measurements per configuration; the column titled "Nmeas" gives the total numbers of sloppy and exact light-quark propagators used on each ensemble. The values NEV and N CG, sl are the numbers of eigenvectors used for the deflation of the light-quark solver, and the conjugate-gradient counts used for the sloppy propagators.
Our calculation uses smeared point-to-all propagators for the up and down quarks (the smearing parameters are given in Sec. III A). The computational cost of generating these propagators was reduced using the all-mode-averaging technique [38,39]. On each configuration, a small number of samples of "exact" correlation functions is combined with a large number of samples of "sloppy" correlation functions in such a way that the expectation value is equal to the exact expectation value, but the variance is reduced significantly due to the large number of sloppy samples [38,39]. The exact correlation functions are generated from light-quark propagators computed with high precision (relative solver residual of 10 −8 ), while the sloppy correlation functions are generated from approximate light-quark propagators. We used the conjugate gradient (CG) solver combined with low-mode deflation, where in the case of the approximate propagators the CG iteration count is fixed to a smaller value, N CG, sl , than needed for the exact propagators. The lowest N EV eigenvectors of the domain-wall operator were computed using Lanczos with Chebyshev-polynomial acceleration. The values of N EV and N CG,sl are also listed in the table. On a given gauge-link configuration, the different samples were obtained by displacing the source locations on a four-dimensional grid, with a randomly chosen overall offset.

B. Bottom-quark propagators
The heavy b quarks are treated with the framework of lattice nonrelativistic QCD (NRQCD) [40,41]. We use the same lattice NRQCD action and parameters as in Ref. [42]. This action includes all quark-bilinear operators through order v 4 in the heavy-heavy power counting, and order Λ 2 /m 2 b in the heavy-light power counting. The bare heavy-quark mass was tuned on the C005 and F004 ensembles such that the spin-averaged bottomonium kinematic mass agrees with experiment, using the lattice spacing determinations of Ref. [43]. We use the same masses also on the other coarse and fine ensembles, respectively, as shown in Table II. The matching coefficients c 1 , c 2 , c 3 were set to their tree-level values (= 1), while for c 4 we use a result computed at one loop in lattice perturbation theory [44]. The gauge links entering the NRQCD action are divided by the mean link u 0L in Landau gauge to achieve tadpole improvement [45,46]. Since we focus on the computation of the energy spectrum in this work, it is sufficient to use the leading-order, tree-level relation between the full-QCD bottom-quark field b and the two-spinor NRQCD quark and antiquark fields ψ and χ when constructing the hadron interpolating fields. In the Dirac gamma matrix basis, and omitting the phase factor that produces the tree-level energy shift, this amounts to and the bottom-quark propagator becomes with the two-spinor NRQCD quark and antiquark propagators G ψ and G χ .

III. INTERPOLATING OPERATORS AND CORRELATION FUNCTIONS
A.bbud four-quark system We are interested in the spectrum of the doubly-bottomed system with quantum numbers I(J P ) = 0(1 + ). The lowest two thresholds in this channel correspond to the meson pairs BB * and B * B * . From Ref. [42] we can see that the Ξ bb N antibaryon-baryon threshold is already much higher than the B * B * threshold: 11.1 GeV for Ξ bb N compared to 10.6 GeV for B * B * . Other possible meson-meson or antibaryon-baryon thresholds with bottomness 2 are even higher.
The J P = 1 + quantum numbers appear in the T g 1 irreducible representation (irrep) of the O h point group [47]. To determine the low-lying spectrum in this irrep we make use of two types of interpolating operators: local and nonlocal. The first three operators are local operators, in which all four (smeared) quark fields are multiplied at the same space-time point and the product is projected to zero momentum (in the following we omit the time coordinate): Operators four and five are nonlocal operators, where each color singlet is projected to zero momentum individually: Above, a, b, c, . . . denote color indices, j, k, l spatial vector indices, and C = γ 0 γ 2 is the charge-conjugation matrix. We expect that operators O 1 to O 3 generate sizable overlap to a tetraquark state. O 1 is an obvious choice in studying this channel. Our reason to include O 2 is a bit more subtle. Since two B * mesons are around 45 MeV heavier than a B and a B * meson, one might expect that O [B * B * ](0) will have less overlap with the ground state than O [BB * ](0) . However, a previous investigation with static quarks (see Ref. [22] and Fig. 3 therein) has determined the wave functions for both BB * and B * B * contributions and found that both spin structures are of similar importance. The inclusion of the color-triplet diquark-antidiquark interpolator O 3 is motivated by the heavy-quark limit, in which the two heavy antiquarks are expected to form a compact color-triplet object, and the tetraquark becomes equivalent to a singly heavy baryon [9][10][11]. The importance of diquark operators was also discussed in Refs. [48,49].
Note that in exotic channels it is typically difficult to properly resolve the ground state due to the coupling to nearby thresholds. An example of this phenomenon is the positive-parity D s spectrum, where the authors of Refs. [29,30] had to include nonlocal interpolating operators to resolve the D s0 (2317) mass puzzle. Previous studies of the I(J P ) = 0(1 + ), bottomness-2 sector [25,26] did not include multi-hadron operators, which might have affected their results. To make our determination of the spectrum more robust, we include operators O 4 and O 5 , which are nonlocal meson-meson scattering operators built from two color-singlets separately projected to zero momentum. We expect them to generate sizable overlap with the nearby first excited state, which will help us isolate the ground state in the multi-exponential fits of the correlation matrices.
To improve the overlap to the low-lying states, we employ standard smearing techniques. The quark fields in O 1 to O 5 are Gaussian-smeared using where ∆ is the nearest-neighbor gauge-covariant spatial Laplacian. For the up and down quarks, the gauge links in ∆ are spatially APE-smeared [50], while the unsmeared gauge links are used for the bottom quark. The smearing parameters are collected in Table III. To determine the spectrum we compute the temporal correlation functions of the interpolating operators O 1 to O 5 , where . . . denotes the path integral expectation value. The corresponding nonperturbative quark-field Wick contractions in a given gauge-field configuration are shown in Fig. 1. Because the calculation of the light-quark propagators is computationally expensive, we reuse existing smeared point-to-all propagators. With such propagators we are limited to operators O 1 , O 2 , O 3 at the source, which have only a single momentum projection, allowing us to remove the summation over x at the source (using translational symmetry). Thus, we do not determine the elements C 44 (t), C 45 (t), C 54 (t) and C 55 (t) of the correlation matrix, where two momentum projections are needed both at the source and at the sink. For a detailed discussion of this approach, see, for example, Refs. [52,53].
The correlation matrix has analytical properties that follow from the symmetries of lattice QCD, in particular time reversal and charge conjugation. These symmetries imply that (C jk (t)) * = C kj (t) and that all C jk (t) are real. Moreover, one can relate C jk (+t) and C jk (−t). We exploit these analytical findings to improve our lattice QCD results, by averaging related correlation functions appropriately and by setting all imaginary parts (which are pure noise) to zero. Since we will compare the energy levels of thebbud four-quark system to the BB * threshold, we need to determine the energies of the B and B * mesons within the same setup. We use the interpolating operators where we also consider nonzero momenta p = 2πn/L, n ∈ Z 3 , to allow the determination of the kinetic masses (needed for the scattering analysis in Sec. VI). The quark fields are smeared with the same parameters as in Table  III. We compute the correlation functions O B (p, t)O † B (p, 0) and O B * (p, t)O † B * (p, 0) . As discussed in Sec. III A, we perform the summation over x at the sink only.

IV. ENERGIES AND KINETIC MASSES OF THE B AND B * MESONS
We determined the energies of the B and B * mesons from single-exponential fits of the two-point functions; the results are listed in Table IV. An example of a corresponding effective-energy plot is shown in Fig. 2. The energy of a state containing n b bottom quarks is shifted by n b times the NRQCD energy shift, which at tree-level would be equal to −m b . This energy shift is not known with high precision, but cancels in energy differences with matching numbers of bottom quarks, including the quantities of interest in the following sections: E n − E B − E B * , where E n is the n-th energy level of thebbud system.
For the scattering analysis in Sec. VI (and to assess the tuning of the b-quark mass), we also need the momentumdependence of the B and B * energies. We find that within statistical uncertainties our results are consistent with the form   up to the highest momenta we computed (p 2 = 3(2π/L) 2 ). The kinetic masses given in Table IV were extracted using with the smallest possible non-vanishing momentum, p 2 = (2π/L) 2 . On the C00078 and C005 ensembles, the spinaveraged kinetic masses m kin,spinav = 1 4 m B, kin + 3 4 m B * , kin agree with the experimental value of 5.28272(2) GeV [54]. On the other ensembles, the lattice results are up to 5% higher, which can be attributed mainly to the following: (i) The tuning of the bare b-quark mass was performed using bottomonium; the results here are affected by the heavier-than-physical light-quark masses and by discretization errors.
(ii) The tuning was performed using a different determination of the lattice spacing (that of Ref. [43], while here we use the lattice spacing determinations of Ref. [32] to convert to physical units).
However, the effect of a possible < ∼ 5% mistuning of the b-quark mass is expected to be even smaller in the energy differences E n − E B − E B * due to partial suppression by heavy-quark symmetry. The hyperfine splittings E B * − E B are of order Λ/m b , and are therefore affected by the same relative error as the b-quark mass, which corresponds to an absolute error of < ∼ 2.5 MeV. Our results for the hyperfine splittings are also shown in Table IV. The value from the physical-pion-mass ensemble C00078 is consistent with the experimental value of 45.3(2) MeV [54]. We find a trend of increasing hyperfine splitting as the light-quark mass is increased.

V. THE LOWEST ENERGY LEVELS OF THEbbud SYSTEM
A. Multi-exponential matrix fitting The spectral decomposition of the correlation matrix (9) reads where |Ω denotes the vacuum state, |n are the energy eigenstates of thebbud system in the T g 1 irrep and for isospin I = 0, and E n are the corresponding energy eigenvalues. Eq. (15) assumes an infinite time extent of the lattice, which is a good approximation in our case. As discussed in Sec. III A, all C jk (t) are real and, thus, so are the overlap factors To extract the energies E n from our numerical results for C jk (t), we perform fully correlated, least-χ 2 , multiexponential fits using a truncated version of (15), where we must choose the time range t min ≤ t ≤ t max such that contributions from higher excited states are negligible.
To enforce the ordering of the E n returned from the fit, for n > 0 we actually use the logarithms of the energy differences, l n = ln(aE n − aE n−1 ), as our fit parameters. Furthermore, we rewrote the overlap factors for n > 0 as Z n = B n Z 0 and used B n as the fit parameters. Note that C jk (t) does not need to be a square matrix. In particular, the multi-exponential matrix fit method allows us to also include correlation matrix elements C jk (t) with j = 4, 5 and k = 1, 2, 3, i.e. the scattering operators O 4 and O 5 . An example of such a (5 × 3) matrix fit is shown in Fig. 3.
On each ensemble, we performed fits to the full matrix as well as to various sub-matrices, for different fit ranges and different numbers of states. The results are listed in Tables VII-XI in Appendix A.
The matrix fits have a large number of degrees of freedom, and the covariance matrix, whose inverse enters in χ 2 , can become poorly determined or singular if the number of statistical samples used to estimate this matrix is not much larger than the number of degrees of freedom. This is particularly problematic when using all-mode-averaging (or, more generally, covariant-approximation-averaging, CAA), where the samples are given by [38,39] (CAA sample) e = (exact sample) e − (sloppy sample) e,0 + 1 N sloppy Here, e labels the exact samples (gauge-link configuration and source location), and the different sloppy samples, computed from quark propagators with reduced solver precision, in our case originate from applying many space-time displacements s to the initial source location e, with s = 0 corresponding to no displacement. As shown in Table I, the number of exact samples, and hence CAA samples, is as low as 80 in the case of the C00078 ensemble. To obtain a meaningful χ 2 even for large numbers of degrees of freedom, we use a "modified CAA" (MCAA) procedure, given by This procedure provides N sloppy (= 32 in our case) times as many samples as the standard CAA procedure, without changing the overall average, and allows robust matrix fits even in the (5 × 3) case. The drawback is that there are autocorrelations between the different choices of s, introduced by the constant (but very small) term (exact sample) e − (sloppy sample) e,0 and by any possible autocorrelations between the different sloppy samples on the same configuration. As a result of these autocorrelations, the uncertainties of parameters obtained from fits based on the MCAA procedure are initially slightly underestimated. We correct for these autocorrelations by rescaling all uncertainties in the fitted E n with a factor estimated for each ensemble using the simple B meson two-point functions. The factor is given by the ratio of uncertainties of the B meson energies obtained from single-exponential fits using the CAA and MCAA procedures. The factors range from 1.08 to 1.27. The uncertainties of all results shown in this paper are already corrected with these factors. As a cross-check of the multi-exponential matrix fitting method, we also determined the spectrum using the variational approach, which involves solving the generalized eigenvalue problem (GEVP) [55][56][57]. We found that the GEVP method, where applicable, gives results consistent with the direct multi-exponential matrix fits. The comparison of the two methods is presented in Appendix B.
Example of a multi-exponential matrix fit for the C005 ensemble. This fit has N = 3, tmin/a = 11, tmax/a = 24, and gives

B. Dependence of the fit results on the choices of interpolating fields
Even though the actual energy levels for a chosen set of quantum numbers are independent of the interpolating operators used in the two-point functions, in practice the numerical results depend on these choices, due to limited statistical precision. In Fig. 4 we present the two lowest energy levels relative to the BB * threshold as determined on ensemble C005 using multi-exponential matrix fitting. The interpolators used are indicated by the five bars below each column. The three black bars at the bottom correspond to the local interpolators O 1 , O 2 , O 3 , while the two red bars at the top correspond to the nonlocal interpolators O 4 and O 5 . A filled bar indicates that an interpolator is included, an empty bar that it is not included. Two things are evident from Fig. 4: first, stable results for the two lowest energy levels are only obtained once the nonlocal two-meson interpolators are included in the analysis; second, once the nonlocal interpolators are included, the estimates of both the ground-state energy and the first excited energy drop significantly. Both observations clearly indicate the importance of the nonlocal interpolators for our study of thebbud system. Note that the first excited energy is close to the BB * threshold, while the ground state energy is significantly below threshold. This is a first indication that the ground state corresponds to a stable tetraquark, while the first excitation corresponds to a meson-meson scattering state. Similar plots for the other four ensembles are collected in Appendix A.

C. Overlap factors
For a given j, the overlap factors Z n j indicate the relative importance of the energy eigenstates |n when expanding the trial state O † j |Ω in terms of energy eigenstates, Therefore, the overlap factors Z n j provide certain information about the composition and quark arrangement of the energy eigenstates |n . In particular, if the overlap factor Z m j for one state |m is significantly larger than all other Z n j , n = m, this might be a sign that |m is quite similar to O † j |Ω . It is convenient to consider rescaled squared overlap factors, which are normalized such that max m (|Z m j | 2 ) = 1 for each trial state O † j |Ω . Here, the indices n and m can take on values from 0 to N − 1, where N is the number of states included in the fit. The results for |Z n j | 2 obtained from our fits are qualitatively similar for all ensembles, and do not strongly depend on the temporal fit range t min ≤ t ≤ t max . In Fig. 5 we show the normalized overlap factors obtained on ensemble C005 using the full 5 × 3 correlation matrix for N = 2 and N = 3, with fit ranges of 11 ≤ t/a ≤ 24 and 10 ≤ t/a ≤ 24, respectively. In particular, for N = 3, one can clearly see that for the three trial states created by the local operators O 1 , O 2 and O 3 , the ground-state overlaps are significantly larger than the overlaps to the first excited state. Vice versa, for the two trial states created by the nonlocal operators O 4 and O 5 , the overlaps to the first excited state are much larger than the ground-state overlaps. This supports our above interpretation concerning the composition of the lowest two energy eigenstates: the ground state |0 seems to be a four-quark bound state and the first excitation |1 a meson-meson scattering state. Finally, it is interesting to note that in the three-exponential fit, O 1 and O 2 appear to produce a large overlap with the second excited state, while the diquark-antidiquark operator O 3 and the nonlocal operators O 4 and O 5 do not. What appears in the fit as the "second excited state" (with a rather high energy and large uncertainty, as shown in Table VIII) FIG. 5. The normalized overlap factors |Z n j | 2 as determined on ensemble C005, indicating the relative contributions of the energy eigenstates |n to the trial state O † j |Ω . The upper row corresponds to a two-exponential fit with 11 ≤ t/a ≤ 24, while the lower row corresponds to a three-exponential fit with 10 ≤ t/a ≤ 24.

D. Final results for the lowest two energy levels on each ensemble
To obtain final results for the energy differences ∆E 0 and ∆E 1 on each ensemble, we select a subset of fits that we deem reliable, based on the discussion in the previous two subsections and based on χ 2 /d.o.f.. All of these fits, which are indicated with filled symbols in Fig. 4 and Figs. 9-12, include at least one of the nonlocal interpolating operators O 4 and O 5 . We then repeat all selected fits for 500 bootstrap samples, where each bootstrap sample consists of randomly drawn gauge-link configurations and source locations (using the same random numbers for the different types of fits to preserve correlations). This produces 500 samples for ∆E 0 and ∆E 1 for each type of fit (and on each ensemble). We then average these values over the different types of fits with equal weights, bootstrap sample by bootstrap sample. Finally, we compute the mean and standard deviation (rescaled again by the MCAA autocorrelation correction factors; see the discussion in Sec. V A) of these new 500 samples. The results for ∆E 0 and ∆E 1 obtained in this way are listed in Table V and are indicated with the horizontal lines and shaded uncertainty bands in Fig. 4

VI. SCATTERING ANALYSIS
The energies E n determined in a spectroscopy computation can be related to the infinite-volume scattering amplitude using Lüscher's method [58] and its generalizations [59][60][61][62][63][64][65][66]; see Ref. [67] for a recent review. Here, we apply Lüscher's method to the lowest two energy levels of thebbud system in the T g 1 irrep, assuming that these energy levels can be described in terms of the elastic S-wave B-B * scattering amplitude (and its analytic continuation below threshold). Thus, we obtain the S-wave scattering amplitude for the two scattering momenta k 0 and k 1 that correspond to E 0 and E 1 . Having only these two points limits the choice of parametrization of the scattering amplitude to a function with two parameters, for which we use the effective-range expansion (ERE). The ERE parametrization then allows us to determine the infinite-volume bound-state energy.

A. Relation between finite-volume energy levels and infinite-volume phase shifts
We define the scattering momentum k n corresponding to the n-th energy level of thebbud system E n through the equation where E B = E B (0) and E B * = E B * (0) are the energies of the single-B meson and single-B * meson states at zero momentum. Solving for k 2 n gives k 2 n = ∆E n (∆E n + 2m B, kin )(∆E n + 2m B * , kin )(∆E n + 2m B, kin + 2m B * , kin ) 4(∆E n + m B, kin + m B * , kin ) 2 , where The mapping between the finite-volume scattering momentum in the rest frame and the infinite-volume scattering amplitude expressed as the S-wave scattering phase shift δ 0 is where Z 00 is the generalized zeta function [58]. The scattering amplitude is given by B. Effective-range expansion and determination of the bound-state pole On each ensemble, we use the lattice QCD results for the energy differences ∆E 0 and ∆E 1 combined with the corresponding results for the B and B * meson energies and kinetic masses to calculate k 2 n for n = 0 and n = 1 using Eq. (23), and we determine the corresponding k n cot δ 0 (k n ) using Eq. (25). We parametrize the scattering amplitude using the effective-range expansion (ERE), and determine the two parameters a 0 (the S-wave scattering length) and r 0 (the S-wave effective range). Bound states correspond to poles in the scattering amplitude (26) below threshold, where −ik > 0. Such a pole occurs for cot δ 0 (k BS ) = i, where k BS is the (imaginary) bound-state momentum. Combining this condition with the ERE then gives where terms of O(|k BS | 4 ) are neglected. We solve Eq. (28) for |k BS | and obtain the binding energy via This approach was previously utilized in Refs. [29,30,68,69]. Determining masses of bound states in this way is equivalent to taking the infinite-volume limit up to exponentially small finite-volume corrections proportional to e −|kBS|L (for more details, see, for example, Refs. [70,71]).

C. Numerical results
In Fig. 6 we show k cot δ 0 (k) = 1/a 0 + r 0 k 2 /2 as a function of the scattering momentum k together with the corresponding lattice data points, for the five different ensembles. One can see that k cot δ 0 (k) is consistent with zero within uncertainties at k = 0, which implies an inverse scattering length also consistent with zero. The results for the inverse scattering length 1/a 0 , the effective range r 0 , the binding momentum |k BS |, and the binding energy E binding are given in Table VI. The binding energies are essentially identical to the finite-volume energy differences ∆E 0 collected in Table V. This supports our interpretation of the ground state as a stable tetraquark. We also applied the consistency check proposed in Ref. [72], but our statistical uncertainties are too large to reach a definitive conclusion.

VII. FIT OF THE PION-MASS DEPENDENCE AND ESTIMATES OF SYSTEMATIC UNCERTAINTIES
Given the statistical uncertainties in our results for E binding (shown in Table VI), we cannot resolve any significant dependence on the lattice spacing or pion mass. We expect lattice discretization errors to be at the level of a few MeV, as discussed further below. Since this is well below our statistical uncertainties, we choose to perform a fit of the pion-mass dependence of our results from all ensembles without including a-dependence. We consider a quadratic pion-mass dependence, corresponding to linear dependence on the light-quark mass, where we use m π,phys = 135 MeV for the physical pion mass in the isospin-symmetric limit. The fit gives and has χ 2 /d.o.f. = 0.27. A plot of the fit function together with the data is shown in Fig. 7. Given that, (i), the fit has excellent quality, (ii), the resulting coefficient c is consistent with zero, and (iii), the result for E binding (m π,phys ) is nearly identical with the result from the C00078 ensemble with m π = 139(1) MeV, we conclude that any remaining systematic uncertainties associated with the extrapolation to m π,phys are negligible. The lattice discretization errors associated with our light-quark and gluon actions and are expected to be at the 1% level for the fine lattices, and at the 2% level for the coarse lattices [32]; multiplying by the QCD scale of Λ ∼ 300 MeV yields 3 MeV to 6 MeV. The NRQCD action introduces additional systematic uncertainties. For our choice of lattice discretization and matching coefficients, the most significant systematic errors in the energy of a heavy-light hadron are expected to be the following: • Four-quark operators, which arise at order α 2 s in the matching to full QCD, are not included in the action. The analysis of Ref. [73] suggests that their effects could be as large as 3 MeV.
• The matching coefficient c 4 of the operator − g 2m b σ · B was computed to one loop. Missing higher-order corrections to this coefficient introduce systematic errors of order • The matching coefficients of the operators of order (Λ/m b ) 2 were computed at tree-level only. The missing radiative corrections to these coefficients introduce systematic errors of order These estimates are appropriate for E B and E B * , which contribute to our calculation of the binding energy via For the tetraquark energy E BS , the power counting is more complicated due to the presence of two bottom quarks. Conservative estimates of the systematic errors can be obtained by replacing the scale Λ in the heavy-light power counting by the binding momentum |k BS | ∼ 800 MeV, which suggests systematic errors of order 10 MeV. It is likely that there is a partial cancellation of the systematic errors in Therefore, we estimate the overall discretization and heavy-quark systematic errors to be not larger than 10 MeV (in future work, the estimates of the heavy-quark errors could be made more precise by numerically investigating the dependence of E binding on the lattice NRQCD matching coefficients). Our final results for the tetraquark binding energy and mass are therefore E binding (m π,phys ) = (−128 ± 24 ± 10) MeV, m tetraquark (m π,phys ) = (10476 ± 24 ± 10) MeV, where m tetraquark is obtained by adding the experimental values of the B and B * masses [54] to E binding .

VIII. CONCLUSIONS
In this work we computed the low-lying spectrum in the bottomness-2 and I(J P ) = 0(1 + ) sector. Using both local and nonlocal interpolating operators, we determined the two lowest energy levels for five different ensembles of lattice gauge-link configurations, including one with approximately physical light-quark masses. We carried out a Lüscher analysis for the first time in this sector and used the effective-range expansion to find the infinite-volume binding energies (but we found these to be nearly identical to the finite-volume binding energies). Our calculation confirms the existence of abbud bound state that is stable under the strong and electromagnetic interactions.
In Fig. 8 we compare our result (34) for the binding energy with several previous determinations. These include direct lattice QCD calculations that also treated the heavyb quarks with NRQCD [25,26], calculations in the Born-Oppenheimer approximation using staticbb potentials (in the presence of u and d valence quarks) computed on the lattice [19,20,22], as well as the studies of Refs. [15] and [11], which are based on a quark model and heavy-quark symmetry relations.
The calculations using staticbb potentials [19,20,22] consistently give a binding energy that is about a factor of 2 smaller than our result (34), but this disagreement might be due the approximations used there, in particular the neglect of 1/m b corrections to the potentials.
The two previous direct lattice QCD calculations [25,26] employed only local four-quark interpolating operators. According to our observations, the lack of nonlocal operators can affect the reliability of the extracted ground-state energy, as the local operators are not well-suited to isolate the lowest BB * threshold state. While the result of Ref. [26] agrees with ours, Ref. [25] gives a significantly larger binding energy. Apart from the lack of nonlocal interpolating operators, another possible source of this discrepancy might be the use of ratios of correlation functions as input to the generalized eigenvalue problem in Ref. [25]. It is interesting to observe that the effective energies shown in Refs. [25,26] approach the ground state from below, corresponding to a decrease in the magnitude of the extracted binding energy as the time separation is increased. Both of these studies used wall sources for the quark fields, while we use Gaussian-smeared sources. These two types of sources can behave quite differently with regard to excited-state contamination [74]. Furthermore, the excited-state spectrum of B-B * and B * -B * scattering states above threshold is very dense, because the changes in the kinetic energy when increasing the back-to-back momenta are suppressed by the heavy-meson masses. For example, on a lattice with L = 6 fm, the energy difference between the threshold and next scattering state is only around 8 MeV. In the context of two-nucleon systems, it has been argued that the dense spectrum can lead to "fake plateaus" at short time separations in the effective energies from ratios of correlation functions [75]; see Ref. [76] for a critical discussion of this issue. Comparison of results for the binding energy of thebbud tetraquark with I(J P ) = 0(1 + ) (black: this work, using lattice NRQCD; blue: previous work using lattice NRQCD [25,26]; red: lattice QCD computations of staticbb potentials and solving the Schrödinger equation [19,20,22]; green: effective field theories and potential models [11,15]).
Even though our work has improved upon previous studies of thebbud system by including nonlocal mesonmeson scattering operators at the sink, our fits still require rather large time separations, leading to large statistical uncertainties. As demonstrated for the case of the H dibaryon in Ref. [53], the results can be vastly improved by including nonlocal operators at both source and sink, and by including additional back-to-back momenta to map out a larger region of the spectrum. This requires more advanced techniques [52,77] for constructing the correlation functions.        Appendix B: Comparison of multi-exponential matrix fitting and solving the GEVP In this section we compare the multi-exponential matrix fitting method, used in the main part of this work, to the variational method [55][56][57]. The latter is based on the generalized eigenvalue problem C jk (t)v n k (t, t 0 ) = λ n (t, t 0 )C jk (t 0 )v n k (t, t 0 ), where C jk (t) must be a hermitian square matrix. As discussed in Sec. III A, our use of point-to-all propagators did not allow us to compute the correlation matrix elements C jk (t) with j, k ∈ {4, 5}, which means that we can only include the operators O 1 to O 3 in the GEVP analysis (in contrast, the multi-exponential matrix fitting method does not require a square matrix and allows us to include O 4 and O 5 at the sink). For large time, the eigenvalues λ n (t, t 0 ) are expected to satisfy where E n is the energy of the nth state. We define the effective energy E eff,n (t) as E eff,n (t) = 1 a ln λ n (t, t 0 ) λ n (t + a, t 0 ) , which for large t should plateau at E n . An example for E eff,n (t) is shown in Fig. 13. We determined E n from constant fits to E eff,n (t) in a suitable range t min ≤ t ≤ t max such that χ 2 /d.o.f. < ∼ 1. Alternatively, one can fit exponential functions Ae −Ent to the eigenvalues λ n (t, t 0 ). We performed such fits as cross-checks and found consistent energies and statistical uncertainties. We compared the multi-exponential matrix-fit method and the GEVP method by determining the lowest two energy levels with both methods in the following way: • We used the same values for t min (larger t min leads to a stronger suppression of excited states; we performed a comparison for several values for t min ).
• We used similar values for t max (the results only weakly depend on t max ; we chose t max as the largest temporal separation where the signal is not lost in statistical noise).
• We checked that the energy levels obtained by solving the GEVP are independent of the parameter t 0 for t 0 /a = 1, 2, . . . , 6 (results shown in the following were obtained with t 0 /a = 3).
This comparison is shown in Table XII for the C005 ensemble. It is reassuring that the results obtained from multiexponential matrix fits and from the GEVP are in excellent agreement (when using the same operator bases). The two methods were also compared extensively in the study of ππ scattering in Ref. [78], were agreement was also found. TABLE XII. Comparison of the results for the two lowestbbud energy levels (relative to the BB * threshold, i.e., ∆En = En − EB − EB * ) from multi-exponential matrix fitting and from the GEVP, for the C005 ensemble.