QCD spin effects in the heavy hybrid potentials and spectra

The spin-dependent operators for heavy quarkonium hybrids have been recently obtained in a nonrelativistic effective field theory approach up to next-to-leading order in the heavy-quark mass expansion. In the effective field theory for hybrids several operators not found in standard quarkonia appear, including an operator suppressed by only one power of the heavy-quark mass. We compute the matching coefficients for these operators in the short heavy-quark-antiquark distance regime, $r\ll 1/\Lambda_{\rm QCD}$, by matching weakly-coupled potential NRQCD to the effective field theory for hybrids. In this regime the perturbative and nonperturbative contributions to the matching coefficients factorize, and the latter can be expressed in terms of purely gluonic correlators whose form we explicitly calculate with the aid of the transformation properties of the gluon fields under discrete symmetries. We detail our previous comparison with direct lattice computations of the charmonium hybrid spectrum, from which the unknown nonperturbative contributions can be obtained, and extend it to data sets with different light-quark masses.


I. INTRODUCTION
One of the long-standing, unconfirmed, predictions of QCD is the existence of hadrons in which gluonic excitations play an analogous role as constituent quarks in traditional hadrons. This kind of states are divided into two classes depending on whether they contain quark degrees of freedom or not. In the case that the state is formed purely by gluonic excitations it is called a glueball, while when the state contains both quark and gluonic degrees of freedom it is called a hybrid. The experimental identification of any of such states has been up until now unsuccessful. In the case of glueballs, this can be understood as owing to the fact that the lowest-lying states, as predicted by lattice QCD calculations [1,2], have quantum numbers coinciding with those of standard isosinglet mesons, and therefore a strong mixing is expected. Glueballs with exotic J P C , such as 0 +− , 2 −+ , or 1 −+ , are expected to appear at rather large masses.
For hybrid states the experimental identification may be simpler since the interplay of quark and gluonic degrees of freedom enlarges the range of possible quantum numbers J P C , including exotic ones among its lowest mass states. Nevertheless, if the quarks forming the hybrid state are light, the hybrids are still expected to appear at the same scale as conventional mesons, Λ QCD , leading again to the expected large mixings if the quantum numbers J P C of the hybrids are not explicitly exotic. On the other hand, hybrids containing heavy quarks, called heavy or quarkonium hybrids, develop a gap of order Λ QCD with respect to the respective states containing only the heavy-quark component, i.e. the standard quarkonium states. Therefore, quarkonium hybrids are expected to be the states including gluonic excitations that are easier to identify experimentally.
It is precisely in the quarkonium spectrum, close and above the open-flavor thresholds, that in the last decade tens of exotic heavy quarkonium-like states have been discovered in experiments at B-factories (BaBar, Belle, and CLEO), τ -charm facilities (CLEO-c and BESIII) and hadron colliders (CDF, D0, LHCb, ATLAS, and CMS). These states are the so-called XYZ mesons. Several interpretations of the XYZ mesons have been proposed. In these interpretations, XYZ mesons are bound states of a heavy-quark-antiquark pair with non-trivial light degrees of freedom. In the case that the light degrees of freedom are light quarks, different tetraquark pictures emerge depending on the spatial arrangement of the light quarks with respect to the heavy quarks. If the light degrees of freedom are gluonic, the picture that emerges is that of a quarkonium hybrid. So far there is no conclusion on which interpretation is the correct one, see Refs. [3][4][5][6][7] for reviews of the experimental and theoretical status of the subject. Quarkonium hybrids are characterized by the separation between the dynamical energy scales of the heavy quarks and the gluonic degrees of freedom. The gluon dynamics is nonperturbative and, therefore, happens at the scale Λ QCD , while the nonrelativistic heavy-quark-antiquark pair bind together in the background potential created by the gluonic excited state at a lower energy scale mv 2 Λ QCD , where m is the heavy-quark mass and v their relative velocity. This separation of energy scales is analogous to that of the electrons and nuclei in molecules, and has led to the observation that quarkonium hybrids can be treated in a framework inspired by the Born-Oppenheimer approximation [8][9][10][11][12][13][14]. In recent papers [15][16][17][18] an effective field theory (EFT) formulation of the Born-Oppenheimer approximation, called the BOEFT, has been developed and used to compute the quarkonium hybrid spectrum. In this paper we will rely on the hierarchy described above, i.e.
Λ QCD mv 2 and work under the further assumption that mv Λ QCD . The advantage of this assumption is that the nonperturbative dynamics can be factored out and its effects encoded in nonperturbative gluonic correlators, allowing for a clear theoretical analysis of the heavy hybrid spin contributions. In the case in which mv ∼ Λ QCD , the potentials will be given by generalized Wilson loops, however their spin structure will be the same.
The spin-dependent operators for the BOEFT have been presented in Ref. [19] up to O(1/m 2 ).
The most interesting feature, also pointed out in Ref. [20], is that quarkonium hybrids, unlike standard quarkonium, receive spin-dependent contributions already at order 1/m. At order 1/m 2 there are spin-dependent operators analogous to those appearing in the case for standard quarkonium as well as three new operators that are unique to quarkonium hybrids. The matching coefficients of these operators, the spin-dependent potentials, are generically characterized as the sum of a perturbative contribution and a nonperturbative one. The perturbative contribution corresponds to the spin-dependent octet potentials and only appears in the operators analogous to those of standard quarkonium. The nonperturbative contributions can be written as a power series in r 2 with coefficients encoding the nonperturbative dynamics of the gluon fields. In this paper, we compute the spin-dependent potentials by matching weakly-coupled potential NRQCD (pNRQCD) [21,22] to the BOEFT and obtain the detailed expressions for the nonperturbative matching coefficients in terms of gluonic correlators. To complete the computation, it is necessary to use discrete symmetries to reduce the pNRQCD two-point functions into the structures matching the ones in the BOEFT. The values of the nonperturbative contributions are unknown, nevertheless our explicit formulas will allow a future direct lattice calculation of these objects. Alternatively, the nonperturbative matching coefficients can be obtained by comparing with lattice calculations of the charmonium hybrid spectrum and the values used to predict the spin-splittings in the bottomonium hybrid sector as shown in Ref. [19]. We provide in this paper a detailed description of the fitting procedure and enlarge the analysis to older lattice data with larger light-quark masses.
The paper is organized as follows: in Sec. II we review the discussion on the relevant scales for quarkonium hybrid systems and summarize weakly-coupled pNRQCD and the BOEFT for hybrid states. In Sec. III we demonstrate the essential calculation steps and present the results for the matching of the spin-dependent potentials, and give explicit formulas for the gluonic correlators.
In Sec. IV we compute the mass shifts in the hybrid spectrum due to the spin-dependent potentials and compare them with the charmonium hybrid spectrum obtained from two different lattice QCD calculations at different light-quark masses and fit the values of the nonperturbative matching coefficients. We use these values to give a prediction for the spin-dependent mass shifts in the bottomonium sector. We give our summary and conclusion in Sec. V. In Appendix A, using discrete symmetries, we obtain the relations between the gluonic correlators that are needed to complete the matching calculation of the spin-dependent potentials. A detailed overview of the matching of the spin-dependent terms of the two-point functions in pNRQCD and the BOEFT is given in Sec. III and Appendix B. Finally, in Appendix C and D we work out the matrix elements of the spin-dependent operators.

