Heavy meson thresholds in Born-Oppenheimer Effective field theory

We consider heavy meson-antimeson pairs and their coupling to quarkonium in the context of nonrelativistic EFTs incorporating the adiabatic expansion. We work out all the leading order couplings of quarkonium to heavy meson-antimeson pairs and obtain their contributions to the masses and widths of quarkonia. We match the new potentials terms to NRQCD. Using the available lattice data for the coupled system of quarkonium and the lowest laying heavy meson-antimeson pair, we extract the mixing potential and use it to compute numerically the contributions of $D\bar{D}(B\bar{B})$ and $D_s\bar{D}_s(B_s\bar{B}_s)$ to the masses and widths of the charmonium (bottomonium) states for $l=0,1,2$ and up to $n=6$ covering the states in threshold region. When a quarkonium state and a heavy meson-antimeson pair are separated by small energy gaps, their interactions can be described by a threshold EFT with contact interactions. We work out the matching between the two EFTs obtaining the couplings of the threshold EFT in terms of the mixing potential and quarkonium wave functions.


I. INTRODUCTION
The discovery during the past two decades of several dozen exotic hadrons, understood as those that cannot be classified as mesons or baryons in the quark model picture, has made evident an important gap in our understanding of the QCD spectrum and therefore of its underlying dynamics.Many of these states have been discovered in the double-heavy sector in experiments at B-factories (BaBar, Belle, and CLEO), τ -charm facilities (CLEO-c and BESIII) and hadron colliders (CDF, D0, LHCb, ATLAS, and CMS), see Ref. [1] for a review on the experimental status.Since the creation of heavy quarks pairs in hadrons is highly suppressed due to their mass being much larger than Λ QCD , the number of heavy quarks in a hadron can be identified by the hadron mass while the light-quark and gluonic content can be identified from other quantum numbers.Therefore, the identification of an exotic state is more straightforward in the heavy quark sectors.In the charmonium and bottomonium sectors these exotic states are commonly labeled as "XYZ" and appear in the mass region of the heavy meson-antimeson pairs also known as heavy meson thresholds.
A crucial tool in understanding the spectrum of double-heavy hadrons is the adiabatic expansion between the dynamics of the heavy quarks and that of the light-degrees of freedom, either light-quarks or gluons.In this picture the double-heavy hadrons are the heavy quark bound states supported by the spectrum of static energies (also known as adiabatic surfaces) associated to the light degrees of freedom.Therefore, the first step to elucidate the quarkonium spectrum in the threshold region is to determine the spectrum of static energies.These have been studied on the lattice and the emerging picture, which we sketch in Fig. 1, is as follows.The ground state corresponds to the standard quarkonium potential, while the first excited static state corresponds to a heavy meson-antimeson pair [2,3].At higher energies additional static energies corresponding to pairs of heavier heavy mesons should also appear.The static energies of the heavy meson pairs appear as horizontal lines at the energy corresponding to the heavy meson-antimeson pair mass with some possible attractive or repulsive behavior for short distances [2].Beyond the quarkonium sector, the static energies of heavy meson-antimeson pairs have also been studied in the lattice for I = 1 in Refs.[2,4,5] and for heavy meson-meson in Refs.[6][7][8][9][10].Excited heavy-quark-antiquark static energies, corresponding to hybrid quarkonium states also appear in the region above the first heavy meson threshold [11][12][13][14].These are repulsive in the short-distance due to the heavy quarks being in a color octet state but in the long-distance become a linear, confining potential, corresponding to the string excitations of the standard quarkonium potential.The spectrum of hybrid quarkonium static energies is formed by multiplets of static energies corresponding to different gluonic states, refereed as gluelumps [15].
The standard and hybrid quarkonium spectra above the first heavy meson threshold contains enough states to account for all the observed neutral heavy exotic states observed so far [16], however the specific assignation of states to the observed states is not yet clear.One important step to clarify the exotic spectrum is to incorporate the mixing between the different static states.The mixing of standard quarkonium with the lowest lying hybrid static energies was considered in Ref. [16].Due to the 1 +− quantum numbers of the associated gluelump, the mixing is through a . Sketch of the spectrum of static energies for a heavy quark-antiquark pair.The static energies corresponding to standard and hybrid quarkonium states are labeled by their D ∞h representation and in the latter case by the quantum numbers of the gluelump in parenthesis.The shapes are obtained from a fit to the lattice data of Ref. [12].The heavy meson-antimeson static energies are drawn as constant lines at the energies given by the spin and isospin averages of the heavy meson and antimeson minus the heavy quark mass corresponding to the sum of the Λ parameter (defined in Eq. ( 4)) for each heavy meson.The three energy levels for the heavy meson-antimeson pairs without and with closed strangeness correspond to the three blocks of static states of heavy meson pairs in Table I.
heavy-quark spin dependent operator which is 1/m Q suppressed, with m Q being the heavy quark mass.However, this is not the case for the next set of hybrid static energies, associated to the 1 −− gluelump, which can mix at leading order with standard quarkonium.We will focus on the mixing of standard quarkonium with heavy meson pairs.As we will discuss, standard quarkonium couples at leading order with the lowest lying heavy meson pairs.Therefore, this coupling is the most important effect of the threshold region degrees of freedom on the quarkonium masses up to the energies where 1 −− gluelump hybrid states start to be relevant.This problem has been studied in Refs.[17][18][19][20][21] using models based on an extension of the Born-Oppenheimer approximation that includes the mixing potentials between the heavy meson pair and quarkonium, which is also refereed as the diabatic approach.In Refs.[17][18][19] the potentials where extracted from the lattice data of Ref. [2] and the coupled channel scattering problem was solved numerically.From the poles of the heavy meson t matrices for specific angular momenta, the bottomonium spectrum was identified.In the case of in Refs.[20,21], a model for the quarkonium heavy meson pair mixing was used.To obtain the spectra a mix approach was used in which the contributions to a particular quarkonium state of the above-lying thresholds are obtained solving the coupled channel Schrödinger equations while the below-lying ones are obtained in perturbation theory.
In this paper we examine the coupling of quarkonium to the heavy meson pairs in the context of an effective field theory (EFT) incorporating the heavy quark mass and adiabatic expansions.The EFT incorporating these two expansions for quarkonium is known as strongly coupled potential nonrelativistic QCD (pNRQCD) [22,23] while its extension to nontrivial light degrees of freedom content has been called Born-Oppenheimer EFT (BOEFT) [24][25][26].The content of this paper can be considered an extension of either EFT, however for easy reference we will consider it as part of the latter.At leading order in the heavy-quark mass expansion the heavy mesons are characterized by the spin, parity and flavor of the light-quark state [27].Combining these for a heavy meson-antimeson pair one arrives to the total spin, parity and charge conjugation of the light-quarks which characterize the heavy meson-antimeson state.We derive all the leading order couplings of heavy-meson-antimeson states to quarkonium.Using these, we compute the contribution of the heavy meson thresholds to the quarkonium self-energy from which we obtain the contributions to the quarkonium masses and decay widths.We obtain the matching expression of the mixing potential in NRQCD considering both the cases when the mixing can be considered a perturbation and not.We compute numerically the quarkonium spectrum up to O(1/m Q ) for S, P and D waves and up to the principal quantum number n = 6 which covers the mass range for which exotic quarkonium states have been discovered.To do so, we use a parametrization of the potentials that combines lattice data and fix order computations.For these states we compute the contributions of the lowest lying heavy meson pairs without and with closed strangeness, using the mixing potential from the lattice data of Ref. [2].Finally, we work out the matching of BOEFT with a threshold EFT containing a quarkonium state, a heavy meson and a heavy antimeson as effective degrees of freedom with contact interactions, as used for instance in Refs.[28][29][30][31][32][33][34][35][36].We discuss what expansion is involved and under what conditions it can be implemented.
We organize the paper as follows.In Sec.II we show how to incorporate heavy meson pairs into BOEFT and compute their contribution to the masses and widths of quarkonium states in perturbation theory.The matching of the new potentials to NRQCD is discussed in Sec.III.In Sec.IV we extract the form of the static potentials for standard quarkonium and the mixing potential with the lowest lying heavy meson pairs from lattice QCD together with inputs from perturbation theory.Using these potentials, in Sec.IV E we compute the contributions of the lowest lying heavy meson thresholds to the quarkonium masses and decay widths.The matching of BOEFT to the threshold EFT is examined in Sec.V. Finally, we provide our conclusions on Sec.VI.

