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

TL;DR

This study uses lattice QCD with multiple ensembles and advanced analysis to confirm a stable doubly-bottom tetraquark state with specific quantum numbers, providing detailed mass and binding energy estimates.

Contribution

It introduces a comprehensive lattice QCD approach including nonlocal interpolators and Lüscher analysis to study the doubly-bottom tetraquark, improving upon previous methods.

Findings

Identifies a stable $ar{b} ar{b} u d$ tetraquark with a binding energy of about -128 MeV.

Provides a mass estimate of approximately 10476 MeV for the tetraquark.

Confirms the tetraquark's stability against strong and electromagnetic decays.

Abstract

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.

Tables

Table 1. Table 1: Gauge-link ensembles Aoki:2010dy ; Blum:2014tka and light-quark propagators used in this work. N s subscript 𝑁 𝑠 N_{s} , N t subscript 𝑁 𝑡 N_{t} : number of lattice sites in spatial and temporal directions; a 𝑎 a : lattice spacing; a ​ m u ; d 𝑎 subscript 𝑚 𝑢 𝑑 am_{u;d} : bare up and down quark mass; a ​ m s 𝑎 subscript 𝑚 𝑠 am_{s} : bare strange quark mass; m π subscript 𝑚 𝜋 m_{\pi} : pion mass. We use all-mode-averaging Blum:2012uh ; Shintani:2014vja with 32 or 64 sloppy (sl) and 1 or 2 exact (ex) measurements per configuration; the column titled “ N meas subscript 𝑁 meas N_{\textrm{meas}} ” gives the total numbers of sloppy and exact light-quark propagators used on each ensemble. The values N EV subscript 𝑁 EV N_{\textrm{EV}} and N CG, sl subscript 𝑁 CG, sl N_{\textrm{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.

Ensemble	$N_s^3\times N_t$	$a$ [fm]	$a{m}_{u;d}$	$a{m}_s$	$m_{\pi }$ [MeV]	$N_{\text{meas}}$	$N_{\text{EV}}$	$N_{\text{CG,sl}}$
C00078	${48}^3\times 96$	$0.1141\left(3\right)$	$0.00078$	$0.0362$	$139\left(1\right)$	$2560$ sl, $80$ ex	500	400
C005	${24}^3\times 64$	$0.1106\left(3\right)$	$0.005$	$0.04$	$340\left(1\right)$	$9952$ sl, $311$ ex	400	100
C01	${24}^3\times 64$	$0.1106\left(3\right)$	$0.01$	$0.04$	$431\left(1\right)$	$9056$ sl, $283$ ex	400	100
F004	${32}^3\times 64$	$0.0828\left(3\right)$	$0.004$	$0.03$	$303\left(1\right)$	$8032$ sl, $251$ ex	400	120
F006	${32}^3\times 64$	$0.0828\left(3\right)$	$0.006$	$0.03$	$360\left(1\right)$	$14144$ sl, $442$ ex	400	120

Table 2. Table 2: Parameters used in the NRQCD action for the bottom quarks.

Ensemble		$a{m}_b$		$u_{0L}$		$c_4$
C00078		2.52		0.8432		1.09389
C005, C01		2.52		0.8439		1.09389
F004, F006		1.85		0.8609		1.07887

Table 3. Table 3: Parameters for the smearing of the quark fields in the interpolating operators. The Gauss smearing is defined in Eq. ( 8 ). A single sweep of APE smearing Albanese:1987ds with parameter α APE subscript 𝛼 APE \alpha_{\textrm{APE}} is defined as in Eq. (8) of Ref.  Bonnet:2000dc , and we apply N APE subscript 𝑁 APE N_{\textrm{APE}} such sweeps.

Ensemble	Up and down quarks					Bottom quarks
	$N_{\text{Gauss}}$	${\sigma }_{\text{Gauss}}$	$N_{\text{APE}}$	${\alpha }_{\text{APE}}$		$N_{\text{Gauss}}$	${\sigma }_{\text{Gauss}}$
C00078	$100$	$7.171$	$25$	$2.5$		$10$	$2.0$
C005, C01	$30$	$4.350$	$25$	$2.5$		$10$	$2.0$
F004, F006	$60$	$5.728$	$25$	$2.5$		$10$	$2.0$

Table 4. Table 4: Energies of the B 𝐵 B and B ∗ superscript 𝐵 B^{*} mesons at rest and kinetic masses computed using one quantum of momentum. Also given are the hyperfine splittings and spin-averaged kinetic masses ( m kin , spinav = 1 4 ​ m B , kin + 3 4 ​ m B ∗ , kin subscript 𝑚 kin spinav 1 4 subscript 𝑚 𝐵 kin 3 4 subscript 𝑚 superscript 𝐵 kin m_{\rm kin,spinav}=\frac{1}{4}m_{B,\,{\rm kin}}+\frac{3}{4}m_{B^{*}\!,\,{\rm kin}} ) in physical units. All uncertainties are statistical only.

Ensemble	$a{E}_B\left(0\right)$	$a{E}_{B^{\ast }}\left(0\right)$	$a{m}_{B,\mathrm{kin}}$	$a{m}_{B^{\ast },\mathrm{kin}}$	$E_{B^{\ast }}\left(0\right)-E_B\left(0\right)$ [MeV]	$m_{\mathrm{kin},\mathrm{spinav}}$ [GeV]
C00078	0.4582(45)0	0.4820(46)	3.03(14)0	3.10(13)0	41.2(5.1)	5.33(22)0
C005	0.4647(14)0	0.4944(16)	3.002(40)	2.993(42)	53.2(1.7)	5.346(73)
C01	0.4742(12)0	0.5061(16)	3.034(38)	3.030(40)	57.0(1.8)	5.409(69)
F004	0.3750(11)0	0.3975(12)	2.323(21)	2.323(25)	51.3(1.8)	5.536(57)
F006	0.37655(87)	0.3985(10)	2.320(20)	2.311(23)	52.3(1.5)	5.513(54)

Table 5. Table 5: Results for the B 𝐵 B and B ∗ superscript 𝐵 B^{*} meson “speed of light” squared, defined in Eq. ( 15 ), for 𝐩 2 = 2 ​ ( 2 ​ π / L ) 2 superscript 𝐩 2 2 superscript 2 𝜋 𝐿 2 \mathbf{p}^{2}=2(2\pi/L)^{2} and 𝐩 2 = 3 ​ ( 2 ​ π / L ) 2 superscript 𝐩 2 3 superscript 2 𝜋 𝐿 2 \mathbf{p}^{2}=3(2\pi/L)^{2} .

Ensemble	$c_B^2\left(2\right)$	$c_B^2\left(3\right)$	$c_{B^{\ast }}^2\left(2\right)$	$c_{B^{\ast }}^2\left(3\right)$
C00078	0.9999(13)0	0.9998(26)	1.0007(12)0	1.0014(24)
C005	1.0019(14)0	1.0036(28)	1.0014(13)0	1.0026(28)
C01	1.0015(12)0	1.0032(25)	1.0000(15)0	1.0002(30)
F004	1.00005(78)	1.0001(16)	1.00063(85)	1.0014(17)
F006	1.00032(64)	1.0007(13)	1.00030(78)	1.0005(16)

Table 6. Table 6: Final results for the lowest two energy levels of the b ¯ ​ b ¯ ​ u ​ d ¯ 𝑏 ¯ 𝑏 𝑢 𝑑 \bar{b}\bar{b}ud system on each ensemble, in lattice units and in MeV. All results are given relative to the B ​ B ∗ 𝐵 superscript 𝐵 BB^{*} threshold, i.e.  Δ ​ E n = E n − E B − E B ∗ Δ subscript 𝐸 𝑛 subscript 𝐸 𝑛 subscript 𝐸 𝐵 subscript 𝐸 superscript 𝐵 \Delta E_{n}=E_{n}-E_{B}-E_{B^{*}} .

Ensemble	$a\Delta E_0$	$a\Delta E_1$	$\Delta E_0$ [MeV]	$\Delta E_1$ [MeV]
C00078	$-0.075\left(16\right)$	$0.001\left(21\right)$	$-129\left(27\right)$	$2\left(36\right)$
C005	$-0.069\left(19\right)$	$0.015\left(30\right)$	$-123\left(34\right)$	$27\left(53\right)$
C01	$-0.061\left(15\right)$	$-0.014\left(22\right)$	$-109\left(26\right)$	$-26\left(39\right)$
F004	$-0.051\left(11\right)$	$-0.011\left(15\right)$	$-122\left(25\right)$	$-25\left(37\right)$
F006	$-0.037\left(12\right)$	$0.025\left(32\right)$	$-88\left(29\right)$	$61\left(76\right)$

Table 7. Table 7: The inverse scattering length 1 / a 0 1 subscript 𝑎 0 1/a_{0} , the effective range r 0 subscript 𝑟 0 r_{0} , the binding momentum | k BS | subscript 𝑘 BS |k_{\rm BS}| , and the binding energy E binding subscript 𝐸 binding E_{\textrm{binding}} for all ensembles.

Ensemble	$1/a_0$ [fm-1]	$r_0$ [fm]	$\left|k_{\mathrm{BS}}\right|$ [MeV]	$E_{\text{binding}}$ [MeV]
C00078	$-0.7\left(11.0\right)$	$0.4\left(1.3\right)$	$824\left(89\right)$	$-129\left(28\right)$
C005	$0.2\left(7.2\right)$	$0.51\left(88\right)$	$808\left(110\right)$	$-123\left(34\right)$
C01	$-1.24\left(77\right)$	$0.35\left(12\right)$	$765\left(91\right)$	$-109\left(26\right)$
F004	$-1.28\left(76\right)$	$0.33\left(10\right)$	$819\left(85\right)$	$-122\left(25\right)$
F006	$-0.9\left(10.9\right)$	$0.4\left(1.8\right)$	$697\left(113\right)$	$-88\left(29\right)$

Table 8. Table 8: Multi-exponential matrix fit results from the C00078 ensemble.

Matrix	$N$	Fit range	${\chi }^2/d.o.f.$	$E_n-E_B-E_{B^{\ast }}$ [MeV]
$(O_1,O_2)\times (O_1,O_2)$	2	7…10	1.05	$-159\left(31\right)$,	$+61\left(77\right)$
$(O_1,O_2)\times (O_1,O_2)$	2	7…17	1.42	$-136\left(27\right)$,	$+100\left(75\right)$
$(O_1,O_3)\times (O_1,O_3)$	2	7…10	0.82	$-144\left(35\right)$,	$+137\left(85\right)$
$(O_1,O_3)\times (O_1,O_3)$	2	7…17	1.63	$-116\left(30\right)$,	$+204\left(84\right)$
$(O_2,O_3)\times (O_2,O_3)$	2	7…10	1.32	$-151\left(38\right)$,	$+104\left(82\right)$
$(O_2,O_3)\times (O_2,O_3)$	2	7…17	1.29	$-106\left(33\right)$,	$+168\left(81\right)$
$(O_1,O_2,O_3)\times (O_1,O_2,O_3)$	2	8…10	0.96	$-181\left(40\right)$,	$+66\left(129\right)$
$(O_1,O_2,O_3)\times (O_1,O_2,O_3)$	2	8…17	1.24	$-125\left(34\right)$,	$+142\left(123\right)$
$(O_1,O_2,O_3)\times (O_1,O_2,O_3)$	3	6…17	1.41	$-93\left(24\right)$,	$+245\left(55\right)$,	$+612\left(328\right)$
$(O_1,O_2,O_3)\times (O_1,O_2,O_3)$	3	6…24	1.28	$-86\left(23\right)$,	$+260\left(55\right)$,	$+651\left(337\right)$
$(O_1,O_2,O_3)\times (O_1,O_2,O_3)$	3	7…10	0.98	$-165\left(31\right)$,	$+63\left(77\right)$,	$+6283\left(3782\right)$
$(O_1,O_4)\times \left(O_1\right)$	1	11…17	0.84	$-169\left(96\right)$
$(O_2,O_5)\times \left(O_2\right)$	1	11…15	0.84	$-133\left(128\right)$
$(O_1,O_2,O_4,O_5)\times (O_1,O_2)$	2	11…13	0.90	$-181\left(110\right)$,	$-104\left(199\right)$
$(O_1,O_3,O_4)\times (O_1,O_3)$	2	9…15	1.81	$-125\left(51\right)$,	$+86\left(98\right)$
$(O_1,O_3,O_4)\times (O_1,O_3)$	2	9…24	1.52	$-129\left(53\right)$,	$+79\left(94\right)$
$(O_1,O_3,O_4)\times (O_1,O_3)$	2	10…16	1.72	$-217\left(85\right)$,	$+327\left(187\right)$
$(O_1,O_3,O_4)\times (O_1,O_3)$	3	5…10	1.29	$-169\left(35\right)$,	$+13\left(28\right)$,	$+516\left(100\right)$
$(O_1,O_3,O_4)\times (O_1,O_3)$	3	7…10	1.34	$-175\left(159\right)$,	$-75\left(143\right)$,	$+253\left(538\right)$
$(O_1,O_3,O_4)\times (O_1,O_3)$	3	5…16	1.55	$-169\left(35\right)$,	$+7\left(28\right)$,	$+479\left(89\right)$
$(O_1,O_3,O_4)\times (O_1,O_3)$	3	5…24	1.47	$-168\left(35\right)$,	$+8\left(27\right)$,	$+483\left(88\right)$
$(O_1,O_3,O_4)\times (O_1,O_3)$	3	7…16	1.55	$-250\left(120\right)$,	$-76\left(58\right)$,	$+312\left(207\right)$
$(O_1,O_3,O_4)\times (O_1,O_3)$	3	7…24	1.46	$-237\left(97\right)$,	$-70\left(52\right)$,	$+372\left(233\right)$
$(O_2,O_3,O_5)\times (O_2,O_3)$	2	8…24	1.30	$-124\left(36\right)$,	$+0\left(62\right)$
$(O_2,O_3,O_5)\times (O_2,O_3)$	2	10…24	1.30	$-182\left(116\right)$,	$+290\left(173\right)$
$(O_2,O_3,O_5)\times (O_2,O_3)$	3	5…10	1.52	$-146\left(36\right)$,	$+5\left(31\right)$,	$+510\left(105\right)$
$(O_2,O_3,O_5)\times (O_2,O_3)$	3	7…10	1.99	$-163\left(55\right)$,	$-19\left(111\right)$,	$+520\left(708\right)$
$(O_2,O_3,O_5)\times (O_2,O_3)$	3	5…16	1.40	$-134\left(36\right)$,	$+2\left(31\right)$,	$+452\left(95\right)$
$(O_2,O_3,O_5)\times (O_2,O_3)$	3	5…24	1.13	$-138\left(37\right)$,	$-2\left(30\right)$,	$+450\left(92\right)$
$(O_1,O_2,O_3,O_4,O_5)\times (O_1,O_2,O_3)$	2	10…16	1.34	$-102\left(60\right)$,	$-54\left(70\right)$
$(O_1,O_2,O_3,O_4,O_5)\times (O_1,O_2,O_3)$	2	10…24	1.28	$-119\left(57\right)$,	$-63\left(71\right)$
$(O_1,O_2,O_3,O_4,O_5)\times (O_1,O_2,O_3)$	3	6…16	2.61	$-122\left(28\right)$,	$-21\left(31\right)$,	$+406\left(254\right)$
$(O_1,O_2,O_3,O_4,O_5)\times (O_1,O_2,O_3)$	3	6…24	2.11	$-143\left(54\right)$,	$-34\left(30\right)$,	$+359\left(225\right)$
$(O_1,O_2,O_3,O_4,O_5)\times (O_1,O_2,O_3)$	3	7…16	1.74	$-155\left(27\right)$,	$-78\left(39\right)$,	$+1736\left(1580\right)$

Table 9. Table 9: Multi-exponential matrix fit results from the C005 ensemble.

Matrix	$N$	Fit range	${\chi }^2/d.o.f.$	$E_n-E_B-E_{B^{\ast }}$ [MeV]
$(O_1,O_2)\times (O_1,O_2)$	2	7…22	0.84	$-58.7\left(8.2\right)$,	$+187\left(15\right)$
$(O_1,O_2)\times (O_1,O_2)$	2	14…24	0.87	$-121\left(57\right)$,	$+535\left(464\right)$
$(O_1,O_3)\times (O_1,O_3)$	2	7…22	0.94	$-59.1\left(8.8\right)$,	$+194\left(15\right)$
$(O_1,O_3)\times (O_1,O_3)$	2	14…24	0.75	$-165\left(70\right)$,	$+571\left(459\right)$
$(O_2,O_3)\times (O_2,O_3)$	2	7…22	1.03	$-60.6\left(9.1\right)$,	$+189\left(15\right)$
$(O_2,O_3)\times (O_2,O_3)$	2	14…24	0.99	$-175\left(65\right)$,	$+602\left(475\right)$
$(O_1,O_2,O_3)\times (O_1,O_2,O_3)$	2	11…24	0.90	$-72\left(22\right)$,	$+282\left(79\right)$
$(O_1,O_2,O_3)\times (O_1,O_2,O_3)$	2	14…24	0.80	$-116\left(58\right)$,	$+628\left(505\right)$
$(O_1,O_2,O_3)\times (O_1,O_2,O_3)$	3	6…24	0.88	$-50.6\left(7.3\right)$,	$+199\left(11\right)$,	$+696\left(74\right)$
$(O_1,O_2,O_3)\times (O_1,O_2,O_3)$	3	8…24	0.90	$-57\left(11\right)$,	$+180\left(20\right)$,	$+362\left(264\right)$
$(O_1,O_2,O_3)\times (O_1,O_2,O_3)$	3	14…24	0.78	$-174\left(93\right)$,	$+550\left(466\right)$,	$+1220\left(1642\right)$
$(O_1,O_4)\times \left(O_1\right)$	1	13…24	0.77	$-74\left(39\right)$
$(O_1,O_4)\times \left(O_1\right)$	2	7…24	0.50	$-63\left(43\right)$,	$+113\left(74\right)$
$(O_2,O_5)\times \left(O_2\right)$	1	13…24	0.90	$-6\left(52\right)$
$(O_2,O_5)\times \left(O_2\right)$	2	7…24	0.75	$-39\left(34\right)$,	$+232\left(84\right)$
$(O_1,O_2,O_4,O_5)\times (O_1,O_2)$	2	12…24	1.07	$-100\left(26\right)$,	$+22\left(45\right)$
$(O_1,O_2,O_4,O_5)\times (O_1,O_2)$	2	13…24	0.89	$-94\left(35\right)$,	$-8\left(66\right)$
$(O_1,O_2,O_4,O_5)\times (O_1,O_2)$	2	14…24	0.85	$-161\left(49\right)$,	$+13\left(110\right)$
$(O_1,O_3,O_4)\times (O_1,O_3)$	2	12…24	1.10	$-92\left(25\right)$,	$+6\left(50\right)$
$(O_1,O_3,O_4)\times (O_1,O_3)$	2	13…24	0.92	$-93\left(35\right)$,	$+11\left(75\right)$
$(O_1,O_3,O_4)\times (O_1,O_3)$	3	7…24	0.85	$-59.2\left(8.9\right)$,	$+9\left(31\right)$,	$+193\left(16\right)$
$(O_1,O_3,O_4)\times (O_1,O_3)$	3	11…24	0.79	$-106\left(28\right)$,	$-38\left(161\right)$,	$+246\left(119\right)$
$(O_2,O_3,O_5)\times (O_2,O_3)$	2	14…24	1.01	$-184\left(57\right)$,	$+9\left(126\right)$
$(O_2,O_3,O_5)\times (O_2,O_3)$	2	15…24	0.92	$-164\left(78\right)$,	$+302\left(276\right)$
$(O_2,O_3,O_5)\times (O_2,O_3)$	3	6…24	1.00	$-49.9\left(7.8\right)$,	$+88\left(28\right)$,	$+204\left(21\right)$
$(O_2,O_3,O_5)\times (O_2,O_3)$	3	8…24	0.98	$-60\left(11\right)$,	$-15\left(50\right)$,	$+178\left(20\right)$
$(O_2,O_3,O_5)\times (O_2,O_3)$	3	10…24	1.07	$-65\left(29\right)$,	$-23\left(141\right)$,	$+219\left(47\right)$
$(O_1,O_2,O_3,O_4,O_5)\times (O_1,O_2,O_3)$	2	11…24	1.25	$-100\left(19\right)$,	$-8\left(28\right)$
$(O_1,O_2,O_3,O_4,O_5)\times (O_1,O_2,O_3)$	2	12…24	1.07	$-93\left(26\right)$,	$+25\left(45\right)$
$(O_1,O_2,O_3,O_4,O_5)\times (O_1,O_2,O_3)$	2	13…24	0.96	$-82\left(36\right)$,	$+9\left(64\right)$
$(O_1,O_2,O_3,O_4,O_5)\times (O_1,O_2,O_3)$	2	14…24	0.91	$-148\left(51\right)$,	$+22\left(110\right)$
$(O_1,O_2,O_3,O_4,O_5)\times (O_1,O_2,O_3)$	3	10…24	1.47	$-89\left(15\right)$,	$+40\left(122\right)$,	$+215\left(67\right)$
$(O_1,O_2,O_3,O_4,O_5)\times (O_1,O_2,O_3)$	3	11…24	1.08	$-97\left(20\right)$,	$+136\left(199\right)$,	$+288\left(160\right)$
$(O_1,O_2,O_3,O_4,O_5)\times (O_1,O_2,O_3)$	3	12…24	1.03	$-90\left(27\right)$,	$-14\left(231\right)$,	$+360\left(199\right)$

Table 10. Table 10: Multi-exponential matrix fit results from the C01 ensemble.

Matrix	$N$	Fit range	${\chi }^2/d.o.f.$	$E_n-E_B-E_{B^{\ast }}$ [MeV]
$(O_1,O_2)\times (O_1,O_2)$	2	7…24	1.11	$-37.6\left(7.2\right)$,	$+203\left(13\right)$
$(O_1,O_2)\times (O_1,O_2)$	2	14…24	1.03	$-104\left(47\right)$,	$-57\left(97\right)$
$(O_1,O_3)\times (O_1,O_3)$	2	7…24	1.04	$-37.7\left(7.5\right)$,	$+210\left(13\right)$
$(O_1,O_3)\times (O_1,O_3)$	2	14…24	1.02	$-108\left(52\right)$,	$-40\left(91\right)$
$(O_2,O_3)\times (O_2,O_3)$	2	6…24	0.94	$-29.0\left(6.8\right)$,	$+218\left(10\right)$
$(O_2,O_3)\times (O_2,O_3)$	2	14…24	0.80	$-84\left(51\right)$,	$-55\left(85\right)$
$(O_1,O_2,O_3)\times (O_1,O_2,O_3)$	2	11…24	1.02	$-58\left(17\right)$,	$+88\left(51\right)$
$(O_1,O_2,O_3)\times (O_1,O_2,O_3)$	2	14…24	0.96	$-104\left(49\right)$,	$-60\left(66\right)$
$(O_1,O_2,O_3)\times (O_1,O_2,O_3)$	3	6…24	0.99	$-31.1\left(6.5\right)$,	$+211\left(10\right)$,	$+825\left(69\right)$
$(O_1,O_2,O_3)\times (O_1,O_2,O_3)$	3	8…24	1.04	$-36.1\left(8.5\right)$,	$+208\left(18\right)$,	$+1281\left(434\right)$
$(O_1,O_2,O_3)\times (O_1,O_2,O_3)$	3	10…24	1.07	$-50\left(13\right)$,	$+156\left(36\right)$,	$+2821\left(3773\right)$
$(O_1,O_4)\times \left(O_1\right)$	1	13…24	1.16	$-96\left(25\right)$
$(O_1,O_4)\times \left(O_1\right)$	2	8…24	1.13	$-114\left(25\right)$,	$+240\left(63\right)$
$(O_2,O_5)\times \left(O_2\right)$	1	13…24	0.85	$-70\left(31\right)$
$(O_2,O_5)\times \left(O_2\right)$	2	8…24	1.18	$-66\left(26\right)$,	$+394\left(98\right)$
$(O_1,O_2,O_4,O_5)\times (O_1,O_2)$	2	13…24	1.24	$-96\left(25\right)$,	$-26\left(41\right)$
$(O_1,O_2,O_4,O_5)\times (O_1,O_2)$	2	14…24	1.10	$-121\left(30\right)$,	$-90\left(46\right)$
$(O_1,O_2,O_4,O_5)\times (O_1,O_2)$	2	15…24	1.17	$-116\left(46\right)$,	$-79\left(84\right)$
$(O_1,O_3,O_4)\times (O_1,O_3)$	2	13…24	1.17	$-92\left(24\right)$,	$-12\left(50\right)$
$(O_1,O_3,O_4)\times (O_1,O_3)$	2	14…24	0.95	$-138\left(39\right)$,	$-41\left(66\right)$
$(O_1,O_3,O_4)\times (O_1,O_3)$	3	7…24	1.19	$-36.2\left(6.9\right)$,	$+30\left(25\right)$,	$+210\left(11\right)$
$(O_2,O_3,O_5)\times (O_2,O_3)$	2	12…24	1.21	$-107\left(22\right)$,	$-56\left(36\right)$
$(O_2,O_3,O_5)\times (O_2,O_3)$	2	13…24	0.94	$-79\left(32\right)$,	$-21\left(55\right)$
$(O_2,O_3,O_5)\times (O_2,O_3)$	3	6…24	1.11	$-30.0\left(7.5\right)$,	$+87\left(26\right)$,	$+244\left(38\right)$
$(O_2,O_3,O_5)\times (O_2,O_3)$	3	8…24	1.10	$-44\left(12\right)$,	$+59\left(32\right)$,	$+322\left(74\right)$
$(O_2,O_3,O_5)\times (O_2,O_3)$	3	10…24	1.01	$-123\left(64\right)$,	$-26\left(27\right)$,	$+241\left(116\right)$
$(O_1,O_2,O_3,O_4,O_5)\times (O_1,O_2,O_3)$	2	13…24	1.14	$-93\left(25\right)$,	$-21\left(43\right)$
$(O_1,O_2,O_3,O_4,O_5)\times (O_1,O_2,O_3)$	2	14…24	1.07	$-118\left(30\right)$,	$-95\left(47\right)$

Table 11. Table 11: Multi-exponential matrix fit results from the F004 ensemble.

Matrix	$N$	Fit range	${\chi }^2/d.o.f.$	$E_n-E_B-E_{B^{\ast }}$ [MeV]
$(O_1,O_2)\times (O_1,O_2)$	2	8…24	0.98	$-31.4\left(8.6\right)$,	$+288\left(17\right)$
$(O_1,O_2)\times (O_1,O_2)$	2	13…24	0.90	$-81\left(21\right)$,	$+133\left(79\right)$
$(O_1,O_3)\times (O_1,O_3)$	2	8…24	0.96	$-37.0\left(9.3\right)$,	$+300\left(17\right)$
$(O_1,O_3)\times (O_1,O_3)$	2	13…24	0.70	$-83\left(25\right)$,	$+113\left(79\right)$
$(O_2,O_3)\times (O_2,O_3)$	2	8…24	0.98	$-36.2\left(9.5\right)$,	$+302\left(17\right)$
$(O_2,O_3)\times (O_2,O_3)$	2	13…24	0.95	$-76\left(28\right)$,	$+145\left(80\right)$
$(O_1,O_2,O_3)\times (O_1,O_2,O_3)$	2	12…24	0.90	$-59\left(17\right)$,	$+218\left(56\right)$
$(O_1,O_2,O_3)\times (O_1,O_2,O_3)$	2	13…24	0.85	$-78\left(21\right)$,	$+132\left(78\right)$
$(O_1,O_4)\times \left(O_1\right)$	1	14…24	1.20	$-118\left(24\right)$
$(O_1,O_4)\times \left(O_1\right)$	2	6…24	0.71	$-111\left(14\right)$,	$+370\left(24\right)$
$(O_2,O_5)\times \left(O_2\right)$	1	14…24	0.82	$-42\left(29\right)$
$(O_2,O_5)\times \left(O_2\right)$	2	5…24	0.82	$-67\left(11\right)$,	$+495\left(20\right)$
$(O_1,O_2,O_4,O_5)\times (O_1,O_2)$	2	14…22	0.99	$-120\left(22\right)$,	$-55\left(34\right)$
$(O_1,O_2,O_4,O_5)\times (O_1,O_2)$	2	15…22	0.88	$-119\left(29\right)$,	$-58\left(37\right)$
$(O_1,O_3,O_4)\times (O_1,O_3)$	2	14…24	0.81	$-135\left(26\right)$,	$-5\left(49\right)$
$(O_1,O_3,O_4)\times (O_1,O_3)$	2	15…24	0.71	$-135\left(34\right)$,	$+12\left(70\right)$
$(O_2,O_3,O_5)\times (O_2,O_3)$	2	14…24	1.10	$-68\left(29\right)$,	$-37\left(51\right)$
$(O_1,O_2,O_3,O_4,O_5)\times (O_1,O_2,O_3)$	2	14…24	1.00	$-119\left(21\right)$,	$-61\left(32\right)$
$(O_1,O_2,O_3,O_4,O_5)\times (O_1,O_2,O_3)$	2	15…24	0.93	$-115\left(28\right)$,	$-55\left(34\right)$

Table 12. Table 12: Multi-exponential matrix fit results from the F006 ensemble.

Matrix	$N$	Fit range	${\chi }^2/d.o.f.$	$E_n-E_B-E_{B^{\ast }}$ [MeV]
$(O_1,O_2)\times (O_1,O_2)$	2	9…24	0.92	$-31.0\left(6.8\right)$,	$+240\left(13\right)$
$(O_1,O_2)\times (O_1,O_2)$	2	13…24	0.90	$-50\left(13\right)$,	$+205\left(40\right)$
$(O_1,O_3)\times (O_1,O_3)$	2	9…24	0.96	$-31.6\left(7.2\right)$,	$+253\left(13\right)$
$(O_1,O_3)\times (O_1,O_3)$	2	13…24	0.70	$-62\left(15\right)$,	$+192\left(41\right)$
$(O_2,O_3)\times (O_2,O_3)$	2	9…22	1.08	$-32.0\left(7.5\right)$,	$+249\left(14\right)$
$(O_2,O_3)\times (O_2,O_3)$	2	13…24	0.88	$-64\left(15\right)$,	$+208\left(41\right)$
$(O_1,O_2,O_3)\times (O_1,O_2,O_3)$	2	13…24	0.91	$-47\left(12\right)$,	$+199\left(40\right)$
$(O_1,O_2,O_3)\times (O_1,O_2,O_3)$	2	15…24	0.86	$-57\left(19\right)$,	$+231\left(80\right)$
$(O_1,O_4)\times \left(O_1\right)$	1	17…24	1.11	$-66\left(28\right)$
$(O_1,O_4)\times \left(O_1\right)$	2	7…24	0.99	$-107\left(10\right)$,	$+350\left(21\right)$
$(O_2,O_5)\times \left(O_2\right)$	1	17…24	0.73	$-30\left(36\right)$
$(O_2,O_5)\times \left(O_2\right)$	2	6…24	0.99	$-63.6\left(8.9\right)$,	$+479\left(18\right)$
$(O_1,O_2,O_4,O_5)\times (O_1,O_2)$	2	17…24	1.29	$-88\left(27\right)$,	$+28\left(63\right)$
$(O_1,O_3,O_4)\times (O_1,O_3)$	2	17…24	1.06	$-84\left(28\right)$,	$+7\left(70\right)$
$(O_1,O_3,O_4)\times (O_1,O_3)$	2	18…24	0.88	$-89\left(34\right)$,	$+74\left(109\right)$
$(O_2,O_3,O_5)\times (O_2,O_3)$	2	17…24	1.21	$-92\left(35\right)$,	$+109\left(87\right)$
$(O_2,O_3,O_5)\times (O_2,O_3)$	2	18…24	0.98	$-89\left(51\right)$,	$+67\left(120\right)$
$(O_1,O_2,O_3,O_4,O_5)\times (O_1,O_2,O_3)$	2	15…24	1.44	$-95\left(17\right)$,	$+24\left(31\right)$
$(O_1,O_2,O_3,O_4,O_5)\times (O_1,O_2,O_3)$	2	16…24	1.21	$-93\left(21\right)$,	$+32\left(44\right)$
$(O_1,O_2,O_3,O_4,O_5)\times (O_1,O_2,O_3)$	2	17…24	1.09	$-83\left(27\right)$,	$+19\left(61\right)$
$(O_1,O_2,O_3,O_4,O_5)\times (O_1,O_2,O_3)$	2	18…24	0.93	$-74\left(34\right)$,	$+66\left(92\right)$

Table 13. Table 13: Comparison of the results for the two lowest b ¯ ​ b ¯ ​ u ​ d ¯ 𝑏 ¯ 𝑏 𝑢 𝑑 \bar{b}\bar{b}ud energy levels (relative to the B ​ B ∗ 𝐵 superscript 𝐵 ∗ BB^{\ast} threshold, i.e., Δ ​ E n = E n − E B − E B ∗ Δ subscript 𝐸 𝑛 subscript 𝐸 𝑛 subscript 𝐸 𝐵 subscript 𝐸 superscript 𝐵 \Delta E_{n}=E_{n}-E_{B}-E_{B^{*}} ) from multi-exponential matrix fitting and from the GEVP, for the C005 ensemble.

Operators	Energy difference	Multi-exponential fitting			GEVP
		Fit range	${\chi }^2/d.o.f.$		Fit range	${\chi }^2/d.o.f.$
$(O_1,O_2)\times (O_1,O_2)$	$\Delta E_0$	7…22	0.84	$-58.7\left(8.2\right)$	7…20	0.70	$-60.5\left(11.4\right)$
	$\Delta E_1$			$+187\left(15\right)$	7…14	0.53	$+183\left(19\right)$
$(O_1,O_2)\times (O_1,O_2)$	$\Delta E_0$	14…24	0.87	$-121\left(57\right)$	14…20	0.61	$-129\left(67\right)$
$(O_1,O_3)\times (O_1,O_3)$	$\Delta E_0$	7…22	0.94	$-59.1\left(8.8\right)$	7…20	0.80	$-62.1\left(12.1\right)$
	$\Delta E_1$			$+194\left(15\right)$	7…14	0.51	$+188\left(19\right)$
$(O_1,O_3)\times (O_1,O_3)$	$\Delta E_0$	14…24	0.75	$-165\left(70\right)$	14…20	0.42	$-160\left(76\right)$
$(O_2,O_3)\times (O_2,O_3)$	$\Delta E_0$	7…22	1.03	$-60.6\left(9.1\right)$	7…20	0.86	$-61.4\left(12.3\right)$
	$\Delta E_1$			$+189\left(15\right)$	7…14	0.56	$+188\left(19\right)$
$(O_2,O_3)\times (O_2,O_3)$	$\Delta E_0$	14…24	0.99	$-175\left(65\right)$	14…20	0.44	$-164\left(79\right)$
$(O_1,O_2,O_3)\times (O_1,O_2,O_3)$	$\Delta E_0$	6…24	0.88	$-50.6\left(7.3\right)$	6…17	1.32	$-53.4\left(10.2\right)$
	$\Delta E_1$			$+199\left(11\right)$	6…14	0.85	$+198\left(14\right)$
$(O_1,O_2,O_3)\times (O_1,O_2,O_3)$	$\Delta E_0$	8…24	0.90	$-57\left(11\right)$	8…17	0.85	$-62.8\left(13.7\right)$
	$\Delta E_1$			$+180\left(20\right)$	8…14	0.57	$+176\left(24\right)$
$(O_1,O_2,O_3)\times (O_1,O_2,O_3)$	$\Delta E_0$	11…24	0.90	$-72\left(22\right)$	11…17	0.73	$-92\left(29\right)$
$(O_1,O_2,O_3)\times (O_1,O_2,O_3)$	$\Delta E_0$	14…24	0.78	$-174\left(93\right)$	14…17	0.43	$-156\left(80\right)$

Equations

$$b=\left(\begin{array}[]{c}\psi\\ \chi\end{array}\right),\quad\bar{b}=\left(\begin{array}[]{cc}\psi^{\dagger},&-\chi^{\dagger}\end{array}\right),$$

$$G_{b}(\textbf{x},t,\textbf{x}^{\prime},t^{\prime})=\Theta(t-t^{\prime})\left(\begin{array}[]{cc}G_{\psi}(\textbf{x},t,\textbf{x}^{\prime},t^{\prime})&0\\ 0&0\end{array}\right)-\Theta(-t+t^{\prime})\left(\begin{array}[]{cc}0&0\\ 0&G_{\chi}(\textbf{x},t,\textbf{x}^{\prime},t^{\prime})\end{array}\right),$$

$$\displaystyle O_{1}=O_{[BB^{\ast}](0)}=\sum_{\textbf{x}}\Big{(}\bar{b}(\textbf{x})\gamma_{5}d(\textbf{x})\Big{)}\Big{(}\bar{b}(\textbf{x})\gamma_{j}u(\textbf{x})\Big{)}-(d\leftrightarrow u)$$

$$\displaystyle O_{2}=O_{[B^{\ast}B^{\ast}](0)}=\epsilon_{jkl}\sum_{\textbf{x}}\Big{(}\bar{b}(\textbf{x})\gamma_{k}d(\textbf{x})\Big{)}\Big{(}\bar{b}(\textbf{x})\gamma_{l}u(\textbf{x})\Big{)}-(d\leftrightarrow u)$$

$$\displaystyle O_{3}=O_{[Dd](0)}=\sum_{\textbf{x}}\Big{(}\epsilon^{abc}\bar{b}(\textbf{x})^{b}\gamma_{j}\mathcal{C}\bar{b}^{c,T}(\textbf{x})\Big{)}\Big{(}\epsilon^{ade}d^{d,T}(\textbf{x})\mathcal{C}\gamma_{5}u^{e}(\textbf{x})\Big{)}-(d\leftrightarrow u).$$

$$\displaystyle O_{4}=O_{B(0)B^{\ast}(0)}=\bigg{(}\sum_{\textbf{x}}\bar{b}(\textbf{x})\gamma_{5}d(\textbf{x})\bigg{)}\bigg{(}\sum_{\textbf{y}}\bar{b}(\textbf{y})\gamma_{j}u(\textbf{y})\bigg{)}-(d\leftrightarrow u)$$

$$\displaystyle O_{5}=O_{B^{\ast}(0)B^{\ast}(0)}=\epsilon_{jkl}\bigg{(}\sum_{\textbf{x}}\bar{b}(\textbf{x})\gamma_{k}d(\textbf{x})\bigg{)}\bigg{(}\sum_{\textbf{y}}\bar{b}(\textbf{y})\gamma_{l}u(\textbf{y})\bigg{)}-(d\leftrightarrow u).$$

$$q_{\mathrm{smeared}}={\left(1+\frac{{\sigma }_{\text{Gauss}}^2}{4N_{\text{Gauss}}}\Delta \right)}^{N_{\text{Gauss}}}q,$$

$$\displaystyle C_{jk}(t)=\Big{\langle}O_{j}(t)O^{\dagger}_{k}(0)\Big{\rangle},$$

$$O_B(\mathbf{p})=\sum _{\text{x}}{e}^{i\mathbf{p}\text{x}}\bar{b}(\text{x}){\gamma }_{5}u(\text{x}),$$

$$O_{B^{\ast }}(\mathbf{p})=\sum _{\text{x}}{e}^{i\mathbf{p}\text{x}}\bar{b}(\text{x}){\gamma }_{j}u(\text{x}),$$

$$E_B(\mathbf{p})$$

$$E_{B^{\ast }}(\mathbf{p})$$

$$m_{\mathrm{kin}}=\frac{{\mathbf{p}}^2-[E(\mathbf{p})-E(0){]}^2}{2[E(\mathbf{p})-E(0)]}$$

$$c^2=\frac{[E(\mathbf{p})-E(0)+m_{\mathrm{kin},1}{]}^2-m_{\mathrm{kin},1}^2}{{\mathbf{p}}^2},$$

$$C_{jk}(t)=\sum _{n=0}^{\infty }⟨\Omega |O_j|n⟩⟨n|O_k^{†}|\Omega ⟩{\text{e}}^{-E_{n}t},$$

$$Z_j^n=⟨\Omega |O_j|n⟩.$$

$$C_{jk}(t)\approx \sum _{n=0}^{N-1}{Z}_j^n{Z}_k^n{\text{e}}^{-E_{n}t},$$

$$(\text{CAA sample}{)}_e=(\text{exact sample}{)}_e-(\text{sloppy sample}{)}_{e,0}+\frac{1}{N_{\text{sloppy}}}\sum _{s=0}^{N_{\text{sloppy}}-1}(\text{sloppy sample}{)}_{e,s}.$$

$$(\text{MCAA sample}{)}_{e,s}=(\text{exact sample}{)}_e-(\text{sloppy sample}{)}_{e,0}+(\text{sloppy sample}{)}_{e,s}.$$

$$O_j^{†}|\Omega ⟩=\sum _{n=0}^{\infty }|n⟩⟨n|{\mathcal{O}}_j^{†}|\Omega ⟩=\sum _{n=0}^{\infty }{Z}_j^n|n⟩.$$

$$|{\tilde{Z}}_j^n{|}^2=\frac{|Z_j^n{|}^2}{{\text{max}}_m(|Z_j^m{|}^2)},$$

$$E_n=E_B+\sqrt{m_{B,\mathrm{kin}}^2+k_n^2}-m_{B,\mathrm{kin}}+E_{B^{\ast }}+\sqrt{m_{B^{\ast },\mathrm{kin}}^2+k_n^2}-m_{B^{\ast },\mathrm{kin}},$$

$$k_n^2=\frac{\Delta E_n(\Delta E_n+2m_{B,\mathrm{kin}})(\Delta E_n+2m_{B^{\ast },\mathrm{kin}})(\Delta E_n+2m_{B,\mathrm{kin}}+2m_{B^{\ast },\mathrm{kin}})}{4(\Delta E_n+m_{B,\mathrm{kin}}+m_{B^{\ast },\mathrm{kin}}{)}^2},$$

$$\Delta E_n=E_n-E_B-E_{B^{\ast }}.$$

$$\cot {\delta }_0(k_n)=\frac{2Z_{00}(1;(k_{n}L/2\pi {)}^2)}{{\pi }^{1/2}{k}_{n}L},$$

$$T_0(k)=\frac{1}{\cot {\delta }_0(k)-i}.$$

$$k\cot {\delta }_0(k)=\frac{1}{a_0}+\frac{1}{2}{r}_0{k}^2+\mathcal{O}(k^4),$$

$$-|k_{\mathrm{BS}}|=\frac{1}{a_0}-\frac{1}{2}{r}_0|k_{\mathrm{BS}}{|}^2,$$

$$E_{\mathrm{binding}}=E_{\mathrm{BS}}-E_B-E_{B^{\ast }}$$

$$=\sqrt{m_{B,\mathrm{kin}}^2-|k_{\mathrm{BS}}{|}^2}-m_{B,\mathrm{kin}}+\sqrt{m_{B^{\ast },\mathrm{kin}}^2-|k_{\mathrm{BS}}{|}^2}-m_{B^{\ast },\mathrm{kin}}.$$

$$E_{\text{binding}}(m_{\pi })=E_{\text{binding}}(m_{\pi ,\text{phys}})+c(m_{\pi }^2-m_{\pi ,\text{phys}}^2),$$

$$E_{\text{binding}}(m_{\pi ,\text{phys}})=(-128\pm 24)\text{MeV},c=(1.5\pm 2.3)\times 10^{-4}{\text{MeV}}^{-2},$$

Full text

††preprint: JLAB-THY-19-2912††preprint: RBRC-1307

I Introduction

Mesons, i.e., hadrons with integer spin, were first envisioned by Gell-Mann and Zweig GellMann:1964nj ; Zweig:1964jf to be built from one, two or more quark-antiquark pairs. However, systems that manifestly contain more than a single quark-antiquark pair were found only relatively recently, primarily in the heavy-quark sector Belle:2011aa ; Olsen:2015zcy ; Lebed:2016hpi ; Esposito:2016noz ; Richard:2016eis ; Olsen:2017bmm . Exotic mesons can be characterized as having $J^{PC}$ quantum numbers that cannot be constructed in the simple quark-antiquark model, or as having a manifestly exotic quark flavor content. In this work, we consider an example for the latter, a $\bar{b}\bar{b}ud$ tetraquark.111In the literature, the term “tetraquark” is somewhat ambiguous. In certain papers it exclusively refers to a diquark-antidiquark structure, while in other papers it is used more generally for arbitrary bound states and resonances with a strong four-quark component, including, e.g., mesonic molecules. Throughout this paper we follow the latter convention. Moreover, the $\bar{b}\bar{b}ud$ system is a tetraquark in a fully rigorous sense, since it contains four net quark flavors.

It can be shown that QCD-stable $\bar{Q}\bar{Q}qq$ tetraquarks must exist in the limit $m_{Q}\to\infty$ Carlson:1987hh ; Manohar:1992nd ; Eichten:2017ffp . In this limit, the two heavy antiquarks form a color-triplet object with a size of order $(\alpha_{s}m_{Q})^{-1}$ and a binding energy of order $\alpha_{s}^{2}m_{Q}$ due to the attractive Coulomb potential at short distances. The doubly-heavy $\bar{Q}\bar{Q}qq$ tetraquarks then become related to singly-heavy $Qqq$ baryons, just like doubly-heavy $\bar{Q}\bar{Q}\bar{q}$ baryons become related to singly-heavy $Q\bar{q}$ mesons Savage:1990di ; Brambilla:2005yk ; Cohen:2006jg ; Mehen:2017nrh . The question is whether the physical bottom quark is heavy enough for $\bar{b}\bar{b}qq$ bound states to exist below the $\bar{b}q$-$\bar{b}q$ two-meson thresholds. Studies based on potential models, effective field theories, and QCD sum rules suggest that this is indeed the case Carlson:1987hh ; Manohar:1992nd ; SilvestreBrac:1993ss ; Brink:1998as ; Vijande:2003ki ; Janc:2004qn ; Vijande:2006jf ; Navarra:2007yw ; Ebert:2007rn ; Zhang:2007mu ; Lee:2009rt ; Karliner:2017qjm ; Eichten:2017ffp ; Wang:2017uld ; Richard:2018yrm ; Park:2018wjk ; Wang:2018atz ; Liu:2019stu . Possible experimental search strategies for bottomness-$2$ tetraquarks are discussed in Refs. Moinester:1995fk ; Ali:2018ifm ; Ali:2018xfq .

Within lattice QCD, $\bar{b}\bar{b}qq$ four-quark systems were explored for the first time using static $\bar{b}$ quarks and the Born-Oppenheimer approximation. A stable $\bar{b}\bar{b}ud$ tetraquark with quantum numbers $I(J^{P})=0(1^{+})$ around $30\ldots 90\,\textrm{MeV}$ below the $BB^{\ast}$ threshold as well as a $\bar{b}\bar{b}ud$ tetraquark resonance with quantum numbers $I(J^{P})=0(1^{-})$ around $15\,\textrm{MeV}$ above the $BB$ threshold were predicted Bicudo:2012qt ; Brown:2012tm ; Bicudo:2015kna ; Bicudo:2016ooe ; Bicudo:2017szl . Effects from the heavy-quark spin were investigated for the stable $I(J^{P})=0(1^{+})$ tetraquark by solving a coupled-channel Schrödinger equation in Ref. Bicudo:2016ooe . Moreover, several flavor combinations were explored and no stable $\bar{b}\bar{b}qq$ tetraquarks with $qq=ss$ and $qq=cc$ were found in this approach Bicudo:2015vta . Recently, the same $\bar{b}\bar{b}qq$ four-quark systems have been investigated with $\bar{b}$ quarks of finite mass treated within nonrelativistic QCD (NRQCD). A stable $\bar{b}\bar{b}ud$ tetraquark with quantum numbers $I(J^{P})=0(1^{+})$ was also seen in two such computations Francis:2016hui ; Junnarkar:2018twb , but there is a quantitative difference by a factor $\approx 2\ldots 3$ in the binding energy between Refs. Francis:2016hui ; Junnarkar:2018twb and Ref. Bicudo:2016ooe , which is not yet understood. Moreover, $\bar{Q}\bar{Q}qq$ systems with further flavor combinations $\bar{Q}\bar{Q}\in\{\bar{b}\bar{b},\bar{b}\bar{c},\bar{c}\bar{c}\}$ and $q\in\{u,d,s,c\}$ have been investigated and some indication has been obtained that systems with $J^{P}=1^{+}$ and $\bar{Q}\bar{Q}qq\in\{\bar{b}\bar{b}us,\bar{b}\bar{b}uc,\bar{b}\bar{b}sc\ ,\ \bar{b}\bar{c}ud\ ,\ \bar{c}\bar{c}ud\}$ are stable as well Francis:2018jyb ; Junnarkar:2018twb .

In this paper we perform a lattice QCD study of the $\bar{b}\bar{b}ud$ four-quark system with quantum numbers $I(J^{P})=0(1^{+})$, using NRQCD $\bar{b}$ quarks and domain-wall light quarks (results obtained at an early stage of this project have been presented in Ref. Peters:2016isf ). We make use of both local interpolating fields (in which the four quarks are jointly projected to zero momentum) and nonlocal interpolating fields (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 systems Mohler:2013rwa ; Lang:2014yfa that including both types of interpolating fields is required to reliably determine ground-state energies in exotic channels. In this way we expand on the works of Refs. Francis:2016hui ; Francis:2018jyb ; Junnarkar:2018twb , where nonlocal interpolating fields were not considered. Having both local and nonlocal interpolating fields allows us to determine the ground-state and first-excited-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.

II Lattice Setup

II.1 Gauge-link configurations and light-quark propagators

We performed the computations presented here using gauge-link configurations generated by the RBC and UKQCD collaborations Aoki:2010dy ; Blum:2014tka with 2+1 flavors of domain-wall fermions Kaplan:1992bt ; Furman:1994ky ; Shamir:1993zy ; Brower:2012vk and the Iwasaki gauge action Iwasaki:1984cj . We use the five ensembles listed in Table 1, which differ in the lattice spacing $a\approx 0.083\,\textrm{fm}\ldots 0.114\,\textrm{fm}$, the lattice size (spatial extent $\approx 2.65\,\textrm{fm}\ldots 5.48\,\textrm{fm}$) and the pion mass $m_{\pi}\approx 139\,\textrm{MeV}\ldots 431\,\textrm{MeV}$. Ensemble C00078 uses the Möbius domain-wall action Brower:2012vk with length of the fifth dimension $N_{5}=24$, while the other ensembles use the Shamir action Shamir:1993zy with $N_{5}=16$. The lattice spacings listed in the Table were determined in Ref. Blum:2014tka .

Our calculation uses smeared point-to-all propagators for the up and down quarks (the smearing parameters are given in Sec. III.1). The computational cost of generating these propagators was reduced using the all-mode-averaging technique Blum:2012uh ; Shintani:2014vja . 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 Blum:2012uh ; Shintani:2014vja . 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_{\textrm{CG, sl}}$, than needed for the exact propagators. The lowest $N_{\textrm{EV}}$ eigenvectors of the domain-wall operator were computed using Lanczos with Chebyshev-polynomial acceleration. The values of $N_{\textrm{EV}}$ and $N_{\textrm{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.

II.2 Bottom-quark propagators

The heavy $b$ quarks are treated with the framework of lattice nonrelativistic QCD (NRQCD) ThackerLepage91 ; Lepage:1992tx . We use the same lattice NRQCD action and parameters as in Ref. Brown:2014ena . This action includes all quark-bilinear operators through order $v^{4}$ in the heavy-heavy power counting, and order $\Lambda^{2}/m_{b}^{2}$ 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. Meinel:2010pv . We use the same masses also on the other coarse and fine ensembles, respectively, as shown in Table 2. 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 Hammant:2011bt . The gauge links entering the NRQCD action are divided by the mean link $u_{0L}$ in Landau gauge to achieve tadpole improvement Lepage:1992xa ; Shakespeare:1998dt .

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 $\psi$ and $\chi$ 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

[TABLE]

and the bottom-quark propagator becomes

[TABLE]

with the two-spinor NRQCD quark and antiquark propagators $G_{\psi}$ and $G_{\chi}$.

III Interpolating operators and correlation functions

III.1 $\bar{b}\bar{b}ud$ 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^{*}$. The lowest three-meson threshold is $BB\pi$, which is approximately 44 MeV above $B^{*}B^{*}$. From Ref. Brown:2014ena we can see that the $\overline{\Xi}_{bb}N$ antibaryon-baryon threshold is already much higher than the $B^{*}B^{*}$ threshold: $11.1\,\textrm{GeV}$ for $\overline{\Xi}_{bb}N$ compared to $10.6\,\textrm{GeV}$ for $B^{*}B^{*}$.

The $J^{P}=1^{+}$ quantum numbers appear in the $T_{1}^{g}$ irreducible representation (irrep) of the $O_{h}$ point group Johnson:1982yq . 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):

[TABLE]

Operators four and five are nonlocal operators, where each color singlet is projected to zero momentum individually:

[TABLE]

Above, $a,b,c,\ldots$ denote color indices, $j,k,l$ spatial vector indices, and $\mathcal{C}=\gamma_{0}\gamma_{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^{\ast}$ mesons are around $45\,\textrm{MeV}$ heavier than a $B$ and a $B^{\ast}$ meson, one might expect that $O_{[B^{\ast}B^{\ast}](0)}$ will have less overlap with the ground state than $O_{[BB^{\ast}](0)}$. However, a previous investigation with static quarks (see Ref. Bicudo:2016ooe 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 Carlson:1987hh ; Manohar:1992nd ; Eichten:2017ffp . The importance of diquark operators was also discussed in Refs. Jaffe:2004ph ; Cheung:2017tnt .

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. Mohler:2013rwa ; Lang:2014yfa 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 Francis:2016hui ; Junnarkar:2018twb 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

[TABLE]

where $\Delta$ is the nearest-neighbor gauge-covariant spatial Laplacian. For the up and down quarks, the gauge links in $\Delta$ are spatially APE-smeared Albanese:1987ds , while the unsmeared gauge links are used for the bottom quark. The smearing parameters are collected in Table 3.

To determine the spectrum we compute the temporal correlation functions of the interpolating operators $O_{1}$ to $O_{5}$,

[TABLE]

where $\langle\ldots\rangle$ 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 $\mathbf{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. Abdel-Rehim:2017dok ; Francis:2018qch .

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))^{\ast}=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.

III.2 $B$ and $B^{\ast}$ meson

Since we will compare the energy levels of the $\bar{b}\bar{b}ud$ four-quark system to the $BB^{\ast}$ threshold, we need to determine the energies of the $B$ and $B^{*}$ mesons within the same setup. We use the interpolating operators

[TABLE]

where we also consider nonzero momenta $\mathbf{p}=2\pi\mathbf{n}/L$, $\mathbf{n}\in\mathbb{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 3. We compute the correlation functions $\langle O_{B}(\mathbf{p},t)O_{B}^{\dagger}(\mathbf{p},0)\rangle$ and $\langle O_{B^{\ast}}(\mathbf{p},t)O_{B^{\ast}}^{\dagger}(\mathbf{p},0)\rangle$. As discussed in Sec. III.1, we perform the summation over $\mathbf{x}$ at the sink only.

IV Energies and kinetic masses of the $B$ and $B^{\ast}$ 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 4. 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 the $\bar{b}\bar{b}ud$ system.

For the scattering analysis in Sec. VI (and to assess the tuning of the $b$-quark mass), we also need the momentum-dependence of the $B$ and $B^{*}$ energies. We find that within statistical uncertainties our results are consistent with the form

[TABLE]

up to the highest momenta we computed ($\mathbf{p}^{2}=3(2\pi/L)^{2}$). The kinetic masses given in Table 4 were extracted using

[TABLE]

with the smallest possible non-vanishing momentum, $\mathbf{p}^{2}=(2\pi/L)^{2}$. To test whether the $B$ and $B^{*}$ energies at higher momenta still satisfy the dispersion relations (12) and (13) with the same parameters, we computed the square of the “speed of light” via

[TABLE]

where $m_{\rm kin,1}$ are the kinetic masses computed with one unit of momentum. The results for $c^{2}$ are listed in Table 5 and are found to be consistent with 1 within their statistical uncertainties of 0.3% or smaller.

On the C00078 and C005 ensembles, the spin-averaged kinetic masses $m_{\rm kin,spinav}=\frac{1}{4}m_{B,\,{\rm kin}}+\frac{3}{4}m_{B^{*}\!,\,{\rm kin}}$ agree with the experimental value of 5.28272(2) GeV Tanabashi:2018oca . 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. Meinel:2010pv , while here we use the lattice spacing determinations of Ref. Blum:2014tka to convert to physical units).

However, the effect of a possible $\lesssim 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 $\Lambda/m_{b}$, and are therefore affected by the same relative error as the $b$-quark mass, which corresponds to an absolute error of $\lesssim 2.5\,\textrm{MeV}$. Our results for the hyperfine splittings are also shown in Table 4. The value from the physical-pion-mass ensemble C00078 is consistent with the experimental value of $45.3(2)\>\textrm{MeV}$ Tanabashi:2018oca . We find a trend of increasing hyperfine splitting as the light-quark mass is increased.

V The lowest energy levels of the $\bar{b}\bar{b}ud$ system

V.1 Multi-exponential matrix fitting

The spectral decomposition of the correlation matrix (9) reads

[TABLE]

where $|\Omega\rangle$ denotes the vacuum state, $|n\rangle$ are the energy eigenstates of the $\bar{b}\bar{b}ud$ system in the $T_{1}^{g}$ irrep and for isospin $I=0$, and $E_{n}$ are the corresponding energy eigenvalues. Equation (16) assumes an infinite time extent of the lattice, which is a good approximation in our case. As discussed in Sec. III.1, all $C_{jk}(t)$ are real and, thus, so are the overlap factors

[TABLE]

To extract the energies $E_{n}$ from our numerical results for $C_{jk}(t)$, we perform fully correlated, least-$\chi^{2}$, multi-exponential fits using a truncated version of (16),

[TABLE]

where we must choose the time range $t_{\textrm{min}}\leq t\leq t_{\textrm{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\times 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 8-12 in Appendix A.

The matrix fits have a large number of degrees of freedom, and the covariance matrix, whose inverse enters in $\chi^{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 Blum:2012uh ; Shintani:2014vja

[TABLE]

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 1, 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 $\chi^{2}$ even for large numbers of degrees of freedom, we use a “modified CAA” (MCAA) procedure, given by

[TABLE]

This procedure provides $N_{\text{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\times 3)$ case. The drawback is that there are autocorrelations between the different choices of $s$, introduced by the constant (but very small) term $(\text{exact sample})_{e}-(\text{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.

The condition numbers of the data covariance matrices range from approximately $10^{6}$ to $10^{15}$, depending on how many interpolating operators are included and on the fit range. We always computed the full inverse using LU decomposition. The minimization of $\chi^{2}$ was performed using the Levenberg-Marquardt algorithm, with initial guesses for the fit parameters obtained as described in Sec. V B of Ref. Alexandrou:2017mpi . The fit results were stable under variations of the initial guesses within reasonable ranges.

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) Luscher:1990ck ; Blossier:2009kd ; Orginos:2015tha . 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.

V.2 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^{\ast}$ 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 the $\bar{b}\bar{b}ud$ system. Note that the first excited energy is close to the $BB^{\ast}$ 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.

V.3 Overlap factors

For a given $j$, the overlap factors $Z_{j}^{n}$ indicate the relative importance of the energy eigenstates $|n\rangle$ when expanding the trial state $O_{j}^{\dagger}|\Omega\rangle$ in terms of energy eigenstates,

[TABLE]

Therefore, the overlap factors $Z_{j}^{n}$ provide certain information about the composition and quark arrangement of the energy eigenstates $|n\rangle$. In particular, if the overlap factor $Z_{j}^{m}$ for one state $|m\rangle$ is significantly larger than all other $Z_{j}^{n}$, $n\neq m$, this might be a sign that $|m\rangle$ is quite similar to $O_{j}^{\dagger}|\Omega\rangle$.

It is convenient to consider rescaled squared overlap factors,

[TABLE]

which are normalized such that $\textrm{max}_{m}(|\tilde{Z}_{j}^{m}|^{2})=1$ for each trial state $O_{j}^{\dagger}|\Omega\rangle$. 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 $|\tilde{Z}_{j}^{n}|^{2}$ obtained from our fits are qualitatively similar for all ensembles, and do not strongly depend on the temporal fit range $t_{\textrm{min}}\leq t\leq t_{\textrm{max}}$.

In Fig. 5 we show the normalized overlap factors obtained on ensemble C005 using the full $5\times 3$ correlation matrix for $N=2$ and $N=3$, with fit ranges of $14\leq t/a\leq 24$ and $12\leq t/a\leq 24$, respectively. One can see that for the trial state created by the diquark-antidiquark operator $O_{3}$, the ground-state overlap is significantly larger than the overlap to the first excited state. Vice versa, for the two trial states created by the nonlocal meson-meson operators $O_{4}$ and $O_{5}$, the overlaps to the first excited state are larger than the ground-state overlaps. This supports our above interpretation concerning the composition of the lowest two energy eigenstates: the ground state $|0\rangle$ seems to be a four-quark bound state and the first excitation $|1\rangle$ a meson-meson scattering state. Finally, it is interesting to note that in the three-exponential fit, the local meson-meson operators $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 (also note that the large values of $|Z_{1}^{2}|^{2}/|Z_{1}^{0}|^{2}=4\pm 5$, $|Z_{2}^{2}|^{2}/|Z_{2}^{0}|^{2}=8\pm 9$ substantially change the normalization of the plots for $O_{1}$ and $O_{2}$ when going from $N=2$ to $N=3$). What appears in the fit as the “second excited state” (with a rather high energy and large uncertainty, as shown in Table 9) is likely an admixture of the dense spectrum of all the scattering states above threshold, which are all created with similar weights by the local operators $O_{1}$ and $O_{2}$, while the nonlocal operators $O_{4}$ and $O_{5}$ mostly create the lowest-lying scattering state.

V.4 Final results for the lowest two energy levels on each ensemble

To obtain final results for the energy differences $\Delta E_{0}$ and $\Delta 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 $\chi^{2}/{\rm 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 $\Delta E_{0}$ and $\Delta 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.1) of these new 500 samples. The results for $\Delta E_{0}$ and $\Delta E_{1}$ obtained in this way are listed in Table 6 and are indicated with the horizontal lines and shaded uncertainty bands in Fig. 4 and Figs. 9-12.

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 Luscher:1990ux and its generalizations Rummukainen:1995vs ; Kim:2005gf ; Christ:2005gi ; Hansen:2012tf ; Leskovec:2012gb ; Gockeler:2012yj ; Briceno:2014oea ; Briceno:2017tce ; see Ref. Briceno:2017max for a recent review. Here, we apply Lüscher’s method to the lowest two energy levels of the $\bar{b}\bar{b}ud$ system in the $T_{1}^{g}$ 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). We neglect higher partial waves and the coupling to other channels, such as $B^{*}B^{*}$. Thus, we obtain the $S$-wave $B$-$B^{*}$ 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.

VI.1 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 the $\bar{b}\bar{b}ud$ system $E_{n}$ through the equation

[TABLE]

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_{n}^{2}$ gives

[TABLE]

where

[TABLE]

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 $\delta_{0}$ is

[TABLE]

where $Z_{00}$ is the generalized zeta function Luscher:1990ux . The scattering amplitude is given by

[TABLE]

VI.2 Effective-range expansion and determination of the bound-state pole

On each ensemble, we use the lattice QCD results for the energy differences $\Delta E_{0}$ and $\Delta E_{1}$ combined with the corresponding results for the $B$ and $B^{*}$ meson energies and kinetic masses to calculate $k_{n}^{2}$ for $n=0$ and $n=1$ using Eq. (24), and we determine the corresponding $k_{n}\cot\delta_{0}(k_{n})$ using Eq. (26). We parametrize the scattering amplitude using the effective-range expansion (ERE),

[TABLE]

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 (27) below threshold, where $-ik>0$. Such a pole occurs for $\cot\delta_{0}(k_{\rm BS})=i$, where $k_{\rm BS}$ is the (imaginary) bound-state momentum. Combining this condition with the ERE then gives

[TABLE]

where terms of $\mathcal{O}(|k_{\rm BS}|^{4})$ are neglected. We solve Eq. (29) for $|k_{\rm BS}|$ and obtain the binding energy via

[TABLE]

This approach was previously utilized in Refs. Beane:2010em ; Mohler:2013rwa ; Lang:2014yfa ; Moir:2016srx . 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^{-|k_{\rm BS}|L}$ (for more details, see, for example, Refs. Sasaki:2006jn ; Beane:2011iw ).

VI.3 Numerical results

In Fig. 6 we show $k\cot\delta_{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\delta_{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_{\rm BS}|$, and the binding energy $E_{\textrm{binding}}$ are given in Table 7. The binding energies are essentially identical to the finite-volume energy differences $\Delta E_{0}$ collected in Table 6. This supports our interpretation of the ground state as a stable tetraquark. We also applied the consistency check proposed in Ref. Iritani:2017rlk , but our statistical uncertainties are too large to reach a definitive conclusion.

As indicated in Fig. 6, the first excited state is close to the inelastic $B^{*}B^{*}$ threshold for some of the data sets, which means that the single-channel analysis performed here is not well justified. A coupled-channel analysis is not feasible with the data we have, but is unlikely to significantly affect the results for the bound state.

VII Fit of the pion-mass dependence and estimates of systematic uncertainties

Given the statistical uncertainties in our results for $E_{\textrm{binding}}$ (shown in Table 7), 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,

[TABLE]

where we use $m_{\pi,\text{phys}}=135\,\textrm{MeV}$ for the physical pion mass in the isospin-symmetric limit. The fit gives

[TABLE]

and has $\chi^{2}/{\rm 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_{\textrm{binding}}(m_{\pi,\text{phys}})$ is nearly identical with the result from the C00078 ensemble with $m_{\pi}=139(1)$ MeV, we conclude that any remaining systematic uncertainties associated with the extrapolation to $m_{\pi,\text{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 Blum:2014tka ; multiplying by the QCD scale of $\Lambda\sim 300\,\textrm{MeV}$ yields $3\,\textrm{MeV}$ to $6\,\textrm{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 $\alpha_{s}^{2}$ in the matching to full QCD, are not included in the action. The analysis of Ref. Blok:1996iz suggests that their effects could be as large as $3\,\textrm{MeV}$.

• •

The matching coefficient $c_{4}$ of the operator $-\frac{g}{2m_{b}}\>{\boldsymbol{\sigma}}\cdot\mathbf{B}$ was computed to one loop. Missing higher-order corrections to this coefficient introduce systematic errors of order

[TABLE]

• •

The matching coefficients of the operators of order $(\Lambda/m_{b})^{2}$ were computed at tree-level only. The missing radiative corrections to these coefficients introduce systematic errors of order

[TABLE]

These estimates are appropriate for $E_{B}$ and $E_{B^{*}}$, which contribute to our calculation of the binding energy via $E_{\rm binding}=E_{\rm BS}-E_{B}-E_{B^{*}}$. For the tetraquark energy $E_{\rm 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 $\Lambda$ in the heavy-light power counting by the binding momentum $|k_{\rm BS}|\sim 800$ MeV, which suggests systematic errors of order 10 MeV. It is likely that there is a partial cancellation of the systematic errors in $E_{\rm binding}=E_{\rm BS}-E_{B}-E_{B^{*}}$. Therefore, we estimate the overall discretization and heavy-quark systematic errors to be not larger than $10\,\textrm{MeV}$ (in future work, the estimates of the heavy-quark errors could be made more precise by numerically investigating the dependence of $E_{\textrm{binding}}$ on the lattice NRQCD matching coefficients). Our final results for the tetraquark binding energy and mass are therefore

[TABLE]

where $m_{\textrm{tetraquark}}$ is obtained by adding the experimental values of the $B$ and $B^{*}$ masses Tanabashi:2018oca to $E_{\textrm{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 a $\bar{b}\bar{b}ud$ bound state that is stable under the strong and electromagnetic interactions.

In Fig. 8 we compare our result (36) for the binding energy with several previous determinations. These include direct lattice QCD calculations that also treated the heavy $\bar{b}$ quarks with NRQCD Francis:2016hui ; Junnarkar:2018twb , calculations in the Born-Oppenheimer approximation using static $\bar{b}\bar{b}$ potentials (in the presence of $u$ and $d$ valence quarks) computed on the lattice Bicudo:2012qt ; Brown:2012tm ; Bicudo:2016ooe , as well as the studies of Refs. Carlson:1987hh ; SilvestreBrac:1993ss ; Brink:1998as ; Vijande:2003ki ; Janc:2004qn ; Vijande:2006jf ; Navarra:2007yw ; Ebert:2007rn ; Zhang:2007mu ; Lee:2009rt ; Karliner:2017qjm ; Eichten:2017ffp ; Wang:2017uld ; Park:2018wjk ; Wang:2018atz ; Liu:2019stu , which are based on quark models, effective field theories, or QCD sum rules.

The calculations using static $\bar{b}\bar{b}$ potentials from lattice QCD Bicudo:2012qt ; Brown:2012tm ; Bicudo:2016ooe consistently give a binding energy that is about a factor of 2 smaller than our result (36), 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 Francis:2016hui ; Junnarkar:2018twb 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. Junnarkar:2018twb agrees with ours, Ref. Francis:2016hui 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. Francis:2016hui . It is interesting to observe that the effective energies shown in Refs. Francis:2016hui ; Junnarkar:2018twb 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 Yamazaki:2017jfh . 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\>{\rm fm}$, the energy difference between the threshold and next scattering state is only around $8\,\textrm{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 Iritani:2016jie ; see Ref. Davoudi:2017ddj for a critical discussion of this issue.

Even though our work has improved upon previous studies of the $\bar{b}\bar{b}ud$ system by including nonlocal meson-meson 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. Francis:2018qch , 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 Peardon:2009gh ; Abdel-Rehim:2017dok for constructing the correlation functions.

Acknowledgements

We thank Antje Peters for collaboration in the early stages of this project, and we thank the RBC and UKQCD collaborations for providing the gauge-field ensembles. L.L. is supported by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics under contract DE-AC05-06OR23177. S.M. is supported by the U.S. Department of Energy, Office of Science, Office of High Energy Physics under Award Number DE-SC0009913, and by the RHIC Physics Fellow Program of the RIKEN BNL Research Center. M.W. acknowledges support by the Heisenberg Programme of the DFG (German Research Foundation), grant WA 3000/3-1. This work was supported in part by the Helmholtz International Center for FAIR within the framework of the LOEWE program launched by the State of Hesse. Calculations on the LOEWE-CSC and on the on the FUCHS-CSC high-performance computer of the Frankfurt University were conducted for this research. We would like to thank HPC-Hessen, funded by the State Ministry of Higher Education, Research and the Arts, for programming advice. This research used resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility operated under Contract No. DE-AC02-05CH11231. This work also used resources at the Texas Advanced Computing Center that are part of the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number ACI-1548562.

Appendix A The two lowest energy levels for all ensembles

In Figs. 9-12 we show the two lowest energy levels for the ensembles C00078, C01, F004, and F006. The style is identical to Fig. 4, where the same energy levels are shown for ensemble C005, and which is discussed in detail in Sec. V.3. The numerical results of the fits for all ensembles are given in Tables 8-12.

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 Luscher:1990ck ; Blossier:2009kd ; Orginos:2015tha . The latter is based on the generalized eigenvalue problem

[TABLE]

where $C_{jk}(t)$ must be a hermitian square matrix. As discussed in Sec. III.1, our use of point-to-all propagators did not allow us to compute the correlation matrix elements $C_{jk}(t)$ with $j,k\in\{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 $\lambda^{n}(t,t_{0})$ are expected to satisfy

[TABLE]

where $E_{n}$ is the energy of the $n$th state. We define the effective energy $E_{\textrm{eff},n}(t)$ as

[TABLE]

which for large $t$ should plateau at $E_{n}$. An example for $E_{\textrm{eff},n}(t)$ is shown in Fig. 13. We determined $E_{n}$ from constant fits to $E_{\textrm{eff},n}(t)$ in a suitable range $t_{\textrm{min}}\leq t\leq t_{\textrm{max}}$ such that $\chi^{2}/\textrm{d.o.f.}\lesssim 1$. Alternatively, one can fit exponential functions $A\textrm{e}^{-E_{n}t}$ to the eigenvalues $\lambda^{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 symmetric correlation matrices: $2\times 2$ correlation matrices with operators $(O_{1},O_{2})$, $(O_{1},O_{3})$ and $(O_{2},O_{3})$ as well as the $3\times 3$ correlation matrix with operators $(O_{1},O_{2},O_{3})$.

• •

We used the same values for $t_{\textrm{min}}$ (larger $t_{\textrm{min}}$ leads to a stronger suppression of excited states; we performed a comparison for several values for $t_{\textrm{min}}$).

• •

We used similar values for $t_{\textrm{max}}$ (the results only weakly depend on $t_{\textrm{max}}$; we chose $t_{\textrm{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,\ldots,6$ (results shown in the following were obtained with $t_{0}/a=3$).

This comparison is shown in Table 13 for the C005 ensemble. It is reassuring that the results obtained from multi-exponential 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 $\pi\pi$ scattering in Ref. Alexandrou:2017mpi , where agreement was also found.