II. SCALES AND EFFECTIVE FIELD THEORY DESCRIPTION
In heavy quarkonium systems there are several well-separated scales typical of nonrelativistic bound states: the heavy-quark mass m (hard scale), the relative momentum between the heavy quarks mv ∼ 1/r (soft scale), where v 1 is the relative velocity and r the relative distance, and the heavy-quark binding energy mv 2 (ultrasoft scale). Additionally, we also encounter the scale of the QCD nonperturbative physics Λ QCD .
Heavy quarkonium hybrids are bound states of a heavy-quark-antiquark pair with a gluonic excitation. In quarkonium hybrids an interesting scale separation pattern appears similar to the one of diatomic molecules bound by electromagnetic interactions. The heavy quarks play the role of the nuclei and the gluons (and the light quarks) play the role of the electrons. In a diatomic molecule the electrons are non-relativistic and their energy levels can be studied in the nuclei static limit due to the latter larger mass. These electronic energy levels are called electronic static energies and are of order m e α 2 , with m e the electron mass and α the fine structure constant. The nuclei vibrational (bound) states occur around the minima of these electronic static energies and have energies smaller than m e α 2 .
In quarkonium hybrids, the light degrees of freedom are relativistic with a typical energy and momentum of order Λ QCD . This implies that the typical size of a hybrid state is of the order of 1/Λ QCD . The scaling of the typical distance of the heavy-quark-antiquark pair, r ∼ 1/(mv), depends on the details of the full inter-quark potential, which has a long-range nonperturbative part and a short-range Coulomb-like interaction. Therefore, it may happen that the heavy-quarkantiquark pair is more closely bound than the light degrees of freedom. This situation is interesting because the hybrid would present a hierarchy between the distance of the quark-antiquark pair and the typical size of the light degrees of freedom that does not exist in the case of diatomic molecules, where the electron cloud and the distance between the nuclei are of the same size. A consequence of this is that while the molecules are characterized by a cylindrical symmetry, the symmetry group for hybrids at leading order in a (multipole) expansion in the distance of the heavyquark-antiquark pair is a much stronger spherical symmetry. This modifies significantly the power counting of the EFT for hybrids with respect to the case of diatomic molecules, leading to new effects. In the following we consider this case with the interquark distance of order r 1/Λ QCD .
As in diatomic molecules, in order for a Born-Oppenheimer picture to emerge it is crucial that the binding energy of the heavy particles, mv 2 , is smaller than the energy scale of the light degrees of freedom. In summary, we will require the following hierarchy of energy scales to hold true: We can then build an EFT to describe quarkonium hybrids by sequentially integrating out the scales above mv 2 [15,17]. In this paper we focus our attention on the spin-dependent terms up to O(1/m 2 ).