II. HEAVY MESON PAIRS IN BOEFT
To construct the Lagrangian for BOEFT one first identifies the quantum numbers of the light degrees of freedom, chiefly the spin κ, parity p and charge conjugation c, that generate the set of static energies we are interested in.For the case of heavy meson-antimeson pairs this means identifying the light-quark states.From the heavy meson heavyquark-spin multiplets we can identify the quantum numbers of the light-quark states.The ground state corresponds to κ p = (1/2) + and is followed by two states of similar mass with κ p = (1/2) − and κ p = (3/2) − .Each heavy mesonantimeson pair is characterized by the spin and parity of the light-quark states forming the two heavy mesons, which we label as κ p1 1 and κ p2 2 .Combining the spin and parity of the light quark and antiquark we arrive to the allowed total spin, parity and charge conjugation, κ pc , of the light-quarks in Table I.Each of these combinations is represented as a field in BOEFT.We will denote these fields as M κ pc (t, r, R), with r and R being the relative coordinate and center of mass of the heavy quarks.One should keep in mind that the fields M κ pc have spin indices corresponding to the light-quarks (in the κ representation) and corresponding to the heavy-quarks (in the (1/2) * ⊗ (1/2) representation).The field M κ pc carries the light-quark flavor quantum numbers of the heavy meson-antimeson it corresponds, this can be the individual light-quark flavors, isospin or chiral symmetry representations.Note that in the latter cases, the field M κ pc would correspond to a sum of heavy meson-antimeson pairs.The light-quark flavor quantum numbers do not affect the construction of the Lagrangian, and to simplify the notation we will not track them in this section.However, it should be keep in mind that the parameters and potentials of the heavy meson-antimeson pair field do depend on the light-quark flavor.The quarkonium field will be denoted as Ψ(t, r, R).The bilinear terms for both fields can be read off Ref. [26], with quarkonium corresponding to the κ pc = 0 ++ case: The expansion of the Hamiltonian densities h Ψ and h κ p up to 1/m Q is as follows with p = −i∇ r and P = −i∇ R .The symmetry group of two static heavy quarks is D ∞h , which is a cylindrical symmetry group for rotations along the r axis.The representations of D ∞h are labeled as Λ σ η , where Λ is the absolute value of the projection into the heavy quark-antiquark axis of the spin κ of the light degrees of freedom and is labeled by capital Greek letters: Σ, Π, ∆, ... corresponding to Λ = 0, 1, 2, ... .η is the CP eigenvalue, denoted by g = +1 and u = −1.Finally, for Λ = 0, there is a reflection symmetry with respect to a plane passing through the r axis.The eigenvalues of the corresponding symmetry operator being labeled as σ = ±1.The potential terms in Eq. ( 2) should be expanded in representations of D ∞h .For quarkonium there is only one possible projection, into the Σ + g representation, and therefore the corresponding projector is just an identity.Furthermore, both the leading order and next-to-leading order potentials are heavy-quark spin independent and therefore consists of a single potential term.
For the heavy meson pairs all the possible projections for each κ pc state are listed in the third column in Table I.The static potential between a heavy meson-antimeson pair has been studied on the lattice in Refs.[2,[4][5][6][7][8][9][10]. The results show that the potentials are mostly flat lines at the energy corresponding to the heavy meson masses except in some cases in the short-distance limit where the potential is attractive or repulsive depending on the specific heavy meson pair.The hidden heavy flavor isospin I = 0 case has been studied in Refs.[2,3].As we will discuss in detail in Sec.IV B, once the mixing with quarkonium is taken into account, the heavy meson-antimeson static potential is completely flat for the range of data in Ref. [2].Therefore, for this work we will assume that the interaction between the heavy mesons is negligible.Thus, we set the heavy meson-antimeson static potential to be the sum of the heavy meson masses at leading order in in the heavy quark mass expansion minus the origin of energies, which is set at the heavy quark masses with 1 κ an identity in the light-quark spin-space, and Λκ p is related to the heavy meson masses [27] as Notice that, due to this identity all the projections into D ∞h representations of the field M κ pc have degenerate static potentials.We assume that the subleading potential V κ pc corresponds to the sum of the spin-dependent O(1/m Q ) operators in the Hamiltonian of each heavy meson corresponding to M κ pc .This is equivalent to assume that there is no significant heavy meson-antimeson interaction at this order either.These spin-dependent operators are the ones that break the degeneracy between different total spin heavy meson states.
The results of Table I are valid for any light-quark flavor content, however, the Λ value does depend on the lightquark flavor of the heavy mesons and therefore so does the position of the corresponding static energy on the spectrum of static energies.Now let us discuss the mixing terms between quarkonium and the heavy meson-antimeson pair.Since the quarkonium field Ψ inherently belongs to a Λ = 0 representation, we should project the heavy meson pair field into the same representation.This can be achieved with the projection vector P κ0 [24,26], which is defined by the eigenvalue equation with Λ = |λ| and S κ the spin operator for the κ representation.From textbooks, as for instance Ref. [37], one can find that where Y κα (r) is a spherical harmonic.Notice, that one can think of Y κα (r) as rank κ irreducible tensor made out of powers of r, with α acting as the spin index.For instance, for κ = 1 the projection vector is just (P 10 ) α = ir α .The phase in Eq. ( 6) is arbitrary and is chosen so the projection vector transforms under time reversal as a spin-κ irreducible tensor, see Appendix A.
The most general leading order operator containing Ψ and M κ pc which is invariant under D ∞h and O(3) transformations and discrete symmetries is as follows Note that summation over repeated spin indices is implicit throughout the paper.More details on the notation for the spin indices can be found in Appendix A. If we look at Table I, we can see that there is at least one coupling of each heavy meson-antimeson pair to quarkonium through the operators in Eq. (7).One can check that the Lagrangian in Eq. ( 7) is the most general set of leading order mixing operators with the following argument.The only objects one could add to the mixing terms in Eq. ( 7) that do not add a heavy-quark mass suppression are scalar matrices build out of r and S κ .As was shown in Ref. [26] the projectors form a basis for these matrices, since As the projection vectors P κλ are orthogonal, the form of Eq. ( 7) is not altered by adding the projectors P κΛ .We can compare our operator for the mixing of quarkonium to the lowest lying heavy meson-antimeson pair, the κ pc = 1 −− term in Eq. ( 7), with the one in Refs.[17][18][19][20][21].In the case of Refs.[17][18][19] the operator coincides except for an i factor needed for invariance under timer reversal.However, the results for the masses and widths should not be affected by this phase.In the case of Refs.[20,21] the mixing operator does not seem to take into account the light-quark spin state that couples with quarkonium.
The contribution of the heavy meson-antimeson pair on the quarkonium masses and widths can be computed in standard perturbation theory.First, let us define the following states: with ψ nl the wave function solution of the Schrödinger equation where n is the principal quantum number and l(l + 1) is the eigenvalue of L 2 For the heavy meson-antimeson pair the wave functions are plane waves labeled by the relative momentum k.The partial wave decomposition of the plane wave is as follows: where j l is a spherical Bessel function.Ignoring the heavy-quark spin, the total angular momentum of the heavy meson pair is L = L Q Q + S κ and the eigenvalue of L 2 is ℓ(ℓ + 1).A heavy meson pair state with momentum k = |k| and ℓ total angular momentum is given by with C a Clebsch-Gordan coefficient.Recall that all repeated spin indices are summed.In order for the state in Eq. ( 14) to have definite parity the sum over l should be understood to run only over even or odd values.Let us compute the expected value of the mixing term in the Lagrangian in Eq. ( 7) for the states in Eqs.(10) and (14). with Now, we compute the self-energy contribution where with ) and µ and m T being the heavy meson-antimeson pair reduced and total masses, respectively.The contribution to the quarkonium state mass corresponds to the real part of Eq. ( 18) where P stands for Cauchy principal value.The contribution to the width of the quarkonium state is obtained from the imaginary part of Eq. (18).We obtain We note that a similar result has been obtained in Ref. [21] 1 .

