Evidence for charm-bottom tetraquarks and the mass dependence of heavy-light tetraquark states from lattice QCD

We continue our study of heavy-light four-quark states and find evidence from lattice QCD for the existence of a strong-interaction-stable $I(J^P)=0(1^+)$ $ud\bar{c}\bar{b}$ tetraquark with mass in the range of 15 to 61 MeV below $\bar{D}B^*$ threshold. Since this range includes the electromagnetic $\bar{D}B\gamma$ decay threshold, current uncertainties do not allow us to determine whether such a state would decay electromagnetically, or only weakly. We also perform a study at fixed pion mass, with NRQCD for the heavy quarks, simulating $qq^\prime \bar{b}^\prime \bar{b}$ and $q q^\prime \bar{b}^\prime\bar{b}^\prime$ tetraquarks with $q,\, q^\prime =ud$ or $\ell s$ and variable, unphysical $m_{b^\prime}$ in order to investigate the heavy mass-dependence of such tetraquark states. We find that the dependence of the binding energy follows a phenomenologically-expected form and that, though NRQCD breaks down before $m_{b^\prime}=m_c$ is reached, the results at higher $m_{b^\prime}$ clearly identify the $ud\bar{b}^\prime \bar{b}$ channel as the most likely to support a strong-interaction-stable tetraquark state at $m_{b^\prime}=m_c$. This observation serves to motivate the direct $ud\bar{c}\bar{b}$ simulation. Throughout we use dynamical $n_f=2+1$ ensembles with pion masses $m_\pi=$415, 299, and 164 MeV reaching down almost to the physical point, a relativistic heavy quark prescription for the charm quark, and NRQCD for the bottom quark(s).


