Shallow Bound States and Hints for Broad Resonances with Quark Content $\bar{b}\bar{c}ud$ in $B$-$\bar{D}$ and $B^*$-$\bar{D}$ Scattering from Lattice QCD

We present the first determination of the energy dependence of the $B$-$\bar{D}$ and $B^*$-$\bar{D}$ isospin-0, $S$-wave scattering amplitudes both below and above the thresholds using lattice QCD, which allows us to investigate rigorously whether mixed bottom-charm $\bar{b}\bar{c}ud$ tetraquarks exist as bound states or resonances. The scattering phase shifts are obtained using L\"uscher's method from the energy spectra in two different volumes. To ensure that no relevant energy level is missed, we use large, symmetric $7 \times 7$ and $8 \times 8$ correlation matrices that include, at both source and sink, $B^{(*)}$-$\bar{D}$ scattering operators with the lowest three or four possible back-to-back momenta in addition to local $\bar{b}\bar{c}ud$ operators. We fit the energy dependence of the extracted scattering phase shifts using effective-range expansions. We observe sharp peaks in the $B^{(*)}$-$\bar{D}$ scattering rates close to the thresholds, which are associated with shallow bound states, either genuine or virtual, a few MeV or less below the $B^{(*)}$-$\bar{D}$ thresholds. In addition, we find hints for resonances with masses of order $100$ MeV above the thresholds and decay widths of order $200$ MeV.