A. Weakly-coupled pNRQCD
The first step in constructing the quarkonium hybrid BOEFT is to integrate out the hard scale which produces the well-known NRQCD [23][24][25]. The next step is to integrate out the soft scale, i.e., expand in small relative distances between the heavy quarks. In the short-distance regime, r 1/Λ QCD , this step can be performed in perturbation theory and one arrives at pNRQCD [21,22,26], which is the starting point of our discussion. The weakly-coupled pNRQCD Lagrangian ignoring light quarks 1 and including the gluon interaction operators from Ref. [27] that will be needed for the present work is S and O are the heavy-quark-antiquark singlet and octet fields respectively, normalized with respect to color as S = S1 c / √ N c and O = O a T a / √ T F . They should be understood as functions of t, the relative coordinates r, and the center of mass coordinate R of the heavy quarks. All the gluon fields in Eq. (1) are multipole-expanded in r and therefore evaluated at R and t: in particular the gluon field strength G µν a ≡ G µν a (R, t), and the covariant derivatives . The momentum and orbital angular momentum of the reduced mass of the heavy-quark-antiquark pair are respectively denoted by p = m dr dt = −i∇ r and L QQ = r × p. The spin vectors of the heavy quark and heavy antiquark are S 1 and S 2 respectively.
The terms with explicit factors of the chromoelectric field E and the chromomagnetic field B are obtained by matching NRQCD to weakly-coupled pNRQCD at tree level. The coefficients c F and c s are matching coefficients of NRQCD (see e.g. Ref. [25]), calculated in perturbation theory, as α s is small at the scale m that characterizes these coefficients. They are equal to 1 at leading order in α s . The ellipsis denotes other spin-independent operators, operators higher order in the multipole expansion or 1/m, and perturbative corrections of higher orders in α s . The Hamiltonian densities h s and h o of the singlet and octet fields respectively read It is useful to organize V o (r) 2 as an expansion in 1/m and separate the spin-independent (SI) and spin-dependent (SD) terms: where S = S 1 + S 2 and S 12 = 12(S 1 ·r)(S 2 ·r) − 4S 1 · S 2 . The octet-field spin-dependent potentials can be found in Ref. [28] 3 . They are given from the tree-level matching of NRQCD to weakly-coupled pNRQCD by where C F = N 2 c −1 Nc T F and C A = 2N c T F are the Casimir factors for the fundamental and adjoint representations of the color gauge group SU (N c ) respectively. We define where T a are the color generators in the fundamental representation. The renormalization scale, ν, is naturally of order mv ∼ 1/r. The matching coefficients f 8 's originate in heavy-quark-antiquark annihilation diagrams. To O(α s ) they read [24,29] At the accuracy of this work, we will use the tree-level expressions of c F and c s , i.e. c F = c s = 1, for the spin-dependent octet potentials in Eqs. (7)-(9).