I. INTRODUCTION
Many theoretical efforts since the formulation of QCD have hypothesized the existence of exotic states containing four or more quarks and/or antiquarks (for a recent review see [1] and references therein). It is only in the last decade that unambiguously exotic states, including the hidden-charm pentaquark states recently discovered at LHCb [2] and at least some of the XYZ states [1], which fail to fit into the standard quark model picture, have begun to be observed experimentally. These experimental results have shown that complicated many-quark structures do exist in nature and the goal for theorists is to investigate why some such multi-quark structures are preferred, and to elucidate the mechanisms underlying their existence. The mechanisms behind the binding of such configurations should help to provide greater insight into the complex phenomena of QCD in the non-perturbative realm.
In this work we use lattice QCD to investigate configurations of two light quarks and two heavy antiquarks in channels expected to be favorable to the formation of bound, exotic tetraquark states. In general, such four-quark bound states have not yet been definitely proven to exist experimentally, but there are indications, both from models and from lattice QCD , that they should exist. A benefit of employing the lattice approach is that unphysical quark masses can be used as input simulation parameters, allowing for an extended investigation of the underlying binding mechanisms.
In a prior work [29] we predicted the existence of udbb and sbb tetraquarks with quantum numbers I(J P ) = 0(1 + ) and 1 2 (1 + ), respectively, using lattice QCD. The focus on these channels was motivated by features of the splittings in the heavy baryon spectrum. The key observation is that a pair of heavy antiquarks in a color 3 c configuration will serve as the source of a nearly static color 3 c field, analogous to that produced by the heavy quark in a singly-heavy baryon. From the pattern of splittings in the heavy baryon sector, it is clear that a strong spin-dependent attraction exists for light quark ud or s pairs in Jaffe's [30] "good-diquark" configuration (color3 c and I = J = 0 or I = 1/2, J = 0, for ud or s, respectively). The strength of this attraction, moreover, increases the lighter the light quark mass. In a doubly-heavy I(J P ) = 0(1 + ) or 1 2 (1 + ) qq QQ tetraquark channel, this good-diquark attraction is available to a localized four-quark state, but not to the lowest-lying asymptotic two-meson state in the same channel, where the spin-dependent interactions of the light quarks with their heavy antiquark partners from the same meson are suppressed by the heavy quark mass. The attractiveQQ 3 c color-Coulomb interaction provides a further contribution to binding in the localized tetraquark configuration not available to two separated mesons. This picture leads to the expectation that bound udbb and sbb tetraquark states should exist in the I(J P ) = 0(1 + ) and 1 2 (1 + ) channels, with binding of order 150 − 200 MeV for the former and a reduced binding for the latter. These semi-quantitative expectations were confirmed by explicit lattice simulations, in which we found states bound by |∆E udbb | = 189(10)(3) and |∆E sbb | = 98(7)(3) MeV relative to the corresponding two-meson thresholds, BB * and B s B * respectively, at the physical mass point. With such binding energies, these states are not only strong-interaction stable, but can, in fact, decay only weakly. Other lattice studies using static prescriptions of the bottom quarks and heavier than physical sea quarks [31][32][33][34][35][36][37][38][39][40] analyzing similar quantities have observed attractive potentials in this channel and some [28,[41][42][43] have indicated binding as well.
Since the discovery of the doubly-charmed Ξ cc baryon at LHCb [44], sum rule calculations and phenomenological models, see e.g. [24,25,45,46], have also led to the identification of these channels as favorable to doubly heavy tetraquark binding. In the model calculations, an important role is also played by the attractive natures of the light-quark spin-dependent interactions and short distance color-Coulomb potential for heavy antiquark pairs in a 3 c color configuration [1].
Assuming the above picture correctly captures the basic physics involved in the binding observed in the udbb and sbb systems, one should see binding which grows as the heavy quark mass is increased since the short-distance, 3 c color-Coulomb attraction should scale as the reduced mass of the heavy anti-diquark system. One should also see contributions to the binding which are, to a first approximation, independent of the heavy quark mass corresponding to the good-light-diquark attraction in the static heavy quark limit. Corrections to this limit should produce binding corrections proportional to the inverse of the heavy quark mass. The increase with decreasing heavy-quark mass of the net residual light-heavy spindependent attraction in the two-meson (vector-pseudoscalar) threshold will also reduce the tetraquark binding relative to this threshold and produce binding corrections proportional to the inverse of the heavy quark mass.
These qualitative expectations can be tested by extending our previous study to include unphysical values of the masses of one or both of the two heavy antiquarks. The only limitation is the NRQCD approach to treating theb, which breaks down before the charm quark mass is reached. As we will see, the results of this variable-heavy-mass study confirm the picture outlined above. Although the NRQCD-based approach does not allow us to push this study down to the charm mass, the pattern of bindings obtained as the variable heavy mass (or masses) is (are) decreased below m b can, nonetheless, be used to identify which channel (or channels) involving one or two charmed antiquarks is (are) most likely to support bound tetraquark states. These considerations lead us to focus our attention, and direct simulation efforts, on the most favorable of these channels, which turns out to be the udcb channel.
In this paper we will first detail the variable heavy mass study outlined above, and then discuss the results of our direct simulations of the udcb channel. Throughout, we will use the same three dynamical, fixed lattice spacing, n f = 2 + 1 PACS-CS ensembles employed in our previous udbb study. These ensembles have pion masses m π = 164, 299 and 415 MeV. A relativistic prescription will be used for the charm quark and, as before, NRQCD for the bottom quark. Evidence is presented for the existence of an I(J P ) = 0(1 + ) udcb tetraquark bound with respect to the lowest free two-meson threshold,DB * , in this channel.
For the variable-heavy-mass study, we will focus our attention on the ensemble with m π = 299 MeV, and study the heavy anti-diquark mass dependence for unphysical tetraquark candidates qq QQ , with q = u, q = d, s and eitherQ =Q orQ =Q.