III. MATCHING TO NRQCD
In this section we obtain V κ pc and V κ mix as NRQCD [38][39][40] correlators.For a more self-contained discussion we also reproduce the result for the quarkonium static potential which was obtained originally in Refs.[41][42][43].The matching of the quarkonium 1/m Q suppressed potential can be found in Ref. [22].
Let us define the following NRQCD operators with ψ a Pauli spinor field that annihilates a heavy quark and χ the one that creates a heavy antiquark.The operators Q κ p contain the light-quark fields.For κ p = (1/2) + , (1/2) − , and (3/2) − light-quark states, suitable operators are as follows: with q(t, x) a light-quark Dirac field.The Q κ p operators are obtained replacing P + by P − in Eqs. ( 23)- (25).Details on the construction of irreducible tensor products, such as the one in Eq. ( 22), can be found in Appendix A. The Wilson line ϕ is defined as where P is the path-ordering operator.
The operators in Eqs. ( 21) and ( 22) interpolate for the quarkonium and heavy meson pair fields, respectively.The matching condition from NRQCD to BOEFT reads as The normalization factors are in general functions of Z = Z(r, p).The light-quark flavor quantum numbers must match in both sides of Eq. ( 28), therefore M κ pc corresponds to a single heavy meson-antimeson pair.For M κ pc fields belonging to isospin or chiral symmetry representations one should just consider the appropriate sums over the light-quark flavor in the left-hand side of Eq. ( 28).Since quarkonium and heavy meson-antimeson pairs mix at leading order one could argue that these states are not the appropriate ones to describe the system and that one should work with a basis of states that diagonalize the Hamiltonian at leading order.However, we know from experience that quarkonium and heavy meson-antimeson pairs are useful states to describe the heavy quark-antiquark spectrum.Therefore, there must be some regime of r in which the mixing can be treated as a perturbation.Thus, let us first assume that the strings in Eqs. ( 21) and ( 22) overlap with well-separated NRQCD static eigenstates.Let us match the NRQCD and BOEFT correlators: We contract the heavy-quark fields in the correlators of the right-hand side of Eqs. ( 29)-( 31) and define the following objects with ϕ C being a Wilson line along the path with the paths C i , i = 1, .., 4 defined in Fig. 2. The right-hand sides of Eqs. ( 32)-( 34) are the traces of a product of color matrices and Eq. ( 32) is a static Wilson loop.
x 2 Wilson line paths appearing in Eqs. ( 32)-( 34).The bold line represent the paths while the black dots stand for the light-quark operators.
The right-hand side of Eqs. ( 29)-( 31) is computed from the BOEFT Lagrangian in Eqs. ( 1) and ( 2) and together with Eqs. ( 32)-( 34) we arrive at where the trace acts on the light-quark spin space and the projectors P κΛ are defined in Eq. ( 8).The mixing potential reads as Now we consider the case when the strings in Eqs. ( 21) and ( 22) overlap with two static states of similar energy with |i, Σ + g ; x 1 , x 2 ⟩ (0) being eigenstates2 of the NRQCD Hamiltonian in the static limit, H (0) , with energies The coefficients a x i (x 1 , x 2 ), are the overlaps of the string x = Ψ, κ pc with the static state i.The BOEFT potentials for Ψ and (P * κ0 ) α M κ pc α can be arranged as matrix with V (0) diagonalizes the BOEFT potential matrix in terms of the mixing angle θ = θ(r) The potentials from Eqs. ( 42) and ( 44) are related as follows V (0) Now we set the following matching condition i.e., we match the NRQCD static eigenstates to the BOEFT eigenstates resulting from diagonalizing the potential matrix in Eq. ( 44).From the matching condition in Eqs. ( 48) and (49) follows that To obtain the NRQCD overlap factors a x i in terms of BOEFT quantities we bracket Eqs. ( 39) and ( 40) with the static states on both sides of the equations.Then the remaining bracket in the left-hand side is obtained using Eqs.( 48) and (49) and Eqs. ( 27) and (28).We find Inserting the expansion into static eigenstates in Eqs. ( 39) and ( 40) into the correlators in left-hand side of Eqs. ( 29) and ( 31) and using Eqs.( 50)- (52), we arrive at the following expressions If, using a nonperturbative technique, the left-hand side of Eqs. ( 53)-( 55) is computed, then one can use the parametrization of the right-hand side to fit the data and obtain V 1Σ + g , V 2Σ + g and θ as has been done in Ref. [2].Finally, inverting Eqs. ( 45)-( 47) we find If we expand Eqs. ( 56)-( 58) for 2V l mix ≪ |V (0) Ψ | and use the result in Eqs. ( 53)-( 55) we recover the matching expression of the first part of this section in Eqs. ( 45)-( 47) for the case of well separated static energies.Therefore the mixing can be considered as a perturbation for mixing angles close to 0 or π/2, corresponding to the condition 2V l mix ≪ |V (0)

IV. COMPUTATION OF THE QUARKONIUM SPECTRA
In this section we compute the charmonium and bottomonium spectra and wave functions and use these results to compute the contribution of the lowest lying heavy meson-antimeson pair, corresponding to the κ pc = 1 −− state in Table I, without and with closed strangeness to the quarkonium states masses and widths.For these two heavy meson-antimeson pairs there is available lattice data for the mixing potential from Ref. [2].To improve the overall accuracy in the determination of the quarkonium spectra in the threshold region and since the heavy meson threshold contributions are sensitive to the energy gap between a quarkonium state and the threshold we compute it up to O(1/m Q ) accuracy.

