Heavy-Quark spin symmetry breaking in the Born-Oppenheimer approximation

Heavy-quark spin symmetry is explicitly broken by the mass splitting between a heavy-light pseudoscalar meson and its vector partner. This fact plays a pivotal role in the physics of states whose mass lies close to the threshold of an open-flavor meson pair, like $X(3872)$. We show that this source of heavy-quark spin symmetry breaking can be systematically included within the diabatic representation of the Born-Oppenheimer approximation. We verify that including all the appropriate coupled channels guarantees conservation of total angular momentum, parity, and charge-conjugation parity. This marks a fundamental step towards a unified, first-principles study of quarkonia and meson molecules with hidden flavor.


Introduction
In the Born-Oppenheimer (BO) approximation for QCD, the spectrum of a heavy quarkantiquark system is determined by potentials which are numerically accessible in lattice QCD [1].Because of this feature, the BO picture is one of the most promising tools for the study of conventional and unconventional heavy quark-antiquark states from first principles; see, for instance, refs.[2,3] and references therein.The idea of using lattice QCD calculations of the BO potentials that include the effects of string breaking to precisely determine the spectrum and decays of heavy quark-antiquark systems has been around for quite some time [4].Although direct lattice measurements of string breaking have been available for many years [5,6], achieving such a description remains a current theoretical challenge.
Until recently, the main difficulty consisted in overcoming the single-channel approximation traditionally associated with BO [2], as string breaking couples quark-antiquark and di-meson channels.Significant theoretical progress in the application of BO with coupled channels has been reported.This includes, for instance, calculations of the bottomonium spectrum in BO using the first calculation of string breaking in lattice QCD [7,8], phenomenological studies of quarkoniumlike mesons and open-flavor di-meson scattering in the diabatic representation of BO [9,10], analyses of exotic heavy hadrons in the diabatic representation of the diquark model [11], and calculations of threshold effects in a BO effective field theory [12].
Nevertheless, there is yet another constraint customarily linked with BO that must be removed in order to carry out precise, ab initio studies of heavy quark-antiquark states: heavy-quark spin symmetry (HQSS).Current BO calculations of the bottomonium spectrum with string breaking [7,8] rely on HQSS, since it is realized for lattice QCD simulations with static quarks.However, because of the mass splitting between a heavy-light pseudoscalar meson B and its vector partner B * , HQSS is broken by the mixing of quarkoniumlike resonances with di-meson pairs.In fact, this HQSS-breaking effect has been used to explain the hadronic decays of Υ(4S) and Υ(10860); see, for instance, refs.[13,14] and references therein.In light of this, one could question the reliability of HQSS-based calculations of coupled quark-antiquark and di-meson channels, especially for energies close to a meson-pair threshold.The situation is even more problematic in the charmoniumlike sector, with a D * -D splitting three times bigger than the B * -B splitting.Phenomenological analyses [15] indicate that the precise location of the di-meson thresholds may be decisive for the understanding of exotic states like X(3872).
The problem of HQSS breaking in the BO approximation has been addressed in the case of ud bb tetraquarks by including the B * -B mass splitting and solving the resulting Schrödinger equations for the coupled BO channels [16].However, the ud bb tetraquark case is made somewhat simpler by the absence of any string breaking, hence the formalism described in ref. [16] is not easily extended to other situations where string breaking plays a prominent role.
In this work, we show that HQSS breaking can be systematically and consistently included in the diabatic representation of the BO approximation.The effect of HQSS breaking is to split the thresholds and replace the single string-breaking transition rate from lattice QCD by a matrix of transition rates between quark-antiquark and di-meson configurations, with coefficients determined from Dirac algebra.The quark-antiquark static energy, di-meson thresholds, and the transition rates form an interaction matrix that depends on the spin configuration of the heavy sources.Linear combinations of the interaction matrix elements for different spin configurations, with coefficients determined from angular momentum theory, form the entries of a potential matrix from which the spectrum of the heavy quark-antiquark system can be calculated.As a byproduct, we verify that the cylindrical symmetry group of a static quark-antiquark pair, D ∞h ⊗ CP , gets naturally promoted to the full O(3) ⊗ C symmetry group of QCD when the kinetic energies of the heavy quarks are introduced.
This result constitutes a substantial amendment to the previous phenomenological analyses in refs.[9,10] by determining the matrix of coupling coefficients that in these studies was arbitrarily fixed ad hoc.It provides a foothold for a detailed, unified study of heavy quark-antiquark systems from ab initio lattice QCD simulations with static quarks.
The contents of this paper are organized as follows.In section 2, we thoroughly review the static approximation for a heavy quark-antiquark pair.In section 3, we calculate the transition rates between quark-antiquark and the various di-meson configurations in terms of the string-breaking transition rate.In section 4, we include violation of HQSS by threshold spin splittings.In section 5, we introduce the kinetic energies of the heavy quarks in the diabatic representation of BO and write down the master Schrödinger equation with quark-antiquark and di-meson channels.Finally, in section 6, we summarize our findings and indicate possible future developments.

