Bound isoscalar axial-vector $bc\bar u\bar d$ tetraquark $T_{bc}$ from lattice QCD using two-meson and diquark-antidiquark variational basis

We report a lattice QCD study of the heavy-light meson-meson interactions with an explicitly exotic flavor content $bc\bar u\bar d$, isospin $I\!=\!0$, and axialvector $J^P=1^+$ quantum numbers in search of possible tetraquark bound states. The calculation is performed at four values of lattice spacing, ranging $\sim$0.058 to $\sim$0.12 fm, and at five different values of valence light quark mass $m_{u/d}$, corresponding to pseudoscalar meson mass $M_{ps}$ of about 0.5, 0.6, 0.7, 1.0, and 3.0 GeV. The energy eigenvalues in the finite-volume are determined through a variational procedure applied to correlation matrices built out of two-meson interpolating operators as well as diquark-antidiquark operators. The continuum limit estimates for $D\bar B^*$ elastic $S$-wave scattering amplitude are extracted from the lowest finite-volume eigenenergies, corresponding to the ground states, using amplitude parametrizations supplemented by a lattice spacing dependence. Light quark mass $m_{u/d}$ dependence of the $D\bar B^*$ scattering length ($a_0$) suggests that at the physical pion mass $a_0^{phys} = +0.57(^{+4}_{-5})(17)$ fm, which clearly points to an attractive interaction between the $D$ and $\bar B^*$ mesons that is strong enough to host a real bound state $T_{bc}$, with a binding energy of $-43(_{-7}^{+6})(_{-24}^{+14})$ MeV with respect to the $D\bar B^*$ threshold. We also find that the strength of the binding decreases with increasing $m_{u/d}$ and the system becomes unbound at a critical light quark mass $m^{*}_{u/d}$ corresponding to $M^{*}_{ps} = 2.73(21)(19)$ GeV.

