Investigation of Doubly Heavy Tetraquark Systems using Lattice QCD
IInvestigation of Doubly Heavy Tetraquark Systemsusing Lattice QCD
Martin Pflaumer ∗ , Luka Leskovec , , Stefan Meinel , Marc Wagner , Goethe-Universit¨at Frankfurt am Main, Institut f¨ur Theoretische Physik, Max-von-Laue-Straße 1,D-60438 Frankfurt am Main, Germany Thomas Je ff erson National Accelerator Facility, Newport News, VA 23606, USA Department of Physics, Old Dominion University, Norfolk, VA 23529, USA Department of Physics, University of Arizona, Tucson, AZ 85721, USA Helmholtz Research Academy Hesse for FAIR, Campus Riedberg, Max-von-Laue-Straße 12, D-60438Frankfurt am Main, GermanyE-mail: [email protected]
Asia-Pacific Symposium for Lattice Field Theory - APLAT 2020, 4 - 7 August 2020
Abstract.
We search for possibly existent bound states in the heavy-light tetraquark channels with quarkcontent ¯ b ¯ bud , ¯ b ¯ bus and ¯ b ¯ cud using lattice QCD. We carry out calculations on several gauge link ensembleswith N f = +
1. Introduction
One of the major challenges in QCD is to understand exotic hadrons. This is mainly motivated byexperimental observations of mesons whose quantum numbers, masses or decays cannot be explained byordinary quark-anti-quark pairs. A prominent example are the charged Z ± b states with masses and decaychannels suggesting the existence of a b ¯ b pair, whereas the non-vanishing electrical charge indicates thepresence of another light quark-antiquark pair [1]. These systems are theoretically extremely challengingto investigate due to several existing decay channels. Doubly-heavy tetraquark systems, which aresimpler to study have quark content ¯ Q ¯ Q (cid:48) qq (cid:48) , where Q , Q (cid:48) ∈ { b , c } are heavy quarks and q , q (cid:48) ∈ { u , d , s } represent light quarks. Previous studies showed that this system forms a stable bound state in the heavyquark limit m Q → ∞ [2–4]. Moreover, many investigations were carried out within quark models,e ff ective field theories and QCD sum rules for physical b quark mass m Q = m b predicting a hadronicallystable state [4–10]. Recently, Born-Oppenheimer investigations of a four-quark system containing aheavy antidiquark ¯ b ¯ b and a light diquark ud based on lattice QCD four-quark potentials predict ahadronically stable tetraquark in the I ( J P ) = + ) channel [11–14, 37] while in the I ( J P ) = − )channel a resonance has been predicted. [16]. More rigorous full lattice QCD studies confirmed the¯ b ¯ bud bound state and considered further heavy-light four-quark systems [17–23]. Here we report on ourfindings for the ¯ b ¯ bud , I ( J P ) = + ) channel [21] and our progress concerning tetraquark systems withquark content ¯ b ¯ bus and ¯ b ¯ cud . a r X i v : . [ h e p - l a t ] S e p . Lattice Setup All computations were performed using gauge link configurations generated by the RBC and UKQCDcollaboration [24, 25] with 2 + ff ering in the latticespacing a ( ≈ .
083 fm . . . .
114 fm), lattice size (spatial extent ≈ .
65 fm . . . .
48 fm) and pion mass( ≈
139 MeV . . .
431 MeV). In this talk, we focus mainly on results obtained on the C005 ensemblehighlighted in Tab. 1. Smeared point-to all propagators were used for all quarks. The heavy b quarkswere treated in the framework of Non-Relativistic QCD (NRQCD) [31, 32], while the charm quarkpropagators were computed using a relativistic heavy quark action as described in [33]. In order toreduce the numerical cost of the computation, the all-mode-averaging technique was applied [34, 35].Ensemble N s × N t a [fm] am u ; d am s m π [MeV]C00078 48 ×
96 0 . . . ×
64 0 . .
005 0 .
04 340(1)C01 24 ×
64 0 . .
01 0 .
04 431(1)F004 32 ×
64 0 . .
004 0 .
03 303(1)F006 32 ×
64 0 . .
006 0 .
03 360(1)
Table 1.
Gauge-link ensembles [24, 25] used in this work. N s , N t : number of lattice sites in spatial andtemporal directions; a : lattice spacing; am u ; d : bare up and down quark mass; am s : bare strange quarkmass; m π : pion mass.
3. Hadronically stable Tetraquark ¯ b ¯ bud in the I ( J P ) = + ) channel In a first step, we studied the energy spectrum of the doubly-bottom four quark system with two lightquarks ud and quantum numbers I ( J P ) = + ) [21]. The two lowest thresholds in this channel are the BB ∗ and the B ∗ B ∗ meson pairs, which di ff er only by 45 MeV.In order to extract the energy spectrum, we considered two types of interpolating operators. On the onehand, we constructed local operators, where all four quarks are located at the same space-time positionand the total momentum is projected to zero. On the other hand, we used so-called non-local operators,which describe two spatially separated mesons, each with definite momentum.We included three local operators, namely two mesonic ones which create a BB ∗ ( O ) and a B ∗ B ∗ ( O )structure, respectively, and a diquark-antidiquark operator ( O ). Additionally, two non-local operatorsdescribing a BB ∗ ( O ) and B ∗ B ∗ ( O ) scattering state were added to the operator basis. The detailedconstruction of the operators is discussed in [21], but for completeness, we summarize them in Tab. 2using the following notation: T ( Γ , Γ ) = (cid:88) x ¯ Q Γ q ( x ) ¯ Q Γ q ( x ) , T ( Γ , Γ ) = (cid:88) x ¯ Q Γ q ( x ) ¯ Q Γ q ( x ) D ( Γ , Γ ) = (cid:88) x ¯ Q a Γ ¯ Q b ( x ) q a Γ q b ( x ) , D ( Γ , Γ ) = (cid:88) x ¯ Q a Γ ¯ Q b ( x ) q b Γ q a ( x ) M ( Γ ) = (cid:88) x ¯ Q Γ q ( x ) , M ( Γ ) = (cid:88) x ¯ Q Γ q ( x ) , M ( Γ ) = (cid:88) x ¯ Q Γ q ( x ) , M ( Γ ) = (cid:88) x ¯ Q Γ q ( x ) . (1)We expect that the local operators will generate a state that predominantly overlaps with the ground state,i.e. describe the stable four-quark state of interest. Additionally the non-local operators are expected tohave sizable overlap to the first excited state, i.e. a two meson state, which helps to isolate the groundstate from the excitations.uark content¯ Q ¯ Q q q quantumnumbers I ( J P ) local operators non-local operators¯ b ¯ bud + ) T ( γ j , γ ) − T ( γ j , γ ) M ( γ j ) M ( γ ) − M ( γ j ) M ( γ ) (cid:15) i jk (cid:16) T ( γ k , γ j ) − T ( γ k , γ j ) (cid:17) (cid:15) i jk (cid:16) M ( γ k ) M ( γ j ) − M ( γ k ) M ( γ j ) (cid:17) D ( γ j C , C γ ) Table 2.
List of operators that were considered for the ¯ b ¯ bud correlation matrix. C = γ γ denotes thecharge conjugation matrix.To determine the energy spectrum we considered the correlation matrix C jk ( t ) = (cid:104)O j ( t ) O † k (0) (cid:105) where O j and O k represent one of the previously introduced interpolating operators. The correspondingschematic representation of the Wick contractions can be found in Fig. 1. As we were using point-to-all T T T T M M T M M TM M T M M T M M M M M M M M Figure 1.
Schematic representation ofWick contractions for the elements ofthe correlation matrix. M and M represent the B and B ∗ mesons inde-pendently projected to zero momentum,while T represents local four-quark op-erators. The black lines represent b quark propagators and the red lines rep-resent light quark propagators.propagators for the light quarks, we were restricted to correlation matrix elements with local operatorsat the source, which means that the resulting correlation matrix is a 5 × C jk ( t ) ≈ N − (cid:88) n = Z nj Z nk e − E n t (2)to the correlation matrix elements where E n is the n -th energy eigenvalue and Z nj = (cid:104) n |O † j | Ω (cid:105) are theoverlaps of the trial states O j | n (cid:105) and the vacuum | Ω (cid:105) . We present the results for the two lowest energylevels for a large number of di ff erent fits in Fig. 2. One observes that the energy levels become onlystable if the non-local operators are included in the operator basis. Furthermore the energy levels aresignificantly lower compared to cases where non-local operators are absent. This shows the importanceof scattering operators for the present study of the ¯ b ¯ bud system. The ground state energy level issignificantly below the BB ∗ threshold, while the first excited level is close to that threshold. This isa first indication that indeed a hadronically stable tetraquark state exists.Additionally, we can gain certain information about the composition of the energy eigenstates | n (cid:105) byconsidering the overlap factors Z nj . For a given operator index j , the overlap factor Z nj introduced in Eq.(2) indicates the relative importance of an energy eigenstate | n (cid:105) , when the trial state O † j | Ω (cid:105) is expandedin terms of energy eigenstates, O † j | Ω (cid:105) = ∞ (cid:88) n = | n (cid:105)(cid:104) n |O † j | Ω (cid:105) = ∞ (cid:88) n = Z nj | n (cid:105) . (3) − − E − E B − E B ∗ [ M e V ] N = , t / a = − : . N = , t / a = − : . N = , t / a = − : . N = , t / a = − : . N = , t / a = − : . N = , t / a = − : . N = , t / a = − : . N = , t / a = − : . N = , t / a = − : . N = , t / a = − : . N = , t / a = − : . N = , t / a = − : . N = , t / a = − : . N = , t / a = − : . N = , t / a = − : . N = , t / a = − : . N = , t / a = − : . N = , t / a = − : . N = , t / a = − : . N = , t / a = − : . N = , t / a = − : . N = , t / a = − : . N = , t / a = − : . Figure 2.
Results for the lowest two¯ b ¯ bud energy levels relative to the BB ∗ threshold, ∆ E n = E n − E B − E B ∗ ,as determined on ensemble C005 fromseveral di ff erent fits. The five bars beloweach column indicate the interpolatorsused, a filled black box indicates a localoperator included, a filled red box ascattering operator included. Aboveeach column, we give the number ofexponentials, the fit range, and the valueof χ / d.o.f. The shaded horizontal bandscorrespond to our final estimates of ∆ E and ∆ E , obtained from a bootstrapaverage of the subset of fits that areshown with filled symbols.If the overlap factor Z mj for a specific | m (cid:105) is dominant compared to all other Z nj with n (cid:44) m , this mightindicate that the state O † j | Ω (cid:105) is quite similar to | m (cid:105) . In Fig. 3 we present the overlap factors obtained bya 3-exponential fit for all five operators. The overlap factors | ˜ Z nj | = | Z nj | / max m ( | Z mj | ) are normalizedsuch that max m ( | Z mj | ) =
1. From Fig. 3, one can see that the trial states created by the non-local meson- | ˜ Z n j | f o r e x p o n e n t i a l s operator O | i | i | i operator O | i | i | i operator O | i | i | i operator O | i | i | i operator O | i | i | i Figure 3.
The normalizedoverlap factors | ˜ Z nj | as de-termined on ensemble C005indicating the relative contri-butions of the energy eigen-states | n (cid:105) to the trial state O † j | Ω (cid:105) .meson scattering operators O and O have a large overlap to the first excited state | (cid:105) , which supportsour assumption that the first excitation is a two-meson scattering state. Additionally, one can see that thetrial state generated by the diquark-antidiquark operator O has significant overlap to the ground state | (cid:105) . This confirms our interpretation that the ground state represents the hadronically stable tetraquark.The previously computed finite volume energy levels E n can be related to the infinite volumescattering amplitude by applying L¨uscher’s method [36]. Here, we used the two lowest energy levelsshown in Fig. 2 and applied L¨uscher’s method to determine the BB ∗ S wave scattering amplitude. First,the finite volume scattering momenta k n defined by E n = E B + (cid:113) m B , kin + k n − m B , kin + E B ∗ + (cid:113) m B ∗ , kin + k n − m B ∗ , kin , with: m kin = p − [ E ( p ) − E (0)] E ( p ) − E (0)] (4)are related to the infinite-volume phase shifts δ ( k n ) via cot ( δ ( k n )) = Z (1; ( k n L / π ) ) π / k n L (5)here Z is the generalized zeta function [36]. The scattering amplitude is given by T ( k ) = δ ( k ) − i (6)and can be parametrized by the e ff ective range expansion (ERE) k cot δ ( k ) = a + r k + O ( k ) , (7)where the two parameters a and r were determined using the finite volume energy levels ∆ E and ∆ E shown in Fig. 2. We illustrate the ERE for ensemble C005 in Fig. 4. − . − . − . . . ak ) − a k c o t δ ( k ) Ensemble C005
Figure 4.
Plot of the e ff ectiverange expansion for ensemble C005.The red line corresponds to the EREparametrization of ak cot( δ ( k )). Be-low the BB ∗ threshold (which is locatedat k = ak cot δ ( k ) + | ak | (again using the EREparameterization for ak cot δ ( k )), whoselowest zero gives the binding momen-tum. The vertical green line indicatesthe inelastic B ∗ B ∗ threshold.Bound states appear as poles of the scattering amplitude below threshold. Combining this polecondition with the ERE yields − | k BS | = a − r | k BS | , (8)where k BS is the bound state scattering momentum. Solving Eq. (8) for | k BS | , the binding energy isobtained via the NRQCD energy-momentum relation (see Sec. VI. in Ref. [21]). We found that theinfinite volume mass of the bound state determined via the pole of the scattering amplitude is essentiallyidentical to the finite volume ground state energy, which confirms the existence of a hadronically stabletetraquark.Note, that all computations were repeated for all five ensembles discussed in Sec. 2, which allowedus to perform an extrapolation to the physical pion mass. In Fig. 5, we present the final results forall ensembles as well as the extrapolation, for which we assume a quadratic pion mass dependence, E binding ( m π ) = E binding ( m π, phys ) + c ( m π − m π, phys ). Our final results for the tetraquark binding energy andmass are E binding ( m π, phys ) = ( − ± ±
10) MeV , m tetraquark ( m π, phys ) = (10476 ± ±
10) MeV . (9)
4. Search for Bound States in the ¯ b ¯ bus and ¯ b ¯ cud sector As discussed in the previous section, we were able to confirm the bound state ¯ b ¯ bud in the 0(1 + ) chan-nel, which according to the literature is certainly the most promising candidate for a hadronically stabletetraquark. We are currently investigating further candidates, which have either a heavier light quark, i.e.an s quark instead of a d quark, or a lighter heavy quark, i.e. a ¯ c quark instead of a ¯ b quark. .
00 0 .
05 0 .
10 0 .
15 0 .
20 0 . m π [GeV ] − − − − − − − − E b i nd i n g [ M e V ] coarsefine Figure 5.
Extrapolation to the physicalpion mass, indicated by the verticaldashed line.In the case of ¯ b ¯ bus , the most promising channel has again angular momentum J = P = + , sowe consider similar quantum numbers as for the ¯ b ¯ bud system, namely I ( J P ) = / + ). Previous latticeinvestigations of this state predict the ground state significantly below B s B ∗ threshold [17, 20].In the ¯ b ¯ cud sector the picture is less clear. First, there are two channels in which a stable tetraquark couldbe expected, the 0(0 + ) and the 0(1 + ) channel. The first study found some indication for a ¯ b ¯ cud boundstate [18, 19]. This was, however, not confirmed by a more recent work [22].In our study of those tetraquarks, we utilize the same techniques as discussed for ¯ b ¯ bud in Sec. 3.In particular, we consider in all three cases local mesonic and diquark-antidiquark operators as well asmeson-meson scattering operators. The employed interpolating operators are listed in Tab. 3, where weuse again the notation given in Eq. (1). We expect that these operators are su ffi cient to extract the groundstates correctly.quark content¯ Q ¯ Q q q quantumnumbers I ( J P ) local operators non-local operators¯ b ¯ bus / + ) T ( γ , γ j ) − T ( γ , γ j ) M ( γ ) M ( γ j ) − M ( γ ) M ( γ j ) (cid:15) i jk (cid:16) T ( γ j , γ k ) − T ( γ j , γ k ) (cid:17) (cid:15) i jk (cid:16) M ( γ j ) M ( γ k ) − M ( γ j ) M ( γ k ) (cid:17) D ( γ j C , C γ )¯ b ¯ cud + ) T ( γ , γ j ) − T ( γ , γ j ) M ( γ ) M ( γ j ) − M ( γ ) M ( γ j ) T ( γ j , γ ) − T ( γ j , γ ) M ( γ j ) M ( γ ) − M ( γ j ) M ( γ ) D ( γ j C , C γ ) − D ( γ j C , C γ )¯ b ¯ cud + ) T ( γ , γ ) − T ( γ , γ ) M ( γ ) M ( γ ) − M ( γ ) M ( γ ) D ( C γ , C γ ) − D ( C γ , C γ ) Table 3.
List of operators considered for the corresponding correlation matrix. Note that for ¯ b ¯ bus weassume an approximate SU(3) flavor symmetry as strange and light quarks are almost massless comparedto the b quark.For ¯ b ¯ bus we found that certain elements of the corresponding correlation matrix are strongly correlated ifall local operators are included in the operator basis. Thus, we restricted the operator basis by generatinga new set of local interpolating operators via linear combination of the original ones. The contribution ofthe operator O j to the new operator O (cid:48) n is determined by the eigenvector component v nj , where (cid:126)v n is theigenvector obtained by solving a 3 × O (cid:48) | n (cid:105) should overlap most strongly with the ground state, while O (cid:48) | n (cid:105) will be more similar to the first excitedstate. The new operators read consequently O (cid:48) n = (cid:88) j = v nj O j , (10)where the eigenvector components can be found in Tab. 4. The scattering operators were not changed,i.e. we used a 4 × b ¯ bus . For ¯ b ¯ cud , 0(0 + ) and ¯ b ¯ cud , 0(1 + ), we used a 3 × v nj j = j = j = n = + . + . + . n = + . − . + . Table 4.
Eigenvector components calculated via a GEP for the 3 x × b ¯ cud , 0(0 + ) we only show resultsfor the lowest energy level, as the first excitation could not be properly isolated from higher excitations.We found in the ¯ b ¯ bus channel a ground state energy level somewhat below the B ∗ B s threshold, whichindicates a bound state. Moreover, the first excited state seems to be close to the threshold and, thus,might be consistent with a meson-meson scattering state. Averaging di ff erent fits on the C005 ensemble,we generated a crude estimate of the binding energy E ¯ b ¯ bus , binding ( m π =
340 MeV) ≈ ( − ±
40) MeV,which is in reasonable agreement with results from Ref. [17, 20].Assuming that the extracted energies for ¯ b ¯ cud shown in Fig. 6, in particular those, where non-localmeson-meson scattering operators were used, reflect the ground state energies, there is no indicationfor a bound state, neither for 0(0 + ) nor for 0(1 + ). In both cases, the ground state seems to be close tothe threshold, which suggests that the lowest energy level is rather a scattering state than a bound four-quark state. This interpretation is further supported by observing a significant decrease of the lowestenergy level when non-local meson-meson scattering operators were considered in the fits. This indicatesa rather large overlap of the ground state to the scattering trial states. In particular for 0(0 + ) this isexpected from lattice QCD studies based on static-static-light-light potentials, where the 0(0 + ) potentialis significantly less attractive than the 0(1 + ) potential (see Ref. [37] , Eqs. (12b) and (16b)). − − E − E B − E B ∗ s [ M e V ] N = , t / a = − : . N = , t / a = − : . N = , t / a = − : . N = , t / a = − : . N = , t / a = − : . N = , t / a = − : . N = , t / a = − : . N = , t / a = − : . N = , t / a = − : . N = , t / a = − : . p r e li m i n a r y Figure 6.
Preliminary results for the lowest energylevels in the ¯ b ¯ bus , 1 / + ) (top left), ¯ b ¯ cud , 0(0 + )(bottom left) and ¯ b ¯ cud , 0(1 + ) (bottom right) sectorsrelative to the relevant threshold determined onensemble C005 from several di ff erent fits. Thebars below each column have the same meaning asdiscussed in Fig. 2. − − E − E B − E D [ M e V ] N = , t / a = − : . N = , t / a = − : . N = , t / a = − : . N = , t / a = − : . N = , t / a = − : . N = , t / a = − : . N = , t / a = − : . N = , t / a = − : . N = , t / a = − : . N = , t / a = − : . N = , t / a = − : . N = , t / a = − : . N = , t / a = − : . p r e li m i n a r y − E − E D − E B ∗ [ M e V ] N = , t / a = − : . N = , t / a = − : . N = , t / a = − : . N = , t / a = − : . N = , t / a = − : . N = , t / a = − : . N = , t / a = − : . N = , t / a = − : . N = , t / a = − : . N = , t / a = − : . N = , t / a = − : . N = , t / a = − : . N = , t / a = − : . N = , t / a = − : . N = , t / a = − : . N = , t / a = − : . N = , t / a = − : . N = , t / a = − : . N = , t / a = − : . p r e li m i n a r y Acknowledgments
We thank Antje Peters for collaboration in the early stages of this project. We thank the RBC andUKQCD collaborations for providing the gauge field ensembles. L.L. acknowledges support fromthe U.S. Department of Energy, O ffi ce of Science, through contracts DE-SC0019229 and DE-AC05-06OR23177 (JLAB). S.M. is supported by the U.S. Department of Energy, O ffi ce of Science, O ffi ceof High Energy Physics under Award Number DE-SC0009913. M.W. acknowledges support by theHeisenberg Programme of the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)- project number 399217702. Calculations on the GOETHE-HLR and on the FUCHS-CSC high-performance computers of the Frankfurt University were conducted for this research. We wouldlike 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 ScientificComputing Center (NERSC), a U.S. Department of Energy O ffi ce of Science User Facility operatedunder Contract No. DE-AC02-05CH11231. This work also used resources at the Texas AdvancedComputing Center that are part of the Extreme Science and Engineering Discovery Environment(XSEDE), which is supported by National Science Foundation grant number ACI-1548562. References [1] A. Bondar et al. [Belle], “Observation of two charged bottomonium-like resonances in Y(5S) decays,” Phys. Rev. Lett. , 122001 (2012) [arXiv:1110.2251 [hep-ex]].[2] J. Carlson, L. Heller and J. Tjon, “Stability of Dimesons,” Phys. Rev. D , 744 (1988)[3] A. V. Manohar and M. B. Wise, “Exotic QQ ¯ q ¯ q states in QCD,” Nucl. Phys. B , 17-33 (1993) [arXiv:hep-ph / Q i Q j ¯ q k ¯ q l ,” Phys. Rev. Lett. , no.20, 202002 (2017) [arXiv:1707.09575 [hep-ph]].[5] M. Karliner and J. L. Rosner, “Discovery of doubly-charmed Ξ cc baryon implies a stable ( bb ¯ u ¯ d ) tetraquark,” Phys. Rev.Lett. , no.20, 202001 (2017) [arXiv:1707.07666 [hep-ph]].[6] Z. G. Wang, “Analysis of the axialvector doubly heavy tetraquark states with QCD sum rules,” Acta Phys. Polon. B ,1781 (2018) [arXiv:1708.04545 [hep-ph]].[7] W. Park, S. Noh and S. H. Lee, “Masses of the doubly heavy tetraquarks in a constituent quark model,” Acta Phys. Polon.B , 1151-1157 (2019) [arXiv:1809.05257 [nucl-th]].[8] B. Wang, Z. W. Liu and X. Liu, “ ¯ B ( ∗ ) ¯ B ( ∗ ) interactions in chiral e ff ective field theory,” Phys. Rev. D , no.3, 036007 (2019)[arXiv:1812.04457 [hep-ph]].[9] M. Z. Liu, T. W. Wu, M. Pavon Valderrama, J. J. Xie and L. S. Geng, “Heavy-quark spin and flavor symmetry partnersof the X(3872) revisited: What can we learn from the one boson exchange model?,” Phys. Rev. D , no.9, 094018(2019) [arXiv:1902.03044 [hep-ph]].[10] E. Braaten, L. P. He and A. Mohapatra, “Masses of Doubly Heavy Tetraquarks with Error Bars,” [arXiv:2006.08650[hep-ph]].[11] P. Bicudo et al. [European Twisted Mass], “Lattice QCD signal for a bottom-bottom tetraquark,” Phys. Rev. D , no.11,114511 (2013) [arXiv:1209.6274 [hep-ph]].[12] P. Bicudo, K. Cichy, A. Peters and M. Wagner, “BB interactions with static bottom quarks from Lattice QCD,” Phys. Rev.D , no. 3, 034501 (2016) [arXiv:1510.03441 [hep-lat]].[13] Z. S. Brown and K. Orginos, “Tetraquark bound states in the heavy-light heavy-light system,” Phys. Rev. D , 114506(2012) [arXiv:1210.1953 [hep-lat]].[14] P. Bicudo, K. Cichy, A. Peters, B. Wagenbach and M. Wagner, “Evidence for the existence of ud ¯ b ¯ b and the non-existenceof ss ¯ b ¯ b and cc ¯ b ¯ b tetraquarks from lattice QCD,” Phys. Rev. D , no. 1, 014507 (2015) [arXiv:1505.00613 [hep-lat]].[15] P. Bicudo, J. Scheunert and M. Wagner, “Including heavy spin e ff ects in the prediction of a ¯ b ¯ bud tetraquark with latticeQCD potentials,” Phys. Rev. D , no. 3, 034502 (2017) [arXiv:1612.02758 [hep-lat]].[16] P. Bicudo, M. Cardoso, A. Peters, M. Pflaumer and M. Wagner, “ ud ¯ b ¯ b tetraquark resonances with lattice QCD potentialsand the Born-Oppenheimer approximation,” Phys. Rev. D , no. 5, 054510 (2017) [arXiv:1704.02383 [hep-lat]].[17] A. Francis, R. J. Hudspith, R. Lewis and K. Maltman, “Lattice Prediction for Deeply Bound Doubly Heavy Tetraquarks,”Phys. Rev. Lett. , no. 14, 142001 (2017) [arXiv:1607.05214 [hep-lat]].[18] A. Francis, R. J. Hudspith, R. Lewis and K. Maltman, “More on heavy tetraquarks in lattice QCD at almost physical pionmass,” EPJ Web Conf. , 05023 (2018) [arXiv:1711.03380 [hep-lat]].[19] A. Francis, R. J. Hudspith, R. Lewis and K. Maltman, “Evidence for charm-bottom tetraquarks and the mass dependenceof heavy-light tetraquark states from lattice QCD,” Phys. Rev. D , no. 5, 054505 (2019) [arXiv:1810.10550 [hep-lat]].[20] P. Junnarkar, N. Mathur and M. Padmanath, “Study of doubly heavy tetraquarks in Lattice QCD,” Phys. Rev. D , no. 3,034507 (2019) [arXiv:1810.12285 [hep-lat]].[21] L. Leskovec, S. Meinel, M. Pflaumer and M. Wagner, “Lattice QCD investigation of a doubly-bottom ¯ b ¯ bud tetraquarkwith quantum numbers I ( J P ) = + ),” Phys. Rev. D , no. 1, 014503 (2019) [arXiv:1904.04197 [hep-lat]].[22] R. J. Hudspith, B. Colquhoun, A. Francis, R. Lewis and K. Maltman, “A lattice investigation of exotic tetraquarkchannels,” [arXiv:2006.14294 [hep-lat]].[23] P. Mohanta and S. Basak, “Construction of bb ¯ u ¯ d tetraquark states on lattice with NRQCD bottom and HISQ up / downquarks,” [arXiv:2008.11146 [hep-lat]].[24] Y. Aoki et al. [RBC and UKQCD], “Continuum limit physics from 2 + ,074508 (2011) [arXiv:1011.0892 [hep-lat]].[25] T. Blum et al. [RBC and UKQCD], “Domain wall QCD with physical quark masses,” Phys. Rev. D , no. 7, 074505(2016) [arXiv:1411.7017 [hep-lat]].[26] D. B. Kaplan, “A Method for simulating chiral fermions on the lattice,” Phys. Lett. B , 342-347 (1992) [arXiv:hep-lat / , 90-106 (1993) [arXiv:hep-lat / , 54-78 (1995)[arXiv:hep-lat / ff and K. Orginos, “The M¨obius domain wall fermion algorithm,” Comput. Phys. Commun. , 1-19(2017) [arXiv:1206.5214 [hep-lat]].[30] Y. Iwasaki and T. Yoshie, “Renormalization Group Improved Action for SU(3) Lattice Gauge Theory and the StringTension,” Phys. Lett. B , 449-452 (1984)[31] B. A. Thacker and G. P. Lepage, “Heavy quark bound states in lattice QCD,” Phys. Rev. D , 196-208 (1991)[32] G. P. Lepage, L. Magnea, C. Nakhleh, U. Magnea and K. Hornbostel, “Improved nonrelativistic QCD for heavy quarkphysics,” Phys. Rev. D , 4052-4067 (1992) [arXiv:hep-lat / , no.9, 094507 (2014) [arXiv:1409.0497 [hep-lat]].34] T. Blum, T. Izubuchi and E. Shintani, “New class of variance-reduction techniques using lattice symmetries,” Phys. Rev.D , no.9, 094503 (2013) [arXiv:1208.4349 [hep-lat]].[35] E. Shintani, R. Arthur, T. Blum, T. Izubuchi, C. Jung and C. Lehner, “Covariant approximation averaging,” Phys. Rev. D , no.11, 114511 (2015) [arXiv:1402.0244 [hep-lat]].[36] M. L¨uscher, “Two particle states on a torus and their relation to the scattering matrix,” Nucl. Phys. B , 531-578 (1991)[37] P. Bicudo, J. Scheunert and M. Wagner, “Including heavy spin e ff ects in the prediction of a ¯ b ¯ bud tetraquark with latticeQCD potentials,” Phys. Rev. D95