Static approximation
For the system of a heavy quark-antiquark pair with strong interactions, the time scale for the motion of the heavy quarks, dictated by the heavy-quark mass m Q , is considerably larger than the time scale for the evolution of the gluon and light-quark fields, dictated by the nonperturbative QCD energy scale Λ QCD .Because of this sharp separation of the time scales involved, it makes sense to expand the Hamiltonian of the heavy quark-antiquark system in powers of Λ QCD /m Q .The leading order contribution in this expansion corresponds to the infinite mass limit for the heavy quarks, m Q → ∞.In this limit, the heavy quark-antiquark pair has no motion and acts as a static source for the light QCD fields.The energy levels of the dynamical fields, that is, gluons and light quarks, are not affected by the spins of the static quarks.This decoupling of the heavy-quark spins is referred to as HQSS.
Hence, let us begin by considering QCD with static quark (Q) and antiquark ( Q) sources at the positions + 1 2 r and − 1 2 r.In this static approximation, r is a constant parameter.The dynamical fields are gluons and light quarks (q) and antiquarks (q).The relevant symmetry group of QCD is O(3) ⊗ C, where O(3) is the group formed by three-dimensional rotations and the parity transformation P , and C is the charge conjugation.However, the static quarks break the symmetry down to the cylindrical subgroup D ∞h ⊗ CP , which consists of: • rotations around the r axis; • a reflection R through a plane containing r; • the combined transformation CP .
There are several terms in the total angular momentum J of a heavy quark-antiquark system.We denote the light QCD angular momentum of the dynamical fields by J light .The total spin of the heavy quarks is S Q Q.The orbital angular momentum of the sources is L. The total angular momentum is then It is also convenient to introduce the static spin defined by The representations of D ∞h ⊗ CP are conventionally labeled as Λ ϵ η , where: • Λ = |J light • r| is the modulus of the projection of the light QCD angular momentum J light on r.The values of Λ are indicated by Greek letters: Σ, Π, ∆ for Λ = 0, 1, 2, respectively, and so on.
• ϵ = ± is the eigenvalue of R in the case Λ = 0, and it is omitted for Λ > 0.
The energy levels of light QCD in the presence of the static sources can be labeled by the light BO quantum numbers Λ ϵ η .One could also define static BO quantum numbers for the system that consists of the light QCD fields and the spins of the heavy quarks, with static spin quantum number s and projection λ = S • r.In the case λ = 0, the value of ϵ is completely determined by the static spin s, the intrinsic parity p of the source, and the light BO quantum numbers J light and ϵ light as ϵ = ϵ light p(−1) J light −s (see appendix A), making this quantum number redundant.Thus, we label the configurations by the static spin projection λ with sign and η, omitting ϵ.We shall denote these static BO quantum numbers that also take into account the spin of the heavy quarks by λ η , such as 0 g , 0 u , ±1 g , ±1 u , ±2 g , ±2 u , and so on.
In this paper, we shall concentrate on light QCD configurations with BO quantum numbers Σ + g , which can be produced by a Q Q source and by di-meson sources.For a Q Q source, that is, the product of Q and Q fields connected by a gauge transporter, Σ + g corresponds to J light = 0, so that S = S Q Q.For a di-meson source, that is, a product of Qq and q Q operators, S coincides with the total spin of the di-meson pair.So, for Q Q and di-meson sources producing light Σ + g configurations, s and λ are just total spin quantum numbers.Note, however, that for a Q Q source the individual spin quantum numbers are 12 , while for a di-meson source the individual spin quantum numbers are integers.
In theory, a Q Q source generates all Σ + g light QCD configurations, labeled by n = 1, 2, 3, . . . in increasing energy levels.However, at any given Q Q distance r = |r|, many of these light QCD configurations will have a very small overlap with the quark-antiquark source, making it difficult to detect them in lattice QCD simulations with finite numerical precision.A practical workaround is to consider the correlation between quark-antiquark and di-meson sources.Following the calculations of the Σ + g BO potentials with string breaking in lattice QCD [5], we truncate the operator expansion by considering only di-meson sources for negative-parity mesons, which yields a good description of the ground n = 1 and the first excited n = 2 Σ + g BO configurations. 1 Here we do not consider di-meson sources for positive-parity mesons, which are presumably required for the calculation of excited Σ + g light QCD configurations with larger n.Neither do we include configurations with light BO quantum numbers Π u , Σ − u , etc., also referred to as hybrid configurations.Hence, the current analysis is restricted to energies well below the lowest threshold for the production of a positive-parity meson (B B1 ) and the energy of the lowest quarkonium hybrid meson.
We shall label the quark-antiquark and di-meson sources producing λ η configurations by their total spin s and projection λ.Note that in the Q Q, B B, and B * B * cases the quantum number η is determined by the total spin of the source as CP = (−1) s+1 for quark-antiquark sources and CP = (−1) s for di-meson sources.On the other hand, one can Table 1.Total spin s for quark-antiquark and di-meson sources producing a light QCD configuration Σ + g that combines with the spins of the heavy quarks to give the specified static BO quantum numbers λ η (λ is the total spin projection and η = g, u stands for CP = +, −).Since Q and Q each have two spin states, there are four sets of static BO quantum numbers λ η that correspond to the light BO configuration Σ + g .Namely, the four spin states of a Q Q pair can be grouped in a triplet corresponding total spin s = 1, with λ = −1, 0, +1 and η = g (that is, CP = +), and a singlet corresponding to s = 0, with λ = 0 and η = u (that is, CP = −).Therefore, the relevant sets of static BO quantum numbers λ η are 0 g , +1 g , −1 g , and 0 u .The corresponding values of the total spin quantum number s for Q Q and B ( * ) B( * ) sources are summarized in Table 1.There are also four sets of static BO quantum numbers λ η that can be produced by di-meson sources but not quark-antiquark sources: +2 g , −2 g , +1 u , and −1 u .Although irrelevant for string breaking, these configurations have also been listed in Table 1 as they will become useful later on.
We proceed to describe the sources with separation r that produce light QCD configurations Σ + g that when combined with the spins of the heavy quarks have specified static BO quantum numbers λ η .A quark-antiquark source that creates 0 g configurations can be constructed as [5,6]) where Q(+ 1 2 r) and Q(− 1 2 r) are the heavy quark and antiquark operators and W(− 1 2 r, + 1 2 r) is the equal-time parallel gauge transporter along the straight line between − 1 2 r and + 1 2 r.The B B source that creates 0 g configurations is (see refs.[5,6]) where q − 1 2 r and q + 1 2 r are light quark and antiquark operators.Notice that we consider only one flavor of light quarks, for simplicity.
The lattice-QCD calculation of the BO potentials with string breaking in ref. [5] was carried out using only Q Q and B B sources in a static BO configuration 0 g , that is, Q 1,0 and B 0,0 .Although these operators are sufficient to generate the first two Σ + g energy levels, they do not exhaust all the possible di-meson sources creating a 0 g configuration.Moreover, one needs B B * η and B * B * sources to create the additional ±1 g and 0 u configurations, which also participate in string breaking; see Table 1.Therefore, we proceed to describe a complete set of quark-antiquark and di-meson sources, including B B * η and B * B * , which is also instrumental for incorporating the B * -B mass splitting later on.
For ease of notation, from here on we shall suppress the argument of the quark operators, in the implicit understanding that Q and q are created at + 1 2 r, while Q and q are created at − 1 2 r.We shall also identify the arbitrary axis ẑ with the r direction, which is legitimate in the static limit since r is a constant.This allows us to simplify expressions involving r such as γ • r → γ 3 , for instance.
One can construct a gauge-invariant quark-antiquark operator creating light BO quantum numbers Σ + g by connecting Q and Q sources via a straight Wilson line W and contracting the Q and Q spins with a 4 × 4 Dirac matrix.Since the static quark and antiquark sources satisfy there are only 4 independent Q Q operators in total.There is one spin-0 operator, and three spin-1 operators, Qγ 3 WQ and Qγ ± WQ with spin projection 0 and ±1, where On the other hand, one can construct 8 independent operators with the quark content and quantum numbers of a ground-state heavy-light meson.There is the B operator, with spin projection 0 and ±1, plus the corresponding B and B * operators that can be obtained by substituting q → Q and Q → q.From the product of B and B * operators with B and B * operators, one can construct 16 B ( * ) B( * ) operators: 1 for B B, 6 for B B * η (3 for each η = g, u), and 9 for B * B * .The 9 B * B * operators can be combined in multiplets with definite total spin, that is, a singlet with s = 0, a triplet with s = 1, and a quintuplet with s = 2, using the Clebsch-Gordan coefficients.These 20 operators, 4 for Q Q plus 16 for B ( * ) B( * ) , constitute a complete basis for generating static BO configurations associated to the ground and first excited Σ + g light QCD configurations.Next, we provide explicit expressions for the subset of operators that create λ η configurations participating in string breaking.Following Table 1, the Q Q and B ( * ) B( * ) sources that create 0 g configurations are The Q Q and B ( * ) B( * ) sources that create ±1 g configurations are The Q Q and B ( * ) B( * ) sources that create 0 u configurations are (2.3c)