II. PHENOMENOLOGY OF HEAVY-LIGHT TETRAQUARKS
We focus on heavy-light tetraquark candidates qq QQ that can be pictured as a combination of a "good" light-diquark qq and heavy anti-diquarkQQ with qq = ud, or s and Q, Q = b or c. The color3 c , J = 0, flavor-antisymmetric good-light-diquark configuration is accessible only when the heavy anti-diquark is in a color 3 c . Assuming no spatial excitation between the heavy antiquarks, the heavy anti-diquark spin is necessarily J h = 1 when Q = Q . The favored tetraquark configuration is then J P = 1 + .
In the limit of infinitely heavy Q, Q , m Q,Q → ∞, the attractive nature of the color-Coulomb potential guarantees a bound ground state of the qq QQ type [47][48][49]. Whether a bound state is realized away from this limit, in particular, when Q, Q are charm or bottom quarks, depends on the details of non-perturbative effects in such systems.
The phenomenological arguments outlined in the previous section, based on observed splitting patterns in the heavy baryon system, were shown in Ref. [29] to suggest the likelihood of the existence of tetraquark bound states of the qq bb -type with tetraquark binding increasing with decreasing light quark mass(es).
The lattice results of Ref. [29] not only confirmed the existence of these bound states, but produced binding energies for physical light quark masses in line with those expected based on the heavy baryon spectrum arguments. Assuming the picture underlying this successful prediction is correct, and the good-diquark contribution to binding indeed increases with decreasing light quark mass, this implies it is imperative to have access to light quark masses as close to physical as possible in lattice simulations of such tetraquark channels. This is a firm prediction of this binding mechanism, one that may seem counterintuitive given the common experience with lattice calculations in other channels, for example, in the study of multi-baryon states, where a decrease in constituent quark masses also decreases the binding energy (for a collection of recent results and presentation of the issues faced see [50] and references therein).
It is possible to further test the qualitative physical picture underlying this understanding of the udbb and sbb tetraquark binding observed in Ref. [29] by studying related systems with variable heavy antiquark mass(es). This study is carried out using the same NRQCD action used previously, in Ref. [51], for the physical bottom quark case. A brief outline of the NRQCD framework, together with details of the implementation of the variable b mass, are provided in Appendix A. The NRQCD heavy-mass parameter, m Q , was tuned by measuring the dispersion relation of the spin-averaged Υ and η b . We also computed static propagators, allowing the extrapolation of m b to ∞ to be carried out for the Q = Q case. For the variable heavy mass study, we focus, to be specific, on the intermediate ensemble, Given the qualitative physical picture outlined above, we expect there to be a contribution to tetraquark binding from the color-Coulomb attraction between the two heavy antiquarks in the color 3 c configuration which scales linearly with the reduced mass of the heavy antidiquark system. There should also be a contribution to the binding which, for a given light-diquark channel, should be independent of the heavy quark masses, reflecting the attractive nature of the good-light-diquark configuration in the infinite heavy quark mass limit. Finally, there should be contributions to the binding resulting from the presence of residual heavy-light interactions, scaling as the inverse of the heavy quark mass(es), in both the tetraquark and two-meson threshold states. Contributions of the former type should scale as 1 m h1 + 1 m h2 for tetraquark states with heavy antiquark masses m h1 and m h2 . Residual interactions between a heavy quark and light-diquark are also present in the heavy baryon systems. Comparing the Σ h − Λ h and Ξ h − Ξ h splittings for h = b, c, one finds heavy baryon residual interactions depending on both the inverse of the heavy quark mass and the type (ud or s) of light-diquark. We thus expect the coefficient of 1 m h1 + 1 m h2 for the residual heavy-light tetraquark interactions to be different for tetraquarks containing ud and s good diquarks. With the ratio of observed charm and bottom vector-pseudoscalar splittings in good agreement with expectations based on the assumption that these scale as the inverse of the relevant heavy quark mass, the contributions to tetraquark binding from residual heavy-light interactions in the corresponding two-meson threshold state can be directly determined from the observed B * − B, B * s − B s , D * − D and D * s − D s splittings, bearing in mind that the correct two-meson threshold must be chosen. Thus, for example, denoting by V and P theb vector and pseudoscalar states, one has that, for tetraquarks of the udbb type, the relevant threshold is B * P for m b < m b but BV for m b > m b . We assume that the observed 1/m h scaling of the bottom and charm vector-pseudoscalar splittings persists for the variable-b-mass V − P splittings. The two-meson threshold contributions for a given physical-to-variable b quark mass ratio, r = m b /m b , are then fixed by the observed charm/bottom meson splittings, and depend on m h1 and m h2 in a manner that varies depending on the relation between these two masses. We use that the vector meson, V , and pseudoscalar meson, P , lie, respectively, 1 4 (m V − m P ) above and 3 4 (m V − m P ) below the spin-average of the V and P masses.
Taking these expectations into account, and writing the results in terms of the physicalto-variable mass ratio r = m b /m b , one expects to obtain a good-quality fit to the binding energies of tetraquarks with at least one unphysical variable-mass antiquark using an expression having the form for the udb b case, where the first term represents the Coulomb binding contribution, the second the good-ud-diquark attraction, the third the residual heavy-light interactions in the tetraquark state and the fourth the two-meson threshold contribution. The numerical value appearing in the fourth term follows from the observed meson splittings. Similarly, for the udb b case, one expects the form to provide a good representation for m b > m b and the form to provide a good representation for m b < m b . The corresponding expectations for the cases involving an s, rather than ud, good-diquark are for sb b , for sb b with m b > m b and For the variable-b-mass study just described, details of the setup and determination of the resulting tetraquark binding energies may be found in Appendices A and B. Fig. 1 and Tab. I display these results for m b running from 6.29 to 0.60 times m b . The restriction to m b ≥ 0.60m b is designed to ensure that the b masses considered are all sufficiently heavy that the NRQCD approximation is reliable. The same 2 × 2 GEVPs used in Ref. [51] are employed for the udb b and sb b channels, while new 3 × 3 GEVPs, described in more detail in the next section, are used for udb b and sb b . The results of a fit to this data using the forms detailed above are shown in Fig. 1 and Tab. II. The success of this fit in describing tetraquark binding energies over a wide range of variable heavy quark masses confirms that the physical picture underlying those fit forms successfully captures the main features responsible for the binding observed in these systems. While the use of NRQCD precludes extending the results of the variable-b-mass study down to the charm mass, the pattern of binding energies shown in Fig. 1 clearly points to the udcb channel as by far the most likely among the four channels, udcb , scb , udcc , or scc , in which one or both of theb antiquarks in udbb is replaced byc, to support a strong-interaction-stable bound state. A naive extrapolation of the results for this channel, moreover, produces a result very nearDB * threshold, strongly motivating direct udcb simulations using a relativistic action for the charm quark. Similar naive extrapolations of the variable-heavy-mass results suggest none of the other three channels is likely to support a strong-interaction-stable bound state. We thus focus, in what follows, on the udcb channel, leaving detailed simulations of the other channels for a subsequent work. (11) -86(10) [51]. (14) 40(5) -116(10) 22(3) TABLE II. Fit results for the constants parametrizing the heavy-quark mass dependence of tetraquark binding energies on the ensemble E M with pion mass m π = 299 MeV.