B. The BOEFT
The final step is to build an EFT, which we call the BOEFT, that describes the heavy-quarkantiquark pair dynamics in the presence of a background gluonic excited state by integrating out the scale Λ QCD . First, we have to identify the degrees of freedom in the BOEFT.
In the short-interquark-distance limit r → 0 and the static limit m → ∞, quarkonium hybrids reduce to gluelumps, which are color-singlet combinations of a local static octet color source coupled to a gluonic field. The gluonic excitations can be characterized by the so-called gluelump operators [15,22]. The Hamiltonian for the gluons at leading order in the 1/m-and multipole expansions, corresponding to the G a µν G aµν term in the Lagrangian in Eq.(1), is given by 3 A contribution to V o S 2 , proportional to the f8's, which originate in quark-antiquark annihilation diagrams, is missing in Ref. [28]. Setting cF = cs = 1 and neglecting the contribution from the quark-antiquark annihilation diagrams in Eqs. (7)-(9) would recover the corresponding expressions in Ref. [28].
We define the gluelump operators, G ia κ , as the Hermitian color-octet operators that generate the eigenstates of H 0 in the presence of a local heavy-quark-antiquark octet source: where a is the color index, κ labels the quantum numbers K P C of the gluonic degrees of freedom and i labels its spin components. The spectrum of the mass eigenvalues, Λ κ , called the gluelump mass, has been computed on the lattice in Refs. [30][31][32].
At the next-to-leading order in the multipole expansion the system is no longer spherically symmetric but acquires instead a cylindrical symmetry 4 around the heavy-quark-antiquark axis.
Therefore it is convenient to work with a basis of states with good transformation properties under D ∞h . Such states can be constructed by projecting the gluelump operators on various directions with respect to the heavy-quark-antiquark axis: where summations over indices i and a are implied. P i κλ is a projector that projects the gluelump operator to an eigenstate of K ·r with eigenvalue λ, where K is the angular momentum operator for the gluonic degrees of freedom andr the unit vector along the heavy-quark-antiquark axis. It is therefore natural to define the degrees of freedom of the BOEFT as the operatorΨ κλ (r, R, t) defined by where P is the momentum operator conjugate to R, and Z is a field renormalization factor, normalized such that the following commutation relations hold: where I is the identity matrix of the spin indices of the heavy quark and antiquark. The BOEFT is obtained by integrating out modes of scale Λ QCD , i.e. the gluonic excitation. The Lagrangian of the BOEFT reads as 4 The symmetry group is changing from O(3) × C to D ∞h , with P replaced by CP .
where the trace is over spin indices of the heavy quark and antiquark, and the ellipsis stands for operators producing transitions to standard quarkonium states and transitions between hybrid states of different κ. The former are beyond the scope of this work 5 and the latter are suppressed at least by 1/Λ QCD since the static energies for different κ are separated by a gap ∼ Λ QCD . The potential V κλλ can be organized into an expansion in 1/m and a sum of spin-dependent (SD) and independent (SI) parts: In Ref. [15] the static potential V (0) κλ (r) was matched to the quark-antiquark hybrid static energies computed on the lattice. In Fig. 1 we show the QCD static energies computed using lattice NRQCD from Ref. [33]: they are plotted as a function of the quark-antiquark distance r and only states with excited glue are presented. The standard quarkonium static energy, without gluonic excitations, would lie below in energy and is not shown. Recently, a new comprehensive lattice study of the hybrid static energies has appeared in Ref. [34].
One of the major features of this spectrum is that in the short-distance region the static energies can be organized in quasi-degenerate multiplets corresponding to the gluelump spectrum. This is a direct consequence of the breaking of spherical symmetry into a cylindrical symmetry once the subleading contributions in the multipole expansion are included. Indeed, at leading order in the multipole expansion V That is, the potential in the short-distance limit only depends on the quantum numbers of the gluelump κ and not on its projection λ.
The lowest gluelump has quantum numbers κ = 1 +− . In Ref. [15] the matrix elements of r m P i κλ were obtained for κ = 1 +− and it was shown to contain off-diagonal terms in λ-λ that lead to coupled Schrödinger equations. The Schrödinger equations were solved numerically and the spectrum and wave functions of hybrid states generated the static energies labeled by Σ − u and Π u were obtained.
The lowest hybrid static energies [33] and gluelump masses [30,31] in units of r 0 ≈ 0.5 fm. The absolute values have been fixed such that the ground state Σ + g static energy (not displayed) is zero at r 0 . The data points at r = 0, labeled with κ = K P C , are the gluelump masses. The gluelump spectrum has been shifted by an arbitrary constant to adjust the 1 +− state with the Π u and Σ − u potentials at short distances, with the dashed lines indicating the expected extrapolation to degeneracy at r = 0. The behavior of the static energies at short distances becomes rather unreliable for some hybrids, especially the higher excited ones. This is largely due to the difficulty in lattice calculations to distinguish between states with the same quantum numbers, which mix. The figure is taken from [31].