String breaking and transition rates
In lattice QCD, string breaking is manifested explicitly as a nonvanishing transition rate between states created by quark-antiquark and di-meson sources.Lattice QCD has been used to calculate the transition rate between Σ + g states created by Q Q and B B sources in a 0 g static BO configuration [5].In this section, we express the transition rates between Q Q and B ( * ) B( * ) sources for different λ η in terms of the string-breaking transition rate g between a gluonic string with light BO quantum numbers Σ + g and a light quark-antiquark pair within the same configuration.
Any di-meson operator can be factorized into products of light-quark operators and heavy-quark operators using the Fierz identity, where the overall minus sign comes from the anticommutation of fermionic operators.Specifically, the di-meson sources in eqs.(2.1b)-(2.1e) that create static BO configurations HQSS implies that the Q and Q spins are conserved.The 4 operators Qγ 3 Q, Qγ + Q, Qγ − Q, and Qγ 5 Q in eqs.(3.1) correspond to 4 independent Q Q spin states.In the quarkantiquark source Q 1,0 in eq.(2.1a), the Q Q spin state is given by the operator Qγ 3 Q.Hence, the static correlation function between Q 1,0 and any of the di-meson sources in eqs.(3.1) proceeds exclusively through the terms in the latter that are proportional to Qγ 3 Q, where s,0 is a shorthand for any of the di-meson sources in eqs.(3.1) and T τ is the Euclidean time-evolution operator from time 0 to τ .Factoring out the Q Q spin given by Qγ 3 Q, the correlator on the right side of eq.(3.2) can be cast symbolically as ⟨0|(qP − γ 3 q)T t W † |0⟩.From it, the string-breaking transition rate g can be obtained by evaluating the time derivative with an appropriate normalization coefficient after some relaxation time for which higher excited light QCD configurations have decayed; see ref. [5].
Thus, it follows from eq. (3.2) that the transition rates between Q Q and B ( * ) B( * ) are all proportional to the string-breaking transition rate g.For the 0 g configuration, the coefficients are simply the numerical factors multiplying Qγ 3 Q inside eqs.(3.1).The transition rates then read where we have specified the static BO quantum numbers η and λ of the configuration as superscripts, the total spin of the quark-antiquark and di-meson sources as Q Q(s) and B ( * ) B( * ) (s), respectively, and the distance r between the sources as an argument.One can proceed in the exact same manner for the ±1 g and 0 u configurations, by considering the corresponding quark-antiquark and di-meson sources.For the ±1 g configurations, the di-meson sources in eqs.(2.2b)-(2.2c)become Since the spin states of the heavy quark-antiquark pair in Q 1,±1 in eq.(2.2a) are given by Qγ ± Q, the transition rates are For the 0 u configuration, the di-meson sources in eqs.(2.3b)-(2.3c)become Since the spin state of the heavy quark-antiquark pair in Q 0,0 in eq.(2.3a) is given by Qγ 5 Q, the transition rates are As explained in more detail in the next section, one can use eqs.(3.3)-(3.7) to verify that, for all four λ η configurations, a unitary change of basis reduces the interaction matrix between Q Q and B ( * ) B( * ) down to a 2 × 2 string-breaking interaction sub-matrix G(r) between W and qP − γ 3 q, plus rows and columns associated to decoupled qP − γ ± q and qP − γ 5 q components.Otherwise said, the relevant string-breaking interaction matrix is the same independently of the spins of the heavy quarks, consistently with HQSS.
The 2 × 2 string-breaking interaction matrix reads where V Q (r) is the Q Q static energy relative to twice the static B meson mass.Notice that the zero in the lower right corner of eq.(3.8) corresponds to approximating the B ( * ) B( * ) static energy by its constant value at large r.The constant can be set to zero by taking the zero of the energy to be the static meson pair threshold.This ansatz provides a good fit to the lattice QCD data in [5]; see, for instance, ref. [7].
It is worth mentioning that the results presented here are straightforwardly extended to the case of dynamical quark fields with N f flavors and SU (N f ) flavor symmetry.The only notable difference is that the string-breaking transition rate with N f light flavors is equal to the single-flavor one times a factor N f .If, alternatively, both light and strange quarks are considered [6], then one needs to substitute eq.(3.8) with an appropriate 3 × 3 string-breaking interaction matrix, like