A. RS' scheme
The quarkonium spectrum at leading order is obtained by solving the Schrödinger equation in Eq. ( 11), for which a heavy-quark (pole) mass, m Q , value must be specified.The second input for the Schrödinger equation is the static potential, V Ψ .Both these objects suffer from renormalon ambiguities when computed in perturbation theory [44].The total energy of a quarkonium system is a physical observable and therefore must be free of these ambiguities.At leading order the total energy is given by Ψ , hence the renormalon ambiguities of the heavy quark mass and the static potential cancel each other.Therefore, it is convenient to work in a scheme in which the renormalons are subtracted from these two quantities.We use the modified renormalon subtraction scheme (RS') of Ref. [45].The subtracted heavy-quark mass and the static potential are defined as follows: All the quantities must be computed to the same order in α s and at the same renormalon subtraction scale ν f .We use the expressions up to O(α 4 s ) that can be found in Ref. [46].We work with ν f = 0.7 GeV and take the heavy quark mass values m RS ′ c = 1.592 (41) GeV and m RS ′ b = 4.949 (41) GeV determined in Ref. [47] and the normalization of the renormalon N m = 0.5626(260) from Ref. [48].The values of α s in the MS scheme are obtained using RunDec at 4-loop accuracy [49,50].

B. Static potentials
In Ref. [2] the static energies of the heavy quark-antiquark pair coupled to a heavy meson-antimeson pair were studied using lattice QCD.The lattice computation was done with n f = 2 degenerate light quarks with masses corresponding to an unphysical pion mass ≈ 640 MeV and a lattice spacing a −1 ≈ 2.37 GeV.Using the data for the ground and first excited states as well as the mixing angle in Eqs. ( 45)- (46) one can obtain the lattice determination of the quarkonium static potential and the quarkonium-heavy-meson pair mixing potential.Similarly in Ref. [3] the static energies were studied in n f = 2 + 1 light-quarks and in addition to the states of Ref. [2] the first strange heavy meson-antimeson-pair was included.Unfortunately, in Ref. [3] the mixing angles are not available and the static potentials cannot be extracted without heavy modeling.We show the original data of Ref. [2] in Fig. 3 and data transformed with Eqs. ( 45)- (46) in Fig. 4. It is interesting to note that the small bump in the first excited state (yellow triangles in Fig. 3), which in the range of r where the bump occurs is dominated by the heavy mesonantimeson component, disappears in the transformed data for the heavy meson-antimeson static potential (yellow triangles in Fig. 4).This seems to indicate that the short-distance heavy meson-antimeson interaction is a result of the mixing with quarkonium.Therefore, since the transformed data for the heavy meson-antimeson static potential is completely flat, our choice for the heavy meson-antimeson static potential in Eq. ( 3) is consistent with the lattice data.Furthermore, the plot of the data for the mixing angle in the right-hand side of Fig. 3 shows that the mixing angle is close to 0 or π/2 except for a narrow region between r ∼ 1.2 − 1.3 fm around the string-breaking distance.Therefore, the mixing potential can be considered a perturbation for most of the range of r. θ(r) Figure 3. Plot of the lattice data of Ref. [2].In the left-hand side we plot the data for the ground and first excited static energies in the quarkonium sector.The open blue circles and the open yellow triangles correspond to the ground and first excited states respectively.Note that, in Ref. [2] the origin of energies is set at the energy of the heavy meson pair for the largest r computed.In the right-hand side we plot the data for the mixing angle.
In the following we focus on finding a parametrization of the quarkonium static potential.The quarkonium potential to be used in the Schrödinger equation is constructed combining the short-distance perturbative expression with the long-distance lattice data.For distances r ≤ 1 GeV −1 the static potential is set to the perturbative expression in the RS' scheme.The convergence of the r-dependence of the perturbative potential is poor at short distances if the renormalization scale is fixed [51] due to the presence of large logarithms.If we set ν ∼ 1/r the large logarithms associated with the soft scale are resummed into the running of α s (ν) and the convergence is improved.Specifically, we set ν = 2/r.
For distances r > 1 GeV −1 the static potential is set to a fit of the lattice data.We parametrize the lattice data with the following function: where the linear coefficient is fixed to σ = 0.21 GeV 2 to reproduce the long-range behavior found in Ref. [52].We constrain the parameters so that the slope at the matching point r m = 1 GeV −1 is equal to that of the perturbative  transformed into the quarkonium and heavy meson static potentials using Eqs.( 45) and ( 46), respectively.The lines correspond to our parametrization of the quarkonium static potential.The dashed orange line is the perturbative potential, shifted by δE offset to match the scale of the lattice data, plotted up to the matching point rm.The continuous red line is the expression in Eq. ( 61) fitted to the lattice data plotted from the matching point onward.The full potential formed by the orange and red lines corresponds to Eq. ( 66).Note that, in Ref. [2] the origin of energies is set at the energy of the heavy meson pair for the largest r computed.
potential.The rest of parameters are obtained by fitting the lattice data, we find Due to possible powerlike divergences in the lattice computations the normalization of the lattice data is unknown.This ambiguity is removed by shifting V L by a constant chosen to ensure that at the matching point the shifted lattice parametrization is continuous with the perturbative potential.Finally, the static potential we use in the Schrödinger equation reads as with δE offset = 0.741 GeV.We plot Eq. ( 66) in Fig. 4.

C. Mixing potential
The lattice determination of the mixing potential is obtained using Eq. ( 47) and the data of Ref. [2].Since the mixing potential is proportional to the difference between the ground and first excited states it is not affected by the ambiguity in the normalization of the energies of the lattice computation.
To parametrize the mixing potential we use a function that interpolates between the short-and long-distance behaviors as introduced in Ref. [53].For distances r ≪ 1/Λ QCD the relative momentum between the heavy quarks can be integrated out and the quarkonium-heavy-meson pair mixing can be studied in weakly coupled pNRQCD [54,55].In this regime the heavy-quark fields can be decomposed into color-singlet and color-octet fields.At leading order in the multipole expansion the quarkonium states overlap with the singlet field while the heavy meson pair states can in principle overlap with both.The transition between the quarkonium state and the octet piece of the heavy meson pair is generated at leading order by a chromoelectric dipolar operator.The r factor in this operator provides the correct dependence on r of the mixing operator for κ pc = 1 −− in Eq. ( 7) and produces a linear dependence of the mixing potential at leading order.Furthermore, the chromoelectric operator has the right quantum numbers to create the light-quark content of the heavy-meson pair.The second contribution corresponds to the overlap of the quarkonium state with the singlet piece of the heavy-meson pair.In this case the leading order transition is generated by three dipolar operators, therefore this contribution to the mixing potential has a r3 dependence at leading order.This second contribution is in principle suppressed respect to the first one in the multipole expansion, however the size of the overlaps of the heavy meson pair state with the singlet and octet fields are unknown.For this reason we keep both contributions in our short-distance description of the mixing potential For distances r ≫ 1/Λ QCD the mixing potential can be expanded in powers of 1/(Λ QCD r) n .If we assume that the mixing potential vanishes in the r → ∞ limit then only n ≥ 0 is allowed.By fitting the data with r > 1 fm we find that the long-distance part is well described by We construct the interpolation by summing the short-and long-distance descriptions multiplied by interpolating functions depending on r and a new r 0 parameter.The interpolating functions are w s = (r 0 /(r + r 0 )) n and w l = (r/(r + r 0 )) n for the short-and long-distance pieces, respectively.The r 0 parameter determines the value of r where both interpolating functions are equal.We pick r 0 = 0.25 fm as it is a reasonable point for the breakdown of the multipole expansion.The full parametrization of the potential is as follows The value of n is set to the smallest value that the short-and long-distance potentials dominate in their respective limits, which in this case is n = 7.The rest of the parameters are fitted to the lattice data, we find It is interesting that the value of c 2 is not as suppressed with respect to the one of c 1 as one would expect from the pNRQCD counting which might indicate that the heavy meson pair has a larger overlap with the singlet field than the octet one in the short-distance regime.This is consistent with the slightly attractive behavior of the first excited static state (in yellow triangles in Fig. 3) in the first few short-distance data points 3 .In Fig. 5 we plot the lattice data for the mixing potential and the parametrization in Eq. ( 69).The normalization of the mixing operator in Ref. [2] and ours in the Lagrangian in Eq. ( 7) for κ = 1 coincide, however the heavy meson-antimeson pair interpolating operator is isospin I = 0 unlike ours, in Eq. (22), which corresponds to a single light-quark flavor.In other words, the heavy meson-antimeson pair operator in Ref. [2] interpolates for a field in BOEFT corresponding to the normalized sum of the charged and neutral heavy meson-antimeson pairs.Since the two lattice light quarks are degenerate we have We will use this mixing potential for all three light-quark flavors, which is an approximation.However, the light-quark mass used in Ref. [2] is in fact closer to the strange mass than the up or down masses, therefore it can be expected that the approximation produces more accurate results for pairs of heavy mesons with strangeness.In Ref. [3] the static energy spectrum of quarkonium coupled to heavy meson pairs without and with strangeness was obtained in lattice QCD for distances between r ∼ 0.25 fm and r ∼ 1.6 fm.Since the mixing angles were not given the mixing potentials can only be extracted by assuming specific forms of the quarkonium and heavy meson pair static potentials.Assuming a Cornell type potential for the former and a constant for the latter one can fit the mixing potential.Since only the long-distance regime is available we fitted mixing potentials ∼ r −1 , ∼ r transformed into the mixing potential using Eq. ( 47).The red line corresponds to the parametrization in Eq. ( 69) fitted to the lattice data.
constant, the latter being the choice in Ref. [3] analysis.In all three cases the fits were of similar quality.Therefore, we conclude that a reliable estimation of the mixing potentials from the data of Ref. [3] is not possible.The mixing potential has been extracted from the lattice data of Ref. [2] in Refs.[17][18][19].However a different procedure was followed.It was argued that the heavy-meson pair to heavy-meson pair correlator of Ref. [2], interpolates not only for the Σ + g representation but also for Π g and Σ − u representations.As a result the authors of Refs.[17][18][19] argue that lattice data should be fitted with a parametrization that takes into account these extra states in the meson pair to meson pair correlator.Lets us note, that fits to the correlators with extra states where also considered in Ref. [2] but where found not to describe the data well.Since the original data on the correlators of Ref. [2] is not available, the authors of Refs.[17][18][19] resample the lattice correlators using the original parametrization and then fit this resampled data with their parametrization containing the extra states.While we agree on the initial point about the extra states in the meson pair to meson pair correlator, we do not think the resampled data can contain information on these extra states as it was produced from the original parametrization.The quarkonium static potential obtained from the resampled data can be found in Fig. 3 of Ref. [17].If we compare it to the one we obtain, in Fig. 4, it can be observed that the shapes are notably different.In our determination, the shape of the static potential is closer to previous lattice determinations of the static potential that did not include the mixing with the threshold, which is the behavior expected away from the string breaking distance as we discussed at the end of Sec.III.A possible explanation for the extra states in the meson pair to meson pair correlator of Ref. [2] not showing up in their fits is that the extra states are degenerate with the Σ + g one, as we do in the Lagrangian in Eq. ( 3).Therefore, in our opinion the most appropriate way of extracting the static and mixing potentials is to use the two state parametrizations of the lattice correlators.We hope that in the future new lattice studies clarify this issue.
In Refs.[20,21] a model for the mixing potential is used.This consists of a Gaussian shape with the maximum at the string breaking distance, i.e. the value of r when ). Comparing with the mixing potential extracted from the lattice data in Fig. 5 we can see that the model misses important features.The maximum of the mixing potential occurs at r ∼ 0.25 fm instead of the string breaking distance and as a result the overall shape is different.Moreover, as the value of the mixing potential at the string breaking distance must be equal to half the avoided crossing separation 4 , the maximum value of the model mixing potential of Refs.[20,21] is much smaller than the one of the potential extracted from the lattice data. 4The avoided crossing distance is (V

D. 1/mQ quarkonium potential
To improve the accuracy in the determination of the quarkonium spectrum we compute the contribution of the 1/m Q suppressed potential [22] using standard time independent perturbation theory.To obtain an expression for this potential we follow an analogous approach to the one in Sec.IV B for the static potential.We construct the potential by combining the short-distance perturbative expression with the available lattice data for long-distances.We use the leading order perturbative result from Ref. [56].It reads as Notice, that the form of the potential depends on the matching scheme.We use the expression for the Wilson loop matching scheme in accordance with the rest of the paper.As in the static potential, we resumme large soft logs by setting ν = 2/r.The 1/m Q suppressed quarkonium potential has been studied in the lattice in Refs.[57][58][59].We use two datasets with simulation parameters β = 5.85, a = 0.123 fm and β = 6.00, a = 0.093 fm in the quenched approximation.The lattice data is plotted in Fig. 6.To parametrize the lattice data we use the following function which interpolates between the dependence on r −2 from the perturbative expression in Eq. ( 74) and the log r dependence obtained from Effective String Theory [52, 53, 60-64]5 .The parameters of Eq. ( 75) are constrained to reproduce the slope of the perturbative potential at the matching point r m = 1 GeV −1 , the rest of the parameters are fitted to the lattice data.The values we obtained are as follows: The ambiguity in the normalization of the lattice data is removed by shifting V (1) L so that the value at the matching point is equal to that of the perturbative expression in Eq. (74).The full expression of the 1/m Q suppressed potential is with δE (1) offset = −0.088GeV.In Fig. 6 we plot the potential in Eq. (81).We should note that the value of δE (1) offset depends noticeably on the specific matching point r m and the order at which we take the perturbative potential in Eq. ( 74).This is a result of the lack of lattice data at shorter distances and the possible existence of renormalon ambiguities in V p.t. (r).

E. Numerical results
We solve numerically the Schrödinger equation for the quarkonium static potential in Eq. ( 66) obtaining the wave functions Ψ nl and eigenenergies nl , with n and l the principal and angular quantum numbers of the quarkonium state, respectively.Using the wave functions we compute the contribution of the 1/m Q suppressed potential in Eq. (81) in standard quantum mechanical perturbation theory.
) [GeV] Figure 6.Plot of the 1/mQ-suppressed quarkonium potential given in Eq. ( 81) with parameters fitted as described in the text.The orange dashed line corresponds to the perturbative part in Eq. ( 74) while the red continuous line corresponds to the parametrization of the lattice data in Eq. ( 75).The blue circles and green triangles corresponds to the lattice data with lattice coupling β = 5.85 and β = 6.00, respectively, from Refs.[57,58].
The results for E (0) nl and E nl can be found in Tables II-IV for bottomonium states and Tables V-VII for charmonium states.The uncertainties in these two quantities are estimated by their difference when changing the matching point between the perturbative expressions and the lattice data fits from r m = 1 GeV −1 to r m = 0.66 GeV −1 .
At the order we are working, isospin and heavy-quark spin contributions can be neglected, thus we use heavy meson masses reflecting these approximations.The values of the heavy meson masses are obtained by first computing the average of the neutral and charge states, when both are available, and then the spin average of the scalar and vector states.For strange heavy mesons only the last step is necessary.The masses of the physical heavy meson states are taken from the PDG [65].We find (87) The heavy meson-antimeson pair contributions to the quarkonium masses are computed using Eq.(19).For this computation we take the quarkonium binding energy up to next-to-leading order.Since the results in this section are all for heavy meson-antimeson pairs with κ pc = 1 −− we will drop this label.On the other hand, since we will compute the contributions for heavy mesons without and with strangeness we will add a label f indicating the light-quark flavor.We use f = q to denote the sum of the neutral and charged heavy meson-antimeson pair contributions.These two contributions are equal in the isospin limit.To denote the contributions of heavy mesons pairs with closed strangeness we use f = s.The vertex form factor a l ′ nl (k) is computed numerically from the expression in Eq. ( 16) using the wave functions form solving Eq. ( 82) and the mixing potential in Eq. ( 73).We sample it for a range of k from 0 GeV to 6 GeV at 10 MeV intervals.A linear interpolation of this data is then used to compute the mass contribution in Eq. (19).The uncertainty of E f l ′ nl is estimated as follows.The uncertainty of k 2 d is obtained by combining in quadrature the uncertainties of E  nl are taken as our central value and its uncertainty, respectively.The contributions of the heavy meson pairs to the quarkonium spectrum is displayed in Tables II-IV for bottomonium and Tables V-VII for charmonium.
The total masses of the quarkonium states are obtained as The uncertainty of M nl is obtained by adding in quadrature the uncertainties of each term in Eq. (88).Our result show the contribution of these two thresholds to the quarkonium masses is comparable to that of the 1/m Q suppressed potential.Contrary to what it could be expected intuitively, the contribution of the thresholds is larger for the lower lying quarkonium states that the excited ones.The underlying reason is in the shape on the  16) 15 (7) 3(2) −0.4(1.5)5037(91) 5 2158(18) 22 (6) 4 mixing potential (see Fig. 5) which peaks, in absolute value, at r ∼ 0.25 fm, while the excited states wave functions extend to far longer ranges.The uncertainty of our results for the quarkonium masses in Tables II-VII is dominated by the uncertainty in the determination of the heavy quark masses in the RS' scheme.This source of uncertainty cancels out in mass differences which are consequently much more accurate.Furthermore, we can reconstruct the spectrum by adding to the experimental mass of a given state our mass differences with respect to the same state.Since our computation does not include spin-dependent contributions it is convenient to consider as a reference an S-wave state, since the spin-averages of these are independent of such contributions.Additionally, we expect the spin-independent O(1/m 2 Q ) contributions to be smaller for higher excited states, in a similar way as the spin-independent O(1/m Q ) contribution we have computed.Therefore, we choose as experimental reference state the 2S doublet, since this is the higher laying S-wave doublet which has been measured for both charmonium and bottomonium.In Table VIII we show the bottomonium and charmonium spectra shifted so the 2S state mass matches the spin-average of the corresponding experimental masses.Both shifts are within the uncertainty of the heavy quark masses.In Figs. 7 and 8 we show all the experimental bottomonium and charmonium states listed in the PDG with definite J P C [65] compared to the shifted spectra.We also display the hybrid quarkonium states expected to appear in the energy range of the figures and with J P C matching those allowed for S, P and D wave quarkonium.However, one should keep in mind that exotic J P C are possible for quarkonium hybrids, including, for instance, heavy-quark spin partners of the states displayed.The mass values of the hybrid quarkonium are taken from Ref. [66] and also shifted to match the experimental value of spin-average mass of the 2S doublet.To do this, we obtain the 2S state mass for the Σ + g static energy from the lattice data of Ref. [12], which is the same source as for the Π u − Σ − u static energy data, used in Ref. [66] to obtain the hybrid spectra.It should be kept in mind that the hybrid quarkonium states displayed in Figs.7 and 8 are not computed to the same accuracy as the conventional quarkonium ones, since they do not include O(1/m Q ) corrections 6nor heavy meson-antimeson pair contributions.
The widths are computed according to Eq. ( 20).We take Eq.(88) as the input mass.As in the mass computation, we compute the uncertainty k 2 d adding up in quadrature the uncertainties of each term and the size of higher order contributions.We create a random Gaussian sample of k 2 d values and compute the decay widths for each value in the set.The average and standard deviation of these computations are assigned as the central value and uncertainty of the widths.The notation for the widths follows the one for the energy contributions, in particular recall that f = q corresponds to the sum of the widths for the neutral and charged heavy meson pairs.The results for the widths are shown in Tables IX-XI for bottomonium states and Tables XII-XIV for charmonium states.We find values of about 5 − 10 MeV for bottomonium and 10 − 50 MeV for charmonium.The threshold contributions to the mass and width of P -wave quarkonia turn out to be slightly larger than S and D wave quarkonia due to it coupling to the @ @ @ n l 0 1 2  We also include the hybrid bottomonium states (blue lines) from Ref. [66] in the mass range and J P C of the figure.Both conventional and hybrid bottomonium spectra are shifted so the 2S state mass matches the experimental spin average one.
heavy meson-antimeson pairs in two partial wave channels instead of one.As it can be seem from these results the uncertainties are large, particularly in comparison to the uncertainties in the mass contributions.This is a result on a strong dependence of the widths on the value of k 2 d .To illustrate this we plot the value of width as a function of the difference between the quarkonium and the heavy meson pair mass for two specific cases in Fig. 9.We can see that a large range of values of the widths can be produced within the uncertainty of the mass difference.