III. LATTICE CORRELATORS AND OPERATORS
The generic form of a lattice QCD correlation function for a particle at rest, i.e. p = 0, in Euclidean time is given by with the interpolating operators O i being chosen to have the quantum numbers of the continuum state to be studied. For example, the simplest local meson operator is where upper (Greek) indices denote Dirac spin and lower (Roman) color. The q and q represent the constituent quark flavors.
In the case of the3 F , J P = 1 + tetraquarks, the relevant free-streaming two-meson thresholds are given by the sums of the lowest-lying pseudoscalar (P) Γ = γ 5 and vector (V) Γ = γ i meson masses. Tab. III provides a list of these thresholds for the cases of interest here. Given the phenomenological picture of the previous sections, a natural first choice of interpolating operator for qq QQ -type tetraquarks is one with a diquark-anti-diquark structure. Using the epsilon identity abc dec = δ ad δ be − δ ae δ bd , the local version of this operator takes the form where C = iγ y γ t is the charge-conjugation matrix. This operator hasQQ color 3 c , spin 1 and light-quark flavor-spin-color (3 F , 0,3 c ). A second possible local operator is one whose discrete structure is meson-meson-like. For tetraquark channels with Q = Q this could be: When Q = Q , a second local flavor antisymmetric, J P = 1 + meson-meson-like combination can also be constructed. Suppressing the spin indices, the two possible "meson-meson" interpolating operators in this case are: With these interpolating operators, there are several options to study the ground state energies of the proposed tetraquarks channels. One is to form the so-called binding correlator, the ratio of tetraquark correlation functions to the product of correlation functions, C P P (t) and C V V (t), of the pseudoscalar and vector mesons making up the corresponding non-interacting two-meson threshold in the channel in question: This definition is beneficial in the present study as the additive mass renormalisation for NRQCD quarks explicitly cancels. The binding correlator behaves, for large t, as ∼ e −∆Et , with ∆E = m 0 −m V −m P , where m 0 is the ground state mass in the channel. Thus, observing an exponentially-increasing binding correlator signals the existence of a ground state lighter than the corresponding two-meson threshold and provides compelling evidence for a bound tetraquark state.
In the binding correlator one might also gain signal through cancellations between the two-meson and tetraquark fluctuations. At the same time, however, forming the binding correlator may introduce a difficult-to-control systematic through contamination of the tetraquark signal with residual excited state effects originating from the two-meson correlators. The potential problem is that ground state saturation might occur at different lattice times for each of the three different correlation functions entering the binding correlator. The binding correlator plateau will then depend on the slowest-plateauing of the three constituent correlators.
To handle (and quantify, if present) this effect, we also compute the individual correlation functions for the tetraquark candidates and two-meson threshold combination(s). This is especially important, since we expect a significantly smaller binding energy for udcb tetraquarks than was observed for the udbb channel. An unambiguous determination of the ground state energies is thus essential. To combine the two mesons we compute the product C P P (t)C V V (t). The resulting mass from a single-exponential fit to this combination for the case of theDB * threshold is given in Tab. IV. Note, however, that in this case the NRQCD additive mass renormalisation does not drop out and the results are given in lattice units with this shift included.