B( * )
s sources is negligible, since it requires the creation of ss which is OZI suppressed.The ansatz in eq.(3.9) with g(r) and g s (r) replaced by constants is identical to the parametrization used in ref. [6] for fitting the lattice QCD interaction matrix with 2 + 1 light-quark flavors.

Threshold spin splittings
To go beyond the static approximation, one has to introduce corrections in increasing powers of 1/m Q that typically depend on the heavy-quark momenta and/or spins.We shall limit ourselves to the first-order corrections, that is, the nonrelativistic kinetic energies of the heavy quarks and heavy-quark spin effects of first order in 1/m Q .In this section we shall focus on the latter, leaving the treatment of the kinetic energy for later.
At first order in 1/m Q , the heavy-quark spin corrections to the static approximation consist of couplings between the angular momentum of the light QCD fields and the spins of the heavy quarks [17].For the Q Q configurations with light BO quantum numbers Σ + g (J light = 0), these 1/m Q couplings vanish.For the B ( * ) B( * ) configurations, one can apply heavy-quark effective theory to calculate these 1/m Q effects for each individual heavy-light meson, which are responsible for the B * -B mass splitting [18][19][20][21].As for the Q Q-B ( * ) B( * ) transition rates, we shall assume that they are reasonably approximated by their HQSS expressions in eqs.(3.3), (3.5), and (3.7) at first order in 1/m Q .Therefore, for the current study, taking into account heavy-quark spin corrections of first order in 1/m Q boils down to splitting the B B, B B * η , and B * B * threshold masses.
Threshold spin splittings are easily included in the interaction matrix between Q Q and B ( * ) B( * ) in each λ η configuration.Let ∆ = m B * − m B be the mass difference between B * and B mesons and m B ( * ) = (m B + 3m B * )/4 the spin-average mass.Taking, for instance, the 0 g configuration, the corrected interaction matrix between Q 1,0 , B 0,0 , B * g 1,0 , B * * 0,0 , and where v is the unit row vector of coefficients from eqs. (3.3), and is the 4 × 4 threshold mass matrix relative to the spin-average B ( * ) B( * ) threshold, Notice that, in the HQSS limit ∆ → 0, the 5 × 5 matrix G g,0 (r) in eq. ( 4.1) can be reduced to the 2 × 2 matrix G(r) in eq.(3.8) plus vanishing rows and columns by a unitary change of basis for the 4 di-meson channels that transforms v into (1, 0, 0, 0).However, when ∆ > 0 the threshold mass matrix M is also affected by this transformation, such that one cannot have both a diagonal M and v = (1, 0, 0, 0) at the same time.This shows that string breaking with threshold spin splittings breaks HQSS.
This HQSS-breaking effect can also be pointed out by calculating the interaction matrices in the ±1 g and 0 u configurations.For the ±1 g configurations, the interaction matrices between Q 1,±1 , B * g 1,±1 , and B * * 2,±1 are 3 × 3 matrices with the off-diagonal entries from eqs. (3.5), For the 0 u configuration, the interaction matrix between Q 0,0 , B * u 1,0 , and B * * 1,0 is a 3 × 3 matrix with the off-diagonal entries from eqs. (3.7), Notice that, in the HQSS limit ∆ → 0, the 3 × 3 interaction matrices G g,±1 (r) and G u,0 (r) in eqs.(4.2) and (4.3) can be reduced to the same 2 × 2 matrix G(r) in eq.(3.8) plus vanishing rows and columns by a unitary change of basis for the 2 di-meson channels.
Note that B B * u and B * B * can also produce λ η configurations that have no overlap with those created by Q Q; see Table 1.Such configurations yield trivial "interaction matrices" which are just the threshold mass matrices for the corresponding di-meson channels.As will be shown later, a formal BO approximation must include the interaction matrices for all the λ η configurations that can be produced by the sources involved.Hence, for the sake of completeness, we list the interaction matrices for the ±2 g and ±1 u configurations as well.For ±2 g , one has the 1 × 1 matrices 5 Diabatic Born-Oppenheimer approximation

