Effective QCD string and doubly heavy baryons

Expressions for the potentials appearing in the nonrelativistic effective field theory description of doubly heavy baryons are known in terms of operator insertions in the Wilson loop. However, their evaluation requires nonperturbative techniques, such as lattice QCD, and the relevant calculations are often not available. We propose a parametrization of these potentials with a minimal model dependence based on an interpolation of short- and long-distance descriptions. The short-distance description is obtained from weakly-coupled pNRQCD and the long-distance one is computed using an effective string theory. The effective string theory coincides with the one for pure gluodynamics with the addition of a fermion field constrained to move on the string. We compute the hyperfine contributions to the doubly heavy baryon spectrum. The unknown parameters are obtained from heavy quark-diquark symmetry or fitted to the available lattice QCD determinations of the hyperfine splittings. Using these parameters we compute the double charm and bottom baryon spectrum including the hyperfine contributions. We compare our results with those of other approaches and find that our results are closer to lattice QCD determinations, in particular for the excited states. Furthermore, we compute the vacuum energy in the effective string theory and show that the fermion field contribution produces the running of the string tension and a change of sign in the L\"uscher term.


I. INTRODUCTION
leading-order parameters can be extracted from the heavy-light meson spectrum using heavy quark-diquark symmetry. In Sec. III, we propose an EST with fermionic degrees of freedom in order to describe the long-distance behavior of the potentials. Based on the D ∞h group, we put forward a mapping from the NRQCD operator insertions in the Wilson loop to EST operators, and use it to compute the potentials. At leading order, they turn out to depend on two parameters only. In Sec. IV A we review the expressions of the doubly heavy baryon hyperfine splittings. In Sec. IV B we model the spin-dependent potentials using suitable interpolations between the known short-distance behavior and the just calculated long-distance one. Then, the remaining unknown parameters of the parametrization of the potentials are fitted to lattice data on the hyperfine splittings of doubly heavy baryons. Using these parameters, we predict the spectrum of doubly heavy baryons including hyperfine contributions in Sec. IV C. We compare our results with other approaches in Sec. V. We close the paper with some conclusions in Sec VI. In Appendix A we calculate the Casimir energy (Lüscher term) in the EST. Finally, in Appendix B we give expressions for the short-distance regime constants as correlators in weakly-coupled pNRQCD.

