The T_bc tetraquarks near the BD threshold

The past two decades have witnessed a dramatic expansion of the known hadron spectrum, including the discovery of numerous exotic states that cannot be accommodated within the conventional $ q\bar{q} $ and $ qqq $ quark model. Among these, manifestly flavor-exotic multiquark states occupy a particularly important position because they cannot mix with ordinary mesons or baryons and therefore provide an especially clean probe of nonperturbative QCD dynamics. The recent observation by LHCb of the doubly charmed tetraquark $ T_{cc}^+ $ [1, 2], with a mass only about 360 keV below the $ D^0 D^{*+} $ threshold, has definitively established that anti-heavy–anti-heavy–light–light systems $ \bar{Q}\bar{Q}'qq' $ can form near-threshold bound structures in nature and has intensified interest in their bottom-sector analogues.

From a theoretical standpoint, $ \bar{Q}\bar{Q}'qq' $ tetraquarks are among the most promising candidates for strong-interaction-stable exotic hadrons. The possibility that four-quark configurations containing heavy quarks may be stable against strong decay was first raised in the pioneering study of Ref. [3], which preceded subsequent model analyses of doubly heavy systems. In the heavy-quark limit $ m_Q \to \infty $, the attractive color-Coulomb interaction between two heavy antiquarks in a color-$ \mathbf{\bar{3}}_c $ configuration guarantees the existence of bound states [3–6], in close analogy with the role of a heavy quark in singly heavy baryons. The binding is reinforced by the attraction between the two light quarks in Jaffe’s “good diquark” configuration ($ \mathbf{\bar{3}}_c $, spin 0) [7]. Both mechanisms are absent, or strongly suppressed, in a configuration of two well-separated heavy-light mesons. These considerations single out the $ I(J^P)=0(1^+) $ $ ud\bar{b}\bar{b} $ and $ I(J^P)=\dfrac{1}{2}(1^+) $ $ \ell s\bar{b}\bar{b} $ channels as the most favorable candidates for binding. Indeed, multiple independent lattice QCD calculations have now firmly established the existence of deeply bound, strong-interaction-stable tetraquarks in both channels [8–12].

Doubly bottom tetraquarks, however, are experimentally challenging because their production requires two $ b\bar{b} $ pairs. The mixed bottom-charm system $ T_{bc} $ is considerably more favorable in this respect: simultaneous $ b\bar{b} $ and $ c\bar{c} $ production is already well established through the copious observation of $ B_c $ mesons, and the estimated production cross section of a $ T_{bc} $ tetraquark at the LHC is approximately forty times larger than that of its doubly bottom counterpart [13]. Whether a bound state actually exists in this channel is, however, a much more delicate question than in the doubly bottom case. Replacing one $ \bar{b} $ by $ \bar{c} $ reduces the reduced mass of the heavy antidiquark and hence the strength of the short-distance color-Coulomb attraction. At the same time, it enhances the residual spin-dependent heavy-light interactions that are known from the heavy-baryon spectrum to oppose binding. The $ T_{bc} $ system thus lies at the boundary between bound and unbound, and its study provides a stringent test of our understanding of the dynamics responsible for exotic binding.

Until very recently, the theoretical situation for $ T_{bc} $ was strikingly inconclusive. Heavy-quark symmetry arguments [6] predict no bound state in either the $ I(J^P)=0(0^+) $ or $ 0(1^+) $ channel. Most nonrelativistic and relativized quark models [14–20] find ground states that are at most very weakly bound or unbound, whereas many chiral quark model analyses [21–23] predict substantial binding of order 130–220 MeV. QCD sum-rule studies are similarly dispersed, ranging from unbound or modestly bound states to very deeply bound ones [24, 25]. A coherent picture has therefore had to come from first-principles lattice QCD.

Direct lattice QCD studies of the $ T_{bc} $ system have themselves evolved considerably and initially produced apparently contradictory results. The first such calculation employed nonrelativistic QCD (NRQCD) for the bottom quark and a relativistic heavy-quark action for the charm quark on three PACS-CS ensembles with pion masses reaching $ m_\pi \simeq 164 $ MeV [26]. That work reported evidence for a strong-interaction-stable $ I(J^P)=0(1^+) $ ground state with a binding energy in the range $ -61\; {\rm{MeV}} \lt \Delta E \lt -15\; {\rm{MeV}} $ relative to the $ D\bar{B}^* $ threshold, placing the predicted mass near the electromagnetic $ D\bar{B}\gamma $ decay threshold. The lightest ensemble, however, had $ m_\pi L = 2.4 $, raising concerns about finite-volume effects, and the plateaus of the wall-local correlators were short and occurred only at late times. A follow-up study [11], performed on a larger-volume ensemble with $ m_\pi \simeq 192 $ MeV and $ m_\pi L = 4.2 $ and employing an improved “box-sink” construction together with an expanded operator basis, did not confirm the earlier signal. The $ I(J^P)=0(1^+) $ $ T_{bc} $ ground state was found to lie marginally above the $ D\bar{B}^* $ threshold, with no evidence for binding in either the $ 0^+ $ or $ 1^+ $ channel. The earlier indication of binding was attributed to late-onset wall-local plateaus and residual finite-volume effects. Subsequently, Ref. [12] studied both the $ 0(0^+) $ and $ 0(1^+) $ channels using domain-wall light quarks, NRQCD bottom quarks, and an anisotropic clover action for the charm. Crucially, that work included nonlocal meson-meson scattering interpolating operators at the sink on five ensembles, including one at the physical pion mass. The authors found finite-volume ground-state energies consistent with the respective lowest two-meson thresholds $ BD $ and $ BD^* $ and concluded that there was no clear evidence of binding in either channel, although a shallow bound state close to threshold could not be ruled out.

Two very recent lattice QCD studies have substantially clarified this picture and, for the first time, provided consistent evidence for binding in the $ T_{bc} $ system. Ref. [27] performed a calculation on four gauge-field ensembles generated by the MILC Collaboration with $ N_f = 2+1+1 $ dynamical quark flavors, employing the Highly Improved Staggered Quark (HISQ) action. The ensembles span lattice spacings from $ \sim 0.058 $ to $ \sim 0.12 $ fm and include five valence light-quark masses, with a correlation-matrix basis constructed from both local diquark–antidiquark operators and local two-meson $ D\bar{B}^* $ and $ B\bar{D}^* $ interpolators. Extracting the $ D\bar{B}^* $ S-wave scattering length via Lüscher’s method and performing continuum and chiral extrapolations, they found a positive scattering length $ a_0^{{\rm{phys}}} = +0.57(^{+4}_{-5})(17) $ fm at the physical pion mass, corresponding to a genuine $ I(J^P)=0(1^+) $ $ T_{bc} $ bound state with a binding energy of $ -43(^{+6}_{-7})(^{+14}_{-24}) $ MeV relative to the $ D\bar{B}^* $ threshold. The authors further observed that the binding strength decreases with increasing $ m_{u/d} $ and that the system becomes unbound above a critical pseudoscalar mass $ M_{{\rm{ps}}}^* = 2.73(21)(19) $ GeV. Complementing this work, Ref. [28] carried out the first lattice-QCD determination of the full energy dependence of the $ B\text{-}\bar{D} $ and $ B^*\text{-}\bar{D} $ isospin-0, S-wave scattering amplitudes both below and above threshold. This study used large $ 7\times 7 $ and $ 8\times 8 $ symmetric correlation matrices with $ B^{(*)}\text{-}\bar{D} $ scattering operators at both source and sink and with nonzero back-to-back momenta. Working at a single lattice spacing and pion mass $ m_\pi \simeq 220 $ MeV on two volumes, and fitting the extracted phase shifts with effective-range expansions, they found near-threshold poles in both the $ J=0 $ and $ J=1 $ channels. These poles correspond to shallow states, either genuine bound states or virtual states, only a few MeV below the respective $ B^{(*)}\text{-}\bar{D} $ thresholds, together with indications of broad resonances of width $ \sim 200 $ MeV located about 100 MeV above threshold. Taken together, these two studies, while differing in the depth of the predicted binding and in the resolution of the genuine-versus-virtual bound-state question, both point toward attractive $ B^{(*)}\bar{D} $ interactions strong enough to produce near-threshold poles and mark a qualitative shift in the theoretical status of the $ T_{bc} $ system.

Despite this recent progress, several important questions remain open. The binding energy of the $ I(J^P)=0(1^+) $ $ T_{bc} $ is quoted as $ -43(^{+6}_{-7})(^{+14}_{-24}) $ MeV in Ref. [27] but as at most a few MeV in Ref. [28]. The genuine versus virtual character of the near-threshold poles, the existence of the $ I(J^P)=0(0^+) $ partner, and the nature of the possible broad resonances above threshold all require further scrutiny. Left-hand-cut effects from pion exchange, coupled-channel dynamics involving $ BD $, $ B^*D $, $ BD^* $, and $ B^*D^* $, discretization effects, and the chiral extrapolation all contribute to the current systematic uncertainties and motivate additional independent calculations with different lattice actions, ensembles, and analysis strategies.

In the present work, we compute the masses, hyperfine splittings, and mean inter-cluster radii of the doubly heavy open-flavor tetraquarks $ T_{bc}^{(0)} $ ($ J^{P}=0^{+} $) and $ T_{bc}^{(1)} $ ($ J^{P}=1^{+} $) within the Born–Oppenheimer dynamical diquark model, using the $ \Sigma_g^+(1S) $ potential extracted from lattice QCD. Comparing the resulting spectrum with the lowest $ B^{(*)}\bar{D}^{(*)} $ thresholds, we find that the central value of the scalar mass straddles the $ B\bar{D} $ threshold to within a few MeV, thereby placing it on the boundary between a genuinely bound tetraquark and a narrow near-threshold tetraquark. The axial-vector eigenvalue consistently lies $ \sim 25 $ MeV above $ B^{*}\bar{D} $, making this state a natural candidate for an S-wave resonance with discovery prospects at LHCb and Belle II. Because the states of interest lie in the immediate vicinity of open two-meson thresholds, we devote particular attention to the systematic limitations of the single-channel constituent description employed here: the restriction to the $ \bar{{\bf{3}}}_c\otimes{\bf{3}}_c $ color configuration, the size of nonadiabatic corrections to the BO approximation, the expected impact of the neglected coupled-channel $ B^{(*)}\bar{D}^{(*)} $ dynamics, and the parametric uncertainties of the model inputs. We formulate our conclusions accordingly.