We present the first determination of the energy dependence of the B-D and B * -D isospin-0, S-wave scattering amplitudes both below and above the thresholds using lattice QCD, which allows us to investigate rigorously whether mixed bottom-charm bcud tetraquarks exist as bound states or resonances.The scattering phase shifts are obtained using Lüscher's method from the energy spectra in two different volumes.To ensure that no relevant energy level is missed, we use large, symmetric 7×7 and 8×8 correlation matrices that include, at both source and sink, B ( * ) -D scattering operators with the lowest three or four possible back-to-back momenta in addition to local bcud operators.We fit the energy dependence of the extracted scattering phase shifts using effective-range expansions.We observe sharp peaks in the B ( * ) -D scattering rates close to the thresholds, which are associated with shallow bound states, either genuine or virtual, a few MeV or less below the B ( * ) -D thresholds.In addition, we find hints for resonances with masses of order 100 MeV above the thresholds and decay widths of order 200 MeV.
The majority of experimentally observed mesons can be understood in the quark model as quark-antiquark pairs.However, mesons, which are hadrons with integer spin, can in principle also be composed of two quarks and two antiquarks.The existence of these socalled tetraquarks had already been proposed in the early history of the quark model and QCD [1][2][3], but clear experimental confirmation was obtained only around a decade ago, for example in form of the observation of the charged Z c and Z b states as reviewed in Refs.[4,5].While the masses and decays of the latter strongly indicate the presence of a cc pair or a bb pair, their non-vanishing electric charge implies additionally a light quark-antiquark pair.Recently, there was another experimental breakthrough in the field, namely the detection of the T cc tetraquark with quark flavors ccud by LHCb [6,7].In contrast to previously observed tetraquarks and tetraquark candidates, its mass is slightly below the lowest meson-meson threshold, making it by far the longest-lived experimentally confirmed tetraquark.Following the observation of this doubly-charm tetraquark, possible next targets for experimental searches could be mixed bottom-charm tetraquarks with flavor content bcud.Their production cross section at the LHC is estimated to be about 40 times larger compared to the doubly-bottom bb ud tetraquark [8].The experimental signatures of a tetraquark are completely different depending on whether its mass is above or below the lowest strong-decay threshold.Thus, reliable theoretical predictions concerning bcud tetraquarks are very important and also urgent.
On the theoretical side, for the lightest bb ud tetraquark with I(J P ) = 0(1 + ) (which is the bottom-quark partner of the previously mentioned T cc ), there is a consensus from recent lattice-QCD calculations that it is deeply bound [9][10][11][12][13][14][15] and will decay through the weak interaction only (see Refs. [8,16,17] for discussions of possible decay modes).For the case of bcud, there is no such consensus.After finding initial hints for a possible QCDstable bcud bound state with I(J P ) = 0(1 + ) from lattice QCD [18], the same authors refined their calculation with larger lattice sizes and other improvements, and the hints disappeared [19].In Ref. [13], some of us also performed lattice-QCD calculations of the bcud energy spectra for both I(J P ) = 0(1 + ) and I(J P ) = 0(0 + ), and we likewise did not find any evidence for QCD-stable bound states (although we could not rule out a shallow bound state).In contrast, another independent group very recently reported an I(J P ) = 0(1 + ), bcud bound state 43 +7 −6
In the following, we present a new lattice-QCD study of the bcud systems with both I(J P ) = 0(1 + ) and I(J P ) = 0(0 + ).This study uses a different lattice setup and substantially more advanced methods compared to previous work, allowing us to apply the Lüscher method to multiple excited states in addition to the ground state and hence to reliably determine the detailed energy dependence of the B-D and B * -D isospin-0, S-wave scattering amplitudes.
In lattice QCD, the low-lying finite-volume energy levels with a given set of quantum numbers (the total spatial momentum, the quark flavor content, and the irreducible representation of the full octahedral group) are extracted from numerical results for imaginary-time two-point correlation functions The operators O i are constructed out of quark and gluon fields such that they excite states with the desired quantum numbers, which resemble the low-lying energy eigenstates of interest.For an infinite (in practice, large) time extent of the lattice, the two-point function is equal to where |Ω⟩ is the vacuum state and the sum is over all eigenstates |n⟩ of the finite-volume QCD Hamiltonian for which the product of overlap matrix elements is nonzero.By analyzing the time dependence of the numerical results for C ij (t), the energies E n can be extracted.Because lattice QCD uses a Monte-Carlo sampling of the Euclidean path integral, the numerical results have statistical uncertainties.Moreover, these uncertainties typically grow exponentially with t.
For multi-quark systems, experience has shown that the simplest possible operator choices in which the quark fields are combined at the same spacetime point ("local" operators) are often insufficient to reliably extract even just the ground state [43].The reason is that all or most of the energy levels resemble multi-hadron states with specific relative momenta, and the spectrum of such states in the case of heavy-quark systems is particularly dense.Among the previous lattice studies of bcud systems, Refs.[18][19][20] used only local four-quark operators with various types of smearing (local, wall, box) applied to each quark.Reference [13] improved upon this by including also two-meson (B-D and B * -D) "scattering" operators, that is, operators with each meson individually projected to a specific momentum (equal to zero only, in this case).These operators were included at the sink only, to avoid having to generate expensive all-toall light-quark propagators.The work presented in the following no longer makes this restriction and is the first lattice-QCD calculation of bcud correlation matrices with B ( * ) -D scattering operators at both source and sink, and also the first to include B ( * ) -D scattering operators with nonzero back-to-back momenta.
Specifically, to study the bcud system with I(J P ) = 0(0 + ), we use seven operators are operators with all four quarks at the same spacetime point (but with Gaussian smearing of the quark fields) and jointly projected to zero total spatial momentum, and O ), and √ 2•2π/L.Two different operators are used for the case with one unit of back-to-back momentum to account for the mixing of S and D partial waves [44].The labels A + 1 and T + 1 refer to the octahedral-group irreps of positive parity that contain the angular momenta J = 0, 4, ... and J = 1, 3, ..., respectively.The explicit definitions of all operators are given in the supplemental material.We compute the symmetric 7 × 7 and 8 × 8 correlation matrices of these operators, using combinations of (Gaussian smeared) point-to-all and stochastic timeslice-to-all propagators [45].
Our calculations were performed for two different lattice sizes, 24 3 × 64 and 32 3 × 64, with a lattice spacing of approximately 0.12 fm and pion mass of approximately 220 MeV in both cases.The charm and bottom valence quarks were implemented using the Fermilab method [46] and lattice NRQCD [47], respectively, with approximately physical kinetic masses.Further details are provided in the supplemental material.We computed the bcud correlation matrices for approximately 1000 gauge configurations on each lattice, with 30 source locations per configuration for the elements computed using (Gaussian smeared) point-source propagators, and 3 random Z 2 × Z 2 sources on 4 timeslices per configuration for the elements computed with (Gaussian smeared) stochastic propagators; we also use color and spin dilution and the one-end trick [45].To extract the bcud finite-volume energy levels from these correlation matrices, we follow the well-established approach of solving the generalized eigenvalue problem (GEVP) [22,48] where we set t 0 /a = 3 and verified that the results do not significantly depend on this choice.We then perform single-exponential fits of the form λ n (t, t 0 ) = A n e −Ent to obtain the energy levels E n ; see the supplemental material for further details.
Our results for the lowest five energy levels of each bcud system are shown as a function of the spatial lattice size N s = L/a in Fig. 1.Also shown are the lowest four noninteracting B ( * ) -D energy levels, calculated as E = E B ( * ) (p 2 ) + E D (p 2 ) with momenta p satisfying the periodic boundary conditions [each component an integer multiple of 2π/L], and with the single-meson energies calculated on the lattice and described by the dispersion relations The values of E B ( * ) (0), m B ( * ) ,kin , E D (0), m D,kin , and m D,4 are provided in the supplemental material.In Fig. 1 we see that the actual bcud energy levels are shifted significantly relative to the noninteracting levels due to the meson-meson interactions in the finite volume, except for the third level in the case of J = 1 (we discuss the reason for this behavior farther below).
To rigorously investigate whether bound states or resonances exist, we map the observed finite-volume energy levels E n to infinite-volume S-wave B ( * ) -D scattering phase shifts δ 0 (k n ) using the Lüscher quantization condition where Z 00 is the generalized zeta function [23] and k n is the scattering momentum associated with energy level E n , calculated from ) with the dispersion relations (2).To ensure that the singlechannel, single-partial-wave approximation is applicable, we only extract the phase shifts for the energy levels below the B * -D * (J = 0) and B-D * (J = 1) thresholds.Furthermore, for J = 1, we observe that the third finitevolume energy level is consistent with the noninteracting |p| = 2π/L energy level that has multiplicity 2 once including both S-wave and D-wave structures, as we did in our operator basis.Because finite-volume interactions for higher partial waves are suppressed, we conclude that this energy level is dominantly D-wave, and we therefore exclude it from the Lüscher analysis.This is further corroborated by the eigenvectors from the GEVP, which show that this state has a non-negligible overlap only with the operator O that was subduced from a Dwave structure.Note that there are similar degeneracies for higher non-interacting levels above the energy region used in our analysis [44].
Our results for the scattering phase shifts, along with effective-range expansion (ERE) fits of the form are shown in Fig. 2 (Left).The numerical values of the fitted ERE parameters are given in the supplemental material.The scattering phase shift is related to the S-wave scattering amplitude and cross section by Poles of T 0 (k) at purely imaginary k correspond to genuine or virtual bound states for Im(k) > 0 or Im(k) < 0, respectively, while poles with Re(k) ̸ = 0 and Im(k) < 0 correspond to resonances.Using our ERE fits, we find genuine bound-state poles as well as resonance poles for both J = 0 and J = 1 at the values of given in Table I (s denotes the Mandelstam variable equal to the square of the center-of-momentum energy).We used our lattice results for the kinetic B ( * ) and D masses to evaluate this expression; to obtain predictions for absolute tetraquark bound state or resonance masses, one simply needs to add the experimental value of the threshold energy, m exp B ( * ) + m exp D .where k is the scattering momentum, δ0(k) is the scattering phase shift, and a = 0.11887(80) fm is the lattice spacing.The data points were obtained directly from the lattice energy levels, and the curves correspond to ERE fits through order k 4 .Also shown are the functions − −(ak) 2 (solid red parabolas) whose intersections with ak cot δ0(k) just below threshold correspond to the shallow bcud bound states we predict, and + −(ak) 2 (dashed red parabolas) whose intersections with ak cot δ0(k) would correspond to virtual bcud bound states.The vertical magenta lines show the positions of two-pion-exchange left-hand branch points.Right: Our results for the product of scattering momentum and B ( * ) -D scattering cross section, which is proportional to the scattering rate, as a function of center-of-momentum energy.
The resonances have masses of order 100 MeV above the B ( * ) -D thresholds and decay widths of order 200 MeV.We caution that the resonance poles lie outside the radius of convergence of the ERE, which is limited by the presence of a left-hand cut associated with twopion t-channel exchange in the scattering process (singlepion exchange would require a D * in the initial or final state, and is therefore not relevant here).The center-ofmomentum energy at which the left-hand cut starts is obtained from the kinematic relations for the Mandelstam variables by expressing s in terms of t and the scattering angle θ * , and then setting t = (2m π ) 2 and θ * = π [49]; this gives √ s cut − √ s th ≈ −18 MeV for both J = 0 and J = 1, corresponding to (ak) 2 cut ≈ −0.019, as indicated with the magenta lines in Fig. 2. While our ERE fit is seen to describe the data very well for real (ak) 2 in the full momentum range, the prediction of resonance poles away from the real axis may be less reliable.We note, however, that instead of expanding k cot δ 0 (k) around k 2 = 0 one could as well expand around the midpoint of the left hand cut and the second threshold.The convergence radius in the complex energy plane would then be significantly larger, while the resulting parametrizations of the lattice data points for k cot δ 0 (k) would be identical to those obtained from the ERE.The reason is that in both cases second-order polynomials in k 2 are fitted to the same data points.This suggests that our results can be trusted in a significantly larger region of the complex energy plane, namely disks of radius 94 MeV around E = 76 MeV for J = 0 and of radius 52 MeV around E = 34 MeV for J = 1.The predicted resonances are still located outside, but quite close to the boundaries of these convergence regions.To further test their stability, we also performed ERE fits through order k 6 .The coefficients of k 6 are found to be consistent with zero within the statistical uncertainties, and the other parameters remain consistent with those from the order-k 4 fit.For J = 0, the resonance pole obtained from the order-k 6 fit is at a similar location.For J = 1, where we have fewer data points, the uncertainties from the k 6 fit are too large to determine the pole locations.
The bound-state poles are extremely close to threshold and therefore well within the region of validity of the ERE.However, as can be seen from the ± −(ak) 2 parabolas in Fig. 2, statistical fluctuations could turn the genuine bound states into virtual bound states, which are not asymptotic states in QCD but would still strongly affect the B ( * ) -D scattering rates near threshold [50], or could, at first glance, lead to the complete disappearance of the poles.To quantify this uncertainty, we generated 10 000 multivariate Gaussian random samples for the ERE fit parameters a 0 , r 0 , and b 0 according to their mean values and covariance matrix.For each sample, we first checked whether it is still consistent with our finite-volume energy spectra.To this end, we calculate the two finite-volume energy spectra predicted by the sample using the Lüscher quantization condition.For some of the samples, the intersection of k cot δ 0 (k) and 2Z 00 (1; (k n L/2π) 2 )/(π 1/2 L) below threshold disappears for at least one of the volumes, leaving only the intersections that match the first and higher excited-state energy levels calculated on the lattice.These samples are therefore inconsistent with the observed spectra and we removed them.For each of the remaining samples, we then checked whether there is a genuine bound state (GBS), a virtual bound state (VBS), or no bound state (NBS), with the following outcome: • J = 0: 88.5% GBS, 11.5% VBS, 0.0% NBS, • J = 1: 97.7% GBS, 2.3% VBS, 0.0% NBS.
The lower and upper limits for ∆m GBS given in Table I correspond to the 16th and 84th percentiles of those random samples for which genuine bound states exist, while the central values correspond to the best-fit points.To further test our prediction of shallow bound states, we performed additional ERE fits of order k 0 and order k 2 using only the three data points closest to threshold, which are within the strict radius of convergence of the ERE.These fits, which are shown in the supplemental material, give consistent results.
Returning to the discussion of our main fits as shown in Fig. 2, we note that, in addition to the shallow-boundstate and the broad-resonance poles, there are poles with purely imaginary k below the left-hand branch point.These poles must be discarded because the direction of the crossing corresponds to a pole residue with an unphysical sign [51] and our parametrization of k cot δ 0 (k) is invalid in that region (bound states this far below threshold are also ruled out by the absence of corresponding finite-volume energy levels).
The scattering rate (probability per time) is equal to the product of flux and cross section, and hence proportional to k σ(k) for nonrelativistic k.These products are shown in Fig. 2 (Right) as a function of the center-ofmomentum energy.We emphasize that the scattering rates only depend on our fit functions for real-valued k 2 that interpolate our data very well, so these predictions are also expected to be very reliable.We observe sharp enhancements in the scattering rates close to the thresholds, related to the shallow bound states or virtual bound states.At higher energies, the scattering rates continue to be enhanced, likely by the broad resonances.The scattering rates are very close to the largest possible value allowed by unitarity, given by |T 0 | 2 = 1, up to several tens of MeV above threshold.
In summary, the substantial improvements made here in determining the bcud finite-volume energy levels allowed us to determine the detailed energy dependence of the B-D and B * -D S-wave scattering amplitudes for the first time using lattice QCD, revealing very interesting strong-interaction phenomena.We found poles for both J = 0 and J = 1 corresponding to shallow bound states, as well as hints for poles corresponding to broad resonances.While further lattice-QCD computations at additional lattice spacings and pion masses will be needed to pin down the exact location and nature of each pole at the physical point, we expect our prediction of shallow bound states, either genuine or virtual, to be quite robust.The possible resonances above threshold are very broad and are therefore presumably difficult to observe at the LHC and future experiments.On the other hand, if the J = 0 pole just below the B-D threshold is confirmed as a genuine bound state, this isoscalar, scalar bcud tetraquark will decay through the weak interaction only and could become the first tetraquark to be observed at the LHC with this feature.If the J = 1 pole just below the B * -D threshold is confirmed as a genuine bound state, it will decay electromagnetically into B Dγ (and also into the J = 0 tetraquark plus a photon, if that tetraquark is confirmed as a genuine bound state).
Acknowledgments: We thank Simone Bacchio, Luka Leskovec, Mikhail Mikhasenko, Akaki Rusetsky, Christopher Thomas, Bira van Kolck, and David Wilson for helpful discussions.We thank the MILC collaboration, and in particular Doug Toussaint, for providing the gauge-link ensembles.