V. THRESHOLD EFT
Let us consider an EFT for a quarkonium state ψ nl close to a heavy meson-antimeson pair threshold as the one considered in Refs.[28][29][30][31][32][33][34][35][36].The heavy meson fields will be represented by H κ with κ the spin of the light-quark state.The fields H κ carry two spin indices, the first one corresponding to the antiquark and the second one to the quark with the order being reversed for the Hermitian conjugates.Therefore, one should read expressions as Tr β where α and β indices correspond to the light-quark and heavy quark spin, respectively.The spin indices will be in the spherical basis unless stated otherwise.Since we work in the 0 + 1 1 + 0 + + 1 + + 2 + + 2 + 2 3 3.00 We also include the hybrid charmonium states (blue lines) from Ref. [66] in the mass range and J P C of the figure.Both conventional and hybrid charmonium spectra are shifted so the 2S state mass matches the experimental spin average one.fields.The bilinear terms in the EFT read as The quarkonium-heavy meson pair couplings at leading order in the heavy quark mass expansion have the following general from and ξ = ξ/|ξ|.To match the common choice in the literature, the quarkonium field ψ nl is chosen to transform under time reversal as a spherical harmonic under complex conjugation.Therefore, a factor i l is needed to match the time reversal transformation of a spin l field and for the whole operator to be invariant under this symmetry7 .To conserve parity only the angular momentum that fulfill p(−1) l+l ′ = 1 are allowed.Likewise requiring charge conjugation invariance leads to the constraint c(−1) l+l ′ = (−1) l+l ′ +κ = 1.
The couplings of quarkonium with l = S, P, D to the lowest lying heavy meson-antimeson pairs (κ p1 1 = κ p2 2 = (1/2) + and κ pc = 1 −− ) up to two derivatives read as The spin indices in Eqs. ( 91)-( 93), explicit or implicit, are in the Cartesian basis.The first two terms can be found in Refs.[28][29][30][31] with different normalizations for the couplings.The equivalence with Refs.[28,31] is for Refs.[29,30] an extra minus sign is needed to account for a different definition of ∇ ⃗ ⃗ .Now we match the threshold EFT to BOEFT.At the tree level the matching can be obtained expanding the fields in eigenstates of the relative motion Hamiltonian.For the quarkonium field for instance with ψ nl defined in Eq. ( 12).In the same way we expand the heavy meson-antimeson pair field.Since in this case the eigenfunctions are plane waves we have Introducing the partial wave expansion of Eq. ( 13) to Eq. ( 96) and after some manipulations we obtain So far we have made no approximation.Now, we note that in the case a quarkonium state mass is close to that of a heavy meson-antimeson pair then the relative momentum between the heavy quarks in the quarkonium state ∼ 1/r is larger than the relative momentum between the heavy mesons k.Therefore, we can expand the spherical Bessel function for kr ≪ 1: Using this expansion in Eq. (97) all the dependence on r factorizes.After a few more manipulations we arrive at We set the following matching condition V = + . . .We can now apply these field expansions into the BOEFT couplings and obtain the tree level matching of the couplings of quarkonium to heavy meson pairs in the threshold EFT.The matching is shown diagrammatically in Fig. 10.We arrive at with It is interesting to consider also the matching of the quarkonium bilinear.At tree level these can be easily obtained using the expansion in Eq. (95).However, in this case one should also consider the self-energy contribution which appears at one loop.The matching is shown diagrammatically in Fig. 11.The key point is to identify the momentum region in the loop diagram in the BOEFT side (left) that matches the heavy meson loop in the threshold EFT side (right).The BOEFT diagram is given in Eq. ( 18) and contains two momentum scales k ∼ k d and k ∼ 1/r.The contribution of the first one matches the heavy meson loop, while the latter gives a contribution to the residual mass of the quarkonium state in the Lagrangian in Eq. ( 89) where E nl is the energy of the quarkonium state computed in BOEFT including the heavy meson-antimeson pair contributions except the one of the nearby heavy meson-antimeson pair considered explicitly in the threshold EFT and To determine if the threshold expansion can be carried out, we compare the values of the |k d | momentum scale of the mesons with the value of ⟨1/r⟩ for the quarkonium states in Table XV.The computation of |k d | is carried out for the transitions between quarkonium and the lowest lying heavy meson pair thresholds without and with closed strangeness.The value of ⟨1/r⟩ is computed including O(1/m Q ) corrections to the quarkonium wave function and considering intermediate states up to n = 8.In principle, these corrections could be sizable for charmonium states since Λ QCD /m c ∼ 18%.However, we find that only the ground states get contributions comparable to this parametric estimate.In practice, for excited states the contributions from lower lying intermediate states almost cancel with those of higher laying states leading to very small contributions.Therefore, we expect the parametric estimate of O(1/m 2 Q ) contributions to ⟨1/r⟩, i.e. (Λ QCD /m c ) ∼ 3.5% and (Λ QCD /m b ) ∼ 0.4%, to be significantly larger than actual contributions.The computation of |k d | on the other hand involves uncertainties inherited from the uncertainty in the determination of the quarkonium masses.Unfortunately, the effect is naturally more significant for states close to threshold.Therefore, at the present accuracy in the determination of the quarkonium masses, we can only rule out the threshold expansion validity for certain states while in no case it can be completely confirmed as a good approximation.For the latter cases we compute the values of the quarkonium-heavy-meson-pair couplings, which we collect in Table XVI.The value of these couplings only depend on the quarkonium bound state through the wave function, and therefore we expect our values to be reliable despite the uncertainty on the quarkonium masses.Nevertheless, due to possible O(1/m Q ) contributions to the mixing potential, the accuracy of the couplings is limited to corrections O(Λ QCD /m Q ), which is reflected in the uncertainty.Next, we compute the contribution to the residual mass of the quarkonium field in the threshold EFT, E κl ′ nl defined in Eq. ( 105).This quantity also only depends on the quarkonium state through the wave function and its uncertainty is assessed as of the size of corrections O(Λ QCD /m Q ).The values of E κl ′ nl can be found in Table XVI for bottomonium and charmonium.