A. Variational analysis
With access to several operators with the same quantum numbers, a variational analysis can be used to determine the ground and excited state energies. In the caseQ =Q , this analysis involves the 2 × 2 matrix WhenQ =Q , there are now two meson-meson interpolating operators, allowing use of the 3 × 3 matrix Given the matrix F (t), one can solve the generalised eigenvalue problem (GEVP) for some reference time t 0 1 , where ν are the (generalised) eigenvectors and λ i (t) the eigenvalues. The solution of the GEVP gives independent eigenvalues corresponding to different states in the system, In this work the aim will be to determine the ground state of the system (λ 1 ). As such, a 2 × 2 or 3 × 3 matrix should suffice, so long as the chosen operators have good overlap with the desired ground state.  [29,51,53]). These configurations [54] use the Iwasaki gauge action with β = 1.9 and non-perturbatively tuned clover coefficient c SW = 1.715. TheDB * threshold is extracted from a single-exponential fit to the product of the two relevantD-and B * -meson correlators. Due to NRQCD's additive mass renormalisation, these values have not been converted into physical units.

IV. NUMERICAL SETUP
Throughout this work the calculations are performed on three n f = 2+1, clover-improved [55], Iwasaki gauge [56], PACS-CS ensembles introduced in [54] and which we label by E H , E M , and E L . The lattice spacing is [52] a −1 =2.194(10) GeV for all three ensembles. We use a partially-quenched strange quark tuned to the (connected) φ-meson [57], which gives a near-physical kaon mass in the chiral limit. The labeling (and pion masses) of the ensembles are consistent with those used in our previous work [29,51,53]. In the valence sector we use Coulomb gauge-fixed wall sources using the FACG algorithm [58], as before. We put sources at multiple time positions and compute propagators for light, strange and charm quarks using a modified deflated SAP-solver [59]. An overview of the lattice parameters and ensemble properties can be found in Tab. IV.
To reliably handle charm quarks on these lattices we employ a relativistic heavy quark (RHQ) action [60][61][62][63], in particular, the "Tsukuba" formulation [61]: The common approach of RHQ actions is to re-interpret the discretisation effects and to re-tune the fully relativistic lattice action by introducing anisotropy in the valence sector to reproduce the correct dispersion relation, i.e., the equivalence between rest and kinetic mass, and physical spectrum. For our RHQ action, the tuning parameters have been previ-ously computed for the ensembles studied here. The values of the parameters κ f , r s , ν s , c E and c B can be found in [64]. Meson masses using quark propagators computed from our implementation of this action are seen to be within ∼ 1% of the experimentally observed values, with splittings between, for example, D and D * mesons equally well-behaved.
For bottom quarks, as noted above, we employ the NRQCD action, as detailed in Appendix A. Overall, the NRQCD action is known to capture the relevant heavy-light quark physics and account for relativistic effects at the few percent level [65][66][67].
A list of the bare quark masses used in this study and the ratios of kinetic masses compared to the physical bottom quark is provided in Tab. V.

