Hyperfine splittings of heavy quarkonium hybrids

In the framework of the Born-Oppenheimer Effective Field Theory, the hyperfine structure of heavy quarkonium hybrids at leading order in the 1/m Q expansion is determined by two potentials. We estimate those potentials by interpolating between the known short distance behavior and the long distance behavior calculated in the QCD Effective String Theory. The long distance behavior depends, at leading order, on two parameters which can be obtained from the long distance behavior of the heavy quarkonium potentials (up to sign ambiguities). The short distance behavior depends, at leading order, on two extra paramentes, which are obtained from a lattice calculation of the lower lying charmonium hybrid multiplets. This allows us to predict the hyperfine splitting both of bottomonium hybrids and of higher multiplets of charmonium hybrids. We carry out a careful error analysis and compare with other approaches.


I. INTRODUCTION
Exotic hadrons (those beyond mesons and baryons) have been a matter of research since the early days of QCD [1]. Among them, the so called hybrids are the closests to standard hadrons, as the only difference resides in their non-trivial gluon content. Hence their flavor structure is the same as in standard hadrons, but their J P C quantum numbers may differ since, in general, the non-trivial gluon content contributes to them. Nevertheless, the experimental confirmation of hybrids is an arduous task. On the one hand, hadrons with exotic quantum numbers, which could be associated to hybrid states, are difficult to produce with current beams. On the other hand, for light hadrons, hybrids with standard quantum numbers may mix with standard hadrons in an arbitrary way. However, when heavy quarks are involved the mixing is suppressed by inverse powers of the heavy quark mass, and hence the identification of hybrids with standard quantum numbers should become simpler, provided we have reliable theoretical predictions for the hybrid spectrum.
An economical approach to calculating the hybrid spectrum is the so-called Born-Oppenheimer effective field theory (BOEFT) [2][3][4][5][6][7]. It exploits the fact that heavy quarks move slowly in heavy hadrons. The effect of the nontrivial gluon (or/and light quark) content is encoded in a series of potentials organized in an 1/m Q expansion, m Q being the heavy quark mass. The leading potential (O(1/m 0 Q )) is heavy quark spin and heavy quark mass independent. It has been used to calculate the spin average spectrum [2][3][4]8], decays to heavy quarkonium [4,[9][10][11], and transitions between heavy quarkonium states [12]. The mixing of heavy quarkonium hybrids with heavy quarkonium starts at order 1/m Q [4]. Spin effects also start at order 1/m Q [4,5], and hence, they are more important than in heavy quarkonium, in which they start at order 1/m 2 Q . In that respect, it is important to have spin effects under good control in order to properly identify possible experimental candidates (see [13,14] for recent reviews).
We calculate here the hyperfine splittings (HFS) for the lower lying charmonium and bottomonium hybrids at leading-order (LO) (O(1/m Q )) in the BOEFT. At this order, the HFS depend on two unknown potentials [5,6]. At short distances the form of these potentials is constrained by the multipole expansion and has been given in ref. [6,15], where one can also find the form of the relevant next-to-leading order (NLO) potentials (O(1/m 2 Q )). The HFS was calculated in that reference using the short distance form of the potentials only. At long distances, the form of the potentials can be estimated using the QCD effective string theory (EST) [16,17]. This has been carried out for heavy quarkonium [18,19] and for the hybrid-quarkonium mixing terms [4]. We provide the results here for the spin-dependent terms of the lower lying static hybrid states (Σ u and Π u ), which turn out to be parameter free (up to signs). We emphasize that the typical distance between heavy quarks in heavy quarkonium hybrids states is of O(1/Λ QCD ), and hence an interpolation between the short and long distance forms of the spin-dependent potentials should provide more reliable estimates than sticking to the short distance form only. We propose a simple interpolation and calculate the HFS with it. We use the charmonium spectrum of ref. [20] to fix the unknown parameters in the short distance form of the potential, and to estimate the interpolation dependence. Then we can predict the HFS of higher multiplets and the HFS of bottomonium hybrids.
We distribute the paper as follows. In Section II, we work out the structure of the spin-dependent terms in a convenient basis. We specify the form of the spindependent potentials at short distances in III A, at long distances in III B, and the interpolation we use in III C. In Section IV, we fix our free parameters using hybrid charmonium lattice data for the lower multiplets and obtain the HFS of higher multiplets. In Section V, we obtain the HFS of bottomonium hybrids. Section VI is devoted to the comparison with other approaches and we close with a discussion and our conclusions in Section VII.

II. THE HYPERFINE SPLITTINGS
General expressions for the BOEFT at NLO have been recently obtained in ref. [7]. The lower lying hybrid potentials correspond to the κ p = 1 + quantum numbers of the light degrees of freedom (LDOF), where κ is the total angular momentum and p the parity. We focus on the heavy quark spin-dependent terms in the Hamiltonian V (1) where S QQ and S κ are the spin operators of the heavy quark-antiquark pair and the total angular momentum operator of the LDOF respectively. P κΛ are the projectors to the irreducible representations of the D h∞ group Λ = 0, 1. They read, In our case, the above potentials act on fields H nm 1 (r, t), n, m = 1, 2, 3, where the first index corresponds to the total angular momentum of the LDOF and the second to the spin of the heavy quarkantiquark pair. In the Cartesian basis, S i 1 jk = −iϵ ijk , we then have, P n ′ n 10 =r n ′r n , P n ′ n 11 = δ n ′ n −r n ′r n Analogously, the heavy quark-antiquark spin opera- We find that only two independent potentials survive.
where we have used that time reversal and Hermiticity imply V sb 1 + 01 (r) = V sb 1 + 10 (r) ∈ R, and we have defined, The structure leading to V hf was first noticed in [4] and the one leading to V hf 2 was already considered in the short distance analysis of Refs. [6,21]. We write the field H nm 1 (r, t) in the |SLJJ M > basis, where S and L are the spin and the orbital angular momentum of the QQ pair respectively, J the total angular momentum of the LDOF plus the orbital angular momentum of the QQ, and J and M the total angular momentum and its third component respectively [4] . For a given J , J = J , J ± 1, L = J, J ± 1 and only S = 1 is affected by (1).
The five dimensional box corresponds to the subspace spanned by where we use the short hand notation ± for J ± 1 or J ± 1 and 0 for J or J . For the terms proportional to −2V hf it reads [5,22], and for terms proportional to −2V hf 2 [23], with, The four dimensional box corresponds to the subspace spanned by ( . For the terms proportional to −2V hf it reads [5,22], 1J M and P 0− 1J M do not exist and the matrices are 7 × 7.
If J = 0, 1J M do not exist and the system is reduced to 3 × 3 matrices for both potentials.
At first order in perturbation theory, only the diagonal terms and the off-diagonal terms corresponding to L = J ± 1, L ′ = J ∓ 1 matter. This is because the 0-th order wave functions have a single component for L = J and two components for L = J ± 1. Then the following formulas for the masses M SJ of the spin multiplet J, which are independent of the potentials V hf and V hf 2 , hold, (12) They were initially derived in ref. [5,22] for V hf 2 = 0, but it is not difficult to see that they are also fulfilled for V hf 2 ̸ = 0 [23].

III. THE POTENTIALS
The potentials V sa 1 + 11 (r) and V sb 1 + 10 (r) can be obtained from suitable operator insertions in the expectation values of Wilson loops [7], and hence, they are amenable to lattice evaluations. However, no such evaluation exists to date. The short distance behavior can be worked out with the help of pNRQCD [24,25]. It has been obtained at NLO in the 1/m Q expansion in Refs. [6,15]. However, the typical distances at which the static hybrid potentials reach the minimum [26][27][28][29], the typical quark-antiquark distance in the bound states [3], and the shape of the wave functions (see plots in Section VI of [4]), indicate that most of the time the quark-antiquark separation is of the order of 1/Λ QCD . It may then be more appropriated incorporating reliable long distance information rather than calculating higher orders at short distances. We shall do that by using the Effective String Theory of QCD (EST) [16,30], which describes well the long distance behavior of the static hybrid potentials [26], as well as those of the 1/m Q and 1/m 2 Q potentials for quarkonium [4,18], calculated in lattice SU (3) pure Yang-Mills theory. We will then stick to the Cornell model's philosophy of using potentials that interpolate between known short and long distance behavior.

A. The short distance behavior
In order to estimate the short distance behavior we use the fact that the 1/m Q potentials are analytic in r in pNRQCD. This implies that, We shall keep the LO terms only. A and B are unknown real constants, can be related to expectation values of operator insertions in Wilson lines [15]. c F = c F (m Q ) is a short distance matching coefficient inherited from NRQCD [31][32][33]. We shall take the next-to-leading logarithmic expres- The corrections to the leading short distance behavior are suppressed by powers of Λ 2 QCD r 2 . We display the potentials (13) in Fig. 1.

B. The long distance behavior
The long distance behavior can be estimated using the QCD effective string theory [16,17], following the mapping given in ref. [18]. It has been obtained in [34], see Appendix A: The parameters gΛ ′ ∼ Λ QCD and gΛ ′′′ ∼ Λ QCD also appear in the spin-dependent potentials for heavy quarkonium. They have been obtained in [4] from the lattice data of ref. [35,36]. Using simple interpolations with the right short and long distance behavior a good fit to data is obtained with the following outcome 1 .
If only long distance data points and the long distance form of the potentials are used, the values of |gΛ ′ | and |gΛ ′′′ | are 40% and 35% larger respectively. κ ≃ 0.187 GeV 2 is the string tension and c F is the same NRQCD matching coefficient that appears in the short distance behavior. The corrections to the leading long distance behavior of the potentials in (14) are suppressed by factors 1/r 2 Λ 2 QCD . We display the potentials (5) at long distances in Fig. 1.

C. The interpolating potentials
We use for the hyperfine potentials simple interpolations with the right short and long distance behavior obtained in Sections III A and III B respectively, r 0 ∼ 1/Λ QCD is the matching scale. It is estimated from the short and long distance behavior of the static hybrid potentials to be r 0 ≃ 3.96 GeV −1 . This figure will be eventually moved in order to estimate the error due to the interpolation dependence. We show the interpolated potentials in Fig. 1.

IV. CHARMONIUM HYBRIDS: FIXING THE SHORT DISTANCE PARAMETERS
The short distance potentials depend on two arbitrary parameters at LO, A and B. We shall fix those parameters by comparing our results to the lattice data of ref. [20] for the lower lying hybrid states (the 1(s/d) 1 (H 1 ), 1p 1 (H 2 ), 1(p/f ) 2 (H 4 ) and 1p 0 (H 3 ) multiplets). The spin average of the multiplets in ref. [20] were higher than in ref. [4]. Since we are going to use the same methodology as in the last reference, we correct the lattice data by the spin av- for the ones with the lowest χ 2 /dof . We adapted the code used in ref. [4] by adding the spin dependent potentials above, and neglecting, for simplicity, the mixing with quarkonium [37]. The charm mass is also taken as in ref. [4], m c = 1.47 GeV. If we neglect the long distance behavior, we find A = 0.0699 GeV, B = 0.0008 GeV 3 with χ 2 /dof = 1.1927. When we include the long distance behavior the fit quality improves considerably, we obtain as a best fit A = 0.1356 GeV, B = −0.0022 GeV 3 with χ 2 /dof = 0, 6439. gΛ ′ and gΛ ′′′ are taken as in (15), r 0 = 3.96 GeV −1 in (16) and all possible sign combinations for the long distance potentials in (14) and (15) are considered. The best fit corresponds to a negative V sa ld (r) and a positive V sb ld (r). The former implies gΛ ′′′ < 0. Reversing the sign of V sa ld (r) (V sb ld (r)) worsens the fit considerably (marginally), see Table I.
We have also explored the dependence of the result on gΛ ′ and gΛ ′′′ according to the possible values given in ref. [4]. The fit has a mild preference for larger values of |gΛ ′ | and |gΛ ′′′ | (χ 2 /dof = 0, 632 versus χ 2 /dof = 0, 643) , so we shall take for now on gΛ ′ = −0.0796 GeV and gΛ ′′′ = 0.3105 GeV, which corresponds to the fit to long distance data only in ref. [4]. These values lead to A = 0.1445 GeV and B = −0.0036 GeV 3 . The change in the spectrum is negligibly small (≲ 1 MeV). The interpolation dependence is estimated by moving r 0 ∈ [2.5 , 5] GeV −1 . The χ 2 /dof marginally improves around r 0 ∼ 3.5 GeV −1 , and considerably deteriorates for r 0 ≤ 3 GeV −1 and r 0 > 3.96 GeV −1 . We shall stick to the default value r 0 = 3.96 GeV −1 and use r 0 ∼ 3.5 GeV −1 to estimate the error due to the interpolation. The changes in the spectrum are of 3.1 MeV in average, with a maximum of 7 MeV. A increases about a 30% and B becomes more than twice its value. See Table  II.
In order to establish the error of A and B due to the input data, we assume a linear dependence of the binding energy on them, which holds at first order in perturbation theory. We obtain A = 0.115 ± 0.034  (17). The short distance parameters A and B correspond to those in Table III  GeV and B = 0.0038 ± 0.0154 GeV 3 . Notice that the value of B is compatible with zero. The spectrum is displayed in Table III. The error due to the uncertainty in gΛ ′ and gΛ ′′′ is negligible. We also neglect the error due to the interpolation, which is not negligible in order to obtain a value for A and B, but it is for the spectrum due to its correlation with r0. We include the error due to higher orders in the 1/m Q expansion, which is about 30 MeV for charm. In ref. [4], it was found that two extra multiplets lie below the 1p 0 , the 2(s/d) 2 ) and the d 2 . For complete-  ness we also display the spectrum of these multiplets including the hyperfine splitting in Table IV  nium, A ′ and B ′ also are, We calculate the spectrum for the central values of these parameters (A ′ = 0.02704 GeV, B ′ = 0.00091 GeV 3 ), which provides us with the central value of the masses, and for the four corners of their 1 σ range, which allows us to estimate the error due to the fit parameters. We take it as the larger difference in either sense. The total error is obtained by adding in quadrature to the latter the error associated to higher orders ∼ Λ 3 QCD /m 2 Q ∼ 3 MeV. The bottom mass is taken as in ref. [4], m b = 4.88 GeV. We display the results in Table V. We have also explored the dependence on r 0 . We have recalculated the central values for r 0 = 3.5 GeV −1 and the corresponding parameters (obtained from charm and converted to bottom according to (17), A ′ = 0.0423 GeV, B ′ = −0.0021 GeV 3 ). Recall that this r 0 choice provides a slightly better fit to charmonium lattice data than our default r 0 = 3.96 GeV −1 . The difference is below 1 MeV for the (s/d) 1 multiplet and around 0.1 MeV for the remaining ones. We can then safely neglect it.

VI. COMPARISON WITH OTHER APPROACHES
We shall compare here our results with those also obtained in BOEFT in Refs. [6,15] and with the lattice QCD results of the HSC collaboration [20,38] (we will not use earlier results at larger pion masses [39]). Recall that the main differences with respect to Refs. [6,15] is that only short distance expressions for the spin-dependent potentials at NLO (8 non-perturbative parameters to fit) are used there, whereas we use LO short distance expressions and LO long distance expressions (2 non-perturbative parameters to fit).
Lattice results for charmonium hybrids [20] are used to fit the non-perturbative parameters. Hence, our errors are slightly larger than the ones from lattice data. Our fits have lower χ 2 /dof than those in Refs. [6,15], which supports the inclusion of long distance potentials. The overall picture for the charmonium hyperfine splittings is similar to the one obtained in those references, see Fig. 2. For spin one states we get similar, slightly larger, smaller and larger errors for the (s/d) 1 (H 1 ), p 1 (H 3 ), (p/f ) 2 (H 4 ) and p 0 (H 3 ) multiplets respectively. For spin zero states, which do not enter in our analysis, the errors are only due to neglected higher orders in the 1/m Q expansion. Note that the error that we assign to them is larger than than the one assigned in Refs. [6,15]. See Fig. 2.
Our predictions for the bottomonium hybrid hyperfine splittings are compatible with the few available states from lattice QCD [38], but we have much smaller errors, see Fig. 3. We get again an overall picture similar to the one of ref. [6,15], with smaller errors for spin one states, and smaller hyperfine splittings in general. However, for the (s/d) 1 multiplet we have about 2σ discrepancies with those references. Nevertheless, both sets of splittings are consistent with lattice data of [38] for this multiplet. The remaining available lattice states do not form any obvious spin multiplet. We have assigned the three lighter (two heavier) ones to the p 1 ((p/f ) 2 ) multiplet. However, the fact that the 1 ++ state is lighter than the 1 −− in the lattice results is in conflict with the BOEFT ones.

VII. DISCUSSION AND CONCLUSIONS
The static potentials we use are taken from [4]. They were obtained from pure SU (3) (n f = 0) lattice data of Refs. [26,27]. More recent pure SU (3) (n f = 0) lattice data exist for the static hybrid potentials, with better resolution at short and intermediate distances [28,29,40]. It would be interesting to incorporate it in future analysis. However, the systematic errors due to using pure SU (3) (n f = 0) rather than QCD with dynamical light quarks (n f = 3) would still be difficult to evaluate. The early study of ref. [41] suggest that they are small. The strategy recently presented in [42] may help to resolve this issue.
The lattice data to which we fit our spin-dependent expressions, ref. [20], is obtained from a 2 + 1 + 1 clover action at a fixed spacial lattice spacing of 0.12 fm, a 3.5 times smaller temporal lattice spacing, and light quark masses corresponding to m π ∼ 240 MeV. In ref. [43] results in the continuum limit for realistic light quark masses are obtained using a HISQ action with 2 + 1 + 1 dynamical quarks, but only for two states in the lowest lying multiplet. The lattice bottomonium data we compare with, ref. [38], uses the same setting as [20], tuning the heavy quark parameters to bottomonium observables, and uses light quark masses corresponding to m π ∼ 391 MeV.
We have focused on the hyperfine splittings. The absolute values of the masses we quote correspond to the spin averages of Ref. [4]. The central values quoted in that reference are lower than those in Ref. [3] but compatible within errors. The differences are due to the choice of normalization (quarkonium spectrum versus RS mass scheme). They are also much lower than the lattice results of Ref. [20], a difference that shrinks if they are compared with earlier lattice data at larger pion mass [39]. This suggest that part of the discrepancy may be due to the quenched lattice data used as an input in Refs. [3,4]. They are also lower than in most models (see [44] for a review). Since the discrepancies usually amount to global shifts, they are not expected to affect the bulk of the hyperfine splitting analysis presented here. However, a small dependence on the input lattice data (unphysical) pion mass was noticed in the short distance analysis of Ref. [15], which may be present in our results as well.
We have shown that the inclusion of long distance contributions calculated in the QCD EST to the LO spin dependent potentials in the BOEFT considerable improves the description of the charmonium hybrids hyperfine splittings obtained in lattice QCD [20]. For the LO spin-dependent potentials the χ 2 /dof moves from 1.193 for short distance contributions only to 0.644 for short and long distance contributions. This figure is much lower than the χ 2 /dof = 0.999 obtained in ref. [15] for NLO spin-dependent potential with short distance contributions only. The fact that long distance contributions are important may be anticipated from the results on the size of charmonium hybrids displayed in Table III of ref. [3], ⟨1/r⟩ ∈ [190, 420] MeV, scales of the order of Λ QCD . Using the QCD EST to describe them has the remarkable feature that it introduces no new unknown parameter, beyond sign ambiguities and the scale r 0 that separates short and long distances in the interpolation. Hence we have two parameter fits to data, rather than the eight parameter fits of ref. [ 2. The spectrum of the lower-lying n (s/d)1 (H1), n p1 (H3), n (p/f )2 (H4) and n p0 (H3) charmonium hybrids computed by adding the LO spin-dependent potentials to the static potentials used in ref. [4] is shown in green boxes. The average mass for each multiplet is shown as a red dashed line. Blue and magenta boxes show the results of the BOEFT with NLO spin-dependent short distance potentials only [6,15] and lattice QCD [20] respectively, adjusted to our spin average. The height of the boxes indicates the uncertainty.
ing to a smaller χ 2 /dof as mentioned above.
Once we have the unknown parameters fixed, we can calculate the hyperfine splittings of higher charmonium hybrid states, of the bottomonium ones, and the error associated to them. This is displayed in Tables III, IV and V and in Figs. 2 and 3. For charmonium hybrids, we get results compatible with ref. [15] with similar errors overall, and, as expected, compatible with ref. [20], the source of our fit, with slightly larger errors. For bottomonium hybrids our hyperfine splittings are compatible with those of Ref. [38], with much smaller errors, but smaller than those in Ref. [15], with similar errors.    3. The spectrum of the lower-lying n (s/d)1 (H1), n p1 (H3), n (p/f )2 (H4) and n p0 (H3) bottommonium hybrids computed by adding the LO spin-dependent potentials to the static potentials used in ref. [4] is shown in green boxes. The average mass for each multiplet is shown as a red dashed line. Blue and magenta boxes show the results of the BOEFT with NLO spin-dependent short distance potentials only [6,15] and lattice QCD [38] respectively, adjusted to our spin average. The height of the boxes indicates the uncertainty.