= +0.57(+4 −5 )(17) fm, which clearly points to an attractive interaction between the D and B * mesons that is strong enough to host a real bound state T bc , with a binding energy of −43( +6 −7 )( +14 −24 ) MeV with respect to the D B * threshold.We also find that the strength of the binding decreases with increasing m u/d and the system becomes unbound at a critical light quark mass m * u/d corresponding to M * ps = 2.73 (21) (19) GeV.
The discovery of a doubly charmed tetraquark1 , T cc , marks an important milestone [1] in spectroscopy of hadrons.Phenomenologically, doubly heavy tetraquarks in the heavy quark limit are long hypothesized to form deeply bound states [2][3][4][5][6][7][8][9][10][11][12] with binding energy O(100 MeV) with respect to the elastic strong decay threshold.While doubly bottom tetraquarks are suitable candidates for such deeply bound states, as predicted by multiple lattice QCD calculations [13][14][15][16][17][18], T cc is found to be 360 keV below the lowest two-meson threshold (D 0 D * + ).A handful of recent experimental developments involving multiple heavy quark production such as the recent discoveries of Ξ cc [19], T cc [1], reports of tri-J/ψ [20], associated J/ψΥ [21], and di-Υ [22] productions, and recent proposals of inclusive search strategies [23,24] augment promising prospects for the doubly heavy hadron sector in the near future.In light of these advancements, a doubly heavy tetraquark with a bottom and a charm quark with a valence quark configuration T bc ≡ bcū d is going to be one of the most sought-after hadron in this decade [25].In this work, using lattice QCD calculations, we show a clear evidence of an attractive interaction between the D and B * mesons that is strong enough to host a real bound state T bc .This finding will further boost the search for such bottom-charm tetraquarks.
The phenomenological picture on deeply bound doubly heavy tetraquarks is based on a compact heavy diquarklight antidiquark interpretation [14,26], whereas the shallow binding energy of T cc could possibly be a reflection of its dominant noncompact molecular nature [8,27].Bottom-charm tetraquarks form an intermediate platform, where there could be complicated interplay between these pictures.A collective and refined knowledge of the low energy spectra in all these three doubly heavy systems (T bb , T bc and T cc ) could culminate in a deeper understanding of strong interaction dynamics across a wide quark mass regime spanning from charm to the bottom quarks.The isoscalar bottom-charm tetraquarks with quantum numbers [I(J P ) = 0(1 + )] have been investigated previously both using lattice [28][29][30] and nonlattice methodologies [5,6,8,9,11,12,26,[31][32][33][34][35][36][37][38][39][40].The predictions from nonlattice approaches are quite scattered from being unbound to deeply bound, whereas the difference in conclusions from the three existing lattice QCD investigations [28][29][30] call for more detailed efforts in this regard.
In this work, we perform a lattice QCD simulation of coupled D B * and BD * two-meson channels 2 that are the relevant lowest two strong decay thresholds, in the order of increasing energies, E D B * = M B * + M D and E BD * = M B + M D * , where M h is the mass of the hadron h.The extracted finite-volume ground state energies are utilized to constrain the continuum extrapolated elastic D B * scattering amplitudes following the Lüscher's finitevolume prescription [41,42] The lattice spacing estimates are measured using the r1 parameter [43].L1 refers to large spatial volume, and S1, S2, and S3 refer to small spatial volumes.
Lattice setup: We use four lattice QCD ensembles (see Table I for relevant details) with N f = 2+1+1 dynamical Highly Improved Staggered Quark (HISQ) fields generated by the MILC collaboration [43].The charm and strange quark masses in the sea are tuned to their respective physical values, whereas the dynamical light quark masses correspond to sea pion masses as listed in Table I.We utilize a partially quenched setup on these configurations with valence quark fields up to the charm quark masses realized using an overlap fermion action as in Refs.[44,45].We employ a nonrelativistic QCD (NRQCD) Hamiltonian [46] for the bottom quark.Following the Fermilab prescription [47], the bare charm [48,49] and bottom [50] quark masses on each ensemble are tuned using the kinetic mass of spin averaged 1S quarkonia {aM QQ kin = 3  4 aM kin (V ) + 1 4 aM kin (P S)} determined on the respective ensembles.The bare strange quark mass is set by equating the lattice estimate for the fictitious pseudoscalar ss meson mass to 688.5 MeV [51].
For the valence m u/d , we investigate five different cases: three unphysical quark masses [corresponding to approximate pseudoscalar meson masses M ps ∼0.5, 0.6, and 1.0 GeV], the strange quark mass [M ps ∼0.7 GeV] and the charm quark mass [M ps ∼3.0 GeV].We evaluate the finite-volume spectrum for all these five quark masses on all four ensembles, investigate the scattering of D and B * mesons in all five cases and then extract the m u/d (otherwise M ps ) dependence of the scattering parameters.
Interpolators and measurements: The finite-volume spectrum is determined from Euclidean two-point correlation functions C ij (t), between interpolating operators O i,j (x, t) with desired quantum numbers, given by Here E n is the energy of the n th state and Z n i = ⟨0| O i |n⟩ is the operator-state overlap between the sink operator O i and state n.We use O and Z to represent the source operator and overlaps to distinguish them from that for the sink as we follow a wall-source to pointsink construction in our C ij evaluations.This is wellestablished procedure in ground state energy determination, despite the non-Hermitian setup in Eq. (1) (see Refs. [14,16,28,29,52,53] for details).We use the following set of linearly independent, yet Fierz related [54], operators, Bilocal two-meson interpolators with nonzero internal meson momenta are also not considered, which would be an important step ahead [55].We also compute two-point correlation functions for B, B * , D, and D * mesons, using standard local quark bilinear interpolators (q Γ Q) with spin structures Γ ∼ γ 5 and γ i for pseudoscalar and vector quantum numbers, respectively.Analysis: The correlation matrices C evaluated for the basis in Eq. ( 2) are analyzed following a variational procedure [56] by solving the generalized eigenvalue problem (GEVP), C(t)v n (t) = λ n (t)C(t 0 )v n (t).The eigenvalues in the large time limit represent the time evolution the low lying eigenenergies E n as lim t→∞ λ n (t) ∼ A n e −E n t .The corresponding eigenvectors v n (t) are related to the operator-state-overlaps Z n i .Eigenenergy extraction proceeds via fitting the eigenvalue correlators, λ n (t), or the ratios R n (t) = λ n (t)/C m1 (t)C m2 (t), with the expected asymptotic exponential behaviour.Here, C mi is the two-point correlation function for the meson m i .R n (t) is empirically known to efficiently mitigate correlated noise between the product of two single hadron correlators and the interacting correlator for the two-hadron system [57].Note that the automatic cancellation of the additive quark mass offset, inherent to NRQCD formulation, is an added advantage in using R n (t) for the fits.The systematics associated with the chosen time interval for fitting are assessed by varying the lower boundary of the time interval, t min , with a fixed upper boundary, t max , chosen considering the noise level.In Figure 1, we present a representative plot showing this t min dependence of the energy splittings (∆E n ) determined from the fits to λ n (t) and R n (t), respectively.The energy differences are evaluated from λ n (t) using the relation ∆E n = E n − M m1 − M m2 , whereas the fits to R n (t) directly yield the respective estimates.We choose the optimal t min values where the two different procedures found to agree asymptotically in time.We also perform additional checks considering an alternative quark smearing with different smearing widths to affirm our energy estimates, see Appendix A. Our final results are based on fitting the ratio correlators R n (t).

Finite-volume eigenenergies:
In Figure 1, we present the finite-volume GEVP eigenenergies, in lattice units, for the isoscalar axialvector bcū d channel.The results shown are for the L 1 ensemble at the five different m u/d values corresponding to M ps ∼ 0.5, 0.6, 0.7, 1.0, and 3.0 GeV.Note that these estimates include the additive offsets related to the NRQCD-based bottom quark dynamics.The non-interacting two-meson energy levels corresponding to D B * and BD * thresholds are indicated as dotted horizontal line segments and those related to B * D * threshold by dashed lines for each M ps .The lowest eigenenergy or the ground state energy is dominated by the Z 0 1 factor corresponding to O 1 , that is related to the D B * threshold and is determined unambiguously by the operator O 1 , see Ref. [58] for details.The most important observation is a clear trend for negative energy shifts in the ground state energies, which can be observed in all the cases, indicating a possible attractive interaction between the D and B * mesons [59].A similar pattern of low lying eigenenergies and ground state negative energyshifts are also observed in other ensembles, see details in Ref. [58].We expect that for our choice of interpolating operators and the accessible values of t, E 0 will be an accurate estimate of E 0 , whereas our setup is unable to accurately estimate excited-state energies.This means the excited eigenenergies presented in Figure 1 may not correspond to the higher lying elastic excitations of the D B * channel.The location of lowest two non-interacting finite-volume levels related to the D B * channel along with the ground state eigenenergies are presented in Appendix B. Hence we focus only on the ground state energies (E 0 ∼ E 0 ) for the rest of the analysis.In Figure 3, we present the ground state energy estimates, in units of E D B * , at various M ps and for all the ensembles.These estimates are evaluated as ) refers to the spin averaged mass of the 1S bottomonium measured on the lattice (experiments).The eigenenergies clearly show a trend of decreasing energy spitting, hence decreasing interaction strength, with increasing M ps .Another interesting feature to note here is the nonzero lattice spacing (a) dependence of the ground state energies on similar volume ensembles (S 1 , S 2 , S 3 ), which we account for through an a dependence in the parametrized amplitude as discussed below.
In Figure 3, we also indicate the branch point location of the left hand cut (lhc) arising out of an off-shell pion exchange process for different M ps by horizontal dashed lines.Recent developments point to the importance of lhc effects on virtual subthreshold poles related to the T cc tetraquark [60].Such effects on bound states are the subject of future studies where one could successfully solve the relevant three particle integral equations.This is beyond the scope of this work and we ignore such effects in our analysis.
D B * scattering amplitude: Assuming these energy splittings in ground states are purely described by elastic scattering in the D B * system, we utilize them to con-strain the associated S-wave scattering amplitude following Lüscher's finite-volume prescription [41,42].For the low energy scattering of D and B * mesons, where other multi-particle thresholds are sufficiently high [61], in the S-wave leading to the total angular momentum and parity J P = 1 + , the scattering phase shifts δ l=0 (k) are related to the finite-volume energy spectrum through kcot[δ 0 ).We follow the procedure outlined in Appendix B of Ref. [62] to constrain the amplitude.A sub-threshold pole in the S-wave scattering amplitude t = (cotδ 0 − i) −1 occurs when kcotδ 0 = ± √ −k 2 for scattering in S-wave.We parametrize the elastic D B * scattering amplitude in terms of the scattering length a 0 in an effective range expansion near the threshold, supplemented by a lattice spacing dependence.This is required to incorporate the cut-off effects observed in the ground state energy estimates.We find that a linear functional form given by kcotδ 0 = A [0] + aA [1] , where A [0] = −1/a 0 , accommodate the a dependence of the kcotδ 0 estimates.We present the fit results for A [0] = −1/a 0 in Figure 4 (circle symbols) as a function of M ps involved.Alternative fitting choices with a leading quadratic dependence or using only data from non-charm M ps are consistent with results in Figure 4, see details in Supplemental material [58].
The sign of A [0] = −1/a 0 determines the fate of the near-threshold pole, if there exists one.A negative (positive) value of A [0] (a 0 ) indicates that the interaction potential is strong enough to form a real bound state [63].After considering all possible systematics, we find that for all the non-charm light quark masses, A [0] is negative, which indicates an attractive interaction strong enough to host a real bound state.On the contrary, at the charm point, despite the unambiguous negative energy shifts in the ground states, the attraction is weak to host any real bound state as suggested by the positive value of kcotδ 0 in the continuum limit.This observation goes in line with the phenomenological expectation for doubly heavy four quark (QQ ′ l 1 l 2 ) systems with m l1 = m l2 that the binding increases with increased relative heaviness of the heavy quarks with respect to its light quark content [14,26,64].Now we investigate the light quark mass (m u/d ) or M ps dependence of the fitted parameters.To this end, we consider three different parametrizations: a linear dependence (f l (M ps ) = α c + α l M ps ) to probe the heavy light quark mass case, a leading M 2 ps dependence (f s (M ps ) = β c +β s M 2 ps ) to assess the chiral behaviour, and a quadratic dependence (f q (M ps ) = θ c + θ l M ps + θ s M 2 ps ) to quantify the associated systematics.In Figure 4, we show the fit results for this M ps dependence in colored bands.The two stars represent A [0] at the physical M ps (equivalently the physical scattering length a phys 0 ) and the critical M ps at which A [0]   dependence.Yet, our fits demonstrate near independence in the fit forms as can be observed from the consistency between the error bands from different fit forms.
Based on the fit form f s (M ps ) in the chiral regime, we find that the scattering length of the D B * system at the physical light quark mass (m phys u/d ), corresponding to M ps = M phys π , to be The asymmetric errors indicate the statistical uncertainties, whereas the second parenthesis quotes the systematic uncertainties with the most dominant contribution arising from the chiral extrapolation fit forms.The positive value of the scattering length at M ps = M phys π , at the level of 3σ uncertainty, is an unambiguous evidence for the strength of the D B * interaction potential to host a real bcū d tetraquark bound state T bc with binding energy with respect to E D B * .The first parenthesis indicates the statistical errors and the second one quantifies various systematic uncertainties added in quadrature.The pseudoscalar meson mass, corresponding to the critical light quark mass, where a 0 diverges, is found to be M * ps = 2.73(21) (19) GeV.This critical point also signifies that QCD dynamics within such exotic systems is such that at a heavy light quark mass the system of quarks perhaps reaches the unitary gas limit, as indicated by the divergent scattering length [65].For M ps ≥ M * ps , the T bc system remains unbound.
Systematic uncertainties: Our lattice setup together with the bare bottom and charm quark mass tuning procedure has been demonstrated to reproduce the 1S hyperfine splittings in quarkonia with uncertainties less than 6 MeV [50,53].We observe the effects of such a mistuning of either of the heavy quark mass on the energy splittings we extract are very small compared to the statistical errors.Additionally, our strategy of evaluating the energy differences and working with mass ratios has also been shown to significantly mitigate the systematic uncertainties related to heavy quark masses [52,53].This is observed to be the case in this study as well, leading to transparent signals for the ground state energy as shown in Figures 1, 1, and 3. Our fitting procedure involves careful and conservative determination of statistical errors, and uncertainties related to the excited-state-contamination and fit-window errors.Additional checks using alternative quark smearing procedures also agree with our energy estimates, see Appendix A. The amplitude determination and followed extrapolations are performed with results from varying the fit-windows to evaluate the uncertainties propagated to our final results.The uncertainties related to the fit forms used in chiral extrapolations are observed to be the most dominant, as is evident from Figure 4. We assume the partially quenched setup involving ensembles with different sea pion masses, we utilize, have negligible effects on the energy splittings we extract for the explicitly exotic T bc tetraquark, similar to what was observed for heavy hadrons in Refs.[66,67].Uncertainty related to scale setting is also found to be negligible in comparison to the statistical uncertainties in the energy splittings.
Summary: We have performed a lattice QCD simulation of coupled D B * -BD * scattering with explicitly exotic flavor bcū d and I(J P ) = 0(1 + ).Following a rigorous extraction of finite-volume eigenenergies and continuum extrapolated elastic D B * scattering amplitudes for the five light quark masses studied, we determine the light quark mass dependence of the elastic D B * scattering length a 0 .We observe unambiguous negative energy shifts between the interacting and non-interacting finitevolume energy levels.Our estimate for a phys 0 (Eq.( 3)) is positive, indicating an attractive interaction between the D and B * mesons, which is strong enough to host a real bound state with binding energy δm T bc = −43( +6 −7 )( +14 −24 ) MeV.We find that the strength of interaction is such that this bcū d tetraquark becomes unbound at M * ps , which is close to the η c meson mass.
In this work, we make several important steps ahead to arrive at robust inference on the nature of interaction between the D and B * mesons.Our main strategy has been to determine the signature of scattering length in D B * interactions at the physical pion mass a phys 0 .Our results indicate that a phys 0 is positive, which suggests that attractive D B * interactions are strong enough to host a real bound state.Further theoretical investigations are desired to reduce the uncertainties in the binding energy of T bc with respect to E D B * .Fully dynamical simulations on several more ensembles, with different volumes and improvized fermion actions, high statistics studies with lighter m u/d , etc. are a few other improvisations that can further constrain the relevant scattering amplitude.Additionally, future works involving Hermitian correlation matrices at rest as well as in moving frames and those using bilocal two-meson interpolators with nonzero relative meson momenta aimed at reliable excited state extraction would be a few important steps ahead [55,62,68,69].We hope that our observations and inferences in this work will motivate more theoretical efforts and experimental searches for such states.
This work is supported by the Department of Atomic Energy, Government of India, under Project Identification Number RTI 4002.M.P. gratefully acknowledges support from the Department of Science and Technology, India, SERB Start-up Research Grant No. SRG/2023/001235.We are thankful to the MILC collaboration and in particular to S. Gottlieb for providing us with the HISQ lattice ensembles.We thank Sara Collins for a careful reading of the manuscript.We thank the authors of Ref. [70] for making the TwoHadronsInBox package utilized in this work.We also thank Gunnar Bali, Parikshit Junnarkar, Alexey Nefediev, Sayantan Sharma, Stephen R. Sharpe, and Tanishk Shrimal for discussions.Computations were carried out on the Cray-XC30 of ILGTI, TIFR.Amplitude analyses were performed on Nandadevi computing cluster at IMSc Chennai.N. M. would also like to thank A. Salve and K. Ghadiali for computational support.
Appendix A: Ground state energy plateau.-Inthis work we have utilized a wall-source point-sink setup to construct the necessary two-point correlation functions.The use of such an asymmetric setup implies the effective energies aE ef f = [ln(C(t)/C(t + δt))]/δt could approach their asymptotic values as rising-from-below, due to the nonpositive definite nature of the coefficients in a spectral decomposition, in contrast to a falling-from-above behaviour in a symmetric setup.In Figure 5, we show the effective mass in wall-source point-sink setup with the brown-circle (R 2 = 0) which rises from below.
To avoid any ambiguity in selecting the plateau regions of effective masses of such correlators, we also employ a wall-source box-sink setup [29], which asymptotically approaches the symmetric limit.In the symmetric limit, the effective masses are expected to follow a conventional falling-from-above feature, modulo the statistical noise.To this end, we vary the smearing radius R to investigate the time dependence of effective mass plateaus in the approach to the symmetric limit.In Figure 5, we present a comparison of the effective energy (top) and effective energy splittings (bottom) determined using different quark sink smearing procedures for the case of M ps ∼ 700 MeV on the finest ensemble.Clearly the rising-from-below behaviour is gradually disappearing in the approach to the symmetric limit.It is also evident that the results at the large time limit from point-sink and box-sink are very much consistent with each other affirming our assessment on effective mass plateau in choosing a fit range.Such a behavior of effective masses with varying smearing radii was also observed in Ref. [29].In the large time limit, where the signal quality is still good, all of sink smearing cases suggest consistent negative energy shifts.This is evident from the large time behaviour of energy splittings presented in the bottom panel of Figure 5, where the correlated statistical noise, not related to the excited state contamination, is suppressed between the numerator and denominator in the ratio correlators R n (t).
The agreement of energy splitting estimates from fits to R n (t) with those evaluated from separate fits to the GEVP eigenvalue correlators λ n (t) and the single-meson correlators C D/ B * at large times (see Figure 1) already rules out the usual concern of accidental partial cancellation of excited state contaminations in R n (t).The consistency at large times between ground state energy plateaus from different sink-smearing radii observed in top panel of Figure 5 further affirms the reliable isolation of the ground state plateau.Note also that the magnitude of such cancellations and the ground state saturation times could be different in different lattice QCD ensembles.All the ground state estimates for noncharm M ps values in our study are determined from the time intervals approximately between 1.5 (2)   Appendix B: Elastic D B * excitations.-Gainingaccess to higher lying elastic excitations in the D B * channel is an important step ahead towards constraining the energy dependence of the amplitude over a long energy range.However, within the wall-smearing setup, all the nonzero momentum excitations are significantly suppressed.This suppression is exact in a free theory, and is empirically confirmed from the early plateauing and from the quality of signals in the interacting theory.While this suppression is advantageous in ground state energy determination (see Refs. [14,16,28,29,52,53] for details), the suppressed coupling to the nonzero momentum excitations implies that the access to higher two-meson elastic excitations with nonzero relative meson momenta are restricted in the wall-smearing setup.This implies other methodologies that facilitate the use of bilocal two-meson interpolators with separately momentum projected mesons are necessary in future studies [55,71,72] 3 .In this respect, it is informative to know the location of the lowest noninteracting level with nonzero relative meson momenta and whether it is close enough to influence the ground state energies in any substantial way.Considering this, in Figure 6 we present the ground state eigenenergies along with the D B * threshold and the next lowest elastic D B * excitation with nonzero relative meson momentum determined using the continuum dispersion relation that is assumed in the finite-volume quantization condition [41,42].Clearly, the location of this first non-interacting elastic excitation is sufficiently high to have any nonnegligible effects on the extracted the ground state energies.In this section, we present the finite-volume GEVP eigenenergies of the isoscalar axialvector bcū d channel that we extract on all four ensembles listed in Table I of the main article, at the five different m u/d values corresponding to M ps ∼ 0.5, 0.6, 0.7, 1.0, and 3.0 GeV.The eigenenergies shown in lattice units include the additive offsets related to the NRQCD-based dynamics of heavy bottom quarks.The non-interacting two-meson energy levels corresponding to D B * and BD * thresholds are indicated as dotted horizontal line segments for each lattice and each M ps .The D * B * threshold in each case is also shown in the figure by dashed lines.Note that the use of wall-smearing setup restricts any direct access to the elastic two-meson excitations with nonzero relative meson momenta.This means although a reliable ground state extraction could be made, the excited eigenenergies may not represent the real elastic excitations.