I. Lattice actions and parameters
Our calculation uses the mixed-action setup that was tested and used successfully by the PNDME collaboration for nucleon-structure computations [55,56].This setup employs gauge configurations generated with 2+1+1 flavors of highly improved staggered (HISQ) sea quarks by the MILC collaboration [57], but uses the clover-improved Wilson action with HYP-smeared gauge links for the valence light quarks.Here we include two ensembles that differ only in the spatial lattice extent; their main properties are given in Table II.We set the bare valence light-quark mass to am l = −0.075and the clover coefficient to c SW = 1.05091 [55,56].We implement the valence charm quarks with the same form of clover-Wilson action and same value of c SW , but with mass parameter tuned according to the Fermilab method to eliminate the main heavy-quark discretization errors [46].That is, the bare mass is tuned such that the spin-averaged kinetic D-meson mass m spinav D,kin = (m D,kin + 3m D * ,kin )/4 matches its experimental value [58]; this condition is satisfied for our final choice am c = 0.6835775 at the 3% level.For the valence bottom quarks, we use order-v 4 lattice NRQCD [47] with tadpole improvement and order-α s corrections to the matching coefficients for the kinetic terms; all parameters are given in Ref. [59].The resulting spin-averaged kinetic B-meson mass is within 4% of the experimental value [58] (see Sec. III below for details).