III. MATCHING OF THE SPIN-DEPENDENT POTENTIALS
We present now the results of the matching for the spin-dependent potentials in Eqs. (19) and (20) for the lowest-lying gluelump (κ = 1 +− ) in the short distance regime 1/r Λ QCD . We will first write down the formulation of the matching for general κ and then focus on the case κ = 1 +− , for which we will demonstrate the essential steps of the calculation and present the final results, and leave the more involved details of calculation in Appendix B.
The matching between weakly-coupled pNRQCD and the BOEFT at the scale Λ QCD is performed by considering the following gauge-invariant two-point Green's function defined in terms of the fields in pNRQCD: where only the repeated color indices a, b and spin indices i, j are summed. In the BOEFT, with Eq. (14), the two-point Green's function is given by The Green's function in pNRQCD (Eq. (22)) has the form where and the gluelump mass Λ κ is related to a gluonic correlator by with φ ab (t, t ) the adjoint static Wilson line defined by In Eq. (24), δV κλλ is obtained from the contributions to Eq. (22) from insertions of singlet-octetand octet-octet-gluon coupling operators from the Lagrangian in Eq. (1). Comparing Eqs. (24) and (23), we obtain Z κ (r, R, p, P ) = 1 and which reduces to Eq. (21) at leading order in 1/m and the multipole expansion. The matching condition is schematically depicted in Fig. 2. In Fig. 2, the left-hand side and the right-hand side correspond to the two-point Green's function computed in the BOEFT and pNRQCD respectively.
Diagram (a) gives the perturbative term P i † κλ V o P i κλ in Eq. (28), which is inherited from the octet potential in Eq. (4), as well as a nonperturbative term Λ κ , the gluelump mass. Diagrams like (b), (c), (d), (e), (f), and (g), with black dots, which denote the singlet-octet-and octet-octet-gluon coupling operators in the pNRQCD Lagrangian Eq. (1), give another nonperturbative contribution δV κλλ . All diagrams in pNRQCD are computed in coordinate space, similar to what is done in Ref. [22].
Now we focus on the case κ = 1 +− . To simplify the notation we will drop the subscript κ for the rest of the manuscript and it should be understood that we are always referring to κ = 1 +− , and so λ takes the values 0, ±1. The spin-dependent potentials in Eqs. (19) and (20) for κ = 1 +− read as follows: where K ij k = i ikj is the angular momentum operator for the spin-1 gluonic excitation. The projectors P i λ read withr = (sin(θ) cos(φ), sin(θ) sin(φ) , cos(θ)) , The 1/m operators in Eq. (29), with coefficients V SK (r) and V SKb (r) 6  The coefficients V i (r) on the right-hand side of Eqs. (29) and (30) have the form V i (r) = is the perturbative octet potential and V np i (r) is the nonperturbative contribution. From the multipole expansion, V np i (r) is a power series in We will work at next-to-leading order in the multipole expansion for the 1/mpotentials and leading order in the multipole expansion for the 1/m 2 -potentials. Therefore, up to the precision we work at, we have 6 The operator with coefficient V SKb (r) contains the tensor r i r j contracted to other vectors. In the case of standard quarkonia, for which the symmetry group is SO(3), it is natural to decompose the tensor r i r j into a sum of a trace part and a traceless symmetric part, each of which corresponds to an irreducible representation of SO(3).
This is done for the operators with coefficients V S 2 (r) and VS 12 a(r), since they also appear in the case of standard quarkonia. In the case of hybrid states here, since the symmetry group is D ∞h instead of SO (3), this decomposition is not of particular relevance and we decided to write the V SKb (r)-operator without substracting the trace part.
In Eqs. (36), (39), and (40), V o SL (r), V o S 2 (r), and V o S 12 (r) are the perturbative tree-level spindependent octet potentials given by Eqs. (7)- (9). The constants V and (g) in Fig. 2 with insertions of spin-dependent operators with a chromomagnetic field or a chromoelectric field in the pNRQCD Lagrangian Eq. (1) , and are expressed as nonperturbative purely gluonic correlators. It should be emphasized that the expressions of V SK (r), V SKb (r), Consider diagram (b) in Fig. 2, with an insertion of the c F -term in Eq. (1). Its contribution to δV λλ is given by where and The color combination being a rotationally invariant tensor, can be written as Therefore, Eq. (44) becomes from which it follows that The detailed derivations of V The correlator (U ss EE ) ijkl in Eq. (49) arises from insertions of two singlet-octet vertices with a chromoelectric field from the pNRQCD Lagrangian Eq. (1). The adjoint Wilson lines connecting the gluelump operators to the chromoelectric fields arise from the two propagators of the octet field in diagram (c) of Fig. 2. Similarly, (U ss BB ) ijkl in Eq. (50) is defined like (U ss EE ) ijkl with the chromoelectric field replaced by the chromomagnetic field. The relevant gluonic correlators that correspond to diagram (d) in Fig. 2 are The correlator (U oo EE ) ijkl bcdef g in Eq. (51) arises from insertions of two octet-octet vertices with a chromoelectric field from the pNRQCD Lagrangian. The three adjoint Wilson lines arise from the three propagators of the octet field in diagram (d) of Fig. 2. The correlator (U oo BDE ) ijklm bcdef g in Eq. (53) arises from insertions of two octet-octet vertices, one with a chromomagnetic field at time t and another with a covariant derivative of the chromoelectric field at time t < t. (U oo BB ) ijkl bcdef g in Eq. (52) and (U oo DEB ) ijklm bcdef g in Eq. (54) are similarly defined. The relevant gluonic correlators that correspond to diagram (e) in Fig. 2 are The relevant gluonic correlators that correspond to diagram (f) in Fig. 2 are The relevant gluonic correlators that correspond to diagram (g) in Fig. 2 are In Eqs. (55) to (63), the correlators arise from insertions of three vertices, each being a singletoctet or an octet-octet vertex, with a chromoelectric field or a chromomagnetic field. Analogous to Eq. (45), we define the color combinations where h abc is as defined below Eq. (43), d abc ≡ 1 2 (h abc + h acb ) and f abc = − i 2 (h abc − h acb ). The tensors defined in Eqs. (64)-(71) have the formÛ ijkl orÛ ijklm , which being rotationally invariant, have the tensor decompositions given bŷ The nonperturbative coefficients V are then given by V np (0) In the derivation of Eqs. which we know the dependence on the heavy-quark flavor, and a nonperturbative purely gluonic correlator, which is independent of the heavy-quark flavor.