V. INDICATIONS OF A BOUND udcb TETRAQUARK IN NATURE
The successful phenomenological description of the heavy-quark mass dependence of J P = 1 + tetraquark binding energies described earlier identifies the udcb channel as the only such channel containing at least one charm antiquark likely to support a strong-interaction-stable bound tetraquark state. Given this insight, in this section, we focus our attention and resources on this channel and present results for the direct calculation of the mass of an I(J P ) = 0(1 + ) udcb tetraquark.
Our lattice results for the ground (red) and first excited (blue) state binding energies, obtained from the 3 × 3 GEVPs, are shown in Fig. 2. Results obtained using the corresponding 2 × 2 GEVPs, in which the local operator withD * B discrete structure has been omitted, are shown for comparison in green, offset slightly in t for presentational clarity.
Results shown in the left-hand panels are those obtained using the binding correlators Eq. (12). For the E M and E L ensembles, one observes a clear rising exponential behavior in both the 3 × 3 and 2 × 2 GEVP results, indicating the presence of a bound udcb tetraquark ground state for these ensembles.
The right-hand panels of Fig. 2, in contrast, focus on the binding energies themselves, which are shown in log-effective form, E eff (t), and given in physical units. The choice of log-effective form is for presentational purposes only; our final binding energies are obtained through fits to the eigenvalues shown in the left-hand panels. The right-hand panels also provide a direct comparison of the eigenvalue-fit results (shown by the shaded bands) with the corresponding log-effective-energies.
The data allow the first two eigenvalues to be clearly resolved for each ensemble, with a signal that can be followed to t/a ∼ 16. The overlaps with the chosen multi-quark interpolating operators therefore seem good. Log-effective energies of the states corresponding to these eigenvalues are shown in the right panels of Fig. 5 in Appendix C. These results confirm that the signal is strong enough to reach the plateau regions of aE eff (t). In the case of the binding correlator, the plateaus are in general reached later and are therefore shorter since the point of signal deterioration is the same in both approaches. This is a sign of unwanted contamination from excited states in the meson denominator of the binding correlator ratio, as discussed in Sec. III.
In all cases the third eigenvalue of the 3 × 3 GEVP (not shown here) is found to lie significantly higher than the first two, a fact confirmed by the eigenvalue plots in the lefthand panels of Fig. 5 in Appendix C.
Note that the results for the ground state eigenvalues obtained from the smaller 2 × 2 GEVP defined in Eq. (13), shown by the green vertical dashes in Fig. 2, agree very well with those obtained from the corresponding 3 × 3 GEVP. In contrast, the 2 × 2 results for the second eigenvalues, denoted by the green diagonal crosses, lie lower than those of the 3 × 3 analysis, and hence correspond to higher energies. This is as expected for analyses employing interpolating fields having a good overlap with the ground state, where both the 2 × 2 and 3 × 3 analyses should provide good access to the ground state energy, but the 2 × 2 analysis, where the second eigenvalue will have to represent, not just the effects of the actual first excited state, but also those of all higher excited states, should yield results corresponding to higher energies than those obtained from the 3 × 3 analysis.
For the heavy pion ensemble, E H , with m π = 415 MeV, the ground state plateaus at theDB * threshold, whereas already for E M , with m π = 299 MeV, the second eigenvalue is consistent with the threshold and the first eigenvalue with a bound tetraquark interpretation. Finally, for E L , with m π = 164 MeV, both the first and second eigenvalues have central values corresponding to states below threshold. In light of the error bars, we suspect that, in this case, the second eigenvalue corresponds to a scattering state affected by the finite lattice volume rather than a genuine bound state. A more detailed investigation of this question must, however, await a future study involving sufficiently high-statistics ensembles with larger volumes for near-physical m π .
To estimate the binding energies we have performed (uncorrelated) single-exponential fits Eq. (16) to the individual eigenvalues and accepted results satisfying χ 2 /d.o.f. 1. We perform fits to the eigenvalues of both the tetraquark correlators and binding correlators. In the former case, we subtract the sum of the threshold two-meson state masses to obtain the binding energy. All accepted results are compared and a result that is representative of both procedures is chosen. The final fit ranges are t/a ∈ [10 : 23] for E L , t/a ∈ [14 : 20] for E M , and t/a ∈ [17 : 21] for E H . Further details of this procedure, as well as the resulting fit stabilities, are given in Appendix C. The final results obtained in this manner are listed in Tab. VI and plotted as the grey shaded bands in the right-hand panels of Fig. 2.
Although further control of the dominant systematic uncertainties is necessary before an accurate prediction can be made, the results of our direct simulation provide good evidence for the existence of a strong-interaction-stable J P = 1 + udcb tetraquark state at physical m π . We take the upper bound of the E M result and lower bound of the E L result as providing  the best assessment of the likely range of binding, leading to the expectation for the binding energy of the udcb tetraquark ground state relative to theDB * threshold. For presentational simplicity in what follows, we will also quote this result in the equivalent compact form ∆E udcb = −38(23) MeV, using the midpoint of the range in Eq. (18) as the central value, and half the width as the error estimate. This strategy for estimating the binding energy implied by our results appears plausible, given the ensembles currently available to us; since we expect both some deepening of the binding in going from E M to E L (as a consequence of the increasing spin-dependent, good-light-diquark attraction with decreasing light quark mass) but also potentially non-negligible finite-volume effects due to the small m π L of the E L ensemble. It thus seems reasonable to expect the true, infinite volume binding for physical m π to lie somewhere between the two bounds noted above. Our finite-volume systematic uncertainty, at present, clearly dwarfs the statistical errors on the individual measurements and only a careful finite-volume study will allow for a more accurate prediction. We are currently in the process of generating ensembles with larger volumes at near-physical m π , and, once this is completed, we will report on the results obtained using these ensembles to more fully quantify the impact of finite-volume effects on the estimated physical-m π binding in a future work. Taking PDG [68] values for the physical D and B * masses, 2 the binding energy estimate obtained above corresponds to a tetraquark mass of ≈ 7154(23) MeV.
For completeness, we quote here also the energy differences between the first excited state (corresponding to the second eigenvalue) andDB * threshold, which are 18(6) MeV for E H , 8(8) MeV for E M , and −26 (7) MeV for E L . These are compatible with our interpretation of the second eigenvalue as corresponding to theDB * threshold for the E M and E L ensembles, with potentially non-negligible finite-volume effects present in the E L case. One should, however, also bear in mind that the limited size of our basis of operators may impact how reliably we can extract the energy of the first excited state, particularly if this state corresponds to an infinite-volume scattering state and, as is likely, finite-volume effects are not negligible for the E L ensemble. A larger basis of operators and finite-volume analysis are thus desirable for a more robust study of the nature of this state.