II. Operators
For the bcud system with I(J P ) = 0(0 + ), we use the seven operators O where the repeated index j is summed over the spatial directions, the repeated color indices a and b are summed over the three colors, V S = L 3 is the spatial lattice volume, and The operators O For the scattering operators, the summations over the back-to-back momentum directions ensure that the operators transform in the A + 1 irrep of the octahedral group that contains J = 0.
For the bcud system with I(J P ) = 0(1 + ), we use the eight operators where D − (q) was defined in Eq. ( 14), is constructed as a colorsinglet contraction of two color-nonsinglet diquarks at the same spacetime point that is then jointly projected to zero momentum.Here, the two heavy quarks are combined to a flavor-symmetric spin-1 diquark and the two light ), and √ 2 • 2π/L.Two different operators are used for the case with one unit of back-to-back momentum to account for the mixing of S and D partial waves [44].Again, all quark fields are smeared, with the same parameters as used for J = 0.

III. D( * ) and B ( * ) dispersion relations
Our results for the D and D * meson energies as a function of spatial momentum squared are shown in Fig. 3.We performed fits to the combined data from the two ensembles using the three-parameter form to allow for different values of m D( * ) ,kin , and m D( * ) ,4 due to discretization errors.Higher powers of p are expected to be negligible for the momentum range we use.The fit results are given in Table III.The spin-averaged kinetic mass agrees with the experimental value [58] within 3%, confirming the successful tuning of the charm-quark mass according to the Fermilab method [46].We also find that the results for m D( * ) ,4 are actually consistent with m D( * ) ,kin within the statistical uncertainties.
For the B and B * mesons, we did not expect a significant difference between m B ( * ) ,kin , and m B ( * ) ,4 due to the high level of improvement of the lattice NRQCD action [59], and we therefore performed two-parameter fits of the form These fits included also a third ensemble of gauge configurations, a12m220L, with the same bare parameters and an even larger volume, and are discussed in more detail in Ref. [61].The data from all three ensembles are well-described jointly by Eq. ( 25) with the parameters given in Table IV.The spin-averaged kinetic B-meson mass is within 4% of the experimental value [58].
IV. bcud energies and B ( * ) -D scattering phase shifts Plots of the effective energies of the five lowest eigenvalues obtained from the GEVP for the bcud correlation matrices are shown in Fig. 4. We determine the energy levels E n by carrying out correlated least-χ 2 fits of the form λ n (t, t 0 ) = A n e −Ent .We perform fits for multiple different ranges t min ≤ t ≤ t max with sufficiently large t min to FIG. 3. Results for aE D (p 2 ) and aE D * (p 2 ) for 0 ≤ p 2 ≤ 4(2π/L) 2 from the two ensembles, along with fits of the form (24).  IV.B and B * meson dispersion-relation parameters in lattice units, obtained from combined fits to the data from the a12m220S, a12m220, and a12m220L ensembles [61].
ensure single-exponential behavior (we use all possible ranges with 8 ≤ t min /a ≤ 10, 15 ≤ t max /a ≤ 19 for J = 0, and 9 ≤ t min /a ≤ 11, 15 ≤ t max /a ≤ 19 for J = 1; these ranges were chosen to yield χ 2 /d.o.f. in the vicinity of 1).Then we obtain the final estimate for E n from a weighted average that takes into account correlations and lowers the weights of fits with χ 2 /d.o.f.> 1, following the FLAG averaging procedure [62] (see also our discussion in Appendix B of Ref. [13]).This weighting guarantees, in particular, that fits with too small t min , for which excited states might not be negligible, are suppressed in the final result.The statistical uncertainties are calculated and propagated to the further analysis using jackknife.
Our results for the bcud finite-volume energies, their values relative to the threshold, the corresponding scattering momenta squared, and the corresponding products of scattering momentum and cotangent of S-wave scattering phase )) of the five lowest eigenvalues λ0(t, t0), ..., λ4(t, t0) of the bcud system with I(J P ) = 0(0 + ) on the a12m220S lattice, obtained by solving the GEVP for the 7 × 7 correlation matrix.For sufficiently large t, these effective energies correspond to the five lowest finite-volume energy levels with the given quantum numbers.Also shown are the strong-decay-threshold energies computed on the same lattice.Note that all absolute energies contain an overall constant shift due to the use of NRQCD and of the Fermilab method.Top right: The corresponding plot for the bcud system with I(J P ) = 0(1 + ), using the 8 × 8 correlation matrix.Bottom left and right: Close-up view of aE eff,0 (t) and aE eff,1 (t) together with the corresponding energy levels obtained by following the FLAG averaging procedure.The bcud finite-volume energies in the T + 1 irrep, their values relative to the B * -D threshold, the corresponding B * -D scattering momenta squared, and the corresponding products of scattering momentum and cotangent of scattering phase shift (all in lattice units).Note that the absolute energies contain an overall constant offset due to the use of NRQCD and of the Fermilab method; this offset cancels in the differences to the threshold.shifts are listed in Tables V (for the A + 1 irrep relevant for J P = 0 + ) and VI (for the T + 1 irrep relevant for J P = 1 + ).As discussed in the main article, some high-lying energy levels are excluded from the phase-shift determination because they lie above inelastic thresholds, and the third energy level in the T + 1 irrep is excluded because it corresponds to a state dominated by a D wave; these cases are labeled "N/A" in the tables.
To demonstrate the importance of our operators and to verify the stability of our results, we also varied the set of operators entering the correlation matrix and the GEVP.For example, for ensemble a12m220 and J = 0, where we use the largest number of energy levels (the five lowest energy levels), we observed the following: • It is important to include local operators: when all three local operators (O 3 ) are omitted, the extracted second, third, and fourth energy levels are collectively shifted upward slightly, although the extracted ground-state energy remains essentially unchanged.