II. DOUBLY HEAVY BARYON POTENTIALS
A. General expressions A general EFT framework to describe any doubly heavy hadron has been presented in Ref.
[1]. The EFT was worked out up to 1/m Q including the terms that depend on the heavy-quark spin and angular momentum. The matching expressions of the potentials in terms of operator insertions in the Wilson loop can also be found in Ref. [1]. This EFT framework has been applied to doubly heavy baryons in Ref. [7] where the spectrum associated to the four lowest lying static energies was obtained. These static energies are characterized by the representation of D ∞h and the quantum numbers of the light-quark operator that interpolates them, in particular the spin κ and parity p. In the present work we will only consider the cases with light-quark interpolating operator with κ p = (1/2) ± . These spin-1/2 cases only have one possible projection into the heavy-quark axis and therefore each correspond to a single D ∞h representation. These are (1/2) g and (1/2) u for the κ p = (1/2) + and κ p = (1/2) − operators, respectively.
(1/2) ± = Q (1/2) ± (t/2, R) . . . Q † (1/2) ± (−t/2, R)P e −ig C 1 +C 2 dz µ A µ (z) , (8) with C 1 and C 2 the upper and lower paths of a rectangular Wilson loop. Note that, unlike the quark-antiquark case, the flow is in the same direction for both paths. The interpolating operators are where α = −1/2, 1/2, and we have used the following3 tensor invariants B. Short-distance potentials The short-distance regime is characterized by r Λ −1 QCD . Since, in this regime the heavy-quark-pair distance and Λ −1 QCD are well-separated scales the matching of NRQCD to the BOEFT for doubly heavy baryons can be done in two steps. First, one integrates out the heavy-quark-pair distance, which can be done in perturbation theory. This produces weakly-coupled potential NRQCD (pNRQCD) for doubly heavy systems presented in Ref. [12]. Then, integrating out the Λ QCD modes one recovers BOEFT. This procedure delivers expressions of the potentials in Eqs. (4)- (7) as an expansion in the heavy-quark-pair distance. An analogous approach was used in Refs. [13,14] to determine short distance expansion of the hybrid quarkonium potentials. All the potentials follow the same general structure in the short-distance regime; a possible nonanalytic term in r produced by integrating out the heavy-quark-pair distance and an expansion in powers of r 2 with nonperturbative coefficients. These nonperturbatice coefficients only depend on the Λ QCD scale and can be expressed as weakly-coupled pNRQCD correlators of light quark and gluon operators. Figure 1. Matching of the Wilson loop for the static potential for doubly heavy baryons, with the expansion in weakly-coupled pNRQCD up to next-to-leading order. The single lines represent the antitriplet fields, the double lines the sextet field, the dotted and the curly lines the light-quark and transverse gluon fields respectively (emissions of longitudinal gluon fields from the triplet and sextet fields and from the vertices are omitted). The crossed circles indicate the insertion of a Q operator and the square or diamond the insertion of a chromoelectric dipole or quadrupole operator, respectively.
The expansion of the static potential in Eq. (4) is given diagrammatically in Fig. 1 and corresponds to the following form with the nonperturbative constants given as pNRQCD correlators in Appendix B.
For the heavy-quark spin and angular-momentum dependent potentials the short-distance expansion of the potentials is given diagrammatically in Fig. 2 and are as follows: (1/2) ± (r) = c F ∆ (1/2) ± + ∆ (1,0) (1/2) ± r 2 + . . . , (1/2) ± r 2 + . . . , Matching of the heavy-quark spin dependent potentials up to next-to-leading order in weakly-coupled pNRQCD. The legend is as in Fig. 1 with the addition of the solid dot, a white-dotted square, and a white-dotted diamond representing the insertion of a leading-order, dipole, and quadrupole heavy-quark spin chromomagnetic couplings, respectively. Further next to leading diagrams can be generated by changing the order of the different internal vertices, and by adding an extra transverse gluon emission to the heavy-quark spin chromomagnetic couplings. The potential of the heavy-quark angular-momentum dependent operator is matched to an analogous expansion.
with the nonperturbative constants given in Appendix B. At leading order both Eqs. (13) and (15) depend on the same correlator and the difference in the contribution to the potential stems from different factors in the coupling of the heavy-quark spin and angular momentum to the chromomagnetic field in the Lagrangian of Eq. (9) in Ref. [12]. The potential of the spin-tensor coupling in the Lagrangian of Eq. (3) vanishes at leading order since, to appear, it requires the insertion in the pNRQCD correlator of operators carrying the dependence on r which are suppressed in the multipole expansion. This type of correlator is also responsible for the next-to-leading order contributions to Eqs. (13) and (15). It is interesting to note, that the next-to-leading coefficient of Eq. (15) can be written as a combination of the next-to-leading coefficients of Eqs. (13) and (14).
In the static and r → 0 limits the heavy-quark pair becomes indistinguishable from a single heavy antiquark. This is the so-called heavy quark-diquark duality [40][41][42][43]. One can use this duality to relate the leading-order coefficients of the expansions of the potentials in Eqs. (12)- (15) to the heavy meson masses. The value of Λ (1/2) + is equal to the leading-nonperturbative contribution to the lowest laying D-or B-meson masses, usually referred to as just Λ.
The value of Λ has been obtained in Refs. [44,45] combining lattice determinations of the heavy-meson masses and perturbative computations of the heavy-quark masses. It is therefore necessary to use values of Λ and the heavy-quark masses computed in the same scheme. We use the values of Ref. [44] in the MRS scheme Following from the heavy quark-diquark duality, the difference Λ (1/2) − − Λ (1/2) + is equal to the mass gap between the ground and first excited heavy-light mesons, up to corrections of order Λ 2 QCD /m Q . The values for this difference are collected in Table I. The values are compatible with the short-distance energy gaps between the static energies (1/2) g and (1/2) u of Refs. [8,9] associated to the light-quark operators (1/2) + and (1/2) − , respectively. )   Table I. Determination of Λ (1/2) − − Λ (1/2) + from D meson mass differences. The masses are taken from the PDG [46]. The uncertainty corresponds only to the experimental uncertainty of the meson masses. The uncertainty in the determination of (Λ (1/2) − − Λ (1/2) + ) due to neglected higher-order terms is expected to be about 30%.
Finally, the value of ∆ (1/2) ± can be related to the hyperfine splittings in D or B mesons [12] with corrections expected to be of order Λ 3 QCD /m 2 Q . The values of ∆ (1/2) ± from Eq. (19) for various heavy-meson masses are found in Table II.
Heavy Mesons ∆ (1/2) ± from the heavy meson masses, taken from the PDG [46], using Eq. (19). We take the renormalization group improved expression for cF at 1 GeV. The uncertainty corresponds only to the experimental uncertainty of the meson masses. The uncertainty in the determination of ∆ (1/2) ± due to neglected higher order terms is expected to be of ∼ 30% for the charm mesons and ∼ 10% for the bottom mesons.