A. Discussion of systematic uncertainties
In this calculation, the systematic error originating from the uncertainty of the lattice spacing in Tab. IV is negligible in comparison to the statistical error of the final result.
As noted above, we believe the dominant source of uncertainty in our result comes from finite-volume effects. The light-quark-mass dependence of the udcb masses covers multiple m π L and appears to support the interpretation of the ground states for the E M and E L ensembles as corresponding to genuine bound tetraquark states, and the expectation that such a bound tetraquark state will therefore also exist at physical m π . However, without additional ensembles with larger lattice volumes, a direct study of the scaling of the binding with m π L, and a full extrapolation to infinite volume and physical m π , is not feasible at present, and must be left for future study.
Aside from general terms ∼ exp(−m π L), finite lattice volume may affect particle energies in two ways. First, a bound state receives corrections proportional to exp(−|p|L), where |p| is the binding momentum, defined via the energy-momentum relation −∆E = E 2 where E 1 and E 2 are the energies of the threshold particles. Second, a state which becomes a scattering state in the infinite-volume limit, but which lies below threshold at finite volume, may receive power-law corrections in 1/L, which for the n = 0 states depend on a 0 /L, where a 0 is the scattering length of the particles that define the two-meson threshold in the channel in question [69,70].
Unlike the case of the udbb and sbb channels, where binding energies much larger than expected finite-volume effects were found [29], finite-volume effects may play a more important qualitative role in the udcb channel, where the expected binding and overall energy scale are lower. Even though a rigorous finite-volume study must be left for future work, we observe that, if the ground state energies for the E M and E L ensembles do, indeed, as the evidence above suggests, correspond to bound udcb tetraquark states, the binding momenta for these states would be (|p|) L = 373(171) MeV and (|p|) M = 245(140) MeV, respectively, implying a strong suppression of the finite-volume effects on the determined bound-state energies. Finite-volume studies with larger lattice volumes, however, remain desirable to more strongly test the bound-state interpretation of the ground states for the E M and E L ensembles.

VI. POSSIBLE EXPERIMENTAL DETECTION
With a predicted binding of between 15 and 61 MeV below theDB * threshold, it is unclear whether the J P = 1 + udcb tetraquark predicted by our results should, like its udbb and sbb analogues, be both strong-and electromagnetic-interaction stable, or be able to decay electromagnetically, toDBγ, whose threshold lies ∼ 45 MeV belowDB * . In either case, such a state, with a mass m udcb 7154(23) MeV, should be much more accessible to current experimental programs than are the corresponding doubly bottom states. 3 Should the binding be greater than ∼ 45 MeV, the udcb tetraquark will decay only weakly, and suffer from the experimental disadvantage of having a large number of exclusive decay modes, each with a relatively small branching fraction. Such decay modes would, however, have the compensating advantage of being accompanied by secondary heavy-meson decays with displaced vertices. The left and center panels of Fig. 3 show examples of such decays. Two-or three-meson decays such as those to (D 0D0 ) or π +D0 D − , might then serve as useful potential experimental signatures. Another three-body mode potentially suitable for detection, produced if the W + materializes as a cs rather than a ud pair, is J/ψD 0 K 0 . In contrast, should the udcb binding energy be less than ∼ 45 MeV, the state should decay electromagnetically essentially 100% of the time toDBγ, as in the rightmost panel of Fig.3.
With (given the predicted binding energy) less than ∼ 20 MeV of phase space available to the electromagnetic decay, such a tetraquark, produced with a reasonable momentum in the lab frame, will decay to produce a mixed-heavy-flavour D + B − orD 0 B 0 meson pair highly collimated in the lab frame.