IV. SPIN SPLITTINGS IN THE HYBRID SPECTRA
We obtain the spin-dependent contributions to the quarkonium hybrid spectrum by applying time-independent perturbation theory to the spin-dependent potentials in Eqs. (29)- (30). We carry out perturbation theory to second order for the terms V in Eqs. (29) and (34) The zeroth-order wave functions are obtained following the procedure described in Ref. [15], by solving the coupled Schrödinger equations involving the potentials V where the components from top to bottom correspond to λ = 0, +1, −1. We define L = L QQ + K,  Table I. Lowest-lying quarkonium hybrid multiplets Multiplet l J P C (s = 0) J P C (s = 1) The angular wave functions v λ l m l are eigenfunctions of L 2 and not of L 2 QQ . As a result, the evaluation of matrix elements of operators involving L QQ is not totally straightforward. The details of the calculation of these matrix elements can be found in Appendix C. We will present the results for the four lowest-lying spin-multiplets shown in Table I. Matrix elements of the spin-dependent operators in Eqs. (29) and (30) for the angular part of the wave functions of the states in Table I are listed in Appendix D. The eight nonperturbative parametersṼ that appear in the spin-dependent potentials Eqs. (34)-(41) are obtained by fitting the spin-splittings to corresponding splittings from the lattice determinations of the charmonium hybrid spectrum. Two sets of lattice data from the Hadron Spectrum Collaboration have been used, one set from Ref. [36] with a pion mass of m π ≈ 400 MeV and a more recent set from Ref. [37] with a pion mass of m π ≈ 240 MeV. We take the values m RS c (1GeV) = 1.477 GeV [38] and α s at 4-loops with three light flavors, α s (2.6 GeV) = 0.26. In the fit the lattice data is weighed by (∆ 2 lattice + ∆ 2 high-order ) −1/2 , where ∆ lattice is the uncertainty of the lattice data and ∆ high-order = (m lattice − m lattice spin-average ) × Λ QCD /m is the estimated error due to higher-order terms in the potential. The    Figure 3, except using the lattice results from Ref. [37] with m π ≈ 240 MeV. decreases with J. This trend is opposite to that of the lattice results. This discrepancy can be reconciled thanks to the nonperturbative contributions, in particular due to the contribution from V np (0) SK , which is only suppressed by 1/m, and has no perturbative counterpart. A consequence of the countervail of the perturbative contribution is a relatively large uncertainty on the full result with respect its absolute value caused by a large nonperturbative contribution. Due to this uncertainty the mass hierarchies among the spin-triplet states of the multiplets H 2 and H 4 are not firmly determined. This is reflected on the change of the mass hierarchies for the central values of the lattice data from Ref. [36] to Ref. [37].
All the dependence on the heavy-quark mass of the V np (j) i 's in Eqs. (48) and (74)-(81) is encoded in the NRQCD matching coefficients c F and c s . At leading order in α s these coefficients are known to be equal to 1 and the dependence on the heavy-quark mass only appears when the next-to-leading order is considered [25]. Hence, at the order we are working, only the heavy-quark Table II. Nonperturbative matching coefficients determined by fitting charmonium hybrid spectrum obtained from the hybrid BOEFT to the lattice spectrum from the Hadron Spectrum Collaboration data of Refs. [36] and [37] with pion masses of m π ≈ 400 MeV and m π ≈ 240 MeV respectively. The matching coefficients are normalized to their parametric natural size. We take the value Λ QCD = 0.5 GeV.
Ref. [36] Ref. [37] V with the renormalization scale set as the heavy-quark mass. Taking this mass dependence into account, we can use the set of nonperturbative parameters to predict the spin contributions in the bottomonium hybrid sector, for which lattice determinations are yet not available due to their larger difficulty compared to the charm sector.
We compute the bottomonium hybrids spectrum by adding the spin-dependent contributions from Eqs. (34)- (41) to the spectrum obtained in Ref. [15]. We show the results thus obtained in Figs. 5 and 6 for the values in the second and third columns of Table II respectively. We use the value of the bottom mass m RS b (1 GeV) = 4.863GeV.