II. OPERATOR-STATE OVERLAPS
In Figure 2, we present the modulus of normalized sink operator-state overlaps | Zn i |, normalized such that its largest value for any given operator O i across all the eigenenergies {n} is unity [74,75].Zn i quantifies the relative relevance of any given operator across all the eigenenergies.The | Zn i | values are presented for all M ps cases on the L 1 ensemble.Each square marker corresponds to the | Zn i | for a given operator O i on to a given eigenenergy n.Each horizontal panel stands for an M ps indicated on the right-hand side, whereas the vertical lines in each horizontal panel part | Zn i | for different operators indicated on the top panel.The x-axis ticks refer to the three finite-volume eigenenergies we have extracted.O 1 , the two-meson operator related to D B * threshold, can be seen to have the largest overlap with the ground state and has significantly small overlaps with the excited eigenenergies.O 2 , the two-meson operator related to BD * threshold, has the largest overlap with the first excited eigenenergy and a very small overlap with the ground state.O 2 also have nonnegligible overlap factors with the second excited eigenenergy indicating BD * -type two-meson Fock component, which decreases with increasing M ps .On the other hand, O 3 , the diquark-antidiquark type operator, have substantial overlap factors with all eigenenergies at the two lightest M ps values, whereas with an increased M ps its largest overlap is with the second excited eigenenergy.Note that O 3 is Fierz related to two-meson interpolators [54], and the large Zn 3 values of O 3 for all n could be related to this underlying connection between two-meson and diquark-antidiquark operators.
A summary from the above observations on overlap factors is as follows.O 1 predominantly determines the ground state, whereas it has significantly small coupling with the excited eigenenergies.One could also evaluate and investigate the normalized source operator-state overlaps from the left eigenvectors of C in Eq. (1) in the main draft, which also leads to the same conclusions.In Figure 3, we show the operator basis dependence as determined for M ps ∼ 700 MeV in the L 1 ensemble, for various operator basis build out of O 1 , O 2 , and O 3 operators as defined in Eq. (2) of the main draft.The digital indexing on the x-axis tick labels refers to various operator basis in the order {O 1 , O 2 , O 3 }, with an overline on the third index as a visual aid within the plot to highlight the diquark-antidiqaurk interpolator.1 (0) indicates an operator is included in (excluded from) the basis.The horizontal lines refer to the D B * , BD * and B * D * thresholds.The gray horizontal bands refer to the two lowest levels in the full basis indicated by 111.A level below the threshold appears only when O 1 is present in the basis.The first excited eigenenergy in the full basis 111 is faithfully reproduced in those bases where O 2 is included.O 3 alone does not precisely determine any level among the GEVP eigenenergies using full basis.Similar observations are also made on other ensembles.In summary, the ground state in the full basis 111 is reliably determined with O 1 and is unaffected by the inclusion of O 2 and O 3 operators.