VII. CONCLUSIONS
The study of heavy tetraquarks of the qq Q Q type continues to yield insight into the binding of exotic hadrons. This paper provides evidence of a strong-interaction-stable I(J P ) = 0(1 + ) udcb tetraquark which we estimate to lie 15 − 61 MeV below the corresponding free, two-meson (DB * ) threshold. The predicted mass range, ≈ 7154(23) MeV, thus straddles the electromagneticDBγ decay threshold, making it impossible, at present, to predict whether it will decay electromagnetically or only weakly. With B c production, and hence the simultaneous production of bb and cc pairs, already established experimentally, these results clearly motivate a search for this state at LHCb.
The results of the direct computation above further strengthen the argument that this class of tetraquarks becomes more strongly bound as the light-quark component becomes lighter. The variable, unphysical heavy-quark mass study also confirms the expectation that such tetraquarks become less bound as the heavy anti-diquark reduced mass decreases.
Complementary to the direct computation of the binding of the I(J P ) = 0(1 + ) udcb tetraquark, the NRQCD study of the udb b , udb b , sb b and sb b channels with variable heavy quark mass, m b , extending down to 0.60 m b confirms our qualitative understanding of the basic physics responsible for the observed tetraquark binding. This study also suggests that strong-interaction-stable udcc , scb , and scc tetraquarks are unlikely to exist, see also [40], though direct computations using a relativistic charm action at almost physical light quark masses have yet to be carried out for these channels.
Our prediction for the udcb binding has an error dominated by what we believe to be a conservative estimate of the finite-volume systematic. While statistics does not appear to be a problem, having more operators and/or crafting sources and sinks that better overlap with our desired ground states should provide longer plateaus. In future, futher investigations of the udcb channel using ensembles with additional near-physical light pion masses and larger spatial volumes will allow us to obtain better quantitative control of the extrapolation to physical m π and finite-volume effects. where U t denotes a gauge link in the temporal direction, is used to compute our heavy quark propagators. Static heavy quark propagators can be realised by setting all the coefficients c i in this action to 0.
The NRQCD action can be organized as an expansion in 1/m Q and naively its regime of validity is limited to the region m Q 1 in lattice units.

Appendix B: Details of the variable-heavy-mass study
This study focuses on unphysical bottom quarks using the ensemble E M of Tab. IV with a single source position. Visible effects due to violations of the condition am Q 1 are e.g. the increase of the hyperfine splitting. While tunings of the parameters c i beyond tree level, as well as the inclusion of additional terms in the action, and/or an increase in the stability parameter n, could be considered to expand the region of validity of the NRQCD approximation, here we choose to work with tree-level values and the form of the NRQCD Hamiltonian detailed in Appendix A. Monitoring the hyperfine splitting, we find that m Q 0.594 m b appears to be the lowest acceptable value for our calculation. From extrapolation of our heavy-quark data, the charm quark lies at around m c 0.33 m b and is thus beyond our reach, in the NRQCD approach.
The binding energies of the udb b , udb b , sb b and sb b tetraquark states on the E M ensemble are extracted by fitting the smallest eigenvalue of the respective GEVPs to a singleexponential ansatz. All fit results with χ 2 /d.o.f. ≤ 1 are kept for further analysis. The fit procedure does not take into account correlations in the data. The final numbers are chosen as representative results with the longest fit window possible.
The resulting data and fits are shown in the left panels of Fig. 4. We observe positive exponential growth of the first eigenvalues determined from the binding correlator with Euclidean time t. This growth becomes steeper as the heavy quark mass is increased in all four tetraquark channels (especially so for the equal-heavy-mass case, where the binding energy is unbounded as m b → ∞). The shaded bands denote the corresponding fits and are seen to describe the data well. We also determine the effective binding energies and show these results (data), together with the corresponding binding energy fits (shaded boxes) in the right panels of Fig. 4. The success of the phenomenological fit forms detailed in Sec. II in reproducing this lattice data confirms our understanding of the basic binding mechanism.
FIG. 5. udcb tetraquark results on E H (top), E M (center) and E L (bottom). The left panels show the GEVP eigenvalues, the right panels the corresponding energies. The red and blue symbols show, respectively, the ground and first excited state results obtained from the 3 × 3 GEVP analyses. The green vertical dashes and green diagonal crosses similarly denote the ground and first excited state results obtained from the 2 × 2 GEVP analyses. The two-mesonDB * thresholds, shown for comparison, are given by the brown bands in the right-hand panels. The grey bands in those same panels depict the final energies obtained from the single exponential fits to the eigenvalues described above. Further details are given in the text. 2 × 2 GEVP results have been offset in t for visual clarity.
are shown as a function of t lef t in the right panels of Fig. 6, and the corresponding results FIG. 6. Fit window dependence for fitting the binding correlator GEVP solutions and those obtained from the particle correlators for the ensembles E H (top), E M (center), and E L (bottom). We obtain our final results from this selection of accepted fits in the stable region that is representative of both procedures. obtained using the binding correlator GEVPs in the corresponding left panels. The multiple points plotted for each t lef t are those corresponding to different choices of ∆t. From this selection of results the final, quoted numbers (corresponding to the grey bands in Fig. 5) are chosen in such a way that they lie in the stable regions of the plots, and are representative of both procedures. The final fit ranges used are t/a ∈ [10 : 23] for E L , t/a ∈ [14 : 20] for E M , obtained from the left panels, and t/a ∈ [17 : 21] for E H chosen from the right panels.