VI. CONCLUSIONS
The first excited state in the spectrum of static energies in the quarkonium sector corresponds to a heavy mesonantimeson pair [2,3].Pairs composed of heavier heavy mesons appear as subsequent excited states, along with quarkonium hybrid static energies.The heavy meson-antimeson static energies are characterized by the total spin of the light quarks and its projection onto the heavy quark-antiquark axis, thus each heavy meson-antimeson pair is associated with more than one static state.EFTs have been built to describe the heavy quark-antiquark bound states supported by the spectrum of static energies.These EFTs incorporate the heavy quark mass expansion and the adiabatic expansion between the heavy and light degrees of freedom.For standard quarkonium the EFT is known  as strongly coupled potential NRQCD [22,23], and its extension to nontrivial light degrees of freedom is called Born-Oppenheimer EFT (BOEFT) [16,[24][25][26]66].
In this paper we have shown how to incorporate into BOEFT the heavy meson-antimeson pairs and have obtained all the leading order operators coupling them to quarkonium.The Lagrangian containing these couplings can be found in Eq. (7).Using these couplings, we have obtained the expressions for the contribution of heavy meson-antimeson pairs to the quarkonium masses and decay widths in perturbation theory.These formulas, in Eqs.(19), ( 20) and ( 16), depend on the total spin of the light-quarks in the threshold state that couples to quarkonium, the mixing potential accompanying the coupling operator, the mass gap between the threshold and the quarkonium state and the wave function of the quarkonium state.
In Sec.III we have discussed the matching of the new potentials, in particular of the mixing one, to NRQCD.The matching has been obtained for both when the mixing potentials can be considered a perturbation and when not.We have also shown that the second case reduces to the first when the separation between the static potentials of the quarkonium and the heavy meson-antimeson pair is larger than the mixing potential, which is the case for most of the range of r except for a small region around the string breaking distance.In Ref. [2] the ground and first excited states for the coupled system of quarkonium with the lowest lying heavy heavy meson-antimeson pair were obtained in lattice QCD.Using this data and Eqs. ( 45)-( 47), the quarkonium and heavy meson pair static potentials as well as the mixing potential can be obtained.It is interesting that the small bump in the excited state at short distances, see Fig. 3, which could be interpreted as a heavy meson-antimeson interaction, disappears completely in the heavy meson-antimeson static potential in Fig. 4.This highlights the importance of taking into account the mixing with heavy quark-(anti)quark states when studying the heavy meson-(anti)meson interactions.
We computed the contribution of the lowest lying heavy meson-antimeson pairs without and with closed strangeness to the masses and widths of the bottomonium and charmonium states with l = S, P, D and n = 1, .., 6 covering the mass range where exotic quarkonium states have been discovered.The quarkonium static potential is obtained combining fix order results in the RS' scheme for the short-distance part of the potential and a fit to lattice data for the medium and long-ranges.To increase the accuracy in the determination of the quarkonium spectrum in the threshold region we compute the quarkonium spectra up to O(1/m Q ).The quarkonium 1/m Q suppressed potential is also parametrized using perturbation theory for the short-distance part and a fit to lattice data [57,58] for the remaining range of r.Our results for the quarkonium masses and the contribution of the lowest lying thresholds are shown in Tables II-IV for bottomonium and Tables V-VII for charmonium.Our result show the contribution of these two thresholds to the quarkonium masses is comparable to that of the 1/m Q suppressed potential.The uncertainty associated to the heavy quark mass can be eliminated by shifting the bottomonium and charmonium spectra to match the experimental mass of a given reference state at the price of not giving a prediction for this reference state.In Table VIII we show such spectra taking as a reference the spin-average of the 2S doublet.The resulting spectra is compared to experimental values in Figs.7 and 8.For the widths we find values of about 5−10 MeV for bottomonium, in Tables IX-XI, and 10 − 50 MeV, in Tables XII-XIV, for charmonium.
Unfortunately, the contributions of the thresholds to a quarkonium state mass and width become very imprecise when the mass gap between them is similar to the uncertainty in the quarkonium mass determination.To improve the accuracy of the threshold contributions, the heavy-quark spin dependent contributions should be taken into account.These can be split between the ones affecting the heavy meson masses, which are of O(1/m Q ), and the ones to the quarkonium masses which are O(1/m 2 Q ).For an accuracy of a few MeV the former ones should be enough in the bottomonium sector but both would be necessary for charmonium states.
In Sec.V, we have discussed the matching of BOEFT with heavy meson-antimeson degrees of freedom to a threshold EFT containing as explicit degrees of freedom a quarkonium state and the heavy mesons of a nearby threshold interacting through contact operators.This is possible when the relative momentum of the heavy mesons is smaller than the inverse of the size of the quarkonium state.The latter being a measure of the relative momentum between the heavy quarks in the quarkonium state.We obtained the matching expression, in terms of the mixing potential and the quarkonium wave function, of the series of couplings with increasing derivatives of a quarkonium state to the heavy meson-antimeson pair, which can be found in Eq. (103).At the current level of accuracy we are only able to rule out for which states and heavy meson thresholds this matching is not valid.Nevertheless, we provide the values of the couplings to the two lowest lying meson-antimeson pairs with the quarkonium states which are not ruled out in Table XVI.This is possible since the values of the couplings do not depend directly on the mass gap between the threshold and the quarkonium state.
Using the threshold EFT and following the analysis from Refs.[68][69][70][71] one can study the heavy meson molecule picture for exotic states close to a heavy meson-antimeson threshold.As pointed out in Ref. [68] the molecular nature of X(3872) can be explained by an accidental fine-tuning of χ 1 (2P ) mass to the DD * threshold which would result in a abnormally large scattering length for the heavy meson-antimeson scattering.Such scenario is compatible with our results, however to confirm it, it would require a high precision determination of the χ 1 (2P ) mass.
Hybrid quarkonium states also are expected to appear in the threshold region.The hybrid states associated to the lowest lying gluelump, with κ pc = 1 +− , appear at 4.000 − 4.150 GeV and 10.690 − 10.790 GeV [16,66,67] in the charmonium and bottomonium sectors, respectively.Meanwhile, the ones associated to the second lowest lying gluelump, with κ pc = 1 −− , appear at 4.5 GeV and 11.14 GeV [72] in the charmonium and bottomonium sectors, respectively.For the former case the mixing with quarkonium is heavy quark mass suppressed [16] and its effects are of O(1/m 2 Q ).However, in the latter case the mixing is not suppressed and its effects can be, in principle, of the same order as the heavy meson-antimeson pair of similar mass.Therefore, its study would be interesting.Furthermore, for an accurate computation of their hybrid quarkonium masses, which is also necessary for a good understanding of the spectrum of quarkonium-like states in the threshold region, it is also necessary the study of the mixing of quarkonium hybrids with the heavy meson-antimeson pairs.For this objective, the formulation of the mixing terms of hybrid quarkonium states with heavy meson-antimeson pairs should be worked out as well as the matching expression as NRQCD correlators.The latter should be computed with lattice QCD or models in order to obtain numerical evaluations.