Introducing the kinetic energy
The other first-order correction in 1/m Q to the static approximation is the nonrelativistic kinetic energy for the relative motion of the heavy quark-antiquark pair, The introduction of the kinetic energy operator gives motion to the heavy sources, which can then be determined by solving a coupled-channel Schrödinger equation.This is the core idea of the BO approximation.In this section, we shall detail the formal derivation of this master Schrödinger equation in the the diabatic representation of BO, which allows one to study heavy quark-antiquark systems using the interaction matrix from ab initio lattice QCD.
As we shall see next, one cannot naively plug the lattice-QCD interaction matrix in eq.(3.8) or one of its corrected forms in eqs.(4.1)-( 4.3) as a potential inside a coupledchannel Schrödinger equation.Instead, what governs the motion of the heavy quarks is the diabatic potential matrix, which is constructed from a combination of the interaction matrices corresponding to different λ η configurations.
Formally, the interaction matrix from lattice QCD can be embedded in a quantummechanical framework using the following procedure.Let us indicate the interaction matrix with quark-antiquark distance r within the static BO configuration λ η as G η,λ (r).Let V η,λ n (r) be the nth eigenvalue of G η,λ (r), that is, the nth static energy level with static BO quantum numbers λ η , and |ζ η,λ n (r)⟩ the corresponding static eigenstate.One can then construct a first-quantized operator, acting as a Hamiltonian for the light QCD fields in presence of a heavy quark-antiquark pair at relative position r.From it, the BO Hamiltonian is constructed as so that a quantum state |Ψ⟩ of the heavy quark-antiquark system corresponds to a solution of the eigenvalue equation with E the energy relative to the spin-average B ( * ) B( * ) threshold.Notice that, whilst the static energy levels V η,λ n (r) depend on the distance r = |r| alone, the static eigenstates |ζ η,λ n (r)⟩ depend on the direction r = r/r as well.Most importantly, λ itself depends on r through its definition λ = S • r.In this case one cannot simply identify the fixed axis ẑ with r, since r is a dynamical variable.
Translating eq. ( 5.3) into a differential equation in r requires an appropriate expansion of the heavy quark-antiquark state |Ψ⟩ in terms of some basis for the light QCD fields.One natural choice is the basis formed by the static eigenstates |ζ η,λ n (r)⟩, called adiabatic basis, since it diagonalizes the light QCD energy operator in eq.(5.2).Another convenient choice is a basis that diagonalizes the kinetic energy operator in eq.(5.1), called diabatic basis [22,23].In general, it is possible to express each element of some diabatic basis as a linear combination of the static eigenstates, where the coefficients A λ,σ n,i (r, r 0 ) form a unitary matrix A(r, r 0 ) which can be obtained by solving a first-order differential equation with boundary conditions at r = r 0 ; see, for instance, refs.[9,24].A conventional choice for the boundary conditions is to match the diabatic basis to the adiabatic one at r = r 0 , which is the same as imposing A(r = r 0 , r 0 ) = I with I for the identity matrix.It is important to realize that for r = r0 the cylindrical symmetry of BO implies where r 0 = |r 0 | and A n,i (r, r 0 ) are the coefficients of a radial unitary matrix A(r, r 0 ) which can be obtained by solving a first-order differential equation in r with boundary conditions A(r = r 0 , r 0 ) = I.Notice that the right side of eq.(5.5) does not depend on r0 .Different choices of the matching point r 0 with the adiabatic basis yield different but equivalent choices of diabatic basis.For the present study, it is most convenient to fix r 0 = 0.The reason is that the diabatic basis in this case is matched to the light QCD eigenstates calculated in the static limit where the two sources are placed at the same point r = 0. Unlike the general case, where a nonvanishing separation vector r sets a preferred direction in space, the static limit with r = 0 is symmetric under O(3) ⊗ C. So, the static eigenstates at r = 0 can be classified by the static spin quantum number s and by the J P C light quantum numbers of the light QCD fields, in addition to the static BO quantum numbers λ η .
As mentioned in section 2, we restrict our analysis to the ground and first-excited Σ + In fact, this choice ensures that the mixing problem is clearly organized in terms of channels with well-defined particle number and spin.Hence, we work with a truncated diabatic basis which consists of 4 Q Q and 16 B ( * ) B( * ) channels labeled by their total spin s and the diabatic BO quantum numbers σ η .
This truncated diabatic basis can be formally written as a set {|ζ η,σ i (r, 0)⟩} labeled by σ η and the additional label i that specifies the quark-antiquark or di-meson channels with total spin s as Q Q(s) or B ( * ) B( * ) (s).Following Table 1, there are five channels in the 0 g configuration, three channels in each of the +1 g and −1 g configurations, and three channels in the 0 u configuration, On top of these, one should take into account the channels in σ η configurations that can be produced by di-meson sources only.Each of the +2 g and −2 g configurations has only one channel, and each of the +1 u and −1 u configurations has two channels, with the expansion coefficients Ψ η,σ i (r) acting as wave function components for the various diabatic channels.The eigenvalue equation for the BO Hamiltonian, eq. ( 5.3), translates into a coupled-channel Schrödinger equation, where η is conserved and are the elements of the diabatic potential matrix V η (r).The dimension of the matrix V η (r) can be determined through eqs.(5.6)-(5.10)by counting the total number of channels from all the σ η configurations sharing the same quantum number η.Thus, the diabatic potential matrix V g (r) within CP = + is a 13 × 13 matrix and the diabatic potential matrix V u (r) within CP = − is a 7 × 7 matrix.