• When only the diquark-antidiquark operator O
A + 1 3 is omitted, the extracted energy levels remain essentially unchanged.
• When the scattering operator with the largest back-to-back momentum √ 3 • 2π/L is omitted, the extracted fourth energy level moves upward, to a position between the fourth and fifth energy levels obtained with the full operator set.This indicates that the number of scattering operators employed in this work is essential to study k cot δ 0 (k) up to the B * D * threshold, and raises some concern about the reliability of the extraction of the fifth energy level with the full operator basis.To test the stability of our determination of the ERE parameters (and hence of our predictions for the bound states and resonances), we performed an additional ERE fit for J = 0 in which we omitted the fifth energy level from ensemble a12m220.This fit gives a/a 0 = −0.022(19), r 0 /(2a) = 1.91(70), b 0 /a 3 = −18.3(5.2), which agrees very well with the full fit (see Table VII below).
V. ERE fit parameters using the full momentum range This section contains our results for the ERE fit parameters for the main fits that use all available scattering momenta and phase shifts from Tables V and VI.The fits were performed in lattice units, where a is the lattice spacing.We used the coefficients of (ak) 0 , (ak) 2 , and (ak) 4 as the fit parameters.Our results for these parameters, along with their correlation matrices, are given in Tables VII and VIII.In addition, we provide the values of 1/a 0 , r 0 , and b 0 in physical units in Table IX.
In our ERE fits, we took into account the uncertainties of the lattice results for both ak cot δ 0 and (ak) 2 (the latter was done by promoting the (ak) 2 values to additional fit parameters).We also took into account the correlations between the different ak cot δ 0 values and the correlations between the different (ak) 2 values.However, we did not include the cross-correlations between (ak) 2 and ak cot δ 0 , which would lead to near-zero eigenvalues of the data covariance matrix whose inverse appears in the definition of χ 2 .To investigate how much this choice affects our results, we performed additional fits in which we also included those cross-correlations, but multiplied all off-diagonal elements of the data covariance matrix by a factor of 0.95 to avoid the issue with near-zero eigenvalues.These fits yield very similar results for both the central values and the uncertainties of the ERE parameters.