A. Motivation
The QCD potentials for heavy quarks can be calculated assuming the heavy quarks to be static color sources. For a heavy-quark-antiquark system, the leading-order (static) potential is the energy of a source in the fundamental representation and a source in the complex conjugate representation separated at a distance r. Since the system must be a color singlet object, a certain gluon configuration must exist between the two sources in order to achieve so. When the distance is larger than the typical QCD scale rΛ QCD 1, a flux tube emerges [47], with a typical radius ∼ Λ −1 QCD . Assuming a constant energy per unit length in the flux tube leads to a linear potential. The flux-tube dynamics can be described by an EST, which matches the lattice QCD calculations very well for the static potential at long distances in the absence of light quarks [19,22,23]. When light quarks are present, the flux-tube configuration is still observed [48] even though it may break due to light quark-antiquark pair creation, a phenomenon known as string breaking [49,50]. Nevertheless a flux-tube like configuration leading to a linear potential remains as an excited state for r beyond the string-breaking scale.
For a baryon with two heavy quarks, we have an analogous situation. The two sources are now in the fundamental representation, and the gluon configuration linking them must also contain a valence light quark. When rΛ QCD 1 we expect a flux tube to emerge from each source and to joint at the point between them where the valence light quark is at each time. Hence, the naive expectation would be to have a potential with the same string tension as in the quark-antiquark system plus a constant contribution ∼ Λ QCD due to the extra energy provided by the link to the valence light quark. Lattice QCD simulations indeed observe a linear potential [8,51]. Hence, we expect an EST to account for the long-distance behavior of the potential as well. Locally, the EST should be the same as the one for the quark-antiquark system, but it must contain some additional degrees of freedom describing the link to the valence light quark. In particular it must keep its transformation properties under D ∞h and flavor. We propose to add a fermion to the usual EST which transforms like the light quark under flavor and the Lorentz group. We write down a reparametrization invariant Lagrangian, and expand it at the desired order in the effective theory expansion.
The metric g ab induced on the string reads with η αβ the Minkowsky metric, and e α a ≡ ∂ξ α /∂x a the Zweibein. The action of the gluonic string is just proportional to the area of the string world sheet with σ the string tension and g = |detg ab |.
The action of a four-dimensional Dirac field constrained on a string is given by with ρ a ≡ γ µ e a µ . The antisymmetrization of the partial derivative is required by Hermiticity. Note that the action in Eq. (23) is invariant under reparametrizations of the string if we choose ψ(x) to transform like a scalar, and under Lorentz symmetry if we choose ψ(x) to transform like a four-dimensional Dirac field but keeping x invariant.
Let us choose the Gauge or string parametrization Expanding the action in Eq. (22) for small string fluctuations we arrive at and for the case of the fermionic action in Eq. (23) we find with l = 1, 2, and a = 0, 3. The fermion field mode expansion is where E n = p 2 n + m 2 l.q. . If we consider both periodic and antiperiodic solutions p n = nπ/r, n ∈ Z. The spinors are defined as with χ +1/2 = (1, 0), χ −1/2 = (0, 1) andχ s = −iσ 2 χ * s . The commutation relations for the creation an annihilation operators are all the other anticommutators vanish.
The field mode expansion in Eq. (28) contains both positive-and negative-parity modes. Since the spinors fulfill the relation u ± (E, −p) = ±γ 0 u ± (E, p) a convenient choice for the transformation of the creation and annihilation operators under parity is One can split the field mode expansion into two components of well-defined parity with the following definitions with with η P the parity eigenvalue The field mode expansion in Eq. (28) can be rewritten in term of the two components of well-defined parity as