Diabatic potential matrix
Next, we show that the diabatic potential matrix V η (r) is completely determined by a combination of the interaction matrices G η,λ (r) with various λ.Let us begin by observing that a spin state with fixed projection onto a general direction r can be expressed as a Wigner rotation of the corresponding spin state with fixed projection along ẑ [25].Thus, one can expand the spin state |s, S • r = λ⟩ in terms of the canonical spin states |s, σ⟩ as with D s σ,λ the Wigner D-matrix element and (φ, θ, ψ) the Euler angles, where θ and φ are identified with the polar angles of r while ψ is left arbitrary.Since the static eigenstates at r = 0 have definite spin s and projection σ along an arbitrary axis ẑ, we can use eqs.(5.5) and (5.14) to express eq.( 5.4) for r 0 = 0 as where s i is the spin of channel i. Plugging eq. ( 5.15) into eq.(5.13), the diabatic potential matrix elements can be expressed as the elements of the interaction matrices G η,λ (r) given in eqs.(4.1)-(4.5).Notice that eq.(5.16) cannot be simply expressed as a matrix equation relating the (σ, σ ′ ) block of the diabatic potential matrix, V η,σ,σ ′ (r), and the interaction matrices G η,λ (r), because the spin quantum numbers s i and s i ′ depend on i and i ′ .The diabatic potential matrix V η (r), which determines the motion of the heavy quarks through the master Schrödinger eq.(5.12), is built by a precise combination of the interaction matrices G η,λ (r) for different values of λ.Note that the block σ = σ ′ = λ of the diabatic potential matrix coincides with G η,λ (r) in the special case r = ẑ: (5.17) It is instructive to verify that the introduction of the heavy-quark motion automatically reestablishes the full O(3)⊗C symmetry of QCD that was broken in the static approximation.As illustrated in appendix B, one can expand the heavy quark-antiquark state in eq.(5.11) in partial waves and observe that the diabatic potential matrix V η (r) only connects partial waves with the same values of the total angular momentum quantum numbers J and M .The Schrödinger eq.(5.12) for wave functions that depend on the vector r gets transformed into a system of coupled radial Schrödinger equations, where u η,J i,l (r) is the reduced radial wave function for the partial-wave channel i with orbital momentum l, spin s i , CP -parity η, and total angular momentum J, and V η,J i,i ′ ,l,l ′ (r) is the partial-wave potential with j 1 j 2 j 3 m 1 m 2 m 3 the Wigner 3-j symbols.Notice that eq.(5.19) does not depend on the total angular momentum projection M , as required by the Wigner-Eckart theorem, so that the same Schrödinger eq.(5.18) (and therefore the same wave function) applies to all the 2J +1 values of M .This suffices to show that the BO Hamiltonian respects the full rotational symmetry group SO(3).Moreover, in appendix C we show that the symmetry under a reflection R of H static (r) ensures that the partial-wave potentials in eq.(5.19) conserve parity.This, combined with symmetry under the combined transformation CP , implies that C-parity is also conserved.We have thus verified that the relevant symmetry group of the BO Hamiltonian is O(3) ⊗ C, so that the heavy quark-antiquark states calculated from it are naturally organized in J P C families.
In practice, plugging the interaction matrices in eqs.(4.1)-(4.5)into eq.( 5.19) yields simple expressions for the partial-wave potentials in terms of the Q Q static energy V Q (r), the B * -B mass difference ∆, and the string-breaking transition rate g(r).Specifically, if one neglects heavy-quark spin corrections of order 1/m 2 Q and above, the Q Q diagonal element of the interaction matrices G η,λ (r) is the same regardless of the value of either η or λ, Combining this with the orthogonality relations of the Wigner 3-j symbols, the sum on the right side of eq.(5.19) gives a simple result for the Q Q potentials, Similarly, if one neglects light-meson exchange contributions and heavy-quark spin corrections of order 1/m 2 Q and above, the di-meson potentials are given by the corresponding threshold mass values relative to the spin-average B ( * ) B( * ) threshold.One has for B B * , and for B * B * .Partial-wave potentials coupling different di-meson channels are zero.Finally, under the reasonable assumption that HQSS remains a good approximation for the transition rates even when threshold spin splittings are included, any partial-wave potential coupling quark-antiquark and di-meson channels is given by a fraction of the string-breaking transition rate g(r), V η,J Q Q(s),B ( * ) B( * ) (s ′ ),l,l ′ (r) = g η,J,l,l ′ ,s,s ′ B ( * ) B( * ) g(r) (5.21) with g η,J,l,l ′ ,s,s ′ B ( * ) B( * ) a scalar coefficient that can be calculated by inserting the interaction matrix elements from eqs. (3.3)-(3.7)into eq.(5.19).For the sake of completeness, we list in appendix D all the nonvanishing values of g η,J,l,l ′ ,s,s ′ B ( * ) B( * ) with J = 0, 1, 2. It is worth emphasizing that, taking the coefficients g η,J,l,l ′ ,s,s ′ B ( * ) B( * ) from appendix D, the coupling potentials in eq. ( 5.21) for S-wave (l = 0) Q Q coupling to P -wave (l ′ = 1) B B, B B * g , B * B * (0), and B * B * (2) in J P C = 1 −− have relative ratios of respectively.These ratios reproduce those in ref. [13] for the amplitudes for production of B B, B * B + c.c., (B * B * ) S=0 , and (B * B * ) S=2 in e + e − annihilation in the limit of exact HQSS (see eq. ( 5) of ref. [13]), up to sign differences due to phase conventions.The coefficients g η,J,l,l ′ ,s,s ′ B ( * ) B( * ) , fully constrained by HQSS and cylindrical symmetry, are always smaller than 1 and differ considerably between one another; see appendix D. This finding is in stark contrast with the naive assumption g η,J,l,l ′ ,s,s ′ B ( * ) B( * ) = 1 for all quantum numbers made in the previous phenomenological analyses in refs.[9,10,15].These analyses, which also use the diabatic representation of BO, did not incorporate the correct symmetries of the static approximation.