VI. ERE fits and bound-state poles using only the low-momentum region
Our ERE fits to only the three data points closest to the thresholds are shown in Fig. 5.In this momentum region, 0th-order fits already describe the data well; for comparison, we also show fits through order k 2 .The resulting ERE parameters, genuine-bound-state masses, and probabilities for the presence of a genuine bound state, a virtual bound state, or no bound state are given in Tables X and XI.

2 • 3 • 8 are
D scattering operators with zero total spatial momentum in which the B and D operators have back-to-back momenta of magnitudes 0, 2π/L, √ 2π/L, and √ 2π/L (L is the spatial lattice size).Similarly, for the bcud system with I(J P ) = 0(1 + ), we use eight operators O B * -D scattering operators in which the B * and D have back-to-back momenta of magnitudes 0, 2π/L (for both O

FIG. 1 .
FIG. 1. Left: The finite-volume energies of the bcud system with I(J P ) = 0(0 + ) as a function of the spatial lattice size.The data points with error bars show the actual finite-volume energy levels; points plotted with a lighter shade of gray are excluded from the Lüscher analysis.The solid blue curves correspond to what would be the noninteracting B-D energy levels, the lowest of which coincides with the strong-decay threshold.The dashed green line shows the B * -D * threshold.Right: The corresponding plot for I(J P ) = 0(1 + ).Here, the solid blue curves correspond to what would be the noninteracting B * -D energy levels, and the dashed green line corresponds to the B-D * threshold.

√√− 1 ] J = 1 FIG. 2 .
FIG. 2. Left: Our results for the functions ak cot δ0(k) for S-wave B-D scattering (top) and S-wave B * -D scattering (bottom),where k is the scattering momentum, δ0(k) is the scattering phase shift, and a = 0.11887(80) fm is the lattice spacing.The data points were obtained directly from the lattice energy levels, and the curves correspond to ERE fits through order k 4 .Also shown are the functions − −(ak) 2 (solid red parabolas) whose intersections with ak cot δ0(k) just below threshold correspond to the shallow bcud bound states we predict, and + −(ak) 2 (dashed red parabolas) whose intersections with ak cot δ0(k) would correspond to virtual bcud bound states.The vertical magenta lines show the positions of two-pion-exchange left-hand branch points.Right: Our results for the product of scattering momentum and B ( * ) -D scattering cross section, which is proportional to the scattering rate, as a function of center-of-momentum energy.