C. Mapping
Our aim is to use the EST introduced in Sec. III B to compute the Wilson loops with operator insertions in Eqs. (4)-(7) which correspond to the potentials in the BOEFT. In order to do so we need a correspondence between NRQCD and EST correlators. This correspondence is defined by a mapping of NRQCD operators to the EST ones with matching symmetry properties. The symmetry transformations which leave a system of two static particles invariant form the group D ∞h , which is the symmetry group of a cylinder. The basic transformations are rotations around the cylinder axis, reflections across a plane including the cylinder axis and parity. The conventional notation for the representations of D ∞h is Λ σ η . Λ is the rotational quantum number, which for integer values is customarily labeled with capital Greek letters, Σ, Π, ∆ . . . for 0 , 1 , 2 . . . . The parity eigenvalue is given as the index η which is labeled as g or u for positive and negative parity, respectively. Finally, σ gives the sign under reflections as + or −; however, it is only written explicitly for the Σ states, because for Λ > 0 rotations around the cylinder axis mix states in this quantum number. An operator belonging to SO(3) ⊗ P representation κ p can be projected into D ∞h representations: the rotational quantum number can take values corresponding to the absolute value of the projections of the spin of the operator into the heavy-quark axis 0 ≤ Λ ≤ |κ| and the reflection eigenvalue corresponds to σ = η(−1) κ . To simplify, we align the heavy-quark-pair axis with the z-axis, i.e., r = (0 , 0 , z), set the heavy-quark positions at z = ±r/2 and the center of mass at R = 0.
Both Dirac and string fermions are spin-1/2 fields and have the same properties under rotations and reflections. Moreover they can only be projected to Λ = 1/2. Therefore, to find the mapping of NRQCD to the EST operators we just need to make sure that the parities coincide with P ± = (1 ± γ 0 )/2. Now, let us focus on the mapping for the chromomagnetic field B, which can be projected into Σ − u and Π u representations. Since we have chosen to align the heavy-quark-pair axis with the z-axis, then B l , l = 1, 2 and B 3 correspond to the Π u and Σ − u representations, respectively. The mapping of the chromomagnetic field into string fluctuations can be found in Ref. [24].
This implies that B l , l = 1, 2 is O(1/r 2 ) and B 3 is O(1/Λ QCD r 3 ). However, mappings into string fermion operators are now possible and in fact provide the leading order contribution to the potentials in Eqs. (5)- (7). This mapping is as follows: with Σ = diag(σ, σ). Note that here both B l , l = 1, 2 and B 3 are O(Λ QCD /r), and hence are more important than the corresponding bosonic operators in Eqs. (41) and (42). Finally, to convert the two-dimensional spin operators in Eqs. (5) and (6) into four-dimensional spin operators, we will use the following prescription D. Long-distance potentials Using the mapping of NRQCD operators in the Wilson loop to EST operators defined by Eqs. (39)-(45) we compute the potentials in Eqs. (4)-(7) as correlators in the EST. For example, let us apply the mapping to the Wilson loop with the insertion of just the light-quark operators in the spatial sides of the loop then the static potential is just Similarly, one can apply the mapping to compute the heavy-quark spin and angular-momentum dependent potentials

Schrödinger equation with the static potential and M
(1) njl the hyperfine contribution. Recall that due to the Pauli principle the heavy-quark spin is s QQ = 0 for l odd and s QQ = 1 for l even. Let us denote the expectation values of the potentials between the radial wave functions as The hyperfine contributions for l = 0 are given by and the splitting is In the case l is an odd number the hyperfine contribution is as follows: which for the cases l = 1, 3 leads to the following splittings For l = 2 the hyperfine contributions are more complicated since they depend on all three potentials in Eq. (3) and the states j = 3/2, 5/2 with = 3/2 and = 5/2 are mixed. For j = 1/2 and 7/2 the contributions are For j = 3/2, 5/2 we have the mixing matrices for = 3/2 and = 5/2 states 1 We diagonalize to obtain the physical states For simplicity we consider the following hyperfine splittings among l = 2 which are linear in the expectation values of the potentials These formulas fix V s1 n2 , V s2 n2 and V l n2 in terms of physical masses. Then, we have the following model-independent predictions