Figure 4 .
Figure 4. Plot of the static potentials.The blue open circles and the yellow triangles correspond to the lattice data of Ref.[2] transformed into the quarkonium and heavy meson static potentials using Eqs.(45) and (46), respectively.The lines correspond to our parametrization of the quarkonium static potential.The dashed orange line is the perturbative potential, shifted by δE offset to match the scale of the lattice data, plotted up to the matching point rm.The continuous red line is the expression in Eq. (61) fitted to the lattice data plotted from the matching point onward.The full potential formed by the orange and red lines corresponds to Eq. (66).Note that, in Ref.[2] the origin of energies is set at the energy of the heavy meson pair for the largest r computed.

Figure 5 .
Figure 5. Quarkonium-heavy meson pair mixing potential.The blue open circles correspond to the lattice data of Ref.[2] transformed into the mixing potential using Eq.(47).The red line corresponds to the parametrization in Eq. (69) fitted to the lattice data.

( 1 )
nl , m RS ′ Q and the size of higher order contributions O(Λ 2 QCD /m 2 c ) ∼ 40 MeV, O(Λ 2 QCD /m 2 b ) ∼ 4 MeV.The mass contribution is computed for a random Gaussian sample of k 2 d with the standard deviation set to its uncertainty.The average and standard deviation of the set of results for E f l ′