1 2 1 3 2 •
are constructed as products of color-singlet B, D and B * , D * operators at the same spacetime point that are then jointly projected to zero momentum by summing over the spatial coordinates.The operator O A + is constructed as a color-singlet contraction of two color-nonsinglet diquarks at the same spacetime point that is then jointly projected to zero momentum.The operators O D "scattering" operators in which the B and D operators are individually momentum-projected and have back-to-back momenta of magnitudes 0, 2π/L, √ 2π/L, and √ 3 • 2π/L (L is the spatial lattice size).

1 7
and k = x, y, z denotes the spatial polarization direction (the operator O T + is shown for k = z only).These operators transform in the T + 1 irrep of the octahedral group that contains J = 1.The operators O products of color-singlet B ( * ) and D( * ) operators at the same spacetime point that are then jointly projected to zero momentum by summing over the spatial coordinates.The operator O T + 1 4

FIG. 5 .
FIG.5.ERE fits using only the three data points closest to the thresholds for S-wave B-D scattering (left) and S-wave B * -D scattering (right); here, a is the lattice spacing.Also shown are the functions − −(ak) 2 (solid red parabolas) whose intersections with ak cot δ0(k) just below threshold correspond to the shallow bcud bound states we predict, and + −(ak) 2 (dashed red parabolas) whose intersections with ak cot δ0(k) would correspond to virtual bcud bound states.The vertical magenta lines show the positions of two-pion-exchange left-hand branch points.