B. Interpolation of the full potentials
We have obtained descriptions of the potentials of the spin and angular-momentum dependent operators in the shortand long-distance regimes in Eqs. (13)-(15) and (48)-(50), respectively. In this section we propose an interpolation between the descriptions of the potentials in these two regions to model the potential in the intermediate distance regime r ∼ 1/Λ QCD . Using this interpolation and the wave functions obtained in Ref. [7], we compute the hyperfine splittings of Sec. IV A in terms of the parameters of the short-and long-distance descriptions. These parameters are then determined by fitting the hyperfine splittings of lattice determinations [30][31][32][33][34][35][36][37][38] of the double charm and bottom baryon spectrum and in the case of the short-distance parameters using heavy quark-diquark symmetry.
The interpolation we propose is constructed by summing the short-and long-distance descriptions multiplied by weight functions depending of r and a new r 0 parameter. The weight functions are w s = r n 0 /(r n + r n 0 ) and w l = r n /(r n + r n 0 ) for the short-and long-distance pieces, respectively. The sum of the weight functions is w s + w l = 1 and the r 0 parameter determines the value of r where both weights are equal. The value of the exponent n is chosen as the minimal value that ensures that the product of the short-and long-distance potentials and the respective weight functions vanishes in the long-and short-distance limits, respectively. For the short-distance potentials we consider the contributions up to next-to-leading order. The resulting interpolated potentials are as follows: with E 1 = (π/r) 2 + m 2 l.q. . Note that for r 0 = 0 we recover the long-distance potentials and for r 0 → ∞ we recover the short-distance potentials.
An accurate determination of m l.q. would require lattice data for the static energies at longer distances than the one currently available. Nevertheless, fitting the long-distance part of the static potential to the lattice data of Refs. [8,9], we find the value m l.q. = 0.226 GeV , for which the contribution to the potentials of the terms proportional to m l.q. is small.
To obtain the unknown parameters in the interpolated potentials in Eqs.
a much better fit than the former. This indicates both that the long distance form is important and that the EST provides a good description of it. We observe that the value of ∆ (1,2) (1/2) + changes significantly, carries large uncertainty, and in the best fits it is compatible with 0. In the case of ∆ (1,0) (1/2) + its value also shows variation, however it is significantly different from zero. Nevertheless, the inclusion of the next-to-leading order terms in the short-distance potentials does not improve the overall quality of the fits. Therefore, it seems that with the current lattice data it is not possible to constrain these next-to-leading order terms in the multipole expansion. The values of Λ f and Λ f stay consistent across the different sets of fits, with Λ f being very stable while Λ f decreasing in absolute value as r 0 gets larger and even changing sign. Our preferred fit is the one with minimal χ 2 d.o.f in Table VII corresponding to r 0 = 0.5 fm. In Fig. 3 we plot the interpolated expressions of the potentials in Eqs. (68)-(70) for the parameter set in the entry for r 0 = 0.5 fm in Table VII and r 0 = 0.3 fm in Table V. In the case of the potentials for κ p = (1/2) − , we plot the potentials with ∆ (0) (1/2) − = 0.075 GeV 2 from the neutral D-meson entry in Table II (1/2) − = 0 GeV 4 and the values of Λ f and Λ f from the entries for r 0 = 0.5 fm in Table VII and r 0 = 0.3 fm in Table V. Although in some cases the potentials in Fig. 3 show significant variation depending on the parameter set used, we will show in the following section that this is not the case for the values of the hyperfine splittings.

