A lattice investigation of exotic tetraquark channels
R. J. Hudspith, B. Colquhoun, A. Francis, R. Lewis, K. Maltman
AA lattice investigation of exotic tetraquark channels
R.J. Hudspith , B. Colquhoun , A. Francis , R. Lewis , and K. Maltman PRISMA + Cluster of Excellence & Institut f¨ur Kernphysik, Johannes Gutenberg-Universit¨at Mainz Department of Physics and Astronomy, York University, Toronto, Ontario, M3J 1P3, Canada Theoretical Physics Department, CERN, CH-1211 Geneva 23, Switzerland and Department of Mathematics and Statistics, York University, Toronto, Ontario M3J1P3, Canada and CSSM, University of Adelaide, Adelaide, SA, 5005, Australia (Dated: June 26, 2020)We perform an n f = 2 + 1 lattice study of a number of channels where pastclaims exist in the literature for the existence of strong-interaction-stable light-heavytetraquarks. We find no evidence for any such deeply-bound states, beyond the J P = 1 + , I = 0 ud ¯ b ¯ b and I = 1 / ls ¯ b ¯ b states already identified in earlier latticestudies. We also describe a number of systematic improvements to our previous latticestudies, including working with larger m π L to better suppress possible finite volumeeffects, employing extended sinks to better control excited-state contamination, andexpanding the number of operators used in the GEVP analyses. Our results also allowus to rule out several phenomenological models which predict significant tetraquarkbinding in channels where no such binding is found. a r X i v : . [ h e p - l a t ] J un I. INTRODUCTION
Results from multiple recent lattice studies [1–5] now rather firmly establish the existenceof an exotic, doubly bottom, I = 0, J P = 1 + ud ¯ b ¯ b tetraquark state, bound with respectto BB ∗ , and hence strong-interaction stable. Though results for the binding energy vary,all correspond to masses below BB threshold, making the state stable with respect to, notonly strong-interaction, but also electromagnetic, decays. The results of Refs. [3, 4] alsopredict a strong-interaction-stable flavour ¯3 F , I = 1 / J P = 1 + strange partner. Refs. [6, 7],using heavy quark symmetry arguments, supplemented by either phenomenological input [6]or phenomenological plus lattice input [7] for leading finite heavy mass corrections, concuron the strong-interaction stability of both the I = 0 and I = 1 / b antiquarks in a colour 3 c and thespin-dependent interaction of light quarks in the ¯3 F , spin 0, colour ¯3 c “good light diquark”configuration, whose attractive character is known phenomenologically from observed heavybaryon splittings [3, 9]. Neither of these interactions is accessible in a state consisting of twowell-separated heavy mesons.While the existence of this ¯3 F of J P = 1 + strong-interaction-stable tetraquark statesseems well established theoretically, detecting these states, and confirming this prediction,remains experimentally challenging. Production rates are likely to be very low, given thattwo b ¯ b pairs must first be produced. While the doubly charmed Ξ cc baryon has now beenobserved by LHCb [10, 11], and double-Υ production has been reported by CMS [12], todate, not even bottom-charm, let alone doubly bottom baryon states have been detected.The recently proposed inclusive search strategy, in which a B c meson originating from adisplaced vertex would signal the presence of doubly bottom hadron production at theLHC [13], though clearly of interest, would still leave open the question of whether sucha signal contained a doubly bottom tetraquark component, or was due entirely to doublybottom baryon production.In this paper we use lattice QCD calculations to investigate whether strong-interaction-stable (hereafter “bound”) tetraquarks exist in other channels likely to be more amenable toexperimental detection. A number of such channels have been investigated in the literatureusing various approaches, including QCD-inspired models [14–42] and QCD sum rules [43–55],and we will study a range of channels where claims for the existence of such bound tetraquarkstates have been made. Our investigations will also serve to test those aspects of the modelsnot constrained by fits to the ordinary meson and baryon spectrum as well as approximationsemployed in implementing the QCD sum rule framework in studies of tetraquark channels.In the case of model studies, the parameters of the models (most typically constituentquark models) are fixed from earlier fits to the ordinary meson and baryon spectrum.Interactions of a constituent quark pair in a colour ¯3 c or a constituent quark and antiquarkin a colour 1 c configuration are, thus, phenomenologically constrained. In tetraquark (andother multi-quark) channels, however, additional colour configurations, which have totallyunconstrained constituent interactions, are also present. In those channels already identifiedby the lattice as supporting deeply bound, doubly bottom tetraquark states, where thephysics of the binding is believed to be dominated by a combination of the colour Coulombattraction in the anti-diquark colour 3 c and the diquark attraction in the good-light-diquark¯3 c configuration, these dominant interactions are phenomenologically constrained, and onewould, therefore, expect the models to successfully predict tetraquark binding in thesechannels. This is, indeed, typically the case. However, even in these channels, differences inthe assumed form of the model interactions, which, by construction, produce no differencesin the ordinary hadron spectrum (since the model parameters have been fit to ensure thespectrum is reproduced) can produce significant differences in an exotic channel. An exampleis provided by the comparison of results for predicted binding in the I = 0, J P = 1 + ud ¯ b ¯ b channel produced by a range of dynamical chiral [21, 23, 25, 26, 29, 30, 36, 41] andnon-chiral[16, 17, 19, 20, 22, 24, 26, 27, 29–32, 37, 40, 42] quark models. The latter typicallyemploy a purely one-gluon-exchange (OGE) form for the spin-spin interaction, while spin-spin interactions for the former are produced by a combination of effective Goldstone-bosonexchange, acting between the light constituents only, and nominal OGE acting between allconstituents. Even if the tetraquark state is entirely dominated by the 3 c anti-diquark, ¯3 c good-light-diquark configuration, this configuration includes both 1 c and 8 c ¯ b -light quarkpairs. The fits of the models to the ordinary meson and baryon spectrum place no constraintson these residual 8 c heavy-light interactions. In addition, the colour dependence of theOGE interaction, and the absence of Goldstone-boson-exchange contributions to the spin-dependent heavy-light interactions in the chiral quark model framework mean one must expectsignificant differences in the residual heavy-light interactions, and hence in the predictionsfor the tetraquark binding, in the two classes of models. This expectation is also borne out inthe literature, where the dynamical non-chiral models of Refs. [16, 17, 19, 20, 22, 24, 27, 29–31, 37, 40, 42] produce binding energies between 54 and 160 MeV, in contrast to the resultsof the dynamical chiral models of Refs. [21, 23, 29, 30, 36, 41], which lie between 214 and497 MeV.Different systematic questions arise for predictions generated using QCD sum rules, wherethe underlying dispersive representations are rigorously valid but practical applications requireapproximations on both the spectral and OPE sides. All tetraquark studies we are aware ofwork with Borel-transformed sum rules and the SVZ spectral ansatz, which consists of a singlelow-lying “pole” (characterized by its mass, M , and coupling, f , to the relevant interpolatingcurrent) and a continuum (approximated using the operator product expansion (OPE)) for s above a “continuum threshold” s [56, 57]. The OPE side is taken as input, and the sumrules used to fix M , f , and s . The OPE, however, involves not just known quantities, likequark masses, α s and (cid:104) ¯ qq (cid:105) , but also unknown higher dimension condensates. Contributionsup to dimensions D = 8 or 10 are retained in recent tetraquark analyses. As stressed in theearliest of the sum rule studies of doubly heavy tetraquarks [43], estimating such higher D contributions (typically in terms of products of known lower-dimension condensates) requiresuse of the factorization approximation. The accuracy of this approximation is not knownin general, but for the I = 1 vector and axial-vector current-current two-point functions,extractions of the two relevant D = 6 four-quark condensates from finite-energy sum ruleanalyses of OPAL [58] and ALEPH [59] hadronic τ decay data found it to be off by as muchas a factor of ∼ I = 0, J P = 0 − channel, which has not one, but threenarrow, low-lying η resonances. Its use in exotic channels, where limited prior knowledgeof the qualitative features of the spectral distribution is likely, thus has the potential toproduce difficult-to-assess systematic uncertainties. One way to investigate this issue is toconsider sum rules for different interpolating operators with the same quantum numbers. Theassociated spectral functions will all receive contributions from all states with these quantumnumbers, with only the contribution strengths operator dependent. Finding compatible polemasses and similar continuum thresholds from all the sum rules would increase confidencein the reliability of the SVZ ansatz. Finding incompatible pole masses would, in contrast,signal either that the approximations used on the OPE sides were unreliable, or that thechannel has more than one narrow state, raising questions about the suitability of the use ofthe SVZ ansatze.The lattice approach is ideally suited to investigating channels with a deeply boundtetraquark ground state, and to testing model and sum rule predictions for the existenceof such states. Provided one employs interpolating operators with reasonable ground stateoverlap, the sizeable gap to the lowest meson-meson threshold increases the likelihood theground state will dominate the corresponding two-point function at moderate Euclideantimes, before the signal is lost in noise. Judicious source and sink choices (discussed in moredetail in Sec. IV below) are also relevant to improving the ground state signal. Weakly boundstates represent more of a challenge because finite volume (FV) effects have the potential toproduce a FV analogue of what will become a meson-meson scattering state in the continuumthat is shifted below the continuum meson-meson threshold. FV effects are expected to besmall for bound states, but not necessarily negligible for continuum scattering states. Theyare typically handled either by extrapolations using simulations at multiple volumes, or byworking with ensembles where the product m π L , with L the lattice length, is large enoughto suppress the dominant “round-the-world” FV effects. A general rule of thumb is that thisrequires m π L >
4. The fact that the good-light-diquark configuration is believed to play animportant role in binding for the tetraquarks so far identified and that, phenomenologically,the associated contributions to binding are expected to grow with decreasing light quarkmass, also puts a premium on working with m π as close to the physical point as possible,subject to maintaining a detectable signal and keeping FV effects under control.For the 32 ×
64 ensembles used in our previous studies of the J P = 1 + , ¯3 F ud ¯ b ¯ b and (cid:96)s ¯ b ¯ b ( (cid:96) = u, d ) channels [3] and J P = 1 + , I = 0 ud ¯ c ¯ b channel [9] the m π L values were6 . m π (cid:39)
415 MeV and 4 . m π (cid:39)
299 MeV, but only 2 . m π (cid:39)
164 MeV. Forthe current study, which focuses on searching for channels with bound, non-molecular (i.e.,non-weakly-bound) tetraquark ground states, we have thus generated a new ensemble withthe same lattice spacing but a larger volume (48 ×
64) and a pion mass, m π (cid:39)
192 MeV,still close to physical, but with m π L = 4 . >
4. With only this single m π , we will be unableto identify channels in which a shallow bound state might exist at physical m π , but not atthe slightly higher m π of our simulation. Channels identified as having a moderate-to-deeplybound tetraquark ground state for m π (cid:39)
192 MeV, will, however, certainly also have such aground state at slightly lower, physical m π .We now turn to the channels to be investigated in this paper. The likelihood of binding isincreased by restricting our attention to those channels where no relative spatial excitationsare required, and where light quark pairs have access to the favourable ¯3 F , spin 0, colour¯3 c , good-light-diquark configuration. The antiquark pair must then have access to the 3 c configuration and, with no relative spatial excitation, have spin 1 if the two antiquarks are thesame. Both spin 0 and 1 are possible if the antiquarks are different. These considerations leadus to investigate channels with either I = 0 or 1 / J P = 1 + and two identical antiquarks, or I = 0 or 1 / J P = 0 + or 1 + and two non-identical antiquarks.In what follows, we will first briefly revisit our previous study [3] of the doubly bottom, J P = 1 + , ¯3 F channels already well established to support bound tetraquark states. Our goalhere is not to re-establish this fact, but rather to illustrate an important improvement tothe methods employed in Ref. [3]. In that analysis (as well as in our analysis of the I = 0, J P = 1 + ud ¯ c ¯ b channel [9]) local sinks were combined with gauge-fixed wall sources. Withthis setup, the resulting ground state effective mass plateaus were short and reached at laterEuclidean times. With Wall-Local effective masses typically plateauing from below, thisleaves open the possibility that the true plateaus had not yet been reached and the resultingbinding energies overestimated. Improving the ground state signal is thus desirable, andwe have succeeded in accomplishing this goal using the extended “Box-Sink” construction,described in Sec. IV below. The results for the bound doubly bottom ¯3 F ground states,at the m π of this study, are thus included primarily to illustrate this improvement, and tomotivate the use of the Box-Sink construction in the other channels considered here, whichare identified and discussed in more detail in Sec. II. For the doubly bottom channels, wherebinding is already clearly established, the ultimate goal is to carry out simulations at notjust the larger volume and single m π of the current study, but also, at this same volume, formultiple m π , in order to provide updated results for the binding, extrapolated to physical m π . The results of this extended update, which is in progress, will be reported elsewhere.The rest of this paper is organized as follows. In Sec. II, we list the additional channels tobe studied and review existing predictions for tetraquark binding in each. Sec. III provides alist of the operators used in these studies and Sec. IV a discussion of our gauge-fixed wallsources, and the extended, Box-Sink construction, introduced to improve our ground statesignals. Sec. V provides details of our lattice setup and calculations, and Sec. VI the resultsof these calculations. Our final conclusions, together with a discussion of the implications ofour results, are given in Sec. VII. II. CHANNELS OF INTEREST AND EXISTING TETRAQUARKPREDICTIONS THEREIN
In this section, we specify the additional channels to be considered in this study, providinga review of existing tetraquark binding predictions for, and the reasons for an interest in,each such choice. The predictions are to be tested against results from the lattice in Sec. VII.The existence of sizeable B c data sets already establishes the existence of significantsimultaneous production of b ¯ b and c ¯ c pairs. A bound tetraquark in one of the two mixed-heavycharm-bottom channels is thus expected to be much more amenable to experimental detectionthan either of the J P = 1 + , ¯3 F doubly bottom analogues discussed above. Replacing one ofthe two ¯ b antiquarks by a ¯ c , however, is expected to both reduce the Coulomb attraction ofthe 3 c heavy antiquark pair and increase the residual spin-dependent heavy-light interactions,which heavy baryon spectrum splittings suggest is likely to further reduce binding. Thebottom-charm states (if any) are thus qualitatively expected to be more weakly bound thantheir doubly bottom counterparts [9].The mixed-heavy I = 0, J P = 1 + ud ¯ c ¯ b channel was investigated previously on thelattice [9]. Evidence of a possible bound ground state with binding between ∼
15 and 61MeV was found, though with rather short, late-time plateaus. The heavy quark symmetryarguments of Refs. [6, 7], which produce binding compatible with that found from the latticefor the doubly bottom I = 0, J P = 1 + ground state, in contrast, predict no binding. Anumber of non-chiral models [17, 19, 20, 34, 37, 42] also find an either unbound or onlyweakly bound ground state, with binding energies between − ±
12 MeV, 217 MeVand 199 MeV, respectively. QCD sum rules, in contrast, produce inconclusive results, withbinding of 60 ±
100 MeV reported in Ref. [45].The J P = 0 + , I = 0, ud ¯ c ¯ b analogue has not yet been investigated on the lattice. Theheavy quark symmetry arguments of Refs. [6, 7] again predict no binding in this channel.Non-chiral models [18–20, 33, 34, 37, 42] predict an either unbound or only weakly boundground state, with Refs.[20, 33] and [34] quoting binding energies of between −
11 MeV(unbound) and 13 MeV (depending on the details of the potential used), 11 ±
13 MeV, and 23MeV, respectively. Significantly larger binding is found in most recent chiral model studies,with Refs. [36, 38, 41], for example, reporting bindings of 136 ±
12 MeV, 196 MeV and 178MeV, respectively. Two recent QCD sum rule studies differ in their pole mass results, withRef. [45], which includes OPE contributions up to D = 8, finding 7 . ± .
10 GeV, whileRef. [48], which includes contributions up to D = 10, finds 6660 ±
150 MeV. The lattercorresponds to a very deep, 485 ±
150 MeV, binding relative to DB threshold.The J P = 0 + and 1 + (cid:96)s ¯ c ¯ b channels have yet to be studied on the lattice. Analogous J P = 1 + , (cid:96)s ¯ b (cid:48) ¯ b states, with variable ¯ b (cid:48) mass as low as ∼ . m b , have, however, been considered,for m π = 299 MeV only, in Ref. [9]. With both ¯ b and ¯ b (cid:48) treated using NRQCD, this studycould not be extended to ¯ b (cid:48) masses as low as m c ; the observed variable ¯ b (cid:48) mass dependence,however, was argued to make a bound (cid:96)s ¯ c ¯ b state extremely unlikely in the J P = 1 + channelat this pion mass. Both the heavy quark symmetry arguments of Refs. [6, 7] and non-chiralmodel of Ref. [42] predict no tetraquark bound state in either this channel or the J P = 0 + analogue. In contrast, the QCD sum rule study of Ref. [54], which includes contributions upto D = 8, finds binding energies of 200 ±
130 and 180 ±
110 MeV for the J P = 0 + and 1 + ground states, respectively. Even deeper J P = 0 + bindings, exceeding 400 MeV, are reportedin the QCD sum rule studies of Refs. [49] and [53], which include OPE contributions up to D = 10.The reported observation by the D0 collaboration [62, 63] of a narrow state, X (5568),decaying to B s π ± , and hence having four distinct flavours, prompted speculation that therelated singly-heavy ud ¯ s ¯ b and ud ¯ s ¯ c channels might support bound tetraquark states. Theinitial argument was based on the expectation that the putative X (5568) should have a U -spin-interchanged SU (3) ud ¯ s ¯ b partner with similar mass, coupled with the observation thatthe threshold, 5773 MeV, for the lowest-lying, two-meson ud ¯ s ¯ b state, BK , lies well above thereported X (5568) mass. While this initial motivation is weakened by the fact that searchesby the LHCb [64], CMS [65], CDF [66] and ATLAS [67] collaborations failed to confirm theD0 result, a number of model and QCD sum rule studies exist predicting bound tetraquarkground states in the I = 0, J P = 0 + and 1 + , ud ¯ s ¯ b and/or ud ¯ s ¯ c channels. The quark colourdelocalization screening model, for example, predicts J P = 0 + and 1 + ud ¯ s ¯ b ground statesbound by 74 and 58 MeV relative to the respective two-meson, BK and B ∗ K , thresholds [39].Another SU (3) chiral quark model study, that of Ref. [35], obtains similar bindings, 70 and 68MeV, for the J P = 0 + and 1 + bottom-strange states. Somewhat smaller bindings, 19 and 16MeV, respectively, are found in the alternate chiral quark model study of Ref. [41]. A boundtetraquark with even lower mass, 5380 ±
170 MeV (corresponding to a binding of 394 ± BK threshold), is also found for the J P = 0 + bottom-strange channel inthe QCD sum rule study Ref. [51]. In fact, the only bottom-strange sector investigationwe are aware of which does not produce a bound tetraquark state is the non-chiral-model, I = 0, J P = 1 + channel study of Ref. [17]. This study also found no bound tetraquark inthe I = 0, J P = 1 + ud ¯ s ¯ c channel. An absence of binding was also found for the related I = 0, J P = 0 + ud ¯ s ¯ c channel in the sum rule study of Ref. [51]. To our knowledge, theonly other ud ¯ s ¯ c channel study is that of Ref. [41], which finds I = 0, J P = 0 + and 1 + statesbound by 15 and 9 MeV, respectively. Any bound state found in these channels would be ofconsiderable phenomenological interest, since it would involve light degrees of freedom in thenovel strangeness +1, colour 3 c configuration about which no phenomenological informationis currently known.While it might appear natural to include the I = 0, J P = 1 + ud ¯ c ¯ c and I = 1 / J P = 1 + (cid:96)s ¯ c ¯ c channels in our study, we have not done so, since evidence already exists that bounddoubly charmed tetraquark states, even if they exist in these channels, will be at most weaklybound. This evidence comes from both the lattice and the heavy-quark-symmetry argumentsof Refs. [6, 7], whose results are compatible with those of lattice studies for the doublybottom bound states. The latter predicts no binding in either of the doubly charmed J P = 1 + channels. On the lattice, in the I = 0, J P = 1 + ud ¯ c ¯ c channel, (i) the results of Ref. [68]clearly establish the absence of binding at m π = 391 MeV, (ii) Ref. [4] (which reaches a lower,257 MeV, pion mass for the coarsest of the three lattices studied) finds an extrapolated,physical- m π binding of 23 ±
11 MeV, sufficiently small for the authors to comment that furtherFV studies are required to reach a firm conclusion concerning binding, and (iii) results for the b (cid:48) mass dependence of the binding energies for I = 0, J P = 1 + ud ¯ b (cid:48) ¯ b (cid:48) states at m π = 299 MeV,reported in Ref. [9], strongly disfavour the possibility that binding survives to physical m π and b (cid:48) masses as low as m c . Only one lattice result, whose binding, 8 ± I = 1 / J P = 1 + (cid:96)s ¯ c ¯ c channel [4]. QCDsum rule studies also find no evidence for binding in either of the non-strange [43, 44, 46, 55],0or strange [44, 46, 55] doubly charmed J P = 1 + channels. Non-chiral models whose I = 0, J P = 1 + ud ¯ b ¯ b bindings are compatible with those found on the lattice, similarly, predicteither unbound [17, 22, 27, 29, 33, 36, 37], or only weakly bound [24, 26, 30], doubly charmed J P = 1 + ground states. While many chiral quark model studies predict significant bindingin the I = 0, J P = 1 + ud ¯ c ¯ c channel [21, 23, 26, 29, 30, 38, 41], the models underlying thesepredictions typically also over-bind the I = 0, J P = 1 + , ud ¯ b ¯ b ground state relative to resultsknown from the lattice. We conclude that the absence of deeply bound tetraquark states inthe two doubly charmed channels is already established, and hence have not included thesechannels in the current study. Additional, larger volume ensembles would be required toreliably investigate such channels, where at-best-weak binding makes FV studies mandatory.In future, we plan to generate additional ensembles with similarly low m π , but larger volumethan considered here, and will revisit the doubly charmed J P = 1 + channels when thesebecome available.Two final channels are studied in this paper. These are the triply-heavy J P = 1 + uc ¯ b ¯ b and sc ¯ b ¯ b channels, considered previously only in the lattice study of Ref. [4], which foundground state masses compatible with the corresponding lowest two-meson thresholds. Thesechannels are, of course, not among those where a bound tetraquark state, if it existed, wouldbe more amenable to current experimental detection than the theoretically established doublybottom states. They are included here only to provide an additional, independent check ofthe results of Ref. [4].We close this section by stressing that deeply bound tetraquark ground states are predictedby some models (especially chiral quark models) and/or some QCD sum rule analyses in allof the I = 0, J P = 0 + and 1 + ud ¯ c ¯ b , I = 1 / J P = 0 + (cid:96)s ¯ c ¯ b , and I = 0, J P = 0 + and 1 + ud ¯ s ¯ b channels. Examples where predicted binding exceeds 100 MeV can be found in (i) the recentsum rule studies of the I = 0, J P = 0 + ud ¯ c ¯ b , I = 1 / J P = 0 + (cid:96)s ¯ c ¯ b , and I = 0, J P = 0 + ud ¯ s ¯ b channels, which produce ground states 485 ±
150 MeV [48], 407 ±
160 MeV [49] or467 ±
150 MeV [53], and 397 ±
170 MeV [52], respectively, below the corresponding lowesttwo-meson thresholds and (ii) recent chiral quark model studies of the ud ¯ c ¯ b sector, whichproduce an I = 0, J P = 0 + ground state 136 ±
12 MeV [36], 196 MeV [38], or 178 MeV [41]below BD threshold and an I = 0, J P = 1 + ground state 171 ±
12 MeV [36], 217 MeV [38]or 199 MeV [41] below B ∗ D threshold. States bound this deeply should be readily amenableto detection in the multi-operator lattice analyses detailed below.1It is also worth emphasizing that chiral and non-chiral quark models, all of which admitparametrizations which successfully reproduce the ordinary meson and baryon spectra, makequalitatively different predictions for the tetraquark channels discussed above. The non-chiralmodels, once one moves beyond the doubly bottom sector, predict only a few, usually at-most-weakly-bound tetraquark candidates. The majority of chiral model studies, in contrast,predict a larger number of much more deeply bound states. The lattice studies below shouldthus allow us to rule out at least one of these classes of models as providing an acceptablephenomenological representation of QCD in the low-energy regime. III. OPERATORS
In this paper we work with the following set of operators having couplings to stateswith two quarks and two antiquarks with quark flavours ψ, φ, θ, and ω , and having either“meson-meson” or “diquark-antidiquark” spin-colour-flavour structure: D (Γ , Γ ) = ( ψ Ta C Γ φ b )(¯ θ a C Γ ¯ ω Tb ) ,E (Γ , Γ ) = ( ψ Ta C Γ φ b )(¯ θ a C Γ ¯ ω Tb − ¯ θ b C Γ ¯ ω Ta ) ,M (Γ , Γ ) = (¯ θ Γ ψ )(¯ ω Γ φ ) , N (Γ , Γ ) = (¯ θ Γ φ )(¯ ω Γ ψ ) ,O (Γ , Γ ) = (¯ ω Γ ψ )(¯ θ Γ φ ) , P (Γ , Γ ) = (¯ ω Γ φ )(¯ θ Γ ψ ) . (1)The Dirac structures Γ and Γ relevant to the channels we consider are specified below.The short-hand terminologies “meson-meson” and “diquark-antidiquark” are meant only tocompactly characterize the spin-colour-flavour structures, and should not be over-interpreted;there is nothing, for example, to prevent one of the “meson-meson” operators from couplingto a state with a spatially extended diquark-antidiquark structure or one of the “diquark-antidiquark” operators from coupling to a meson-meson scattering state. An illustration isprovided by the deeply bound ground state in the doubly bottom, I = 0, J P = 1 + channel,which was found to couple strongly to both the “meson-meson” and “diquark-antidiquark”operators considered in Ref. [3].The lattice explicitly breaks rotational symmetry, breaking the continuum rotation groupdown to the little group O h . The operators listed above lie in either the A or T irrep, whichmap directly to the continuum J = 0 and J = 1 quantum numbers respectively. For theoperators that are produced from the product T (cid:78) T = A (cid:76) E (cid:76) T (cid:76) T [69] we pick2 Type ( ψφθω ) I ( J ) P Diquark-Antidiquark Dimeson udcb/udsb/udsc + D ( γ , γ i ) , D ( γ t γ , γ i γ t ) M ( γ , γ i ) − N ( γ , γ i ) M ( I, γ i γ ) − N ( I, γ i γ ) E ( γ , γ i ) , E ( γ t γ , γ i γ t ) O ( γ , γ i ) − P ( γ , γ i ) O ( I, γ i γ ) − P ( I, γ i γ ) (cid:15) ijk M ( γ j , γ k ) udbb + D ( γ , γ i ) , D ( γ t γ , γ i γ t ) M ( γ , γ i ) − N ( γ , γ i ) M ( I, γ i γ ) − N ( I, γ i γ ) lsbb/ucbb/scbb (1) + D ( γ , γ i ) , D ( γ t γ , γ i γ t ) M ( γ , γ i ) , M ( I, γ i γ ) N ( γ , γ i ) , N ( I, γ i γ ) (cid:15) ijk M ( γ j , γ k ) uscb (1) + D ( γ , γ i ) , D ( γ t γ , γ i γ t ) M ( γ , γ i ) , M ( I, γ i γ ) N ( γ , γ i ) , N ( I, γ i γ ) E ( γ , γ i ) , E ( γ t γ , γ i γ t ) O ( γ , γ i ) , O ( I, γ i γ ) (cid:15) ijk M ( γ j , γ k ) TABLE I: J P = 1 + operators used in this work. Type ( ψφθω ) I ( J ) P Diquark-Antidiquark Dimeson udcb/udsb/udsc + E ( γ , γ ) , E ( γ t γ , γ t γ ) M ( γ , γ ) − N ( γ , γ ) M ( I, I ) − N ( I, I ) M ( γ i , γ i ) uscb (0) + E ( γ , γ ) , E ( γ t γ , γ t γ ) M ( γ , γ ) , M ( I, I ) N ( γ , γ ) , N ( I, I ) M ( γ i , γ i ) TABLE II: J P = 0 + operators used in this work.out the scalar and vector components respectively [5].The operator combinations which in principle couple to states in the J P = 0 + and 1 + channels considered below are listed in Table I. While the diquark-antidiquark operators, D and E , were included in initial explorations of the ud ¯ s ¯ c and ud ¯ s ¯ b channels, they were3found to couple only weakly to the corresponding ground states and hence not included inthe final analyses of these channels. They were also omitted from the final version of theanalyses of the ud ¯ c ¯ b channels where using the set of meson-meson operators was found togive more statistically precise results and produce no significant change to the energies ofthe lowest-lying states.Throughout this work we will obtain our ground states from a generalised eigenvalueproblem (GEVP) approach [70–72]. We use the solutions of the GEVP (sample by sample)to extract the “optimised correlator” C i ( t ) = (cid:88) j,k V ij ( τ ) † C jk ( t ) V ki ( τ ) (2)where V is the matrix with columns made of the eigenvector solutions of C ij ( t ) v j ( t ) = λ i C ij ( t + t ) v j ( t ) . (3) τ (the “diagonalization time” [73]) in Eq. (2) is to be chosen large enough to produceimproved projection onto the ground state. We then perform a correlated single-exponentialfit to the resulting optimised correlators to obtain our final levels. In the rest of this work wewill quote results obtained using t /a = 2 and τ /a = 4, values for which the ground statemasses are found to display good stability upon variation of these parameters. Larger valuesof t and τ tend to make the solution statistically less precise, and eventually unstable. IV. WALL SOURCES AND THE BOX-SINK CONSTRUCTION
Throughout this work we will use Coulomb-gauge-fixed wall sources [74, 75] (fixed tohigh-precision using the FACG algorithm of [76]). These sources have some benefits, as wellas some peculiarities, which are discussed below.A gauge-fixed wall source is a time-slice of point sources on a Coulomb gauge-fixedbackground. Here we assume the gauge condition has been applied over the whole volume.The elements of the source do not need to be stochastically drawn as the gauge fixing conditionconnects colors appropriately; this means that these sources can be used for baryons [77–79].Gauge-fixed wall source results can be considerably more statistically precise, especially atlighter pion masses, than point sources at the same cost due to the volume-sampling. Likeother wall sources, these automatically zero-momentum project, and if momentum is needed4it must be put into the source explicitly as a momentum source or introduced via partialtwisting [80, 81], as was done for the NRQCD tuning (see App.A).Use of a guage-fixed wall source is much like that of a point source, though with somepotential complications for channels like those studied here where two-meson thresholds exist.When combined with a local sink, one typically finds a negative sign for the amplitude of(at least) the first excited state, and hence a “Wall-Local” correlator whose effective massapproaches the ground state mass from below. If one cannot measure the signal to largeenough t to ensure the ground state plateau has been reached, the ground state mass maythen be under-estimated, and lead one to either conclude that a bound state exists whenin fact one does not or over-estimate the binding in the case one does exist. This issuedoes not arise for point or stochastic wall sources, where the effective masses approach theground-state mass from above and so offer a rigorous upper bound on that mass.Often it useful to compute “Wall-Wall” correlators [82], defined by contracting objects inthe usual way, with propagators that are summed over the spatial sites of the sink time-slice, S W ( t ) = (cid:88) x S ( x, t ) . (4)The resulting “Wall-Wall” correlation functions are symmetric under the exchange of operatorsat the source and sink and have effective masses which approach the ground state fromabove. They are, however, very statistically noisy. Since Wall-Local correlator effectivemasses empirically approach their ground state plateaus from below (i.e., have negativeexcited state contamination) while those of Wall-Wall correlators approach theirs from above(i.e., have positive excited state contamination), sink combinations intermediate betweenthe two are expected to exist with small excited state contamination. The goal is to find, ifpossible, combinations of this type which also suffer from only a limited loss in statisticalprecision. One possibility would be to use some form of sink smearing, though the productsof link matrices needed with the propagator solution would make that computationally costly.Another would be to simply take advantage of the gauge condition and construct propagatorsby summing over points lying inside a sphere around each reference sink point x , S B ( x, t ) = 1 N (cid:88) r ≤ R S ( x + r, t ) . (5)With this construction, varying R between 0 and its maximum value, 3( L/ , producesa continuous interpolation between the Wall-Local and Wall-Wall cases. We refer to sinks5 t / a a ( m e ff + ∆ ) R = 1728 ( Wall - Wall )R = 49 ( Wall - Box )R = 20 ( Wall - Box )R = 0 ( Wall - Local ) t / a a ( m e ff + ∆ ) R = 49 ( Wall - Box )R = 20 ( Wall - Box )R = 0 ( Wall - Local ) FIG. 1: An illustrative comparison of Box-Sink B c meson correlation function results fordifferent Box-Sink radii. R is in lattice units. The right panel shows a zoom into theplateau region.constructed in this manner as “Box-Sinks”. As the Box-Sink sum can be performed atthe level of contractions and doesn’t require extra inversions, the construction of Wall-Boxcorrelators is reasonably cheap to implement. In channels where the ground state is a compacthadronic object, one would expect an optimal Wall-Box improvement of the ground stateplateau for Box-Sink radii roughly matching the physical ground-state hadron size. Operators Flavor R Flavor R Flavor R Single Mesons ¯ ud su cu cs
49 ¯ bu
49 ¯ bs bc ud ¯ s ¯ c ud ¯ s ¯ b ud ¯ c ¯ b us ¯ c ¯ b uc ¯ b ¯ b sc ¯ b ¯ b TABLE III: Values of the Box-Sink summation radii (in lattice units) used in this work.Values were chosen so the corresponding effective mass plateaus begin reasonably early.One should bear in mind that, in contrast to what happens in the Wall-Wall case, Wall-Boxand Wall-Local correlation matrices are not symmetric [83]. The GEVP method, however,doesn’t necessarily need this symmetry to obtain the underlying real eigenvalues. Carefulmonitoring of the imaginary parts of the eigenvalues is needed and in practice the lack of6symmetry was not found to be an issue.An illustration of the utility of the Box-Sink construction is provided, for the case of the B c meson channel, in Fig. 1. We see that intermediate values of R , indeed, exist whichextend the ground state plateau to significantly lower t than is the case for the Wall-Local orWall-Wall combinations. Results for all Box-Sink radii must, of course, approach the sameground state mass at sufficiently large t . The point here is that the improvement in theground state plateau occurs for small enough R that the onset of the significant increasein noise that occurs as R grows towards the Wall-Wall limit is avoided in the optimal R plateau region. The figure also illustrates that a simple Wall-Local correlator can havesignificant excited state contamination and only approach the plateau at very large times.This is less of a problem in single-meson channels, like the B c channel, where one can followthe signal out to large enough t , but becomes an issue in channels with larger signal-to-noiseproblems, such as the tetraquark channels we study here, where we are typically only able tofollow our signal out to t/a ≈
20 due to an exponential signal to noise problem. Such issueshave been discussed in the literature for a long time in the context of the nucleon and ρ meson, see, e.g., Refs. [77, 78, 84].A further illustration of the utility of the Box-Sink construction, now for tetraquarksignals, is provided by the effective mass plots of the optimised ground-state correlators forthe J P = 1 + , I = 0 ud ¯ b ¯ b and I = 1 / (cid:96)s ¯ b ¯ b channels at various Box-Sink radii R , shownin Fig. 2. While the errors grow, as anticipated, as R increases, the ground state signalsplateau much earlier, and in the region of much lower noise, than was the case for thecorresponding Wall-Local plateaus shown in Ref. [3]. This is especially true of the optimised R region, R (cid:39)
36 or 49. While the signals degrade quickly for larger t/a (and larger R ), itis clear that a “window of consistency” exists in the region of the Box-Sink plateaus, for bothchannels. It is also clear that, while the Wall-Local ( R = 0) results reach the same plateaus,they do so only at much later t/a , rather close to the region where the signals are lost.The fact that enhanced excited-state contamination is seen in the Wall-Local results ofboth the B c and tetraquark channels suggests this behavior may be common for Wall-Localcorrelators. Fig. 2 illustrates, in more detail, the potential danger this creates. In the absenceof the longer, more reliable Box-Sink plateaus, it would be easy, with only the Wall-Localresults, to mistakenly identify the start of the Wall-Local plateaus as occurring around t/a = 14, leading to an overestimate of the binding in both channels. The utility of the7 t/a a ( m e ff + ∆ ) BB*R = 49R = 36R = 25R = 16R = 0 (a) ud ¯ b ¯ b t/a a ( m e ff + ∆ ) BB s *R = 49R = 36R = 25R = 16R = 0 (b) (cid:96)s ¯ b ¯ b FIG. 2: Box-Sink effective mass plots for the lowest GEVP state in the J P = 1 + , I = 0 ud ¯ b ¯ b and I = 1 / (cid:96)s ¯ b ¯ b channels. Also shown is the lowest-lying two-meson threshold at rest.Box-Sink construction in avoiding this problem is clear. Given the results shown in Fig. 2,we expect our earlier Wall-Local-correlator-based estimates for the physical point bindings inthese channels [3] to be reduced once the ongoing improved, multiple- m π Box-Sink analysisis completed. For now, the results for these channels, at the single m π of the current study,serve to establish the Box-Sink construction as an improved method for identifying channelssupporting deeply-bound tetraquark states and quantifying the binding therein, and motivateour use of it in the studies of the other channels, reported below. V. LATTICE DETAILS
For all channels investigated in this work we have used a single n f = 2 + 1, 48 × β and C SW , thelattice-spacing (for which we use a − = 2 . η b and Υ dispersionrelations within a lattice irrep than is possible with only Fourier modes, and thus improvethe statistical resolution. Details of our action and of the new tuning procedure we haveused are provided in App.A Almost 1000 u, s, c and b propagators have been generated foreach of the tetraquark channels considered here.The new ensemble has some attractive features. In particular, its pion mass, ≈
192 M eV ,is reasonably close to physical while still maintaining m π L >
4. As for the previous PACS-CSensembles, we have a strange sea quark heavier than physical and thus have to use a partially-quenched strange (using the κ s -value of [95]) to get a respectable (507 MeV) kaon mass. TheTsukuba tuning also gives close to physical D and D ∗ meson masses. The lattice-valuedmasses we obtain are presented in Table VI in App. C. Note that, for NRQCD, where thereis an additive mass renormalisation, ∆, only mass differences have physical meaning. Wewill present results in this lattice-valued form for the majority of this work.As we use NRQCD for the b quarks, channels with one or more b quarks do not havea rigorous continuum limit. We do, however, include terms that should cancel up to O ( a )discretisation effects in our NRQCD Hamiltonian [96]. The c quarks have been tuned to havea flat dispersion relation and so are expected to have a very mild approach to the continuum.The light-quarks are O ( a )-improved non-perturbatively with the C SW term. In all, we expectdiscretisation effects to be relatively small.Table IV provides a summary of the new ensemble parameters, the number of configura-tions, N conf , and the number of sources per configuration, N src , used in this study.All vector correlators were averaged over the spatial index i , and all correlators weresymmetrised with their backward-propagating states in order to improve statistical resolution.9 V κ l κ sea s κ P Qs N conf N src m π L ×
64 0 . . . TABLE IV: Details of the ensemble used in this work
VI. RESULTS
In this section, we present the measured spectra (GEVP eigenvalues) and associatedtwo-meson thresholds for the channels discussed above. The eigenvalues are obtained fromcorrelated, single-exponential fits to the corresponding optimised correlators. The windowsof t used in these fits are tabulated in App. C and the corresponding effective-mass plots,showing results for the ground states, as well as those excited states for which some signalcould be resolved, collected in App. D. Reliable, early-onset ground-state plateaus are evidentin all of the channels considered. Additional low-lying excited states are also typically welldetermined. Results for a number of higher-lying, less well determined excited states, obtainedfrom fits to signals lost at relatively early times and/or with relatively short fit windows arealso included in the tables and spectrum plots. These serve to identify, qualitatively, theparts of higher-lying two-meson spectrum which couple most strongly to the operators underconsideration. + + a ( m + ∆ ) DB DB*BD*D*B*D *(2300)B I (a) ud ¯ c ¯ b + + a ( m + ∆ ) DB s BD s DB s *D s B*BD s *B s D*D*B s *B*D s * (b) (cid:96)s ¯ c ¯ b FIG. 3: Lowest-lying GEVP eigenvalues and associated two-meson thresholds for the I = 0, J P = 0 + and 1 + ud ¯ c ¯ b (left) and J P = 0 + and 1 + (cid:96)s ¯ c ¯ b (right) channels.0Fig. 3a shows the spectrum results for the I = 0, J P = 0 + and 1 + ud ¯ c ¯ b channels. Theground states in both cases are found to be roughly consistent with the lowest two-mesonthreshold, with no further states near threshold. There is thus no sign of an extra, boundtetraquark state in either channel.The results for the J P = 0 + and 1 + (cid:96)s ¯ c ¯ b channels are shown in Fig. 3b. Here we see theground states in both channels pushed marginally below the corresponding lowest two-mesonthresholds, DB s and DB ∗ s , respectively. With no other states in the vicinity, however, wesuspect this reflects either a residual systematic (from either the GEVP procedure or ourfitting to the optimised ground-state correlators) or a small ( <
10 MeV) FV shift, and thatthe ground states in these channels correspond to the lowest two-meson thresholds. Withonly a single volume, the possibility of weakly bound ground states cannot, however, beunambiguously ruled out. Interestingly, the three highest excited states in the J P = 1 + channel all lie just below higher-lying two-meson thresholds, with splittings consistent withthose between these thresholds. This “pushing down” of the eigenvalues is of the order of afew MeV at most. This pattern is not, however, observed for the J P = 0 + channel, wherethe third excited state lies much higher than the third excited two-meson threshold, D ∗ B ∗ s . + + a ( m + ∆ ) KB KB*BK *(892)K *(892)B*K *(700)B I (a) ud ¯ s ¯ b + + a m KD KD*DK *(892)K *(892)D*K *(700)D *(2300) (b) ud ¯ s ¯ c FIG. 4: Lowest-lying GEVP eigenvalues and associated two-meson thresholds for the I = 0, J P = 0 + and 1 + ud ¯ s ¯ b (left) and ud ¯ s ¯ c (right) channelsFig. 4a shows our results for the I = 0, J P = 0 + and 1 + ud ¯ s ¯ b channels. The ground statesagain come out consistent with the corresponding lowest two-meson thresholds, with no signof an extra state below, or even near, these thresholds, suggesting the observed ground states1 + a ( m + ∆ ) B c B*BB c *B*B c * (a) uc ¯ b ¯ b + a ( m + ∆ ) B c B s *B s B c *B s *B c * (b) sc ¯ b ¯ b FIG. 5: Lowest-lying GEVP eigenvalues and associated two-meson thresholds for the J P = 1 + uc ¯ b ¯ b and sc ¯ b ¯ b channels.once more correspond to the two-meson threshold states. Results for the analogous I = 0, J P = 0 + and 1 + ud ¯ s ¯ c channels are shown in Fig. 4b. The ground states lie slightly above thelowest two-meson thresholds, with, once more, no sign of an additional, bound state belowthreshold in either channel.Finally, in Fig. 5 we show our results for the triply-heavy J P = 1 + , uc ¯ b ¯ b and sc ¯ b ¯ b channels.The ground states are again consistent with the corresponding two-meson thresholds ( BB ∗ c and B s B ∗ c , respectively), confirming the earlier results for these channels reported in Ref. [4]. VII. DISCUSSION AND CONCLUSIONS
In this work we have used lattice simulations to investigate a number of exotic tetraquarkchannels in which model and/or QCD sum rule studies and/or heavy-quark symmetryarguments predict bound, strong-interaction-stable tetraquark states to exist. We havealso revisited our earlier lattice studies of the doubly bottom, J P = 1 + , ¯3 F ud ¯ b ¯ b and (cid:96)s ¯ b ¯ b and bottom-charm, J P = 1 + , I = 0 ud ¯ c ¯ b channels [3, 9], introducing a number ofsignificant improvements over those earlier works. In particular, we have (i) generated anew ensemble with close-to-physical ( (cid:39)
192 MeV) pion mass and significantly larger volume( m π L = 4 . >
4) to better control possible finite volume effects; (ii) added additionaloperators to improve our ability to resolve the tower of states contributing to each correlator,2and (iii) made major improvements to our ground-state signals via the Box-Sink construction.We have also, as described in more detail in App. A, re-tuned our NRQCD bare mass to amuch more precise and reliable value, using partially-twisted boundary conditions. All ofthese represent dramatic improvements over our previous work.The only channels in which we find evidence for deeply bound, strong-interaction-stabletetraquark states are the doubly bottom J P = 1 + , ¯3 F , ud ¯ b ¯ b and (cid:96)s ¯ b ¯ b channels alreadyidentified as supporting such bound states in earlier lattice studies, including our own. Ournew results for these channels, albeit so far only at a single m π (cid:39)
192 MeV, display significantlyimproved ground-state plateaus, largely due to the improved, Box-Sink construction. Withthe bindings at this m π reduced compared to those implied for the same m π by the fits tothe m π dependences reported in Ref. [3], we expect to find similarly reduced physical- m π results once we complete our ongoing multiple- m π analysis of these channels and are ableto perform the physical-point extrapolation. With the binding in these channels expectedto increase with decreasing m π , the final physical-point bindings will be deeper than thoseof the current study. The current, single- m π results, thus provide further confirmationof the conclusions reached in previous lattice studies regarding the existence of a ¯3 F ofstrong-interaction-stable, J P = 1 + , doubly bottom tetraquark states. Once the multiple- m π analysis has been completed, the improved ground state signals produced by the Box-Sinkconstruction should also significantly improve the reliability of the extrapolated physical- m π results.In contrast to the doubly bottom, ¯3 F , J P = 1 + channels, no bound, non-molecular groundstates are found in any of the other ten channels considered here. States with bindings > ∼ −
20 MeV, as predicted by many model and QCD sum rule studies, should have beeneasily detectable in these analyses.Among the channels where current results rule out the possibility of a bound, non-molecular ground state is the I = 0, J P = 1 + ud ¯ c ¯ b channel. This conclusion stands incontrast to that reached in our earlier study of this channel [9], where indications for possiblemodest tetraquark binding were found, albeit based on rather short, late-time plateau signals,and with worries about possible FV effects for the ensemble with m π = 164 MeV, where m π L = 2 .
4. The desirability of studying this channel at larger volumes, in order to morestrongly test the bound state interpretation of the ground state results, was stressed alreadyin Ref. [9], and the current study provides precisely such a larger-volume extension, one which,3moreover, has the advantage of the significantly improved Box-Sink analysis methodology.The indications of binding seen in the earlier study do not survive the improved, larger-volume, Box-Sink analysis: the ground-state mass is found to lie, clearly slightly above, ratherthan below, DB ∗ threshold. The very good associated ground-state effective-mass plateau isshown in the right-hand panel of Fig. 8 in App. D. We conclude that there is, in fact, nonon-molecular bound tetraquark state in this channel, and that the late-time behavior ofthe earlier Wall-Local results was likely affected by the late-Wall-Local plateauing problemdiscussed in detail for the B c and doubly bottom ¯3 F , J P = 1 + channels in Sec. IV above.Turning to a comparison to the results of other approaches, we note first that thepredictions of the heavy-quark-symmetry approach [6, 7] for the binding energies in the I = 0, J P = 1 + ud ¯ b ¯ b and J P = 1 + (cid:96)s ¯ b ¯ b channels, as well as for an absence of binding in the J P = 0 + and 1 + , I = 0 ud ¯ c ¯ b and J P = 0 + and 1 + (cid:96)s ¯ c ¯ b channels, are all compatible withcurrent lattice results.We turn next to tests of QCD sum rule and model predictions. As stressed by the authorsof the various model tetraquark studies, with the parameters of the models already fullydetermined in earlier non-exotic sector studies, the predictions for tetraquark binding areparameter free, and hence allow for highly non-trivial tests of the underlying models. In caseswhere multiple studies of the same model exist, we focus, where possible, on those in whichnumerical convergence of the prediction for the ground state energy has been demonstrated,taking these to supercede earlier, less complete studies, whether by the same or differentauthors.The earliest of the chiral quark model (ChQM) studies, Ref. [21], employs a modelwith pseudoscalar nonet exchange potentials acting between light quarks only and colourCoulomb and pairwise linear confinement potentials acting between all quarks. The resulting I = 0, J P = 1 + ud ¯ b ¯ b binding, 497 MeV, far exceeds that allowed by lattice results, thusruling out this implementation of the ChQM approach. Apart from Ref. [39], all otherChQM studies we are aware of include octet (or nonet) pseudoscalar exchange potentials andscalar exchange potentials, both acting between light quarks only, pairwise or multi-bodyconfinement potentials, acting between all quarks, and, in addition to the colour Coulomb partof the effective OGE interaction, the associated colour-spin-dependent hyperfine interaction,again acting between all quarks. Specific implementations differ, especially with regard to theform of the effective scalar exchange(s). In Ref. [28] u, d and s quarks are assumed to interact4via the exchange of a full nonet of scalar mesons. Other implementations include only theexchange of the σ , either in “SU(3)” form, where the exchange is restricted to u, d and s quarks [23, 25, 29, 31, 35, 36, 38], or “SU(2)” form, where it is restricted further to u and d quarks only [30, 35, 41]. A final ChQM implementation, used to study the I = 0, J P = 0 + and 1 + ud ¯ s ¯ b channels [39], is the quark delocalization color screening model (QDCSM),which retains pseudoscalar exchange between light quarks and OGE between all quarks, butincludes no scalar meson exchange, instead employing a modified confinement form thatallows two-center cluster configurations with quark delocalization between the clusters. Thescalar-nonet implementation is ruled out by its 32 MeV prediction for the I = 0, J P = 1 + ud ¯ b ¯ b binding, which lies well below the range, ∼ −
180 MeV, allowed by current latticeresults. The SU(2) implementation is also ruled out, this time by significant over-bindings relative to lattice results: 318 [41] or 404 MeV [30] for the I = 0, J P = 1 + ud ¯ b ¯ b channeland 178 and 199 MeV [41] for the I = 0, J P = 0 + and 1 + ud ¯ c ¯ b channels (which our resultsshow to be at most weakly unbound). The SU(3) implementation is similarly ruled out, withpredicted I = 0, J P = 1+ ud ¯ b ¯ b bindings of 341 [23], 214 [29], 217 [31], 278 [36] or 359 [38]MeV, I = 0, J P = 0 + and 1 + ud ¯ c ¯ b bindings of 136 ±
12 and 171 ±
12 MeV [36] or 178 and217 MeV [38], and I = 0, J P = 0 + and 1 + ud ¯ s ¯ b bindings of 70 and 68 MeV [35] (our latticeresults again show the ground states to be at most weakly bound in these channels). Finally,the QCDSM implementation predicts bindings of 74 and 58 MeV for the I = 0, J P = 0 + and 1 + ud ¯ s ¯ b channels, and hence is also ruled out by the absence of such deep binding in ourlattice results. We conclude that all implementations of the ChQM approach are ruled outby current lattice results, and hence that the ChQM framework does not provide a reliablephenomenological representation of QCD in the low-energy regime.The results of non-chiral models fare much better when tested against lattice results.Multiple non-relativistic quark model studies of tetraquark channels exist using the AL1 [20]and Bhaduri (or BCN) [97] potentials. These potentials are characterized by pairwise linearconfinement, colour Coulomb and regularized OGE hyperfine interactions. The AL1 modelis a variation of the BCN model in which the strict OGE relation between the strengths ofthe colour Coulomb and hyperfine interactions has been relaxed. Ref. [20] also investigatesthree additional models (the AL2, AP1 and AP2 models), characterized by alterations tothe radial dependence of the AL1 confinement potential and/or the regularization of theAL1 OGE potentials. Predictions for binding in the I = 0, J P = 1 + ud ¯ b ¯ b ; J P = 1 + (cid:96)s ¯ b ¯ b ;5 I = 0, J P = 0 + and 1 + , ud ¯ c ¯ b ; and I = 0, J P = 1 + ud ¯ c ¯ c channels exist for all these models,as well as for the Godfrey-Isgur-Capstick (GIC) model [98, 99], a “relativized” variant of thisclass of model. BCN results may be found in Refs. [17, 19, 20, 24, 26, 29, 30], AL1 results inRefs. [20, 24, 34, 40], AL2, AP1 and AP2 results in Ref. [20], and GIC results in Ref. [42].Predictions for binding in the I = 0, J P = 1 + ud ¯ b ¯ b ; J P = 1 + (cid:96)s ¯ b ¯ b ; I = 0, J P = 0 + ud ¯ c ¯ b ;and I = 0, J P = 1 + ud ¯ c ¯ b channels, obtained from the non-relativistic BCN, AL1, AL2, AP1and AP2 models, range from 131 to 167 MeV [19, 20, 24, 26, 29], 40 to 61 MeV [20], 0 to 23MeV [20, 34] and 0 to 23 MeV [20, 34], respectively. Taking the ∼ −
35 MeV variationswith model choice for these channels seen in Ref. [20] as a measure of the uncertaintiesassociated with the specifics of the implementations of the interactions in this class of model,one sees that the predictions are compatible within errors with current lattice results. Futurelattice studies, using larger-volume ensembles to better quantify possible residual FV effects,will be useful for sharpening these tests for the J P = 0 + and 1 + , I = 0 ud ¯ c ¯ b channels, wherethe results of this paper show no evidence for binding while the most recent AL1 modelstudy [34] finds modest 23 MeV bindings for both. The relativized GIC model, in contrast,predicts only 54 MeV binding in the I = 0, J P = 1 + ud ¯ b ¯ b channel and no binding in the J P = 1 + (cid:96)s ¯ b ¯ b channel [42], and hence is ruled out by lattice results.Another model approach ruled out by lattice results, specifically those for the I = 0 ud ¯ c ¯ b channels, is that of the LAMP model. In this approach, an ansatz for a relativisticconfinement potential is employed in the multiquark sector. Taking the BD interactionstudied in Ref. [32] to be specific, linearly rising mean-field confinement potentials for the B and D are merged into a one-body, two-well, W-shaped confinement potential for the BD system by truncating the individual B and D wells at the midpoint between theircenters. Such a potential allows for the delocalization of an initially single-well light-quarkwavefunction into the second of the two potential wells. The optimal well-center separationand degree of delocalization are determined by a variational minimization. The result, inwhich the effect of the OGE hyperfine interaction has been neglected, is a BD binding of 155MeV, clearly incompatible with the lattice results above. This rules out the LAMP ansatzfor the effective confining interaction.Turning to tests of QCD sum rule results, we see that recent predictions for very deeptetraquark binding (by 569 ±
260 MeV [48] in the I = 0, J P = 1 + ud ¯ b ¯ b channel, 477 ± J P = 1 + (cid:96)s ¯ b ¯ b channel, 485 ±
150 MeV [48] in the I = 0, J P = 0 + ud ¯ c ¯ b ±
160 MeV [49], 467 ±
150 MeV [53] and 200 ±
130 MeV [54] in the J P = 0 + (cid:96)s ¯ c ¯ b channel, 180 ±
110 MeV [54] in the J P = 1 + (cid:96)s ¯ c ¯ b channel, and 394 ±
170 MeV [52] in the I = 0, J P = 0 + ud ¯ s ¯ b channel) are not compatible, within quoted uncertainties, with currentlattice results. Other sum sum rule predictions (400 ±
300 MeV [43, 44], 80 ±
80 MeV [46]and 240 ±
150 MeV [55] in the I = 0, J P = 1 + ud ¯ b ¯ b channel, 290 ±
300 MeV [44], 140 ± ±
160 MeV [55] in the J P = 1 + (cid:96)s ¯ b ¯ b channel, 5 ±
100 and 60 ±
100 MeV[45] in the I = 0, J P = 0 + and J P = 1 + ud ¯ c ¯ b channels), are compatible with lattice resultswithin errors, though, with the exception of the I = 0, J P = 1 + ud ¯ b ¯ b prediction of Ref. [46],with significantly higher central values. It is not clear whether the above disagreements withlattice results are due to shortcomings associated with the use of the SVZ ansatz, inaccuraciesin the factorization approximation for some of the higher dimension condensates, a need forupdates to OPE input (e.g., the quark condensate, where precise results from FLAG arenow available [100]), or other aspects of the sum rule implementations. On this last point itis worth noting the discrepancy between the predictions, 80 ±
80 MeV [46] and 569 ± I = 0, J P = 1 + ud ¯ b ¯ b channel obtained from analyses employingthe same interpolating operator, essentially identical OPE input, and both truncating theOPE at the same dimension, D = 10. The predictions, 142 ±
80 MeV [46] and 477 ± J P = 1 + (cid:96)s ¯ b ¯ b channel also come from analyses using the sameinterpolating operator, essentially identical OPE input and truncation of the OPE at thesame dimension, D = 10. A similar discrepancy exists in the I = 0, J P = 0 + ud ¯ c ¯ b channelwhere very different ground state masses, 7140 ±
100 MeV [45] and 6660 ±
150 MeV [48], arepredicted by analyses using the same interpolating operator and similar OPE input, differingonly in the dimension at which the OPE was truncated ( D = 8 in Ref. [45] and D = 10 inRef. [48]). These observation make clear that significant, yet-to-be-fully-quantified systematicuncertainties exist for applications of the conventional implementation of the QCD sum ruleframework to the multiquark sector. The targets provided by current (and future) latticeresults for tetraquark binding (or lack thereof) in the channels investigated in this papershould help in further investigations of these issues, and have the potential to aid in thedevelopment of improved sum rule implementation strategies.We close with a reminder of improvements/extensions to the present study either alreadyin progress or planned for the near future. The first is the completion of the extension of the m π (cid:39)
192 MeV Box-Sink analysis reported here to additional 48 ×
64 ensembles having both7lighter and heavier m π . This will allow us to update and sharpen our previous physical-pointextrapolations for the bindings of the doubly bottom, J P = 1 + , ¯3 F states, and further testnon-chiral model predictions. This work is ongoing. The next step after that is to generate anew set of ensembles with even larger volume covering a similar range of m π . This will allowus to carry out quantitative FV studies, and extend the current analysis to channels, such asthe doubly charmed channels, where shallow binding remains a possibility and experimentaldetection would be less challenging than in the doubly bottom sector. ACKNOWLEDGEMENTS
BC, RL and KM are supported by grants from the Natural Sciences and EngineeringCouncil of Canada, RJH by the European Research Council (ERC) under the EuropeanUnions Horizon 2020 research and innovation programme through grant agreement 771971-SIMDAMA. All computations were carried out using a crucial allocation from ComputeCanada on the Niagara supercomputer at Scinet.
Appendix A: Tuning NRQCD
We use the tree-level ( c i = 1), tadpole-improved, NRQCD action of [96], defined by thesymmetric time evolution [101], G ( x, t + 1) = (cid:18) − δH (cid:19) (cid:18) − H n (cid:19) n U t ( x, t ) † (cid:18) − H n (cid:19) n (cid:18) − δH (cid:19) G ( t ) , (A1)where H = − c M ∆ ,H I = (cid:18) − c M − c nM (cid:19) (cid:0) ∆ (cid:1) + c i M ( ˜∆ · ˜ E − ˜ E · ˜∆) + c ∆ M H D = − c M σ · (cid:16) ˜∆ × ˜ E − ˜ E × ˜∆ (cid:17) − c M σ · ˜ BδH = H I + H D . (A2)A tilde indicates that an improved version has been used (see e.g. [102] for the tadpole-improved derivative terms) and M is a bare parameter that we tune to recover the spin-averaged kinetic mass of the connected η b and Υ. We use the fourth root of the plaquette forthe tadpole improvement and a value of n = 4 for the stability parameter.8For the tuning we use Coulomb gauge-fixed wall-sources and we implement the momenta bysolving one of the b propagators on a partially-twisted boundary [80, 81]. This is practicallyimplemented by applying a U (1) phase to the gauge field [103], U µ ( x ) → e i πθ µ /L µ U µ ( x ) (A3)which yields p µ = πL µ θ µ . The benefit of this approach is that one can stay within the sameirrep for arbitrary momenta, smoothly interpolating discretisation and finite-volume effects.We choose to twist equally in all directions to minimise the rotation-breaking hypercubicartifacts and the effects introduced by the twisting. The momenta we use is listed in Table Vand the effective masses for these twists is shown in Fig. 6. ( θ, θ, θ ) ( ap ) ( θ, θ, θ ) ( ap ) (0 , ,
0) 0 ( √ , √ , √ ) 0 . , , ) 0 . √ , √ , √ ) 0 . √ , √ , √ ) 0 . √ , √ , √ ) 0 . √ , √ , √ ) 0 . √ , √ , √ ) 0 . , , ) 0 . , ,
1) 0 . TABLE V: Twist angles, chosen so the ( ap ) values are evenly spaced between (0 , ,
0) and(1 , , A (1 + ( ap ) B + C (( ap ) ) ) exp (cid:20) − t (cid:18) ( aM ) + 12 aM K ( ap ) (cid:19)(cid:21) , (A4)where aM is the mass of the ground state and aM K is the kinetic mass. This enables us todescribe our data simply by five fit parameters.We note that, as in Ref. [104], the kinetic mass ordering is inverted compared to nature.This situation can be improved by the inclusion of yet-higher-order (1 /M ) terms in theNRQCD Hamiltonian [105], specifically the parameter c . However, we follow the approachof [104] and tune to the spin average of these kinetic masses instead. Using HPQCD’sestimate [104] for the spin-averaged mass, 9 . a − = 2 .
194 GeV [86]for our lattice spacing, we find the physical b quark corresponds to a bare mass aM = 1 .
20 30 40 50 60 70 t/a a m e ff θ = 1θ = 8/9θ = 7/9θ = 6/9θ = 5/9θ = 4/9θ = 3/9θ = 2/9θ = 1/9θ = 0 (a) Effective masses Number of momenta in fit a M K (b) FIG. 6: (Left) Effective masses for each twist angle, and corresponding fit for aM = 1 . aM M K [ G e V ] η b YPhysical η Physical Spin AveragePhysical Y
FIG. 7: Kinetic mass dependence on the bare mass aM . The fit is correlated andconstrained to have the same slope for the η b , Υ and spin average.This tuning is obtained by measuring the kinetic masses for several bare masses and fittingtheir dependence linearly, as is shown in Fig. 7.When we compute our heavy-light mesons we pack the backward and forward propagating0NRQCD propagators into one as the NRQCD propagator has only 2 Dirac indices whereasthe full light, strange, and charm propagators have four, S ( x ) = S fwd ( x ) 00 S bwd ( x ) (A5)We also apply anti-periodic boundary conditions in time to the NRQCD propagators. Appendix B: Meson masses
As we use Coulomb gauge-fixed wall sources our amplitudes are partially related tophysical matrix elements, but are polluted by the wall source [77]. If we consider for examplethe pseudoscalar J P = 0 − meson channel, we can use any of the following 8 correlationfunctions at large asymptotic times, C W LP P ( t ) = P W P L m P (cid:0) e − m P t + e − m P ( L t − t ) (cid:1) ,C W LAA ( t ) = A W A L m P (cid:0) e − m P t + e − m P ( L t − t ) (cid:1) ,C W LAP ( t ) = A W P L m P (cid:0) e − m P t − e − m P ( L t − t ) (cid:1) ,C W LP A ( t ) = P W A L m P (cid:0) e − m P t − e − m P ( L t − t ) (cid:1) ,C W WP P ( t ) = P W P W m P (cid:0) e − m P t + e − m P ( L t − t ) (cid:1) ,C W WAA ( t ) = A W A W m P (cid:0) e − m P t + e − m P ( L t − t ) (cid:1) ,C W WAP ( t ) = A W P W m P (cid:0) e − m P t − e − m P ( L t − t ) (cid:1) ,C W WP A ( t ) = P W A W m P (cid:0) e − m P t − e − m P ( L t − t ) (cid:1) . (B1)Here we use P to denote the operator with a γ -insertion and A to denote the temporal axialcurrent operator. Using these sources means, for example, that C W LP A is not equal to C W LAP .One can, however, clearly fit all 8 correlators simultaneously with 5 fit parameters. One cansimilarly analyse correlators involving the V i = γ i and T it = γ i γ t currents for J P = 1 − statesor I and γ t for J P = 0 + states. The same fit form can be used when the sink is non-localtoo, although neither matrix element will correspond to a physical one. We use these fits todetermine the masses of mesons not containing a b quark. This fit is not possible for mesonscontaining b quarks as C AA = C P P at our level of heavy-light approximation.1
Appendix C: Tables of results
In this work we use simple, local operators of the form M = ( ¯ ψ Γ φ ) (C1)to determine the masses of ordinary meson states. The results, in lattice units, and with theadditive offset in the case of mesons containing a b quark, are given in Tab. VI. Where thereare several Γ insertions listed, a simultaneous fit of the form of Eq. B1 was used, otherwise acorrelated single-cosh fit was performed on the one channel. As some of the particles listed inthe table are resonances, the ground-state masses obtained from the associated local-operatorcorrelators should not be directly compared to the masses of these resonances in nature. Tosimplify notation, we have, nonetheless, used the name of the lowest-lying resonance in thechannel to label the channel itself in these cases, trusting this will not lead to any confusion.Throughout this work we use the shorthand name B I for the scalar, J P = 0 + B meson.
State Γ am State Γ a ( m + ∆) π γ , γ t γ B γ K γ , γ t γ B ∗ γ i K ∗ (700) I, γ t B I I K ∗ (892) γ i , γ i γ t B s γ D γ , γ t γ B ∗ s γ i D ∗ γ i , γ i γ t B c γ D ∗ (2300) I, γ t B ∗ c γ i TABLE VI: Table of results for the ground-state effective masses associated with the simple,local meson operators used in this work.Tables VII through XII present the results for the lowest resolvable GEVP eigenvaluesin the ten tetraquark channels specified in the main text, obtained from correlated, single-exponential fits to the plateau regions of the corresponding optimised correlators. Entriesin column 1 specify the quantum numbers of the state. Those in columns 2, 3 and 4 givethe corresponding eigenvalue (in lattice units), the fit window interval (in lattice units) andthe associated χ /dof . Fits for higher-lying excited states having reasonable χ /dof have2been included even in cases where the fit window is relatively short since such results helpidentify those parts of the higher-lying two-meson spectrum having the largest couplings tothe operators used in constructing the optimised correlators. J P a ( m + ∆) fit (low,high) χ /dof + + TABLE VII: Eigenvalues and fit parameters for the I = 0, J P = 0 + and 1 + ud ¯ c ¯ b channels J P a ( m + ∆) fit (low,high) χ /dof + + TABLE VIII: Eigenvalues and fit parameters for the J P = 0 + and 1 + (cid:96)s ¯ c ¯ b channels3 J P a ( m + ∆) fit (low,high) χ /dof + + TABLE IX: Eigenvalues and fit parameters for the I = 0, J P = 0 + and 1 + ud ¯ s ¯ b channels J P a ( m + ∆) fit (low,high) χ /dof + + TABLE X: Eigenvalues and fit parameters for the I = 0, J P = 0 + and 1 + ud ¯ s ¯ c channels4 J P a ( m + ∆) fit (low,high) χ /dof + TABLE XI: Eigenvalues and fit parameters for the I = 0, J P = 1 + uc ¯ b ¯ b channel J P a ( m + ∆) fit (low,high) χ /dof + TABLE XII: Eigenvalues and fit parameters for the I = 0, J P = 1 + sc ¯ b ¯ b channel Appendix D: Effective mass plots
If we consider our correlators C ( t ) to behave like single exponentials at large t (a reasonableexpectation for heavy states) we can define an effective mass via am eff ( t ) = tanh − (cid:18) C ( t − − C ( t + 1) C ( t −
1) + C ( t + 1) (cid:19) . (D1)This has the desirable property that tanh − ( x ) is defined for all x in [ − , − definitions, log (cid:16) C ( t ) C ( t +1) (cid:17) or cosh − (cid:16) C ( t +1)+ C ( t − C ( t ) (cid:17) , where fluctuations can lead tonegative values of the ratios appearing in the arguments of the log or cosh − , causing thesedefinitions to break down.Figs. 8, 9, 10 and 11 show the resulting J P = 0 + and 1 + effective masses for the ud ¯ c ¯ b , (cid:96)s ¯ c ¯ b , ud ¯ s ¯ b and ud ¯ s ¯ c channels, respectively. Fig. 12, similarly, shows the effective masses forthe J P = 1 + uc ¯ b ¯ b and sc ¯ b ¯ b channels.5 t/a a ( m + ∆ ) e ff t/a a ( m + ∆ ) e ff FIG. 8: Effective mass plots for the I = 0, J P = 0 + (left) and 1 + (right) ud ¯ c ¯ b channels. t/a a ( m e ff + ∆ ) t/a a ( m + ∆ ) e ff FIG. 9: Effective mass plots for the J P = 0 + (left) and 1 + (right) (cid:96)s ¯ c ¯ b channels. t/a a ( m + ∆ ) e ff t/a a ( m + ∆ ) e ff FIG. 10: Effective mass plots for the I = 0, J P = 0 + (left) and 1 + (right) ud ¯ s ¯ b channels.6 t/a a m e ff t/a a m e ff FIG. 11: Effective mass plots for the I = 0, J P = 0 + (left) and 1 + (right) ud ¯ s ¯ c channels. t/a a ( m + ∆ ) e ff t/a a ( m + ∆ ) e ff FIG. 12: Effective mass plots for the J P = 1 + uc ¯ b ¯ b (left) and sc ¯ b ¯ b (right) channels. [1] P. Bicudo, K. Cichy, A. Peters, and M. Wagner, Phys. Rev. D , 034501 (2016),arXiv:1510.03441 [hep-lat].[2] P. Bicudo, J. Scheunert, and M. Wagner, Phys. Rev. D , 034502 (2017), arXiv:1612.02758[hep-lat].[3] A. Francis, R. J. Hudspith, R. Lewis, and K. Maltman, Phys. Rev. Lett. , 142001 (2017),arXiv:1607.05214 [hep-lat].[4] P. Junnarkar, N. Mathur, and M. Padmanath, Phys. Rev. D , 034507 (2019),arXiv:1810.12285 [hep-lat]. [5] L. Leskovec, S. Meinel, M. Pflaumer, and M. Wagner, Phys. Rev. D100 , 014503 (2019),arXiv:1904.04197 [hep-lat].[6] E. J. Eichten and C. Quigg, Phys. Rev. Lett. , 202002 (2017), arXiv:1707.09575 [hep-ph].[7] E. Braaten, L.-P. He, and A. Mohapatra, (2020), arXiv:2006.08650 [hep-ph].[8] Other arguments favoring the possibility of bound doubly bottom tetraquarks, may be foundin Refs. [106] and [107].[9] A. Francis, R. J. Hudspith, R. Lewis, and K. Maltman, Phys. Rev. D , 054505 (2019),arXiv:1810.10550 [hep-lat].[10] R. Aaij et al. (LHCb), Phys. Rev. Lett. , 112001 (2017), arXiv:1707.01621 [hep-ex].[11] R. Aaij et al. (LHCb), Phys. Rev. Lett. , 162002 (2018), arXiv:1807.01919 [hep-ex].[12] V. Khachatryan et al. (CMS), JHEP , 013 (2017), arXiv:1610.07095 [hep-ex].[13] T. Gershon and A. Poluektov, JHEP , 019 (2019), arXiv:1810.06657 [hep-ph].[14] J. Ader, J. Richard, and P. Taxil, Phys. Rev. D , 2370 (1982).[15] L. Heller and J. Tjon, Phys. Rev. D , 969 (1987).[16] J. Carlson, L. Heller, and J. Tjon, Phys. Rev. D , 744 (1988).[17] S. Zouzou, B. Silvestre-Brac, C. Gignoux, and J. Richard, Z. Phys. C , 457 (1986).[18] H. J. Lipkin, Phys. Lett. B , 242 (1986).[19] B. Silvestre-Brac and C. Semay, Z. Phys. C , 273 (1993).[20] C. Semay and B. Silvestre-Brac, Z. Phys. C , 271 (1994).[21] S. Pepin, F. Stancu, M. Genovese, and J. Richard, Phys. Lett. B , 119 (1997), arXiv:hep-ph/9609348.[22] D. Brink and F. Stancu, Phys. Rev. D , 6778 (1998).[23] J. Vijande, F. Fernandez, A. Valcarce, and B. Silvestre-Brac, Eur. Phys. J. A , 383 (2004),arXiv:hep-ph/0310007.[24] D. Janc and M. Rosina, Few Body Syst. , 175 (2004), arXiv:hep-ph/0405208.[25] J. Vijande, A. Valcarce, and K. Tsushima, Phys. Rev. D , 054018 (2006), arXiv:hep-ph/0608316.[26] J. Vijande, E. Weissman, A. Valcarce, and N. Barnea, Phys. Rev. D , 094027 (2007),arXiv:0710.2516 [hep-ph].[27] D. Ebert, R. Faustov, V. Galkin, and W. Lucha, Phys. Rev. D , 114015 (2007),arXiv:0706.3853 [hep-ph]. [28] M. Zhang, H. Zhang, and Z. Zhang, Commun. Theor. Phys. , 437 (2008), arXiv:0711.1029[nucl-th].[29] J. Vijande, A. Valcarce, and N. Barnea, Phys. Rev. D , 074010 (2009), arXiv:0903.2949[hep-ph].[30] Y. Yang, C. Deng, J. Ping, and T. Goldman, Phys. Rev. D , 114023 (2009).[31] T. Carames, A. Valcarce, and J. Vijande, Phys. Lett. B , 291 (2011).[32] R. R. Silbar and T. Goldman, Int. J. Mod. Phys. E , 1450091 (2014), arXiv:1304.5480[nucl-th].[33] M. Karliner and J. L. Rosner, Phys. Rev. Lett. , 202001 (2017), arXiv:1707.07666 [hep-ph].[34] T. F. Carams, J. Vijande, and A. Valcarce, Phys. Rev. D , 014006 (2019), arXiv:1812.08991[hep-ph].[35] X. Chen and J. Ping, Phys. Rev. D , 054022 (2018), arXiv:1807.10505 [hep-ph].[36] C. Deng, H. Chen, and J. Ping, Eur. Phys. J. A , 9 (2020), arXiv:1811.06462 [hep-ph].[37] W. Park, S. Noh, and S. H. Lee, Nucl. Phys. A , 1 (2019), arXiv:1809.05257 [nucl-th].[38] G. Yang, J. Ping, and J. Segovia, Phys. Rev. D , 014001 (2020), arXiv:1911.00215[hep-ph].[39] H. Huang and J. Ping, Eur. Phys. J. C , 556 (2019), arXiv:1902.05778 [hep-ph].[40] E. Hernndez, J. Vijande, A. Valcarce, and J.-M. Richard, Phys. Lett. B , 135073 (2020),arXiv:1910.13394 [hep-ph].[41] Y. Tan, W. Lu, and J. Ping, (2020), arXiv:2004.02106 [hep-ph].[42] Q.-F. L, D.-Y. Chen, and Y.-B. Dong, (2020), arXiv:2006.08087 [hep-ph].[43] F. S. Navarra, M. Nielsen, and S. H. Lee, Phys. Lett. B , 166 (2007), arXiv:hep-ph/0703071.[44] M.-L. Du, W. Chen, X.-L. Chen, and S.-L. Zhu, Phys. Rev. D , 014003 (2013),arXiv:1209.5134 [hep-ph].[45] W. Chen, T. Steele, and S.-L. Zhu, Phys. Rev. D , 054037 (2014), arXiv:1310.8337 [hep-ph].[46] Z.-G. Wang, Acta Phys. Polon. B , 1781 (2018), arXiv:1708.04545 [hep-ph].[47] W. Chen, H.-X. Chen, X. Liu, T. Steele, and S.-L. Zhu, Phys. Rev. D , 114005 (2017),arXiv:1705.10088 [hep-ph].[48] S. Agaev, K. Azizi, B. Barsbay, and H. Sundu, Phys. Rev. D , 033002 (2019),arXiv:1809.07791 [hep-ph].[49] S. Agaev, K. Azizi, B. Barsbay, and H. Sundu, (2019), arXiv:1912.07656 [hep-ph]. [50] S. Agaev, K. Azizi, and H. Sundu, Nucl. Phys. B , 114890 (2020), arXiv:1905.07591[hep-ph].[51] S. Agaev, K. Azizi, and H. Sundu, Phys. Rev. D , 094020 (2019), arXiv:1907.04017[hep-ph].[52] H. Sundu, S. Agaev, and K. Azizi, Eur. Phys. J. C , 753 (2019), arXiv:1903.05931 [hep-ph].[53] S. Agaev, K. Azizi, B. Barsbay, and H. Sundu, (2020), arXiv:2002.04553 [hep-ph].[54] Q.-N. Wang and W. Chen, (2020), arXiv:2002.04243 [hep-ph].[55] L. Tang, B.-D. Wan, K. Maltman, and C.-F. Qiao, (2019), arXiv:1911.10951 [hep-ph].[56] M. A. Shifman, A. Vainshtein, and V. I. Zakharov, Nucl. Phys. B , 385 (1979).[57] M. A. Shifman, A. Vainshtein, and V. I. Zakharov, Nucl. Phys. B , 448 (1979).[58] K. Ackerstaff et al. (OPAL), Eur. Phys. J. C , 571 (1999), arXiv:hep-ex/9808019.[59] M. Davier, A. Hcker, B. Malaescu, C.-Z. Yuan, and Z. Zhang, Eur. Phys. J. C , 2803(2014), arXiv:1312.1501 [hep-ex].[60] D. Boito, M. Golterman, M. Jamin, A. Mahdavi, K. Maltman, J. Osborne, and S. Peris,Phys. Rev. D , 093015 (2012), arXiv:1203.3146 [hep-ph].[61] D. Boito, M. Golterman, K. Maltman, J. Osborne, and S. Peris, Phys. Rev. D , 034003(2015), arXiv:1410.3528 [hep-ph].[62] V. Abazov et al. (D0), Phys. Rev. Lett. , 022003 (2016), arXiv:1602.07588 [hep-ex].[63] V. M. Abazov et al. (D0), Phys. Rev. D , 092004 (2018), arXiv:1712.10176 [hep-ex].[64] R. Aaij et al. (LHCb), Phys. Rev. Lett. , 152003 (2016), [Addendum: Phys.Rev.Lett. 118,109904 (2017)], arXiv:1608.00435 [hep-ex].[65] A. M. Sirunyan et al. (CMS), Phys. Rev. Lett. , 202005 (2018), arXiv:1712.06144 [hep-ex].[66] T. Aaltonen et al. (CDF), Phys. Rev. Lett. , 202006 (2018), arXiv:1712.09620 [hep-ex].[67] M. Aaboud et al. (ATLAS), Phys. Rev. Lett. , 202007 (2018), arXiv:1802.01840 [hep-ex].[68] G. K. C. Cheung, C. E. Thomas, J. J. Dudek, and R. G. Edwards (Hadron Spectrum), JHEP , 033 (2017), arXiv:1709.01417 [hep-lat].[69] D. C. Moore and G. T. Fleming, Phys. Rev. D74 , 054504 (2006), arXiv:hep-lat/0607004[hep-lat].[70] C. Michael and I. Teasdale, Nucl. Phys. B , 433 (1983).[71] M. Luscher and U. Wolff, Nucl. Phys. B , 222 (1990). [72] B. Blossier, M. Della Morte, G. von Hippel, T. Mendes, and R. Sommer, JHEP , 094(2009), arXiv:0902.1265 [hep-lat].[73] B. Hrz and A. Hanlon, Phys. Rev. Lett. , 142002 (2019), arXiv:1905.04277 [hep-lat].[74] A. Billoire, E. Marinari, and G. Parisi, Phys. Lett. B , 160 (1985).[75] R. Gupta, G. Guralnik, G. W. Kilcup, and S. R. Sharpe, Phys. Rev. D , 2003 (1991).[76] R. Hudspith (RBC, UKQCD), Comput. Phys. Commun. , 115 (2015), arXiv:1405.5812[hep-lat].[77] D. Daniel, R. Gupta, G. W. Kilcup, A. Patel, and S. R. Sharpe, Phys. Rev. D , 3130(1992), arXiv:hep-lat/9204011.[78] S. Aoki, M. Fukugita, S. Hashimoto, Y. Iwasaki, K. Kanaya, Y. Kuramashi, H. Mino,M. Okawa, A. Ukawa, and T. Yoshie (JLQCD), Nucl. Phys. B Proc. Suppl. , 354 (1996),arXiv:hep-lat/9510013.[79] D. Antonio et al. (RBC, UKQCD), Phys. Rev. D , 114501 (2007), arXiv:hep-lat/0612005.[80] P. F. Bedaque, Phys. Lett. B , 82 (2004), arXiv:nucl-th/0402051.[81] C. Sachrajda and G. Villadoro, Phys. Lett. B , 73 (2005), arXiv:hep-lat/0411033.[82] Y. Aoki et al. (RBC, UKQCD), Phys. Rev. D , 074508 (2011), arXiv:1011.0892 [hep-lat].[83] See App. B for more discussion of how to obtain physical amplitudes from these correlators.[84] T. Bhattacharya, R. Gupta, G. Kilcup, and S. R. Sharpe, Phys. Rev. D , 6486 (1996),arXiv:hep-lat/9512021.[85] S. Aoki et al. (PACS-CS), Phys. Rev. D , 074503 (2010), arXiv:0911.2561 [hep-lat].[86] Y. Namekawa et al. (PACS-CS), Phys. Rev. D , 094512 (2013), arXiv:1301.4743 [hep-lat].[87] http://luscher.web.cern.ch/luscher/openQCD/. [88] We implemented this heavy-quark action using AVX/FMA2 vector intrinsics directly inopenQCD. Typical charm-quark propagator inversions are comparable to those of our strange-quark propagators.[89] M. Luscher, Comput. Phys. Commun. , 209 (2004), arXiv:hep-lat/0310048.[90] M. Luscher, JHEP , 011 (2007), arXiv:0710.5417 [hep-lat].[91] M. Luscher, JHEP , 081 (2007), arXiv:0706.2298 [hep-lat].[92] S. Aoki, Y. Kuramashi, and S.-i. Tominaga, Prog. Theor. Phys. , 383 (2003), arXiv:hep-lat/0107009.[93] S. Aoki, Y. Kayaba, and Y. Kuramashi, Nucl. Phys. B , 271 (2004), arXiv:hep-lat/0309161. [94] Y. Namekawa et al. (PACS-CS), Phys. Rev. D , 074505 (2011), arXiv:1104.4600 [hep-lat].[95] D. Mohler, C. Lang, L. Leskovec, S. Prelovsek, and R. Woloshyn, Phys. Rev. Lett. ,222001 (2013), arXiv:1308.3175 [hep-lat].[96] G. Lepage, L. Magnea, C. Nakhleh, U. Magnea, and K. Hornbostel, Phys. Rev. D , 4052(1992), arXiv:hep-lat/9205007.[97] R. Bhaduri, L. Cohler, and Y. Nogami, Nuovo Cim. A , 376 (1981).[98] S. Godfrey and N. Isgur, Phys. Rev. D , 189 (1985).[99] S. Capstick and N. Isgur, AIP Conf. Proc. , 267 (1985).[100] S. Aoki et al. (Flavour Lattice Averaging Group), Eur. Phys. J. C , 113 (2020),arXiv:1902.08191 [hep-lat].[101] R. Lewis and R. Woloshyn, Phys. Rev. D , 014502 (2009), arXiv:0806.4783 [hep-lat].[102] S. Groote and J. Shigemitsu, Phys. Rev. D , 014508 (2000), arXiv:hep-lat/0001021.[103] O. H. Nguyen, K.-I. Ishikawa, A. Ukawa, and N. Ukita, JHEP , 122 (2011), arXiv:1102.3652[hep-lat].[104] R. Dowdall et al. (HPQCD), Phys. Rev. D , 054509 (2012), arXiv:1110.6887 [hep-lat].[105] A. V. Manohar, Phys. Rev. D , 230 (1997), arXiv:hep-ph/9701294.[106] A. V. Manohar and M. B. Wise, Nucl. Phys. B , 17 (1993), arXiv:hep-ph/9212236.[107] A. Czarnecki, B. Leng, and M. Voloshin, Phys. Lett. B778