TABLE II
[57]e main properties of the two gauge-link ensembles used in this work.Here, Ns and Nt are the numbers of lattice sites in spatial and temporal directions, a is the lattice spacing from the r1 scale[57], m (sea) π

TABLE V .
The bcud finite-volume energies in the A + 1 irrep, their values relative to the B-D threshold, the corresponding B-D scattering momenta squared, and the corresponding products of scattering momentum and cotangent of scattering phase shift (all in lattice units).Note that the absolute energies contain an overall constant offset due to the use of NRQCD and of the Fermilab method; this offset cancels in the differences to the threshold.

TABLE VII .
Fit results for the J = 0 ERE parameters and their correlation matrix.

TABLE IX .
The ERE parameters in physical units.

TABLE X .
Results for the inverse scattering lengths and genuine-bound-state masses (relative to the B ( * ) -D thresholds) from the order-k 0 ERE fits to the three data points closest to the threshold.Also shown are the fractions of the random EREparameter samples that yield a genuine bound state (GBS), a virtual bound state (VBS), or no bound state (NBS).

TABLE XI .
Results for the ERE parameters and genuine-bound-state masses (relative to the B ( * ) -D thresholds) from the order-k 2 ERE fits to the three data points closest to the threshold.Also shown are the fractions of the random ERE-parameter samples that yield a genuine bound state (GBS), a virtual bound state (VBS), or no bound state (NBS).