C. Doubly heavy baryon spectra
Now we compute the spectrum of double charm and bottom baryons including the hyperfine contributions using the interpolated potentials in Eqs. (68)- (70). For the states associated to the (1/2) g static energy, we take values of the parameters from the entry r 0 = 0.5 fm in Table VII and for comparison the entry r 0 = 0.3 fm in Table V. The results can be found in Tables VIII and IX for double charm and double  (1/2) − = 0 GeV 4 and take the values of Λ f and Λ f from the r 0 = 0.5 fm entry in Table VII and for comparison the entry r 0 = 0.3 fm in Table V. The results can be found in Tables X and XI for double charm and bottom baryons respectively. The results for the spectra for the two sets of parameters are very close.
We plot the spectra for double charm and bottom baryons in Figs. 4 and 5, respectively, for the parameters of the entry r 0 = 0.5 fm in Table VII. Let us discuss the uncertainties of our results. The leading order masses, M (0) , have uncertainties associated to the r 0 =0.3 Table V r 0 =0.5  values of the heavy-quark masses and Λ (1/2) + , in Eqs. (16)- (18), as well as the uncertainty in the parametrization of the static potentials which was estimated as 10 MeV in Ref. [7]. Adding these uncertainties in quadrature we obtain δM  [8,9] that we used to obtain M (0) in Ref. [7]. We expect the contribution due to the unphysical light-quark mass to be almost independent of r. This is supported by the calculations of the charmonium spectrum (with respect to the η c mass) at m π ∼ 400 MeV [52] and m π ∼ 240 MeV [53], in which almost no difference is observed for the masses of the states below threshold. 2 Hence, it will just produce an overall shift to the static energies computed on the lattice. However, in the computation of Ref. [7] the static energies were rescaled in order for the ground state static energy to be given in the short distance by the expression in Eq. (12). Therefore any additive constant contribution to the static energies produces no change in our results.
The hyperfine contribution, M (1) , has uncertainties associated to the statistical errors of the values of the parameters and interpolation of the potentials. The former ones are displayed in parentheses in Tables VIII-XI and are   for double bottom baryons except for a few cases in Tables VIII and X for double charm states where the difference is larger. Finally, one should consider the size of higher-order contributions to the doubly heavy baryon masses. The most important is the contribution form heavy-quark-spin and angular-momentum independent 1/m Q suppressed potential of O(Λ 2 QCD /m Q ), which we estimate as ∼ 64 MeV and ∼ 19 MeV for double charm and double bottom baryons, respectively. However, in the case of the hyperfine splittings the previous contribution cancels out and the higher order corrections correspond to the 1/m 2 Q suppressed potentials of O(Λ 3 QCD /m 2 Q ), which we take as ∼ 14 MeV and ∼ 1 MeV for double charm and double bottom baryons, respectively.
As an example, in the following we show the value of the masses for the double charm ground state doublet, often refereed as Ξ cc [(1/2) + ] and Ξ * cc [(3/2) + ], adding the different uncertainties in quadrature: and for the double bottom ground state doublet: In the hyperfine splittings most of the uncertainties cancel out and hence our results have higher precision The figures above are compatible with all lattice determinations we are aware of, see Table VI of Ref. [7] and Table XII for doubly charmed and doubly bottom baryons respectively.  Table IX. Hyperfine contributions to the double bottom baryons for the (1/2)g static energy for two sets of parameters of the hyperfine potentials. All masses in GeV units.
Let us finally note that the hyperfine splittings of the (1/2) u states are entirely predicted from the long-distance parameters Λ f and Λ f , obtained from fits to the hyperfine splittings of the (1/2) g states, and the only short-distance parameter, ∆ (1/2) − , obtained from the D-meson spectrum. It is then interesting to compare them with the lattice results of Ref. [35]. For the (1/2) u ground state doublet (1/2 − , 3/2 − ), we obtain from Table X 49 (21) MeV for the hyperfine splitting, which agrees well with the 41(21) MeV of [35]. This is a nontrivial test of the EST we use, since the (1/2) + potentials differ from the (1/2) − at long distances in a very particular way [see (48) (1/2) − = 0 GeV 4 ). All masses in GeV units.