Figure 7 .
Figure 7.Comparison of the experimental bottomonium spectrum (black dots) with the spectrum we have obtained (red lines).We also include the hybrid bottomonium states (blue lines) from Ref.[66] in the mass range and J P C of the figure.Both conventional and hybrid bottomonium spectra are shifted so the 2S state mass matches the experimental spin average one.

Figure 8 .
Figure 8.Comparison of the experimental charmonium spectrum (black dots) with the spectrum we have obtained (red lines).We also include the hybrid charmonium states (blue lines) from Ref.[66] in the mass range and J P C of the figure.Both conventional and hybrid charmonium spectra are shifted so the 2S state mass matches the experimental spin average one.

Figure 9 .
Figure 9. Plots of the decay width dependence with the mass difference between the quarkonium and the two meson masses.In the left we plot the width of the 5S bottomonium state decaying to In the right we plot the width of the 1D bottomonium state decaying to DD.From the plots we can see that for variations of the mass difference of the order of its uncertainty the values of the widths can change by large amounts.

Figure 10 .
Figure 10.Matching of the quarkonium-heavy meson pair vertex between BOEFT (left) and the threshold EFT (right).The double continuous lines represent the quarkonium state, the double dashed line represents the heavy meson pair field in BOEFT and the single dashed lines correspond to the heavy mesons in the threshold EFT.The dots in the right-hand side stand for the vertices with extra derivatives.

Figure 11 .
Figure 11.Matching of the quarkonium two-point function between BOEFT (left) and the threshold EFT (right).The double continuous lines represent the quarkonium, the double dashed line represent the heavy meson pair field in BOEFT while the single dashed lines correspond to the heavy mesons in the threshold EFT.

Table I .
Total spin, parity, charge conjugation, and D ∞h representations of the light quark-antiquark pair combinations of the three lightest light-quark states forming heavy mesons.Each block of states, separated by a single horizontal line, corresponds to degenerate or nearly degenerate states.

Table VII .
Spectrum of D-wave charmonium states.All entries in MeV.

Table VIII .
[65]omonium (left)and charmonium (right) spectra with the origin of energies adjusted to the experimental value of the spin average of the respective 2S states.The experimental spin-average mass is marked by (e) and taken from the PDG[65].

Table IX .
Widths of S-wave bottomonium states.All entries in MeV units.

Table XI .
Widths of D-wave bottomonium states.All entries in MeV units.b.t.stands for below threshold.nonrelativisticregime, we treat the antiparticle fields (denoted with a bar) as independent fields from the particle

Table XIV .
Widths of D-wave charmonium states.All entries in MeV units.

Table XVI .
Couplings g (l ′ ,d) of the charmonium (left) and bottomonium (right) states with quantum numbers (l, n) to heavy mesons in the l ′ partial wave.We display only the cases for which |k d | < ⟨1/r⟩ within uncertainty as displayed in TableXV.All dimension-full entries are in GeV unless indicated otherwise.Note that the couplings for B B and D D pairs apply to both neutral and charged meson pair cases.