The remainder of this paper is organized as follows. In Sec. II we summarize the dynamical diquark model in its Born–Oppenheimer formulation, recalling the role of the $ \Sigma_g^+(1S) $ potential and the construction of the heavy-antidiquark–light-diquark Hamiltonian. We also describe the numerical solution of the radial Schrödinger equation and specify the input parameters. Section III presents our results for the masses, hyperfine splittings, and spatial observables of the $ T_{bc}^{(0,1)} $ states, together with a graphical comparison to the relevant $ B^{(*)}\bar{D}^{(*)} $ thresholds and a discussion of the resulting phenomenology, including the prospects for experimental discovery. We summarize our findings and outline directions for future work in Sec. IV.

The Born-Oppenheimer (BO) approximation was developed in the early days of quantum mechanics [29] as a tool for studying molecular systems. It exploits the large separation of time scales between the slow motion of heavy atomic nuclei and the rapid response of the surrounding electrons. Because the nuclei are much heavier than the electrons, they can be treated as approximately static sources of the electromagnetic field at each instant. The electrons adjust almost instantaneously to the nuclear configuration, and their energy—together with the mutual Coulomb repulsion of the nuclei—defines an effective potential whose eigenvalues determine the vibrational and rotational spectrum of the molecule [30].

An analogous separation of scales arises in QCD whenever a hadronic system contains heavy color sources. The BO approximation for heavy quarkonium $ (Q\bar{Q}) $ was first formulated in Ref. [31], where the role of the electrons is played by the gluon field and any light quarks present in the system, while the role of the nuclei is played by the heavy quark and antiquark. The large ratio $ m_Q/\Lambda_{{\rm{QCD}}} $ ensures that the gluonic and light-quark degrees of freedom adjust adiabatically to the slowly varying heavy-quark configuration, generating a set of static BO potentials $ V_\Gamma(r) $ labeled by the quantum numbers Γ of the gluon field. The heavy-quark dynamics is then obtained by solving the Schrödinger equation in these potentials, in direct analogy with the molecular case. This picture has since been formulated as a systematic effective field theory (EFT), the Born-Oppenheimer EFT [32, 33], and applied to a wide range of exotic hadrons, including the $ \chi_{c1}(3872) $ and $ T_{cc}^+(3875) $ [34], as well as hidden-charm pentaquarks [35]. Systematic studies of tetraquarks in the BO framework and in the large-N expansion can be found in Ref. [36].

When the BO approximation is applied to four-quark systems, at least two sufficiently heavy color sources are required, together with degrees of freedom corresponding to light quarks and gluons. The slow sources—here the heavy antiquarks $ \bar{b} $ and $ \bar{c} $, or, equivalently, the composite clusters built from them—move adiabatically in a confining potential generated by the fast gluonic and light-quark fields [37]. Because the typical excitation energy of the gluon field and the light-quark cloud is of order $ \Lambda_{{\rm{QCD}}}\simeq 200 $–300 MeV, as confirmed by lattice-QCD studies of the BO potentials [38], while the relative motion of the heavy clusters occurs on a parametrically slower time scale, the BO expansion provides a controlled framework for computing the tetraquark spectrum. The dominant systematic uncertainty of the approach is generically of order $ {\cal{O}}(\Lambda_{{\rm{QCD}}}) $ and primarily affects the fine-structure splittings, while the gross features of the spectrum and the spatial properties of the states are expected to be robust.

For the $ T_{bc} $ system considered in this work, both heavy antiquarks are unambiguously “slow” in the BO sense, with $ m_c,m_b\gg\Lambda_{{\rm{QCD}}} $. Thus, the adiabatic approximation is on an even firmer footing than in the charm-strange case studied recently in Ref. [39]. The light-quark pair $ (u,d) $ and the gluon field are treated as the fast degrees of freedom, whose response to the slowly varying positions of the heavy antiquarks is encoded in an effective BO potential between the color clusters that contain those heavy antiquarks.

The size of the nonadiabatic corrections neglected in this treatment can be estimated as follows. The leading corrections arise from mixing of the $ \Sigma_g^+ $ potential with excited BO potentials, induced by the kinetic operator of the heavy clusters. They are parametrically suppressed by the ratio of the characteristic energy of the slow degrees of freedom to the gluonic excitation gap. The former is the level spacing of the radial spectrum in the $ \Sigma_g^+(1S) $ potential, of order a few hundred MeV for the reduced masses relevant here. The latter separates $ \Sigma_g^+ $ from the first excited potentials ($ \Pi_u $, $ \Sigma_u^- $), which lattice QCD places at roughly 1 GeV at the relevant inter-cluster separations [31, 38]. Since the mixing enters the energies at second order, the resulting shifts are expected to be at the level of a few percent of the radial level spacing, i.e. of order 10 MeV. Corrections of this size affect fine-structure quantities but not the gross features of the spectrum. A physically distinct—and, in the present context, more important—limitation is that the adiabatic $ \Sigma_g^+ $ potential does not describe the rearrangement of the light quarks into color-singlet heavy-light mesons, i.e. the coupling of the δ–$ \bar\delta' $ configuration to the open $ B^{(*)}\bar{D}^{(*)} $ continua. For eigenvalues located within a few tens of MeV of these thresholds, as found below, this coupling can shift the energies by an amount comparable to their displacement from threshold; we return to this point quantitatively in Sec. III Statements about the near-threshold character of individual states must therefore be understood within these combined $ {\cal{O}}(10\; \text{MeV}) $ systematic limitations of the single-channel adiabatic description, which are addressed explicitly below.

This picture naturally matches onto the dynamical diquark model, which we now describe.

A compact tetraquark can be described as a bound state of a colored diquark and a colored antidiquark, with the diquark (antidiquark) treated as an effective degree of freedom. A diquark is a colored bound quark–quark pair $ (qq) $, and an antidiquark is a colored bound antiquark–antiquark pair $ (\bar{q}\bar{q}) $; these can be combined into a color singlet. More detailed discussions of diquark degrees of freedom can be found in the reviews of Refs. [40, 41]. Diquark–antidiquark potential models of this type have been applied broadly across the tetraquark landscape, from all-light [42] and singly charmed [43] systems to doubly heavy configurations, and unified quark–diquark effective-mass formalisms calibrated on the baryon spectrum have recently been developed for heavy multiquarks [44]. Since a quark transforms in the fundamental representation 3 and an antiquark in $ \bar{{\bf{3}}} $, combining two quarks gives $ {\bf{3}}\otimes{\bf{3}}=\bar{{\bf{3}}}\oplus{\bf{6}} $, while combining two antiquarks gives $ \bar{{\bf{3}}}\otimes\bar{{\bf{3}}}={\bf{3}}\oplus\bar{{\bf{6}}} $. A color-singlet tetraquark $ |qq\bar{q}\bar{q}\rangle $ can therefore be formed either from the $ \bar{{\bf{3}}}\otimes{\bf{3}} $ channel or from the $ \bar{{\bf{6}}}\otimes{\bf{6}} $ channel.

In the present work, we restrict ourselves to the $ \bar{{\bf{3}}}\otimes{\bf{3}} $ configuration, and it is important to state explicitly the motivation for, and the limitations of, this choice. At short distances, one-gluon exchange between two quarks is attractive in the $ \bar{{\bf{3}}}_c $ channel and repulsive in the $ {\bf{6}}_c $ channel; thus, only the antitriplet supports compact, quasi-bound diquark clusters. It must be stressed, however, that this argument concerns the interaction inside each cluster and does not by itself disfavor the sextet configuration. A $ {\bf{6}}_c $ diquark and a $ \bar{{\bf{6}}}_c $ antidiquark can be coupled to a color singlet, and the static potential between them is attractive. Indeed, Casimir scaling—verified on the lattice at the few-percent level for source separations up to 1 fm [45]—makes it stronger than the $ \bar{{\bf{3}}}_c $–$ {\bf{3}}_c $ potential by a factor of $ C_6/C_3 = 5/2 $, making both the color-Coulomb well deeper and the confining rise steeper by the same factor. The relative energies of the two color channels are therefore determined by a genuine competition: the repulsive intra-cluster interaction, which prevents the formation of a compact $ {\bf{6}}_c $ cluster and thereby raises the configuration in energy, competes with the enhanced inter-cluster attraction, which lowers it at short distances but penalizes it at the separations $ \langle r\rangle\simeq 0.45 $ fm relevant here. We examine this competition quantitatively in Sec. III D and find that, within the present framework, the net effect places the sextet configuration several hundred MeV above the antitriplet configuration, so that its admixture acts as a second-order correction rather than as a competing leading configuration. In addition, the dynamical diquark model is based on a specific production mechanism, in which a compact $ \bar{{\bf{3}}}_c $ diquark and $ {\bf{3}}_c $ antidiquark are formed essentially at a point and separate rapidly before hadronization [46]; the lattice $ \Sigma_g^+ $ potential employed below is precisely the static energy of such a color source–antisource pair. The $ \bar{{\bf{3}}}\otimes{\bf{3}} $ restriction is therefore internally consistent with the dynamical assumptions of the model, rather than an independent approximation. Nevertheless, it must be acknowledged that, in a general four-quark Hamiltonian, the $ \bar{{\bf{3}}}\otimes{\bf{3}} $ and $ \bar{{\bf{6}}}\otimes{\bf{6}} $ channels mix, and static considerations alone do not determine the size of the induced sextet admixture. Coupled-color-channel quark-model calculations indicate that this admixture can be non-negligible, particularly in doubly heavy systems in which the heavy pair is not constrained by the Pauli principle. In particular, the $ bc\bar{q}\bar{q} $ analysis of Ref. [18], employing a Bhaduri-type potential and retaining the full color space, finds that the physical eigenstates contain sizable components beyond a single $ \bar{{\bf{3}}}\otimes{\bf{3}} $ configuration and that binding conclusions are sensitive to the retained color basis; similar observations have been made in chiral quark-model studies [21–23]. Since enlarging the variational space can only lower the ground-state energy, the eigenvalues obtained below within the pure $ \bar{{\bf{3}}}\otimes{\bf{3}} $ sector should be regarded as upper limits within the model, with the neglected color mixing acting in the direction of additional binding. A quantitative assessment of this effect requires a genuine four-body calculation and lies outside the Born–Oppenheimer diquark framework adopted here. Nevertheless, a semiquantitative estimate of its size within the BO framework is presented in Sec. III D, and we take this restriction into account when interpreting the near-threshold results of Sec. III C.