V. COMPARISON WITH MODELS
There is a substantial amount of literature regarding doubly heavy baryons in different approaches; various quark models [54][55][56][57][58][59][60][61][62][63][64][65][66][67], Bethe-Salpeter equations [68][69][70], Born-Oppenheimer approximation with model potential [71,72], semiempirical mass formulas [73][74][75], QCD sum rules [76,77], Faddeev equations [78], and bag models [79]. In this section we compare our results with a selected set of model computations and other approaches (see [65,80,81] for further comparisons). In Table XIII we have collected the masses of the ground state doublet in the double charm baryon sector from different approaches. The values of the Ξ cc mass are in good agreement for about 3/4 of the references, including our own value. Considering the uncertainties only a few works show very significant differences. The values for Ξ * cc show more dispersion with only half of the references being compatible with our own value. On the other hand, the splitting between the two masses is compatible with our value for only 1/4 of the references. This is in contrast with lattice QCD calculations, which are compatible with our current result for the hyperfine splitting (76) (see Table VI of ref. [7]).
The masses of the ground state doublet in the double bottom baryon sector are shown in Table XIV. In this case the differences are a lot more significant. For both the Ξ bb and Ξ * bb only Refs. [58,62,74,79] are compatible with our results and in general there is more dispersion among the values of the different model approaches. Although the values of the hyperfine splittings present less variation in absolute values, in relative terms the variation is also larger than in the double charm baryon sector. Moreover, very few values are compatible with ours. This is due to our small uncertainty for the splitting produced by the cancellation of uncertainties associated to various parameters (1/2) + (3/2) + (5/2) + (7/2) + (9/2) + (1/2) -(3/2) -(5/2) -(7/2) -(9/2) - remaining only the uncertainty on higher order contributions, which are small for double bottom baryons. We note that no reference has compatible results with ours for both for the Ξ bb and Ξ * bb masses and the hyperfine splitting. This is in contrast with the good agreement we get with the available lattice results (see Table XII).
The spectrum of doubly heavy baryons beyond the ground state doublet has also been studied in Refs. [58,59,63,65,67,70]. In Fig. 6 we compare our spectra with the ones in Ref. [59,63,67] obtained with a quark model with a relativistic light quark, a nonrelativistic quark model, and a relativistic quark model with a diquark core respectively. The spectra of Ref. [58] is derived from a similar quark model as in Ref. [59], but the values are shifted down by about a 100 MeV. Ref. [70] uses the Bethe-Salpeter equation in a diquark picture and presents a limited number of states in the spin-symmetry limit. The results of Ref. [65] do not include the Pauli principle for the heavy-quark wave functions and we do not consider it beyond the ground state. From Fig. 6 (a) we can see that for double charm baryons the pattern of states beyond the ground state doublet does not agree with ours in none of the cases or among the quark model approaches themselves. For all displayed model spectra the excited states lie (much) lower than ours. This is in contrast to the overall agreement found with lattice calculations in [1]. For double bottom baryons [see Fig. 6 (b)] the discrepancies reach the ground state doublet, as the results of Ref. [59], and to a lesser extend the ones of Ref. [63], lie higher than ours. However, there is agreement for the first excited (odd-parity) doublet, except for Ref. [59]. For higher states the discrepancies persist, except for the odd-parity states of Ref. [67], which are compatible with ours.

VI. CONCLUSIONS
An EFT describing doubly heavy hadrons was put forward in Ref.
[1]. It is built upon the nonrelativistic expansion of the heavy quarks and the adiabatic expansion between the dynamics of the heavy quarks and the light degrees of freedom corresponding to the gluons and light quarks. The EFT was constructed in the single hadron sector up to the heavy-quark spin and angular momentum terms suppressed by 1/m Q . Expressions of the potentials as operator insertions in the Wilson loop were obtained by matching the EFT to NRQCD. The computation of the Wilson loop with operator insertions cannot be done using perturbative techniques and should be carried out (ideally) in lattice QCD or other nonperturbative approaches (see for instance [83] for an AdS/CFT inspired proposal).
In Ref. [7] this EFT framework was applied to doubly heavy baryons. Using the lattice data of Refs. [8,9] for the static energies the leading-order spectrum of doubly charm and bottom baryons was computed for the four lowest lying static states. However, since there are no available lattice determinations of the potentials of the heavy-quark spin and angular-momentum operators, the computation of the hyperfine contributions to the doubly heavy baryon masses was not possible.
In this paper we presented a parametrization of the 1/m Q suppressed heavy-quark spin and angular-momentum operators with a minimal amount of modeling, the general idea of which can be extended to other potentials for doubly heavy hadrons, such as double charm tetraquark, T + cc , recently discovered by the LHCb Collaboration [84]. This parametrization of the potentials is based in their description in short-and long-distance regimes. In the short-distance regime, defined as r 1/Λ QCD , the Wilson-loop expressions of the potentials can be expanded in the multipole expansion. This can be done using weakly-coupled pNRQCD, which is the EFT that incorporates the multipole expansion systematically, for two heavy quarks [12]. This produces short-distance expressions of the potentials as an expansion in powers of r 2 with coefficients that encode the nonperturbative dynamics of the light Ref. [63] Ref. [67] (1/2) + (3/2) + (5/2) + (7/2) + (1/2) -(3/2) -(5/2) -