V. CONCLUSIONS
The spin-dependent operators for heavy quarkonium hybrids up to order 1/m 2 were presented in Ref. [19]. The most prominent feature is the appearance of two spin-dependent operators already at order 1/m, unlike standard quarkonia, in which case the spin-dependent operators appear at order 1/m 2 . These operators, in Eq. (29), couple the total spin of the heavy-quark-antiquark pair with the spin of the gluonic degrees of freedom that generates the hybrid state. At order 1/m 2 , we have the spin-orbit, total spin squared, and tensor spin operators familiar in the studies of standard quarkonia. In addition, three new operators appear at order 1/m 2 , which involve the 10  In the short heavy-quark-antiquark distance regime, r 1/Λ QCD , the matching coefficients, i.e.  Table II. We found that it is possible to reproduce the lattice data of the charmonium hybrid spectrum with nonperturbative matching coefficients of natural sizes. The values of the pion mass utilized in Refs. [36] and [37] are m π ≈ 400 MeV and m π ≈ 240 MeV respectively. The variation of the values of the nonperturbative matching coefficients obtained by the fits for the two lattice data sets can be tentatively attributed to the light-quark mass dependence of the gluon correlators, in particular for the matching coefficient of the leading spin-dependent operator V np (0) SK . Finally, we have taken advantage of the fact that the gluonic correlators are independent of the heavy-quark flavor to compute the mass spin-splittings for the bottomonium hybrid spectrum. The results are shown in Figs. 6 and 5. The bottomonium hybrid spectrum including spin-dependent contributions has not yet been computed on the lattice 8 . Calculations of the bottomonium hybrid spectrum on the lattice are difficult due to the widely-separated scales of the system, i.e. the bottom-quark mass being much larger than Λ QCD . Therefore, precision calculations of the bottomonium hybrid spectrum on the lattice would require both large volume and small lattice spacing, which is computationally challenging. On the other hand, an EFT approach can take advantage of the same wide separation of scales and is, like lattice QCD, a model-independent approach rooted in QCD. Combining both approaches opens a promising path towards the understanding of exotic quarkonia.
The impact of this calculation is manifold. First, as observed above, the spin dependence of the operators for quarkonium hybrids is significantly different from that for standard quarkonia. This has an important impact on the phenomenological calculation. Second, we have obtained for the first time the expressions of the nonperturbative contributions to the spin-dependent potential for quarkonium hybrids in terms of gauge-invariant correlators depending only on the gluonic degrees of freedom, which are suitable for computation on the lattice or evaluation in QCD vacuum models.
The technology to calculate these correlators in lattice QCD already exists [34,[40][41][42][43] and could be readily applied. Third, we emphasize that the obtained nonperturbative correlators depend only on the gluonic degrees of freedom and not on the heavy-quark flavor. This allows us to extract the unknown nonperturbative parameters in the spin-dependent potential from the charmonium hybrid spectrum and use them for bottomonium hybrids. Finally, since the BOEFT can be generalized to considering also light quarks as the light degrees of freedom [18], the spin dependent operators will likely have similar characteristics also in that case. Therefore, we supply the full list of matrix elements of the spin-dependent operators in Appendix D to facilitate applications to the spectrum of XYZ states for phenomenologists and model builders. Since most of the phenomenological applications for XYZ states up to date either do not contain such spin-dependent terms or construct them in a way inspired by the traditional quarkonium case, we believe that this result can prove to be very useful.
The next step forward in the BOEFT framework will be to release the assumption mv Λ QCD we used in this paper, and work out the spin-dependent corrections using only the hierarchy Λ QCD mv 2 underlying the Born-Oppenheimer approximation [20]. In this case the spin decomposition of the potential will be the same as obtained here but the actual form of the r-dependent potentials will be given in terms of generalized Wilson loops. This is analogous to the computation of the spin-dependent potential for traditional quarkonia as generalized Wilson loops in strongly-coupled pNRQCD [44][45][46] which was later used by lattice groups to obtain the form of the nonperturbative spin-dependent potentials [47][48][49] and can be addressed with the technology developed in [9,33,34].
Its contribution to δV λλ is given by which with Eqs. (A28) and (A29) is simplified to where (Û EE ) ijkl is defined in Eq. (64).
With the tensor decomposition Eq. (72) and Eq. (A58) and using the commutation relation Consider diagrams (c) and (d) in Fig. 2, with insertion of two c F -vertices. Its contribution to δV λλ is given by which with Eqs. (A30), (A31), and (66) is simplified to where Similar to Eqs. (72) and (A58), rotational invariance and parity imply that (Û BB c ) ijkl , also has tensor decomposition of the formÛ ijkl =Ũ I (δ ij δ kl + δ ik δ jl ) +Ũ III δ il δ jk . In Eq. (B9), the terms S j 1 S k 1 and S j 2 S k 2 can be rewritten using σ i σ j = i ijk σ k + δ ij , the first term of which gives zero when contracted with (Û BB c ) ijkl , since (Û BB c ) ijkl = (Û BB c ) ikjl . Therefore, after applying the tensor decomposition of (Û BB c ) ijkl and (Û BB b ) ijkl , the spin-dependent terms in Eq. (B9) are given by which gives Eqs. (80) and (81).
Consider diagram (d) in Fig. 2 with insertion of one c F -vertex and one r i r j D i E j -vertex. Its contribution to δV λλ is given by which with the help of Eqs. (A32), (A33), (67), and (68) becomes Consider diagrams (e), (f), and (g) in Fig. 2, with insertion of one c F -vertex and two r · E-vertices.
Its contribution to δV λλ is given by where each term is the sum of the three possible diagrams with the insertion of c F -vertex in a different location (δV g ) BEE = − ic F 8mr i † λr m λ r k r l (h bcd S j 1 − h bdc S j 2 )(U ooo BEE ) ijklm bcdef ghpq d ef g d hpq , We can eliminate the last term on the right-hand side of Eq. (B32) using the relation Therefore, we have where V L i QQ ,θ j = iφ irj 0 + i cot(θ)θ iφj , (C3) from which one can obtain L i QQ ,r j ± = ± r j 0r i ± + cot(θ)θ irj ± .
Next, we consider the V SLb -operator, r †i λ L i QQ S l + S i L l QQ r l λ .