The dynamical diquark model, originally proposed in Ref. [46], postulates that the confinement mechanism itself is responsible for the binding of exotic states. In a hard production process, a diquark–antidiquark pair is formed essentially instantaneously and, as a consequence of the production kinematics, rapidly separates. As colored objects, the diquark δ and antidiquark $ \bar{\delta} $ cannot move arbitrarily far apart; instead, a color flux tube (or string) develops between them, and the system evolves in a manner closely analogous to heavy quarkonium. The model was subsequently refined and extended in Refs. [47–51], where δ and $ \bar{\delta} $ are treated as quasi-bound hadronic subcomponents whose relative motion is quantized in the BO potentials introduced above. The distinct quantum numbers of the gluonic field give rise to towers of states labelled by the BO quantum numbers $ \Lambda_\eta^\epsilon $ familiar from diatomic molecular physics, with the lowest multiplet supported by the ground-state potential $ \Sigma_g^+ $.

The dynamical diquark model assumes that the diquarks are pointlike and carry no internal orbital excitation, so that the entire quantum-number content of a given multiplet is encoded in the diquark and antidiquark spins, together with the relative orbital angular momentum L and gluonic excitation. For the lowest S-wave multiplet $ \Sigma_g^+(1S) $, in the basis of good diquark-spin eigenvalues, the states of the $ \delta\text{-}\bar{\delta} $ system are [47, 49]

where the inner subscripts denote the spins of δ and $ \bar{\delta} $, and the outer subscript denotes the total spin J, which, for S-wave states, coincides with the total quark spin S. These six basis states span the full $ \Sigma_g^+(1S) $ multiplet. When the two quarks inside δ (or the two antiquarks inside $ \bar{\delta} $) are identical, Fermi statistics, together with the antisymmetry of the color-$ \bar{{\bf{3}}}_c $ configuration and the symmetry of the S-wave spatial wave function, require the spin wave function to be symmetric. Thus, only the $ s_\delta=1 $ (or $ s_{\bar{\delta}}=1 $) configuration remains [51]. This restriction applies, for example, to $ cc $ and $ \bar{b}\bar{b} $ clusters, but not to nonidentical-flavor clusters such as $ [ud] $ or $ [\bar{b}\bar{c}] $, for which both spin-0 and spin-1 configurations are allowed.