Practical applications
As a last remark, we illustrate how the diabatic BO approximation can be used in practical calculations of the quarkoniumlike spectrum.
Since J, P , and C are all conserved quantum numbers, one can concentrate on a specific J P C configuration and consider only the corresponding orbital angular momentum channels when calculating the spectrum.From eq. (5.18), one sees that the BO Hamiltonian in a particular J P C configuration can be represented as a matrix H J P C acting on the reduced radial wave functions.It can be expressed as where the entries of the radial potential matrix V J P C (r) have the form (5.22) Note that η, which determines the relevant quark-antiquark and di-meson channels, corresponds to the product of the given P and C, and that the values of l and l ′ are constrained by conservation of total angular momentum and parity.The eigenvalues and eigenvectors of H J P C can be calculated numerically.From them, the masses of quarkoniumlike bound the 5 × 5 di-meson potential sub-matrix, and the 2 × 5 mixing potential sub-matrix.
In the limit ∆ → 0, there is a unitary transformation that turns the 7 × 7 potential matrix in eqs.(5.25) The last two channels are just trivial free-wave channels decoupled from the others.The other five channels form the 5 × 5 potential matrix with HQSS, with V 0 −+ HQSS (r) the 2 × 2 potential matrix in eq. ( 5.24) and the 3 × 3 potential matrix for J P C = 2 −+ with HQSS.It is clear from eqs. (5.24), (5.26), and (5.27) that, in the limit ∆ → 0, the Q Q S-wave (l = 0) and D-wave (l = 2) channels in the J P C = 1 −− configuration are completely decoupled from each other.This also shows that, in this HQSS limit, the spectrum of quarkoniumlike states with Q Q in S-wave is the same regardless of whether the Q Q spin is 0 or 1 (J P C = 0 −+ or 1 −− ), as expected.
These results should be compared with the BO calculation using lattice-QCD BO potentials with string breaking in ref. [7], where the spectrum of Υ and η b states is also calculated from a coupled-channel radial Schrödinger equation for quark-antiquark and di-meson channels.One can notice that the radial potential matrix of this reference, which may be easily read from the left side of eq.(50) in ref. [7], is identical to our V 0 −+ HQSS (r) in eq. ( 5.24) by substituting m Q → 2µ M , V Q (r) → V QQ (r), g(r) → V mix (r), and noticing that V M M,∥ (r) in eq.(50) of ref. [7] is actually 0 (see eq. (59) in ref. [7]).
Therefore, our results reproduce the Schrödinger equation used in ref. [7] in the HQSS limit ∆ → 0. When threshold spin splittings ∆ > 0 are included, however, we observe breaking of HQSS through the difference between the radial potential matrices in the cases of Q Q spin 0 and 1, V 0 −+ (r) and V 1 −− (r).This HQSS-breaking effect is straightforwardly incorporated into the lattice-QCD potentials with string breaking using our eqs.(5.23) and (5.25) for the specific 0 −+ and 1 −− cases, or eqs.(5.19) and (5.22) for any J P C in general.The same cannot be said of the results in refs.[7,8], for which no clear prescription to go beyond the HQSS approximation is provided.
Notice that, since the heavy-quark kinetic energy and the B * -B mass splitting are both first-order corrections in 1/m Q to the static approximation, there is no formal justification for including the former without the latter in a BO approximation.From a phenomenological point of view, this HQSS-breaking effect is kinematically suppressed for quarkonium states with masses well below the B B threshold, where di-meson channels hardly play any role, but it is relevant in all other cases.This is especially important for exotic candidates, whose masses often lie close to some open-flavor meson-pair threshold.
For these reasons, the diabatic BO approximation studied here should be regarded as a necessary improvement of the scheme used in refs.[7,8], which allows one to perform detailed studies for energies close and above threshold.