IV. RESULTS ON SCATTERING AMPLITUDE
In Table I, we tabulate the results from different amplitude fits that were performed.In Figure 4, we present the quality of these fits by comparing the fit results with the data points (see the figure caption for details).On the left of Figure 5, we present the light quark mass dependence in the chiral regime determined from continuum extrapolated elastic D B * scattering amplitudes following a linear and quadratic lattice spacing dependence.We present the final estimate (black star) from the linear fit form considering the presence of heavy quarks in our system, whereas the difference in the fit results are accounted in the systematics quoted in the main draft.On the right of Figure 5, we present a comparison of the light quark mass dependence in the chiral regime between the fit involving all M ps datasets and the fit involving the lightest four M ps datasets.Either fitting procedures can be seen to be consistent with our final estimate in the chiral limit is shown by black star.  5. Left: Comparison of light quark mass dependence of scattering amplitudes in the chiral regime determined from a linear (red band) and quadratic (green band) dependence of kcotδ0 on the lattice spacing.Right: Comparison of light quark mass dependence of scattering amplitude kcotδ0 in the chiral regime determined using the results from all five light quark masses (red band) and the results from four light light quark masses (black).The results in the physical limit in either cases can be seen to be consistent with the main result indicated by the star symbol.

FIG. 3 .
FIG. 3. The ground state energies in units of E D B * on all ensembles (see Table I for color-symbol conventions) for allMps values (different vertical panels).
where ∆E n is the estimate from fit to R n (t), M B * = M B * − 0.5M bb lat +0.5M bb phys accounts for the NRQCD additive offset, and M bb lat (M bb phys

FIG. 4 .
FIG. 4. Continuum extrapolated kcotδ0 or A [0] = −1/a0 estimates of the D B * system as a function of M 2 ps in units of E D B * .The dotted vertical line close to the y-axis indicates Mps = M phys π .The two star symbols represent the amplitude at Mps = M phys π and the critical Mps = M * ps above which the system becomes unbound.

49 FIG. 5 .
FIG.5.Comparison of effective energy (top) and effective energy splitting (bottom) for the ground state as determined using three different smearing radii applied on the quark fields at the sink timeslice.The legend indicates the smearing radius squared in units of the lattice spacing[29].The blue horizontal band indicates the final fit estimate for the energy and energy splitting.The results presented are for the case Mps ∼ 700 MeV on the finest ensemble.

FIG. 6 .
FIG.6.The ground state energy eigenvalues in the background of lowest two non-interacting D B * finite-volume levels units of E D B * on all ensembles (see TableIfor color-symbol conventions) for all Mps values (different vertical panels).

D
FIG.1.Low-lying finite-volume eigenenergies for isoscalar axialvector bcū d channel on the four ensembles studied.The near degenerate eigenlevels are slightly shifted horizontally for clarity.

FIG. 3 .
FIG.3.Operator dependence of the low lying eigenenergies of the L1 ensemble and Mps ∼700 MeV for all possible operator basis that can be built out of the three operators utilized in this work.

FIG. 4 .
FIG. 4. Left: kcotδ0, in units of the elastic threshold E D B * , versus a (lattice spacing) for all Mps values.We follow the marker/color coding in Table I of the main draft for the data points referring to the simulated data.The colored/gray bands indicate the fit results with linear and quadratic lattice spacing dependence, respectively.Right: kcotδ0 versus k 2 for all Mps values studied in units of the elastic threshold E D B * .The dashed orange (cyan) curve indicates the constraint for the existence of a sub-threshold pole in the scattering amplitude.The horizontal bands are the continuum extrapolated estimates of kcotδ0 for the respective Mps.The black dashed vertical lines in the plots on the right indicate the location of the branch point associated with the left hand cut arising from the D Bπ channel.

TABLE I .
. The light quark mass m u/d dependence of the extracted amplitudes suggests a binding energy of −43( +6 −7 )( +14 −24 ) MeV for the bcū d tetraquark pole with respect to E D B * at the physical point m phys u/d .Relevant details of the lattice QCD ensembles used.
). O 1 and O 2 are two-meson operators of the type D B * and BD * , respectively.O 3 is a diquark-antidiquark type operator.Here Γ k = Cγ k with C = iγ y γ t being the charge changes its sign or above which the system becomes unbound.It is indeed desired to have more points in the intermediate mass regime between the charm and the strange quark masses to further constrain the fm [t min ] to 2.3(2) fm [t max ].The consistent ground state saturation times across different ensembles with different specifications further imply the reliability of our ground state saturation, despite our asymmetric setup.
Similar patterns of overlap factors are also observed for other ensembles, all of which indicate that O 1 predominantly determines the ground state.The two excited eigenenergies have strong two-meson and diquark-antidiquark Fock components in the two lightest M ps values.FIG.2. Modulus of normalized sink operator-state overlaps | Zn i | for an eigenenergy indicated by n = 0, 1, 2 and an operator represented by Oi, where i = 1, 2, 3 on the L1 ensemble.The errors in the normalized overlap factors are smaller than the symbols and hence are suppressed.

TABLE I .
Results from amplitude fits at different light quark mass cases indicated in terms of Mps in the first column.The amplitude is approximated to be determined by the scattering length, with a linear or quadratic lattice spacing dependence as discussed in the main draft.The optimized parameter values in the table are presented in units of energy of the D B * threshold, E D B * .The A[0]parameter in either parametrization is negative of the inverse scattering length in the continuum limit.