Our results
Ref. [59] Ref. [63] Ref. [67] (1/2) + (3/2) + (5/2) + (7/2) + (1/2) -(3/2) -(5/2) -(7/2) -(9/2) -  Tables VIII -XI. degrees of freedom, 3 which we show in Sec. II B and Appendix B. At leading order in the multipole expansion only one coefficient is necessary and it can be determined in a model-independent way using the heavy quark-diquark duality from the heavy-mesons masses. The long-distance regime is characterized by r 1/Λ QCD . In the case of a heavy-quark-antiquark pair it is known from lattice QCD that in this regime a gluonic flux tube connecting the two heavy quarks emerges. It is well-known that an Effective String Theory (EST) [17][18][19] reproduces accurately the lattice determinations [24,26]. In Sec. III we propose an extension of the EST to include the presence of a fermion constrained to move on the string. We obtain a mapping of the NRQCD operators inserted in the Wilson loop to operators in the EST based on imposing the same transformation properties under D ∞h and flavor. Using this mapping we can translate the Wilson-loop expressions for the potentials to EST correlators and evaluate them. This procedure yields long-distance expressions of the potentials depending on two unknown coefficients of the EST. Additionally, we compute the vacuum energy in the EST with fermions in Appendix A and show that (i) the string tension runs with the square of the mass of the fermion and (ii) the sign of the Lüscher term changes. These features can in principle be checked by lattice calculations of the ground state energy of two static quarks separated at a large distance with an additional light quark. 4 The final parametrization of the potentials is obtained by interpolating between the short-and long-distance descriptions. We choose the most simple interpolation that ensures that the correct short-and long-distance behavior are recovered in the corresponding limits. Nevertheless, an extra parameter is introduced in the definition of the interpolation. The hyperfine contributions to doubly heavy baryons can be computed using these parametrizations of the heavy-quark spin and angular-momentum dependent potentials. The values of the remaining unknown parameters are determined by fitting the hyperfine splittings obtained in lattice QCD in Refs. [30][31][32][33][34][35][36][37][38] for several S-, P -and , D-wave multiplets. This guarantees that all our inputs are from QCD, and the modeling is reduced to the choice of interpolation, provided that the EST we use is the correct effective theory at long distances. Using the parameters thus determined, we make predictions for the spectrum of double charm and bottom baryons including the hyperfine contributions. Our results are summarized in Tables VIII-XI and in Figs. 4 and 5.
Finally in Sec. V we compared our results with previous model approaches and sum rules determinations. We observe a huge dispersion of results. In the absence of lattice calculations for many states, specially for double bottom baryons, our EFT approach offers a framework in which modeling is minimal and errors can be reliably quantified, unlike in most models. Since lattice determinations of the potentials for doubly heavy hadrons are difficult, in particular when unquenched simulations are required, the procedure outlined in the paper and in Ref. [29], to obtain reliable parametrizations of the potentials can be of significant utility in future studies of doubly heavy hadrons. In turn, this motivates further development of the EST to cases with multiple light quarks.
where we have used the equation of motion of the fermion field to simplify the fermionic term. Now, let us compute the expected value of the Hamiltonian in the vacuum 0|H|0 = π r ζ(−1) − 2n f ∞ n=1 nπ r 2 + m 2 l.q. .

(A2)
We have obtained the standard Lüscher term for the bosonic string plus a new contribution from the fermion field. The new contribution is ill defined in an analogous way as the Lüscher term is. The latter can easily be defined through the analytic continuation of the Riemann zeta function but the nonvanishing m l.q. in the former requires extra care.
In the nonrelativistic case (which implies a cutoff in the sum so that nπ r m l.q. ) there are no quantum fluctuation contributions of the fermion to the vacuum energy, In the general case, we will use dimensional regularization. Let us first write  2 k 2 + nπ Note that the 1/ pole can be absorbed in a redefinition of the string tension µ-independence leads to µ dσ(µ) dµ = n f π m 2 l.q. , which implies that the string tension runs in such a way that it decreases at large distances. For this to be consistent within an EFT framework, we need m 2 l.q. σ(µ), so that the replacement in Eq. (A8) generates higher-order terms elsewhere. Since σ(µ) ∼ Λ 2 QCD , we need m l.q. Λ QCD , which may be achieved by implementing chiral symmetry linearly in the fermion fields on the world sheet. Numerically, we find m 2 l.q. ∼ 0.051 GeV 2 whereas σ(µ) ∼ 0.21 GeV 2 .