Overview
We have revisited the static approximation for a heavy quark-antiquark pair.In addition to the Q Q and B B sources, we have also included B B * and B * B * sources.We have shown that the transition rates between quark-antiquark and di-meson configurations are given by fractions of the string-breaking transition rate, with coefficients completely determined from Dirac algebra.
We have then introduced the mass splitting between a heavy-light pseudoscalar meson B and its vector partner B * .We have illustrated the pattern of HQSS breaking deriving from string breaking with threshold spin splittings.We have further shown that this symmetry breaking is easily accounted for in the interaction matrix between quark-antiquark and di-meson configurations.
Finally, we have built a systematic treatment of the BO approximation at first order in 1/m Q in the diabatic representation.We have shown that the dynamics of the heavy quark-antiquark system, including the heavy quark-antiquark interaction, string breaking, and HQSS breaking, is governed by the diabatic potential matrix, which is completely determined from an interaction matrix accessible in lattice QCD.We have also shown, without any ad hoc assumption, that this construction of the BO approximation leads to a natural restoration of the O(3) ⊗ C symmetry of QCD from the cylindrical symmetry D ∞h ⊗ CP of the static approximation.
This result provides an essential stepping stone towards precise, first-principles studies of quarkoniumlike mesons and open-flavor di-meson scattering with string breaking.It is worth highlighting that the inclusion of threshold spin splittings is particularly relevant to the understanding of any heavy quark-antiquark state in which meson pairs play a significant role, and to the understanding of near-threshold states like X(3872) in particular.
For the charmoniumlike sector, it may be necessary to consider even the small D + -D 0 and D * + -D * 0 mass differences.For instance, they may be instrumental in explaining the isospin-violating branching fractions of X(3872) [26,27].Note that such isospin mass splittings can be easily incorporated into the diabatic BO approximation; see ref. [15].
For practical applications in the future, it may be necessary to include also other relevant effects, like heavy-quark spin corrections to the quark-antiquark potential of order 1/m 2 Q and light-meson exchange interactions between two open-flavor mesons.The former are well known [28][29][30] and can be simply added to the Q Q partial-wave potentials.As for the latter, we strongly encourage a theoretical effort for their incorporation in the diabatic BO approximation.
Therefore, a static BO state with static spin projection λ = 0 is an eigenstate of R with eigenvalue ϵ.

B Conservation of total angular momentum
To verify that total angular momentum is conserved when the orbital angular momentum of the sources is reintroduced, let us expand eq.(5.11) in partial waves as |Ψ⟩ = η,J,M,i,l dr u η,J,M i,l (r) |r⟩ |ζ η,J,M i,l (r, 0)⟩ with |ζ η,J,M i,l (r, 0)⟩ the partial-wave channel i with orbital momentum l, spin s i , CP -parity η, total angular momentum J, and projection M , and u η,J,M i,l (r) the corresponding reduced radial wave function.

2 1 1 construct
two independent combinations of B and B * , (B B * ∓ B * B)/ √ 2 with CP = ±.In this last case, we shall specify the corresponding quantum number η of the meson pair as B B * η whenever necessary.Hence, we shall indicate the Q Q sources with s = 0, 1 as Q s,λ , the B B source with s = 0 as B 0,0 , the B B * η sources with s = 1 as B * η 1,λ , and the B * B * sources with s = 0, 1, 2 as B * * s,λ .We will sometimes use the shorthand B ( * ) B( * ) to indicate all the various meson pairs, B B, B B * η , and B * B * .

. 9 )
where g s (r) is the string-breaking transition rate by strange quarks and δ s = m Bs − m B is the mass difference between static B s and B mesons.Notice that we have assumed that the B ( * ) B( * ) and B ( * ) sB( * )s static energies are well-approximated by their constant values at large r.We have also assumed that the transition rate between B ( * ) B( * ) and B ( * ) s

g
light BO configurations generated by a Q Q source and by di-meson sources for negativeparity mesons.The Q Q operators for r = 0 generate J P C light = 0 ++ light ; see eqs.(2.1a), (2.2a), and (2.3a).On the other hand, the B ( * ) B( * ) operators for r = 0 generate J P C light = 1 −− light and 0 −+ light ; see eqs.(3.1), (3.4), and (3.6).The weaker cylindrical symmetry for r ̸ = 0 allows the 0 ++ light component of the Q Q configurations to mix with the 1 −− light (Λ = 0) component of the B ( * ) B( * ) configurations through string breaking; see section 3.However, such mixing is forbidden by the stronger O(3) ⊗ C symmetry for r = 0. 2 Paradoxically, the absence of mixing for r = 0 is what makes r 0 = 0 the ideal choice of diabatic basis for studying string breaking.