The $ T_{bc} $ system of interest is an open-flavor tetraquark, built from a light diquark $ \delta=[ud] $ and an antidiquark $ \bar{\delta}'=[\bar{b}\bar{c}] $ that are not charge-conjugate partners. Several consequences follow immediately. First, because the two heavy antiquarks carry different flavors, the Pauli principle does not constrain the $ \bar{b}\bar{c} $ cluster, and both the scalar ($ s_{\bar{\delta}'}=0 $) and axial-vector ($ s_{\bar{\delta}'}=1 $) configurations are allowed. This situation contrasts sharply with the doubly bottom $ ud\bar{b}\bar{b} $ system, in which only the $ s_{\bar{\delta}}=1 $ configuration survives in the S-wave. Second, because δ and $ \bar{\delta}' $ are no longer C-conjugate partners, C parity is not a good quantum number for individual $ T_{bc} $ states, and the symmetric and antisymmetric combinations appearing in the $ X_1 $ and Z states of Eq. (1) do not correspond to distinct physical eigenstates [39]. The labels $ J^{PC} $ in Eq. (1) should therefore be understood only as a convenient shorthand for the underlying spin couplings.

The appropriate basis for the open-flavor system is the direct-product basis $ |s_\delta,s_{\bar{\delta}'}\rangle_J $, with $ s_\delta,s_{\bar{\delta}'}\in\{0,1\} $ and J denoting the total spin of the tetraquark. The full $ \Sigma_g^+(1S) $ multiplet is then spanned by the six states

yielding two $ J^P=0^+ $, three $ J^P=1^+ $, and one $ J^P=2^+ $ states, all with $ P=(-1)^L=+1 $.

A crucial structural feature of the $ T_{bc} $ system is that Fermi statistics for the S-wave, color-$ \bar{{\bf{3}}}_c $ light diquark constrains the isospin of the state through the spin of the light diquark. Since the combined color $ \otimes $ space $ \otimes $ flavor $ \otimes $ spin wave function must be totally antisymmetric, with the color component antisymmetric and the $ L=0 $ spatial component symmetric, the product of the flavor and spin components must be antisymmetric. This leads to the fixed correlation

with the combinations $ (I,s_\delta)=(0,1) $ and $ (1,0) $ forbidden. Consequently, the physical Hilbert space of the light diquark is two-dimensional rather than four-dimensional: I and $ s_\delta $ are not independent quantum numbers but are strictly correlated.

Imposing Eq. (3) on the six states of Eq. (2), the $ \Sigma_g^+(1S) $ multiplet splits into two physically distinct sectors. The isoscalar sector is built on the good light diquark ($ s_\delta=0 $), which can combine with either a scalar or an axial-vector heavy antidiquark to give

These are the only two $ I=0 $ members of the multiplet. In particular, the state $ |0_\delta,1_{\bar{\delta}'}\rangle_{J=1} $ appearing in Eq. (2) is identical to $ T_{bc}^{(1)} $, whereas $ |1_\delta,0_{\bar{\delta}'}\rangle_{J=1} $ has $ s_\delta=1 $ and therefore belongs to the isovector sector rather than to the $ I=0 $ sector. Consequently, the isovector sector, built on the bad light diquark ($ s_\delta=1 $), contains only configurations with $ s_{\bar{\delta}'}=1 $, namely

together with a fourth isovector state,

which is constructed from the scalar heavy antidiquark. The two isovector $ 1^{+} $ states are therefore not degenerate: as shown in Sec. II C, they are separated by exactly the same combination of the heavy-antidiquark mass gap and chromomagnetic coupling that splits the two isoscalar states. Thus, the $ I=0 $ sector contains two states, and the $ I=1 $ sector contains four, for a total of six, consistent with the dimensionality of the $ \Sigma_g^+(1S) $ multiplet.

The scalar state $ T_{bc}^{(0)} $ is the natural dynamical-diquark candidate for the near-threshold $ I(J^P)=0(0^+) $ pole found in the B-$ \bar{D} $ scattering analysis of Ref. [28], whereas the axial-vector state $ T_{bc}^{(1)} $ is the counterpart of the $ T_{bc} $ bound state reported in lattice studies in Refs. [26–28]. The isoscalar states are the primary focus of the present work, because lattice-QCD evidence exists for them and they lie in the immediate vicinity of the lowest open-flavor thresholds. The isovector members of the multiplet are nevertheless computed and discussed in Sec. III D, where they are found to lie substantially higher in mass—above every $ B^{(*)}\bar{D}^{(*)} $ threshold—and are therefore expected to be broad, consistent with the absence of any lattice-QCD signal for them to date.

Following the dynamical-diquark analyses of open-flavor systems in Refs. [39, 49, 50], and specializing to the S-wave $ \Sigma_g^+(1S) $ multiplet relevant to the $ T_{bc} $ ground states, the effective Hamiltonian for the δ–$ \bar{\delta}' $ system is written as

with

Here $ V_{\Sigma_g^+}(r) $ is the BO potential between the color-$ \bar{{\bf{3}}}_c $ diquark $ \delta=[ud] $ and the color-$ {\bf{3}}_c $ antidiquark $ \bar{\delta}'=[\bar{b}\bar{c}] $, obtained from lattice-QCD determinations of the static potential in the fundamental representation [31, 38]. Because we restrict the analysis to $ L=0 $, all spin-orbit and tensor contributions vanish identically, and the remaining spin-dependent dynamics is contained in $ H_\kappa $.

Because $ T_{bc} $ is an open-flavor system with different quark content in the two clusters, the spin-spin Hamiltonian contains two independent short-range color-magnetic couplings:

where $ \kappa_{ud} $ and $ \kappa_{\bar{b}\bar{c}} $ describe the color-magnetic interactions within the light and heavy clusters, respectively. Their expectation values in the diquark-spin basis $ |s_\delta,s_{\bar{\delta}'}\rangle_J $ are obtained from

The matrix elements are applied separately to each cluster. Together with the isospin–spin lock of Eq. (3), they automatically generate the full isoscalar–isovector splitting of the $ \Sigma_g^+(1S) $ multiplet without requiring any additional isospin-dependent operator.

The physical mass of a $ T_{bc} $ tetraquark state $ |n;s_\delta,s_{\bar{\delta}'};J\rangle $ with radial quantum number n, diquark spins $ (s_\delta,s_{\bar{\delta}'}) $, and total spin J is given by

where $ E_n $ denotes the eigenvalue of the radial Schrödinger equation,

where $ L=0 $ is imposed throughout. Eq. (15) is solved numerically on a discretized radial grid with Dirichlet boundary conditions $ \psi(0)=\psi(\infty)=0 $. We use a shooting algorithm that matches the logarithmic derivative of the wave function, with an adaptive step size to ensure numerical stability in the strongly confining regime.

Before proceeding, we clarify the terminology. The quantity M in Eq. (14) is the eigenvalue of the single-channel Hamiltonian in Eq. (10), in which the confining $ \Sigma_g^+ $ potential supports a purely discrete spectrum. It is not the position of a pole of the physical $ B^{(*)}\bar{D}^{(*)} $ scattering amplitude. A rigorous distinction among a bound state, a virtual state, and a resonance requires analytic continuation of a scattering amplitude across the relevant thresholds, which the present framework does not provide. Accordingly, in the following comparisons of M with two-meson thresholds, the terms “bound-state candidate” and “resonance candidate” indicate only the position of the single-channel eigenvalue relative to the threshold in question, along with the interpretation suggested by that position. Whether the corresponding pole of the full amplitude is genuinely bound, virtual, or resonant—and on which Riemann sheet it resides—can be determined only through a coupled-channel scattering analysis, which lies beyond the scope of this work.

We characterize the spatial structure of the δ–$ \bar{\delta}' $ wave functions using two complementary observables: the mean separation $ \langle r\rangle $ and the inverse-radius scale $ \langle 1/r\rangle^{-1} $. These quantities probe different length scales of the bound state and together provide a more discriminating description than a single radius. The expectation value $ \langle r\rangle $ is sensitive to the long-distance behavior of the wave function and therefore reflects the extent of the confining color flux tube connecting the diquark and antidiquark. In contrast, $ \langle 1/r\rangle^{-1} $ is weighted toward small separations and thus probes the short-distance region where the attractive part of the interaction dominates. This distinction is particularly informative for a Cornell-type potential of the form $ V(r) = -\alpha/r + \sigma r $, in which the color-Coulomb term $ -\alpha/r $ governs the short-range dynamics, whereas the linear term $ \sigma r $ is responsible for long-range confinement. In this framework, $ \langle 1/r\rangle^{-1} $ provides a direct measure of the typical strength of the short-range attraction experienced by the cluster pair, while $ \langle r\rangle $ characterizes the overall size of the confined system. The combination of these observables is especially useful for distinguishing a compact diquark–antidiquark configuration from a loosely bound hadronic molecule. Molecular states generally exhibit both a large mean separation and a large $ \langle 1/r\rangle^{-1} $, with both quantities typically of order 1 fm or more, whereas compact $ \delta\bar{\delta}' $ states have $ \langle 1/r\rangle^{-1} $ well below the hadronic scale, indicating genuine short-range overlap between the clusters.

Combining Eqs. (12), (13), and (14) with the state definitions in Eqs. (4)–(8) yields closed-form expressions for the physical masses of the isoscalar ground states $ T_{bc}^{(0)} $ and $ T_{bc}^{(1)} $:

where $ m_{[\bar{b}\bar{c}]} $ and $ m_{\{\bar{b}\bar{c}\}} $ denote, respectively, the antisymmetric (scalar) and symmetric (axial-vector) heavy-antidiquark masses, and $ E_{0}^{(0,1)} $ denotes the BO ground-state eigenvalue computed using the corresponding reduced mass $ \mu^{(0,1)} $. The corresponding isovector masses are

whereas the fourth isovector state in Eq. (9), which is built from the scalar heavy antidiquark, has

The threefold degeneracy of the $ \widetilde{T}_{bc}^{(J)} $ states arises directly because all three are based on the same $ |1_{[ud]},1_{\{\bar{b}\bar{c}\}}\rangle $ internal configuration and differ only in the coupling of the two spin-1 clusters to total $ J=0,1,2 $. Since $ H_\kappa $ of Eq. (12) acts only on the internal spins of the individual clusters and is insensitive to how the cluster spins are subsequently combined, it assigns the same eigenvalue to all three isovector members. This degeneracy is lifted only by spin-orbit, tensor, and higher-order interactions, which are not included in the present analysis. By contrast, the state $ \widetilde{T}_{bc}^{\prime(1)} $ is not degenerate with them: subtracting Eq. (19) from Eq. (18) yields $ \bigl(m_{\{\bar{b}\bar{c}\}}-m_{[\bar{b}\bar{c}]}\bigr) + \bigl(E_{0}^{(1)}-$ $E_{0}^{(0)}\bigr) + 2\kappa_{\bar{b}\bar{c}} $, i.e. precisely the isoscalar hyperfine splitting of Eq. (20). Thus, the two $ I=1 $, $ J^P=1^+ $ levels reproduce the $ 0^+ $–$ 1^+ $ pattern of the isoscalar sector. Since both carry the same conserved quantum numbers, interactions beyond $ H_\kappa $ will generally mix them; such mixing is not included here.

Two important mass splittings follow immediately from Eqs. (16)–(18). Within the isoscalar sector, the axial-vector–scalar splitting receives a kinematic contribution from the heavy-antidiquark mass gap and a chromomagnetic contribution from the heavy-antidiquark spin-spin coupling.

With the inputs adopted in Sec. III, the kinematic gap dominates ($ m_{\{\bar{b}\bar{c}\}}-m_{[\bar{b}\bar{c}]}\simeq 40 $ MeV), the BO contribution $ E_{0}^{(1)}-E_{0}^{(0)} $ is at the sub-MeV level because the two reduced masses are very similar, and the chromomagnetic term $ 2\kappa_{\bar{b}\bar{c}}\in[20,40] $ MeV. Therefore, the slope $ \partial\Delta_{HF}/\partial\kappa_{\bar{b}\bar{c}}=2 $ is the only quantity that responds to variations in the chromomagnetic coupling at fixed kinematic input. A measurement of $ \Delta_{HF} $ determines $ \kappa_{\bar{b}\bar{c}} $ only after the heavy-antidiquark mass gap has been independently fixed. The isoscalar–isovector gap, however, depends on which $ J^P $ channel is compared. Taking the two $ J^P=1^+ $ states, one finds

which is governed solely by the light-diquark spin–spin coupling and directly reflects the energy required to promote the light $ [ud] $ pair from the good to the bad configuration; the heavy-antidiquark mass and BO contributions cancel because both states are built on the symmetric $ \{\bar{b}\bar{c}\} $ cluster. By contrast, comparing the two $ J^P=0^+ $ states gives

This is because the isoscalar $ 0^+ $ state is built on a scalar heavy antidiquark, whereas the isovector $ 0^+ $ state is built on an axial-vector one; thus, the gap includes the antisymmetric–symmetric antidiquark mass difference as well as the spin-spin energy cost of both cluster transitions. For the parameter ranges adopted in Sec. III, Eqs. (21) and (22) give isoscalar–isovector gaps of $ \simeq 2\kappa_{ud}\approx 206 $ MeV and $ \simeq 40+2\kappa_{ud}+2\kappa_{\bar{b}\bar{c}}\approx 266 $–286 MeV, respectively. These values indicate that the isoscalar channels are by far the most favorable for the formation of strong-interaction-stable $ T_{bc} $ tetraquarks. The diquark masses $ m_{[ud]} $ and $ m_{[\bar{b}\bar{c}]} $, along with the spin-spin couplings $ \kappa_{ud} $ and $ \kappa_{\bar{b}\bar{c}} $, are taken as external inputs from QCD sum rules and meson phenomenology, as detailed in Sec. III

We now present the numerical results of the dynamical diquark model for the $ T_{bc} $ system in the $ \Sigma_g^+(1S) $ multiplet. For the light diquark, we adopt $ m_{[ud]} = 0.64 \pm 0.06 $ GeV from the QCD sum-rule analysis [52]. For the heavy antidiquark, we use $ m_{[\bar{b}\bar{c}]} = 6.38 $ GeV and $ m_{\{\bar{b}\bar{c}\}} = 6.42 $ GeV, obtained using the Regge-trajectory approach [53]. For the effective spin–spin coupling of the light diquark, we take $ \kappa_{ud} = 103 $ MeV, as determined from the $ \Sigma_Q $–$ \Lambda_Q $ baryon splittings in the diquark–antidiquark analysis of Ref. [54]. To assess the sensitivity of the spectrum to the heavy-antidiquark spin–spin interaction, we consider three representative parameter sets, labeled Set I, Set II, and Set III, with $ \kappa_{\bar{b}\bar{c}} = 10 $, 15, and 20 MeV, respectively:

All three values lie well within the phenomenological range reported in the QCD Laplace sum-rules analysis [55]. The light-diquark mass, the heavy-antidiquark masses, and the light-diquark spin–spin coupling $ \kappa_{ud} $ are held fixed across all three sets. The BO potential $ V_{\Sigma_g^+}(r) $ is taken from the lattice-QCD parametrizations of Refs. [31, 38], and the radial Schrödinger equation is solved numerically as described in Sec. II. For ease of reference, Table 1 summarizes all model inputs, their adopted values, sources, quoted uncertainties, and roles in the calculation.

The inputs collected in Table 1 originate from four distinct frameworks—QCD sum rules, Regge phenomenology, baryon spectroscopy, and lattice QCD—so it is necessary to explain why this is unavoidable and in what sense the resulting set remains internally consistent. The Hamiltonian in Eq. (10) requires three logically distinct types of input. The first is the static energy of a color-$ \bar{{\bf{3}}}_c $ source in the presence of a color-$ {\bf{3}}_c $ antisource. This is a gauge-invariant, well-defined quantity, and lattice QCD is the only first-principles source for it. The second is the masses of the two clusters. These are not observables: a diquark is a colored, gauge-variant object with no asymptotic state. Therefore, they can be defined only within a model or through a correlator supplemented by gauge fixing or a Wilson line, as in the QCD sum-rule and Regge analyses employed here. The third is the short-range chromomagnetic couplings, which are defined operationally through their action on measured hadron splittings. No single calculation supplies all three inputs, and the same allocation of inputs is used in every application of the dynamical diquark model, in which the BO potentials are taken from lattice QCD while the diquark masses and couplings κ are taken from spectroscopy [47–50]. Consistency in this context therefore cannot require a common source. Instead, it requires a common definition: each input must be precisely the quantity appearing in Eq. (10)—the mass of a compact color-triplet cluster treated as a static source, and a contact spin-spin coupling between the constituents of that cluster. The determinations of Refs. [52–55] all satisfy this requirement, and it is in this sense that they may be combined.

A second, more technical point governs how much of the residual scheme dependence is physical. The spectrum generated by Eqs. (14)–(15) depends on the cluster masses and the BO potential only through the combination $ m_\delta + m_{\bar{\delta}'} + V_{\Sigma_g^+}(r) $: the substitution $ V_{\Sigma_g^+}(r)\to V_{\Sigma_g^+}(r)+\Delta $ together with $ m_\delta+m_{\bar{\delta}'}\to m_\delta+m_{\bar{\delta}'}-\Delta $ leaves every eigenvalue and every wave function unchanged. Constituent cluster masses are therefore meaningful only relative to the additive convention adopted for the static energy, and the lattice $ \Sigma_g^+ $ potential carries precisely such a scheme-dependent self-energy constant. A mismatch between the conventions of two determinations consequently appears as a common shift of the entire spectrum rather than as a distortion of it. This is the technical reason why absolute masses carry a substantially larger systematic uncertainty than splittings. It also implies that quantities defined as differences—the chromomagnetic couplings and the cluster mass gaps—are transferable between frameworks, whereas absolute cluster masses are not. We accordingly cross-check the former and assign a conservative uncertainty to the latter.

The least constrained input, $ \kappa_{\bar b\bar c} $, can be checked against a determination entirely independent of the sum rules of Ref. [55]. Because the color factor $ \langle{\boldsymbol{\lambda}}_i\cdot{\boldsymbol{\lambda}}_j\rangle $ equals $ -8/3 $ for a $ \bar{{\bf{3}}}_c $ diquark but $ -16/3 $ for a color-singlet meson, the chromomagnetic coupling inside a heavy antidiquark is, apart from differences in the wave function at the origin, one half of its mesonic counterpart.

Lattice QCD predicts $ M(B_c^{*})-M(B_c)=54(3) $ MeV [56], whereas quark-model estimates range up to approximately 70 MeV, giving $ \kappa_{\bar b\bar c}\simeq 13 $–18 MeV. The interval 10–20 MeV examined in Sets I–III therefore brackets this independent estimate, with Set II close to its midpoint. Thus, the spectral spread reported below is representative of the underlying uncertainty in this coupling rather than an arbitrary choice.

One source-dependent ambiguity, however, directly affects a physical quantity, and we quantify it here. It concerns the allocation of the heavy-antidiquark spin dependence, which, in the parameter set of Table 1, enters in two places: through the two distinct Regge masses $ m_{[\bar b\bar c]} $ and $ m_{\{\bar b\bar c\}} $, whose 40 MeV gap is itself of chromomagnetic origin, and through the term $ 2\kappa_{\bar b\bar c}({\bf{s}}_{\bar b}\cdot{\bf{s}}_{\bar c}) $ of Eq. (12). Two limiting conventions bracket the possibilities. In convention (A), used for our central results, the two contributions are treated as physically distinct: the Regge gap measures the spin splitting of the isolated $ \bar b\bar c $ cluster, whereas $ \kappa_{\bar b\bar c} $ describes the residual chromomagnetic interaction inside the tetraquark, where the cluster is embedded in the color field of the light diquark and its wave function at the origin is modified. In convention (B), which is normally adopted in dynamical-diquark analyses [47, 49], the cluster is assigned a single spin-averaged mass.

and all spin dependence is generated by $ H_\kappa $; the Regge gap is then reproduced for $ 2\kappa_{\bar b\bar c}=40 $ MeV, i.e. $ \kappa_{\bar b\bar c}=20 $ MeV, which is precisely the value adopted in Set III. Convention (A) maximizes the hyperfine splitting, whereas convention (B) minimizes it; thus, the physical result is bracketed by these two conventions. Table 2 compares the resulting spectra, obtained from Eqs. (16)–(17) with the same radial eigenvalue. The absolute masses differ by at most $ +30 $ MeV in the $ 0^{+} $ channel and $ -10 $ MeV in the $ 1^{+} $ channel, well within the parametric uncertainty established below. Meanwhile, $ \Delta_{HF} $ ranges from $ 2\kappa_{\bar b\bar c}=20 $–40 MeV in convention (B) to 59–79 MeV in convention (A). We quote the latter as our central values but regard $ \Delta_{HF}\simeq 20 $–80 MeV as the interval allowed by the present inputs. The qualitative picture of Sec. III C is only mildly affected: in both conventions, the axial-vector state lies above $ B^{*}\bar{D} $ and well below $ B\bar{D}^{*} $, whereas the scalar moves from being essentially at the $ B\bar{D} $ threshold in convention (A) to 24–39 MeV above it in convention (B). Since this displacement is smaller than the uncertainty in the absolute masses, both conventions are compatible with a scalar located at or slightly above the $ B\bar{D} $ threshold, and neither establishes binding.

Before presenting the spectrum, we assess the total parametric uncertainty associated with Table 1, because the treatment of uncertainties determines how the near-threshold comparisons below should be interpreted. The dominant source is the light-diquark mass. Propagating $ \delta m_{[ud]}=\pm 60 $ MeV shifts the absolute masses of both $ T_{bc}^{(0)} $ and $ T_{bc}^{(1)} $ by up to $ \pm 60 $ MeV; the direct constituent shift is partially compensated by the induced change in the BO eigenvalue through the reduced mass. The heavy-antidiquark masses carry uncertainties of a few tens of MeV, typical of Regge determinations. The lattice-based parametrization of the BO potential contributes an estimated $ {\cal{O}}(10) $ MeV through the offset and string tension entering $ E_0 $, and the neglected color mixing and coupled-channel effects discussed in Secs. II B and III C act at a comparable scale. A realistic total uncertainty on the absolute masses is therefore of order $ \pm 60 $–70 MeV, substantially larger than the spread induced by $ \kappa_{\bar b\bar c} $ alone. Crucially, however, these dominant uncertainties are strongly correlated between the two isoscalar channels: $ m_{[ud]} $, $ \kappa_{ud} $, and the BO potential enter $ M(T_{bc}^{(0)}) $ and $ M(T_{bc}^{(1)}) $ identically and therefore cancel in the hyperfine splitting $ \Delta_{HF} $, whose uncertainty is instead controlled by the heavy-antidiquark mass gap and by $ \kappa_{\bar b\bar c} $. The robust outputs of the calculation are thus the mass differences, the ordering of the states relative to one another and to the pattern of thresholds, and the parametric dependences. The absolute proximity of any individual eigenvalue to a given threshold at the few-MeV level is a statement based on central values and should not be overinterpreted.

A related point, which deserves to be stated explicitly, concerns the comparison with physical thresholds. The framework adopted here describes only the relative motion of the δ–$ \bar\delta' $ pair in the BO potential. The $ B^{(*)} $ and $ \bar D^{(*)} $ mesons are quark–antiquark color singlets that lie outside this Hilbert space, so the model neither computes nor is calibrated to the physical meson masses. The thresholds quoted in Sec. III C are fixed by the experimental PDG masses, whereas the tetraquark masses are built from independently determined constituent inputs. Statements at the few-MeV level about the position of an eigenvalue relative to a threshold therefore compare quantities of different origin and are meaningful only at the level of central values, within the parametric uncertainty budget established above. This limitation is inherited by all constituent-model analyses of near-threshold states and is made explicit here so that the reader can weigh the conclusions accordingly.

The mass spectra and characteristic spatial scales of the doubly heavy open-flavor tetraquark states $ T_{bc}^{(0)} $ ($ J^{P}=0^{+} $) and $ T_{bc}^{(1)} $ ($ J^{P}=1^{+} $), obtained for the three parameter sets, are summarized in Tables 3 and 4.

Across the three parameter sets, the scalar state is predicted within the narrow range $ M(T_{bc}^{(0)}) = 7.143\text{–}7.158 $ GeV, whereas the axial-vector state lies at $ M(T_{bc}^{(1)}) = 7.217\text{–}7.222 $ GeV. The total spread induced by varying $ \kappa_{\bar{b}\bar{c}} $ over the range 10–20 MeV is only $ \Delta M \simeq 15 $ MeV in the $ 0^{+} $ channel and $ \Delta M \simeq 5 $ MeV in the $ 1^{+} $ channel. This mild sensitivity is expected on general grounds. In the dynamical diquark picture, the tetraquark mass is dominated by the constituent masses $ m_{[\bar{b}\bar{c}]} $, $ m_{\{\bar{b}\bar{c}\}} $, and $ m_{[ud]} $ and by the BO potential $ V_{\Sigma_g^+}(r) $, all of which are kept fixed. The parameter $ \kappa_{\bar{b}\bar{c}} $ enters only through the heavy-antidiquark chromomagnetic term $ \kappa_{\bar{b}\bar{c}}\, \vec{S}_{\bar b} \cdot \vec{S}_{\bar c} $, which is diagonal in the heavy-antidiquark spin basis.

The two trends—the decrease in $ M(T_{bc}^{(0)}) $ and the increase in $ M(T_{bc}^{(1)}) $ as $ \kappa_{\bar{b}\bar{c}} $ grows—are therefore physically transparent. The scalar $ T_{bc}^{(0)} $ is built on the spin-0 antisymmetric antidiquark $ [\bar{b}\bar{c}] $, for which $ \langle \vec{S}_{\bar b} \cdot \vec{S}_{\bar c}\rangle = -3/4 $, whereas the axial-vector $ T_{bc}^{(1)} $ contains the spin-1 symmetric antidiquark $ \{\bar{b}\bar{c}\} $, for which $ \langle \vec{S}_{\bar b} \cdot \vec{S}_{\bar c}\rangle = +1/4 $. An increase in $ \kappa_{\bar{b}\bar{c}} $ thus lowers the $ 0^{+} $ mass by $ 3\kappa_{\bar{b}\bar{c}}/2 $ and raises the $ 1^{+} $ mass by $ \kappa_{\bar{b}\bar{c}}/2 $ relative to the spin-independent baseline, in quantitative agreement with the results in Tables 3–4. The asymmetry between the two shifts—a factor of three stronger in the scalar channel—directly reflects the different Casimirs of the two antidiquark configurations and explains why Set III exhibits the largest hyperfine splitting while leaving $ \langle r\rangle $ essentially unchanged. Defining $ \Delta_{HF}\equiv M(T_{bc}^{(1)})-M(T_{bc}^{(0)}) $, we obtain $ \Delta_{HF} \simeq 59 $, 69, and 79 MeV for Sets I, II, and III, respectively. The dominant contribution to the absolute value of $ \Delta_{HF} $ is the $ \sim 40 $ MeV mass gap between the symmetric and antisymmetric heavy-antidiquark configurations, $ m_{\{\bar{b}\bar{c}\}}-m_{[\bar{b}\bar{c}]} $, which is independent of $ \kappa_{\bar{b}\bar{c}} $. The chromomagnetic coupling adds the term $ 2\kappa_{\bar{b}\bar{c}} $, so that each 5 MeV increment in $ \kappa_{\bar{b}\bar{c}} $ generates the expected 10 MeV shift in $ \Delta_{HF} $ to within a fraction of an MeV, in agreement with the analytic slope $ \partial\Delta_{HF}/\partial\kappa_{\bar{b}\bar{c}}=2 $ derived from Eq. (20). Once the heavy-antidiquark mass gap has been fixed by independent input, this linear response to $ \kappa_{\bar{b}\bar{c}} $ provides a direct means by which future experimental data can constrain the chromomagnetic coupling in a regime that has so far been inaccessible.

The two radial observables reported in the tables encode complementary information about the internal structure. The mean inter-cluster separation $ \langle r\rangle $ characterizes the average size of the heavy-antidiquark–light-diquark configuration, weighted by the long-distance tail of the radial probability density $ r^{2}|\psi(r)|^{2} $. The inverse mean radius $ \langle 1/r\rangle^{-1} $ is most sensitive to the short-distance region of the wave function, where the Coulomb-like one-gluon exchange dominates. Both quantities are independent of the parameter set at the displayed precision, reflecting the fact that $ \kappa_{\bar{b}\bar{c}}\langle \vec{S}\cdot\vec{S}\rangle\sim {\cal{O}}(20\; \text{MeV}) $ is a genuinely small perturbation compared with the typical kinetic and potential scales of the $ \Sigma_g^+(1S) $ state. Moreover, in the BO formulation adopted here, the chromomagnetic interaction enters the radial Schrödinger equation only as an additive constant, leaving the wave function itself unaffected at this level.

At the displayed precision, the two channels differ by one unit in the last quoted digit.

However, this difference should not be assigned physical significance, as we now show. The only quantity that distinguishes the two radial problems is the reduced mass that enters the radial Schrödinger equation,

where the antisymmetric and symmetric heavy-antidiquark masses, $ m_{[\bar{b}\bar{c}]}=6.38 $ GeV and $ m_{\{\bar{b}\bar{c}\}}=6.42 $ GeV, respectively, differ by 40 MeV. Because the light diquark is more than an order of magnitude lighter than the heavy antidiquark, this 40 MeV gap corresponds to a fractional change in the reduced mass of only $ \Delta\mu/\mu \simeq 5\times 10^{-4} $. Moreover, both components of the BO potential imply a size that decreases monotonically with increasing μ—as $ a_0\sim 1/(\mu\alpha) $ in the Coulombic limit and as $ \langle r\rangle \sim (\mu\sigma)^{-1/3} $ in the linear limit—so the axial-vector state, which has the slightly larger reduced mass, is expected to be very slightly more compact than the scalar state, by a relative amount of order $ 10^{-4} $, i.e. by approximately $ 10^{-4} $ fm. This is three orders of magnitude below the two-digit precision at which the radii are quoted in Tables 3–4. The one-unit difference in the last digit is therefore a rounding artifact of the numerical grid and has neither physical magnitude nor physical sign. The meaningful statement is that, within the accuracy of the present calculation, the two states share the same $ \Sigma_g^+(1S) $ radial wave function: the chromomagnetic term enters Eq. (15) only as an additive constant, and the heavy-antidiquark mass gap is negligible on the scale of the reduced mass.

The two states share an essentially common shape, as quantified by the dimensionless ratio

which, by Jensen’s inequality ¹ applied to the convex function $ 1/r $, is bounded above by unity (saturated only for a delta-like distribution). Both values lie well between the pure Coulombic ($ \sim 2/3 $) and pure harmonic-oscillator ($ \sim 0.85 $) regimes, indicating that, within the present model, the short-range gluon exchange and the long-range confining string contribute comparably to the binding in both $ J^{P} $ sectors—a pattern familiar from heavy quarkonium phenomenology and found here to extend to the diquark–antidiquark description of four-quark configurations.

The compactness of the wave functions— $ \langle r\rangle\simeq 0.45 $–0.46 fm, significantly smaller than the typical confinement scale ($ \sim 1.0 $ fm) and well below the deuteron-like scale of loosely bound molecular candidates—is consistent with an interpretation of $ T_{bc}^{(0)} $ and $ T_{bc}^{(1)} $ as compact diquark–antidiquark configurations rather than extended hadronic molecules. We emphasize, however, that this conclusion is not model-independent: the radii are computed from the wave function generated within the assumed δ–$ \bar\delta' $ framework, which imposes a compact configuration from the outset, and they cannot by themselves provide an unambiguous discriminator between compact and molecular structures. An extended molecular component, if present in the physical state through the near-threshold coupling discussed below, would lie outside the Hilbert space of the present description. The radii quoted here should therefore be interpreted as characterizing the compact core of the state within the model, not as a measurement of the physical size of the state.

We further note that within a molecular interpretation the $ J^{P} $ ordering of the radii would be governed by the binding-energy hierarchy of the $ B^{(*)}\bar D^{(*)} $ thresholds rather than by the heavy-antidiquark mass splitting, and would be expected to display a substantially larger relative spread than the essentially exact degeneracy of the two radial profiles found here. The near-coincidence of the two radial scales between the scalar and axial-vector channels, together with their insensitivity to the parameter set, thus reinforces the picture of a common $ \Sigma_g^+(1S) $ spatial profile shared by the entire $ T_{bc} $ multiplet, perturbed only marginally by the heavy-antidiquark mass splitting and the chromomagnetic coupling.

The lowest meson–meson thresholds with the same flavor and $ J^{P} $ quantum numbers as the $ T_{bc} $ system, computed using PDG-averaged meson masses [57], are

The last of these, which carries $ J^{P}=0^{+},1^{+},2^{+} $, plays no role in the isoscalar ground states discussed in this subsection, because they lie far below it. It is, however, the relevant reference for the isovector members of the multiplet considered in Sec. III D. A compact visual summary of the resulting spectroscopic situation is provided in Fig. 1, in which the predicted masses are displayed alongside these thresholds.

Figure 1.  Predicted masses for the doubly heavy tetraquarks $ T_{bc}^{(0)} $ ($ J^{P}=0^{+} $, left panel) and $ T_{bc}^{(1)} $ ($ J^{P}=1^{+} $, right panel), obtained using the three parameter sets employed in this work (solid color-coded lines: Set I, blue; Set II, green; Set III, orange), compared with the lowest open-flavor meson–meson thresholds with the same quantum numbers (red dashed lines). All masses are given in MeV.

The figure shows two qualitatively distinct dynamical situations. In the scalar sector (left panel), the three central values bracket the $ B\bar D $ threshold at 7149.0 MeV, with displacements of $ +9 $, $ +2 $, and $ -6 $ MeV for Sets I, II, and III, respectively. These displacements are much smaller than the $ {\cal{O}}(60\; \text{MeV}) $ parametric uncertainty on the absolute masses established above. Thus, the calculation supports the conclusion that the scalar eigenvalue is compatible with the $ B\bar D $ threshold, not that its position relative to threshold is determined. At the central input values, $ T_{bc}^{(0)} $ effectively coincides with the two-meson threshold, and whether the physical state would emerge as a genuinely bound state or as a narrow near-threshold structure cannot be decided within the present single-channel framework.

Should $ M(T_{bc}^{(0)}) \lt M_B + M_{\bar D} $, as Set III already suggests, S-wave hadronic decay into $ B\bar D $ is kinematically forbidden. The remaining decay modes would be the weak $ b\to c\ell\bar\nu $ and $ b\to c\bar ud $ transitions, the radiative $ T_{bc}^{(0)}\to B\bar D\gamma $ channel, and possible three-body modes through the off-shell tails of $ B^{*} $ and $ \bar D^{*} $. For the analogous $ T_{cc}^{+}(3875) $ system, the latter mechanism is in fact the dominant decay channel, with $ T_{cc}^{+}\to D^{0}D^{0}\pi^{+} $ proceeding through the off-shell tail of the $ D^{*+} $. For $ T_{bc}^{(0)} $, the corresponding $ B^{*} $ and $ \bar D^{*} $ tails lie $ \sim 45 $ and $ \sim 140 $ MeV above the elastic threshold, respectively, and are therefore kinematically more suppressed than the $ D^{*+} $ tail in the charm system. If the resulting radiative and three-body modes could be neglected, the remaining weak width $ \Gamma\sim {\cal{O}}(10^{-10}) $ MeV would constitute a distinctive signature of a strong-interaction-stable tetraquark. We stress that this scenario is conditional on $ M(T_{bc}^{(0)}) $ lying below threshold, which the present calculation cannot establish, and that a dedicated estimate of the radiative and three-body channels would be required before any quantitative statement about the width can be made. The near-degeneracy visible in Fig. 1 therefore strongly motivates a dedicated coupled-channel analysis that explicitly incorporates the $ B\bar D $ continuum, since the physical pole may migrate from the real axis onto the second Riemann sheet depending on the residual $ B\bar D $ interaction.

In the axial-vector sector (right panel), by contrast, the three central values form a tight cluster around 7220 MeV. They lie above $ B^{*}\bar D $ by 23–28 MeV and well below $ B\bar D^{*} $ by $ \sim 70 $ MeV for all parameter sets. Within the single-channel framework, this places $ T_{bc}^{(1)} $ as a natural candidate for an S-wave resonance in the $ B^{*}\bar D $ channel, with the caveat, stated in Sec. II C, that an eigenvalue above a two-meson threshold does not by itself establish a resonance. Whether the corresponding pole of the scattering amplitude survives as a resonance, and at what position and width, can only be determined from the amplitude itself. Should such a resonance exist near the eigenvalue found here, the $ \sim 25 $ MeV excess above threshold would be large enough to open the $ B^{*}\bar D $ decay channel but small enough that phase-space suppression and the proximity of the $ B^{*}\bar D $ cusp would visibly distort the invariant-mass line shape. This situation is reminiscent of $ T_{cc}^{+}(3875) $ relative to $ D^{*+}D^{0} $, so a Flatté-type or explicitly coupled-channel parametrization would be needed in the experimental analysis.

The expected size of the neglected coupled-channel effects can be estimated on general grounds. The δ–$ \bar\delta' $ configuration couples to the open and nearby closed $ B^{(*)}\bar D^{(*)} $ channels through light-quark rearrangement, and this coupling generates a self-energy whose real part shifts the bare eigenvalue. For a compact state coupled in S-wave to a two-meson continuum, the magnitude of this shift is set by the same coupling that controls the hadronic width of an above-threshold state and the level repulsion from a below-threshold channel. For hadronic couplings of typical strength and the near-threshold momenta relevant here ($ k \lesssim 300 $ MeV), both dimensional analysis and explicit coupled-channel studies of analogous doubly-heavy systems—including the $ T_{cc}^{+}(3875) $, for which the compact-core pole is drawn toward the $ D^{*+}D^{0} $ threshold by the channel coupling—indicate shifts of order 10–30 MeV. These shifts are predominantly attractive, i.e. directed toward or below the nearest threshold. Three consequences follow. First, for the scalar state, a downward shift of this size would move an eigenvalue lying a few MeV above $ B\bar D $ below threshold, reinforcing the conclusion that the bound-versus-unbound character of $ T_{bc}^{(0)} $ cannot be decided at the present level of analysis. Second, for the axial-vector state, an attractive shift of 10–30 MeV would move the pole from $ \sim 25 $ MeV above $ B^{*}\bar D $ into the immediate vicinity of that threshold, substantially reducing—and possibly closing—the gap to the lattice determinations discussed below. Third, since these shifts are comparable to the displacements from threshold on which the qualitative classification of the states rests, the coupled-channel dynamics is not a small correction in this system. The classifications offered here must therefore be regarded as provisional until a calculation that explicitly incorporates the $ B^{(*)}\bar D^{(*)} $ continua becomes available.

The contrast between the two channels—a quasi-bound scalar versus a near-threshold resonant axial-vector—is a direct consequence of the interplay between the 60–80 MeV hyperfine splitting found here and the heavy-quark-symmetry-breaking pattern of the surrounding open-flavor spectrum. Because the $ D^{*} $–D hyperfine gap ($ \sim 140 $ MeV) is substantially larger than the $ B^{*} $–B gap ($ \sim 45 $ MeV), exciting the charmed meson costs markedly more energy than exciting the bottom meson. This pushes $ B\bar D^{*} $ far above $ B^{*}\bar D $ and opens a wide window in which $ T_{bc}^{(1)} $ can be trapped as a narrow resonance. Fig. 1 thereby encapsulates, at a glance, both the internal tetraquark dynamics and the external heavy-quark-symmetry-breaking pattern that shapes the accessible decay channels. Any future observation that inverts this hierarchy would point to physics beyond the diquark–antidiquark picture adopted in this work.

Our results can be compared with a growing body of theoretical estimates obtained using a variety of independent methods. In QCD sum-rule analyses, Ref. [24] finds the $ T_{bc} $ states in the range 7.0–7.2 GeV. Constituent quark-model calculations, including the relativistic diquark-antidiquark treatment of Ref. [16] and the chromomagnetic-interaction analyses of Refs. [6, 58], cluster around 7.1–7.3 GeV. Ref. [58] notes that the $ I(J^P)=0(1^+) $ $ T_{bc} $ state lies close to, but possibly above, the $ B^*\bar{D} $ threshold. A similar conclusion is reached in Ref. [17] within a constituent quark model with parameters fitted to the $ \Xi_{cc}^{++} $ mass. More recently, the improved chromomagnetic-interaction model of Ref. [59] and the quark-model study of Ref. [18] find $ T_{bc} $ ground-state masses in the range 7.15–7.25 GeV. Our values fall comfortably within this band, while the hyperfine splitting $ \Delta_{{\rm{HF}}}\sim 60 $–80 MeV is somewhat smaller than the $ \sim 100 $ MeV typical of pure chromomagnetic estimates. This difference can be traced to the inclusion of the spatial wave-function dependence of the spin-spin operator and to the moderate values of $ \kappa_{\bar{b}\bar{c}} $ adopted in the present work.

A more demanding benchmark is provided by recent lattice-QCD results, which we now compare systematically. Table 5 summarizes, for each available lattice calculation of the isoscalar $ T_{bc} $ system, the interpolating-operator basis, the channels studied, and the quantitative outcome, alongside the eigenvalues obtained in this work. The table makes explicit that the residual differences among the lattice studies themselves—binding by tens of MeV in Ref. [27] versus at most a few MeV in Ref. [28]—are correlated with the operator content (local versus nonlocal scattering operators), the treatment of the left-hand cut from one-pion exchange, and the chiral and continuum extrapolations. They are therefore of the same nature as the systematic effects that separate those results from the present constituent calculation.

The two-meson-plus-diquark variational analysis of Ref. [27] reports a binding energy of $ -43(^{+6}_{-7})(^{+14}_{-24}) $ MeV below $ B^{*}\bar{D} $ for the $ I(J^P)=0(1^+) $ ground state, corresponding to a physical mass of $ \simeq 7151 $ MeV when added to the experimental $ B^{*}\bar{D} $ threshold of 7194.4 MeV. By contrast, the $ B^{(*)}\bar{D} $-scattering analysis of Ref. [28] locates the near-threshold pole at most a few MeV below $ B^{*}\bar{D} $, i.e. near 7190 MeV. Our axial-vector eigenvalue $ M(T_{bc}^{(1)})\simeq 7220 $ MeV therefore lies $ \sim 70 $ MeV above the lattice binding-energy value of Ref. [27] and $ \sim 30 $ MeV above the scattering-pole value of Ref. [28], placing it on the opposite side of the $ B^{*}\bar{D} $ threshold in both cases. This $ {\cal{O}}(30 $–$ 70) $ MeV gap lies within the systematic uncertainty generically expected for a constituent BO description. In such a description, the dominant systematic uncertainty, of order $ \Lambda_{{\rm{QCD}}} $, primarily affects fine-structure quantities that determine the position of the pole relative to threshold. We therefore view our result as the constituent δ–$ \bar{\delta}' $ value, while a comparison with lattice results at better than $ {\cal{O}}(\Lambda_{{\rm{QCD}}}) $ accuracy lies beyond the scope of the present analysis. In this connection, we note that the single-channel $ \Sigma_g^+(1S) $ description adopted here omits two physical effects that are explicitly included in the lattice analyses and are natural candidates for explaining the residual gap: long-range one-pion exchange between the $ B^{(*)} $ and $ \bar D^{(*)} $ clusters, which generates a left-hand cut in the elastic amplitude, and explicit coupled-channel $ B^{(*)}\bar D^{(*)} $ dynamics, which can generically pull the physical pole downward toward—or below—the nearest two-meson threshold. A quantitative incorporation of either ingredient lies beyond the present single-channel BO treatment.

From an experimental standpoint, a resonant structure in the $ B\bar D $ invariant-mass distribution near 7.15 GeV, displaced from threshold by at most a handful of MeV, would be a natural candidate for $ T_{bc}^{(0)} $, whereas a structure in $ B^{*}\bar D $ around 7.22 GeV, roughly 25 MeV above threshold, would be the expected signature of $ T_{bc}^{(1)} $. The absence of a nearby eigenvalue in the $ B\bar D^{*} $ channel further suggests that the axial-vector state, if observed, could exhibit a pronounced branching asymmetry in favor of $ B^{*}\bar D $ over $ B\bar D^{*} $. Such an asymmetry would provide an additional discriminator against competing molecular interpretations; a quantitative statement, however, would require explicit decay-width calculations that are not performed here. With the high-luminosity runs of LHCb and the upcoming Belle II dataset, $ bc\bar q\bar q $ states in the 7.1–7.3 GeV window are within reach, and confirmation or refutation of the hierarchy and near-threshold character suggested here would significantly sharpen our understanding of multiquark dynamics in the doubly-heavy sector.

The two isoscalar states discussed thus far exhaust neither the $ \Sigma_g^+(1S) $ multiplet nor the color structures available to a $ bc\bar q\bar q $ system. In this subsection, we complete the spectrum of this multiplet and estimate the effect of the $ \bar{{\bf{6}}}\otimes{\bf{6}} $ configuration that was set aside in Sec. II B.

The four isovector members of the multiplet are obtained from Eqs. (18) and (19) without any additional numerical input. Comparison with Eqs. (16)–(17) yields the exact relations

Consequently, the entire isovector sector is obtained from the isoscalar sector by a uniform upward shift of $ 2\kappa_{ud} \approx 206 $ MeV, corresponding to the energy required to promote the light $ [ud] $ pair from the good to the bad diquark configuration. This is a structural prediction of the isospin–spin lock in Eq. (3): the same operator that generates the good–bad diquark splitting also determines the isoscalar–isovector gap, without introducing any ad hoc isospin-dependent interaction. The resulting spectrum of the full multiplet is presented in Table 6.

Three features of Table 6 deserve comment. First, every isovector state lies above all four open-flavor thresholds of Eqs. (34–37), including $ B^{*}\bar D^{*} $. Consequently, all S-wave decay channels allowed by the quantum numbers are open, with sizable phase space, and these states are expected to be broad. This expectation is consistent with the absence of any lattice-QCD signal for isovector $ T_{bc} $ states. It also explains why the isoscalar channels—the only channels in which the light pair can occupy the good-diquark configuration—are those in which a narrow, potentially strong-interaction-stable structure can be expected. Second, the three states $ \widetilde{T}_{bc}^{(0,1,2)} $ are exactly degenerate at the level of $ H_\kappa $, as explained in Sec. II C; a physical splitting of this triplet requires spin-orbit and tensor forces, which vanish identically for $ L=0 $ in the Hamiltonian adopted here and would first be generated by mixing with excited BO potentials and by relativistic corrections. The degeneracy pattern of Table 6 is therefore a prediction of the leading BO approximation, not a statement about the physical spectrum at the 10 MeV level. Third, since $ \widetilde{T}_{bc}^{\prime(1)} $ and $ \widetilde{T}_{bc}^{(1)} $ carry identical conserved quantum numbers, the physical $ I=1 $, $ J^P=1^+ $ eigenstates will be mixtures of the two internal spin configurations once interactions beyond $ H_\kappa $ are included; the entries of Table 6 for these two states should be read accordingly.

Beyond the $ \Sigma_g^+(1S) $ multiplet, the same framework supports a tower of radial and orbital excitations. The first radially excited S-wave level of the $ \Sigma_g^+ $ potential lies roughly 0.6 GeV above the ground state for the relevant reduced masses. The excited BO potentials $ \Pi_u $ and $ \Sigma_u^- $ generate hybrid-like multiplets in the same region, so the 7.5–8.0 GeV window is expected to be densely populated. Since all of these states lie far above the open-flavor thresholds, and since the corresponding level orderings are considerably more sensitive to the BO potentials and to the neglected spin-dependent forces than the ground-state doublet, we do not attempt a quantitative treatment here.

As emphasized in Sec. II B, the argument for discarding the sextet channel on the grounds that the interaction between two (anti)quarks in the $ {\bf{6}}_c $ ($ \bar{{\bf{6}}}_c $) representation is repulsive is incomplete. The interaction between a $ {\bf{6}}_c $ cluster and a $ \bar{{\bf{6}}}_c $ cluster is attractive and, by Casimir scaling, stronger than the $ \bar{{\bf{3}}}_c $–$ {\bf{3}}_c $ interaction. Both effects can be assessed, at least semiquantitatively, within the present framework.

For the intercluster dynamics, Casimir scaling of the static energy gives

where $ V_0 $ denotes the additive self-energy constant, which is not Casimir-scalable and, as discussed in Sec. III A, is in any case degenerate with the cluster masses. Lattice measurements of static potentials in higher representations support Eq. (40) at the few-percent level for separations up to about 1 fm [45], which includes the region $ \langle r\rangle\simeq 0.45 $ fm relevant here. Solving Eq. (15) with the scaled potential and the same reduced mass, we find that the S-wave eigenvalue of the sextet channel lies 0.1–0.3 GeV above the antitriplet one. The spread reflects the Coulomb-to-linear partitioning of the parametrization: the factor $ 5/2 $ deepens the short-distance well but also increases the slope of the confining term, and at the separations relevant here, the latter effect dominates. The associated wave function is correspondingly more compact, with $ \langle r\rangle \simeq 0.3 $ fm. The intra-cluster energy penalty must also be included. One-gluon exchange between two quarks in the $ {\bf{6}}_c $ channel carries a color factor equal to $ -1/2 $ times that of the $ \bar{{\bf{3}}}_c $ channel, so the interaction is repulsive and no compact $ {\bf{6}}_c $ cluster is formed. Strictly speaking, the cluster description itself ceases to apply, and the energy of the configuration is raised further relative to the estimate above.

Taken together, these considerations place the diagonal energy of the $ \bar{{\bf{6}}}\otimes{\bf{6}} $ configuration several hundred MeV above that of the $ \bar{{\bf{3}}}\otimes{\bf{3}} $ configuration. Its effect on the ground state is therefore of second order,

i.e. a downward shift of approximately 10 MeV, reaching a few tens of MeV only for off-diagonal matrix elements at the upper end of the plausible range for a color-rearrangement interaction. Three consequences follow. First, the sextet channel does not constitute a competing leading configuration, so the restriction to $ \bar{{\bf{3}}}\otimes{\bf{3}} $ adopted here is justified as a leading approximation. Second, by the variational argument of Sec. II B, the neglected mixing acts in the direction of additional binding, so the eigenvalues of Tables 3–4 remain upper limits within the model. Third, the magnitude of the induced shift, $ {\cal{O}}(10) $ MeV, is comparable to that of the neglected coupled-channel effects estimated in Sec. III C and is already included in the uncertainty budget of Sec. III A. We stress that Eqs. (40) and (41) are estimates within a static, two-cluster picture; a genuine assessment requires a four-body calculation retaining the full color space, such as those performed in Refs. [18, 21–23], in which sizable components beyond a single $ \bar{{\bf{3}}}\otimes{\bf{3}} $ configuration are indeed found.

In this work, we carried out a systematic study of the doubly-heavy open-flavor tetraquark states $ T_{bc}^{(0)} $ ($ J^{P}=0^{+} $) and $ T_{bc}^{(1)} $ ($ J^{P}=1^{+} $) within the dynamical diquark model. We treated the $ T_{bc} $ system as a heavy antidiquark–light diquark pair propagating in the $ \Sigma_g^+(1S) $ BO potential extracted from lattice QCD. The light-diquark mass, the heavy-antidiquark masses of both spin configurations, and the light-diquark spin–spin coupling were fixed from independent phenomenological inputs. The heavy-antidiquark chromomagnetic coupling $ \kappa_{\bar{b}\bar{c}} $ was varied over the range 10–20 MeV to probe the sensitivity of the spectrum to the least constrained parameter of the model. A full parametric uncertainty budget for all inputs was presented in Sec. III A, together with an explicit discussion of why the inputs must originate from different frameworks, of the degeneracy between the constituent cluster masses and the additive constant of the static potential, and of the two limiting conventions for apportioning the heavy-antidiquark spin dependence, which bracket the hyperfine splitting between 20 and 80 MeV. We also completed the $ \Sigma_g^+(1S) $ multiplet by computing its isovector members and estimated the size of the neglected color mixing, non-adiabatic corrections, and coupled-channel dynamics.

The principal numerical findings can be summarized as follows. The scalar eigenvalue lies in the narrow window $ M(T_{bc}^{(0)}) = 7.143 $–7.158 GeV, while the axial-vector eigenvalue lies in $ M(T_{bc}^{(1)}) = 7.217 $–7.222 GeV, with $ \kappa_{\bar b\bar c} $-induced spreads of only $ \sim 15 $ and $ \sim 5 $ MeV, respectively. The overall parametric uncertainty in the absolute masses is of order $ \pm 60 $–70 MeV and is dominated by the light-diquark mass. However, this uncertainty is largely common to the two channels and therefore cancels in mass differences. The associated hyperfine splitting varies linearly with $ \kappa_{\bar{b}\bar{c}} $ in the range $ \Delta_{HF}\simeq 59 $–79 MeV. Most of this splitting $ (\sim 40\; \text{MeV}) $ arises from the kinematic mass gap between the symmetric and antisymmetric heavy-antidiquark configurations, while the slope $ \partial\Delta_{HF}/\partial\kappa_{\bar{b}\bar{c}}=2 $ follows analytically from the $ \langle\vec{S}_{\bar b} \cdot \vec{S}_{\bar c}\rangle $ eigenvalues of the two heavy-antidiquark configurations. Once the heavy-antidiquark mass gap has been independently fixed, this linear response gives the framework genuine predictive power: a future measurement of $ \Delta_{HF} $ would then fix $ \kappa_{\bar{b}\bar{c}} $ in a regime that is currently inaccessible to both lattice QCD and phenomenology. The spatial observables, the mean inter-cluster separation $ \langle r\rangle\simeq 0.45 $–0.46 fm and the inverse mean radius $ \langle 1/r\rangle^{-1}\simeq 0.33 $–0.34 fm, were found to be essentially independent of both $ \kappa_{\bar{b}\bar{c}} $, which enters the radial equation only as an additive constant, and $ J^{P} $. The 40 MeV gap between the symmetric and antisymmetric heavy-antidiquark configurations changes the reduced mass by only $ 5\times 10^{-4} $, so the two channels share the same radial profile to a precision far beyond that at which the radii are quoted. This indicates that the spatial wave function is governed almost entirely by the heavy-antidiquark–light-diquark reduced mass and that the chromomagnetic interaction acts as a small perturbation on an otherwise nearly universal radial profile. The ratio $ \langle 1/r\rangle^{-1}/\langle r\rangle\simeq 0.73 $–0.74 places the $ \Sigma_g^+(1S) $ wave function squarely between the sharply localized and long-tailed limits, indicating that both short-range gluon exchange and the long-range confining string contribute appreciably to the binding. This pattern is familiar from heavy quarkonia and is now seen to extend to genuinely four-quark configurations.

Comparison with the relevant open-flavor thresholds reveals two qualitatively different dynamical situations. The central value of the scalar $ T_{bc}^{(0)} $ mass lies essentially on the $ B\bar D $ threshold, with displacements ranging from $ +9 $ MeV (Set I) to $ -6 $ MeV (Set III). Whether the physical state would ultimately emerge as a genuinely bound, weakly decaying tetraquark—in close analogy with the role played by the $ D^{*+}D^{0} $ threshold for $ T_{cc}^{+}(3875) $—or as a narrow near-threshold structure cannot be determined within the present framework. The answer depends not only on the precise value of $ \kappa_{\bar{b}\bar{c}} $ but also on the parametric uncertainties of the inputs and on the neglected coupled-channel dynamics, whose estimated $ {\cal{O}}(10 $–$ 30) $ MeV effects are comparable to the displacement from threshold. A definitive answer will require a coupled-channel analysis that incorporates the $ B\bar D $ continuum explicitly in order to track the migration of the physical pole between the first and second Riemann sheets. The axial-vector $ T_{bc}^{(1)} $ eigenvalue, in contrast, is consistently found 23–28 MeV above $ B^{*}\bar D $ but $ \sim 70 $ MeV below $ B\bar D^{*} $. This suggests an interpretation as an S-wave resonance candidate in the $ B^{*}\bar D $ channel, with a line shape that would be strongly distorted by the nearby cusp. Here too, the identification of the eigenvalue with a physical pole position is subject to the coupled-channel caveats discussed in Sec. III C. The contrast between the two channels—a quasi-bound scalar versus a near-threshold resonant axial-vector—is not accidental within the model. Rather, it is a direct consequence of the interplay between the $ \sim 70 $ MeV tetraquark hyperfine splitting and the heavy-quark-symmetry-breaking pattern of the surrounding open-flavor spectrum, in which the $ \sim 140 $ MeV $ D^{*} $–D gap substantially exceeds the $ \sim 45 $ MeV $ B^{*} $–B gap.

Completing the $ \Sigma_g^+(1S) $ multiplet, we find that its four isovector members follow from the isoscalar ones by a rigid shift of $ 2\kappa_{ud}\approx 206 $ MeV. This is a structural consequence of the isospin–spin lock rather than of an independent isospin-dependent force. All four isovector members lie above every $ B^{(*)}\bar D^{(*)} $ threshold, so they are expected to be broad and are unlikely to be observed as distinct structures. We have also revisited the restriction to the $ \bar{{\bf{3}}}\otimes{\bf{3}} $ color configuration. The interaction between a $ {\bf{6}}_c $ cluster and a $ \bar{{\bf{6}}}_c $ cluster is attractive and, by Casimir scaling, stronger than the $ \bar{{\bf{3}}}_c $–$ {\bf{3}}_c $ interaction by a factor of $ 5/2 $. Thus, the sextet channel cannot be dismissed on the basis of the intra-cluster repulsion alone. Solving the radial problem in the Casimir-scaled potential, however, places the sextet configuration 0.1–0.3 GeV higher, because the steeper confining rise outweighs the deeper short-distance well at the separations relevant here, and the intra-cluster repulsion raises it further. The resulting second-order admixture shifts the ground state downward by $ {\cal{O}}(10) $ MeV, comparable to the neglected coupled-channel effects and already covered by the quoted systematic uncertainty.

From a broader perspective, the present analysis illustrates how the dynamical diquark framework, supplemented with lattice-QCD BO potentials and a minimal set of independently constrained inputs, can deliver definite and falsifiable central-value estimates for exotic states in the doubly-heavy sector with very modest parametric freedom, provided the systematic limitations described above are kept in view. Several natural extensions merit further investigation. The inclusion of explicit coupled-channel $ B^{(*)}\bar D^{(*)} $ dynamics is essential to convert the near-threshold scalar result into a quantitative statement about binding versus virtuality. The tower of orbitally and radially excited states built on higher BO potentials—most notably $ \Sigma_u^- $ and $ \Pi_u $, which generate hybrid-like multiplets—would extend the spectrum to the 7.5–8.0 GeV region and provide additional tests of the diquark organization. Electromagnetic and weak transition rates, together with production cross sections in $ pp $ collisions and $ B_c $ decays, would supply the dynamical observables needed to discriminate the compact-tetraquark interpretation from molecular or hadrocharmonium alternatives. Finally, an analogous treatment of the $ T_{bb} $ and $ T_{cc} $ sectors, using the same BO potential and the corresponding heavy-antidiquark masses, would test the universality of the underlying picture and clarify the systematic trends predicted by heavy-quark symmetry.

On the experimental side, $ T_{bc} $ states in the 7.1–7.3 GeV window are within reach of the high-luminosity runs of LHCb and the upcoming Belle II data set. A resonant structure observed in the $ B\bar D $ invariant-mass distribution near 7.15 GeV, displaced from threshold by no more than a few MeV, would be a natural candidate for $ T_{bc}^{(0)} $. A structure in $ B^{*}\bar D $ around 7.22 GeV, roughly 25 MeV above threshold, would constitute the expected signature of $ T_{bc}^{(1)} $. A pronounced branching asymmetry in favor of $ B^{*}\bar D $ over $ B\bar D^{*} $ would further reinforce the compact-tetraquark interpretation against molecular alternatives. The discovery—or non-observation—of these states with the mass hierarchy and near-threshold character estimated here would significantly sharpen our understanding of multiquark dynamics in the doubly-heavy sector and, by pinning down $ \kappa_{\bar{b}\bar{c}} $, provide a long-sought empirical anchor for the chromomagnetic interaction between two distinct heavy antiquarks.

This work was supported by the Scientific Research Projects Coordination Unit of Ondokuz Mayis University under project number BAP01-2025-5761.