Effective field theory for double heavy baryons at strong coupling

We present an effective field theory for doubly heavy baryons that goes beyond the compact heavy diquark approximation. The heavy-quark distance $r$ is only restricted to $m_Q\gg 1/r \gg E_{bin}$, where $m_Q$ is the mass of the heavy quark and $E_{bin}$ the typical binding energy. This means that the size of the heavy diquark can be as large as the typical size of a light hadron. We start from Non-Relativistic QCD, and build the effective field theory at next-to-leading in the $1/m_Q$ expansion. At leading order the effective field theory reduces to the Born-Oppenheimer approximation. The BO potentials are obtained from available lattice QCD data. The spectrum for double charm baryons below threshold is compatible with most of the lattice QCD results. We present for the first time the full spin averaged double bottom baryon spectrum below threshold based on QCD. We also present model-independent formulas for the spin splittings.


I. INTRODUCTION
The recent discovery of the Ξ ++ cc baryon by the LHCb Collaboration [1,2], together with the expectation that other states can be confirmed or discovered in the near future, has revitalized the interest of the theoretical community on double heavy baryons. Earlier SELEX claims [3,4] on the discovery of Ξ + cc appear to clash with LHCb searches [5,6], as well as earlier ones by BABAR [7] and BELLE [8].
From the theoretical side, a QCD based approach to double heavy baryons was already considered in the early days of Heavy Quark Effective Theory [9]. The key observation was that a QQ state at short distances has an attractive channel in the 3 * representation. Then, if the heavy quark masses are large enough, the QQ would form a compact hard core and the lowest-lying excitations would be given by the typical hadronic scale Λ QCD . The spectrum would then be analogous to the one of heavy-light mesons. Predictions for the hyperfine splitting were put forward within this approach [9], sometimes referred to as heavy quark-diquark duality. This framework was put in a more solid theoretical basis in Ref. [10] by working out a lower energy effective theory (EFT) for the QQ system, similar to potential nonrelativistic QCD (pNRQCD) [11,12]. The quark-diquark duality assumption could then be quantified: it would hold when the typical binding energies of the QQ systems, E bin , are much larger than the typical hadronic scale (E bin Λ QCD ). In practice, however, the duality hypothesis does not hold for double charm nor for double bottom baryons. Indeed, it is well known that E bin Λ QCD does not hold for charmonium and bottomonium (QQ), the attractive channel of which is twice stronger than the one of their QQ counterpart. The EFT built in Ref. [10] is valid whenever the typical size of the system is smaller than typical hadronic size, r 1/Λ QCD , and in particular it is still correct for E bin ∼ Λ QCD . In this case, the excitations due to the internal QQ dynamics compete with the excitations of the light degrees of freedom (which include a light quark) surrounding the QQ compact core, and the spectrum of the lower-lying states will not mimic the one of the heavy-light mesons anymore. For charmonium and bottomonium the hypothesis that E bin ∼ Λ QCD is only reasonable for the ground state and the gross description of excited states clearly requires the introduction of a confining potential, in addition to the Coulomb-like potential that arises from the hypothesis r 1/Λ QCD . Therefore, for doubly heavy charm or bottom baryons, for which the attraction of the Coulomb-like potential is twice weaker than that for quarkonium, the assumption E bin ∼ Λ QCD seems unlikely to hold.
It is the aim of this paper to build an EFT for QQq systems in which the hypothesis r 1/Λ QCD is released, along the lines suggested in Ref. [13]. In fact, this paper may be considered a concrete example of a more general formalism developed in an accompanying one [14]. This EFT is built upon the heavy quark mass expansion, m Q Λ QCD , and an adiabatic expansion between the dynamics of the heavy quarks, and the light degrees of freedom, the gluons and light quarks, Λ QCD E bin . Under these assumptions, the heavy quark mass and the typical hadronic scale can be integrated out producing an EFT that at leading order (LO) consist of a set of wave function fields for the QQ system with the quantum numbers of the light degrees of freedom, in addition to the ones of the QQ, interacting through a number of Born-Oppenheimer (BO) potentials. Since the BO potentials cannot be calculated in perturbation theory, we shall use available lattice data for them. Four different BO potentials turn out to be relevant for describing the spectrum of double charm and double bottom baryons below the first heavy-meson-heavy-baryon threshold. At LO, the BO potentials are flavor independent. They neither depend on the heavy quark mass, nor on light quark ones m q , if Λ QCD m q is assumed. Hence, the LO Lagrangian enjoys heavy quark spin symmetry and chiral symmetry. We calculate the spectrum using the lattice data of Refs. [15,16] as the input for the BO potential. We also work out the EFT at next-to-leading order (NLO) in the 1/m Q expansion for the terms depending on the heavy quark spin and angular momentum. Heavy quark spin symmetry is violated at this order. We put forward model-independent formulas for the spin splittings. In particular, we make a prediction for the spin partner of Ξ ++ cc . Another QCD-based approach to double heavy baryons is lattice QCD. The study of double heavy baryons on the lattice is quite challenging due to the wide spread of the characteristic scales. The light quark and gluon dynamics occurs at low energies and requires of large lattices for accurate simulations, while the heavy quarks necessitate small lattice spacings. The combination of both requirements results in computationally demanding simulations. To reduce the computational cost, early studies relied on the quenched approximation and were carried out in lattice nonrelativistic QCD (NRQCD) [17][18][19][20]. For doubly charmed baryons, relativistic actions were latter used in Refs. [21,22] and full QCD simulations in Refs. [23][24][25][26][27][28][29][30][31][32][33][34][35]. The latter, however, were limited to the lowest-lying spin 1/2 and 3/2. The spectrum of a wider array of j η P states was obtained in Refs. [36,37]. In the bottom sector unquenched computations have been carried out in Refs. [24,32,38,39] but still using nonrelativistic bottom quarks. Lattice NRQCD, expanded about the static limit, is also necessary to compute the matching coefficients of the EFT presented here, namely the BO potentials. This can be done by using the expressions of the matching coefficients as operator insertions in the Wilson loop that we present in an accompanying paper [14]. An example of this are the static energies computed in Refs. [15,16] that we use as input.
We have organized the paper as follows. In Sec. II we construct the EFT at NLO. In Sec. III we focus on the LO Lagrangian. We find suitable parametrizations of the lattice data of Refs. [15,16] that fulfill expected short and long distance constraints, and calculate the spectrum. In Sec. IV, we discuss the hyperfine splittings and produce a number of model-independent results. We compare our findings with lattice QCD results in Sec. V. Sec. VI is devoted to the conclusions. In Appendix A we derive the coupled Schrödinger equations for the κ p = (3/2) − states, which are a mixture of the (1/2) u and (3/2) u static energies, and in Appendix B we collect the plots of the double heavy baryon radial wave functions.

II. NONRELATIVISTIC EFT FOR DOUBLE HEAVY BARYONS
Double heavy baryons are composed of two distinct components: a heavy quark pair in a 3 * color state and a light quark. The heavy quark mass is larger than the typical hadronic scale, m Q Λ QCD , and therefore the natural starting point to study double heavy baryons is NRQCD [40][41][42]. At leading order in the 1/m Q expansion the heavy quarks are static and the spectrum of the theory is given by the so-called static energies. The static energies are the energy eigenvalues of the static eigenstates which are characterized by the following set of quantum numbers: the flavor of the light quark, the heavy quark pair relative distance r, and the representations of D ∞h . The latter is a cylindrical symmetry group, also encountered in diatomic molecules. The representations of D ∞h are customarily written as Λ η , with Λ the absolute value of the projection of the light quark state angular momentum on the axis joining the two heavy quarks,r, and η = ±1 is the parity eigenvalue, denoted by g = +1 and u = −1 1 . The static energies are nonperturbative objects and have to be computed in lattice QCD. Preliminary calculations for the static energies of two heavy quarks and a light quark state in the isospin limit were presented in Refs. [15,16]. The lattice data for the four lowest lying states is shown in Fig. 1. The lowest lying static energy corresponds to the representation (1/2) g followed by three very close states corresponding to the representations (1/2) u , (3/2) u and (1/2) u . The prime is used to indicate an excited state with the same representation as a lower lying one. In the short-distance limit the symmetry group is enlarged from D ∞h to O(3) and the states can be labeled by their spin (κ) and parity (p). In this limit the κ p of the three lowest lying states found in Refs. [15,16] correspond to (1/2) + , (3/2) − , and (1/2) − . In the short-distance limit the two heavy quarks act as a single heavy antiquark; this is sometimes refereed to as quark-diquark symmetry [9,[43][44][45]. Therefore, the double heavy baryon system in the short-distance limit is equivalent to a heavy-light meson and in fact the spectrum found in Refs. [15,16] is consistent with the D and B meson spectra within a few tenths of MeV in this limit. Projecting the κ p states into the heavy quark axis one can obtain states in representations of D ∞h . We show the correspondence in Table I. The most significant feature in Table I is that (3/2) − projects both to (1/2) u and (3/2) u , thus we expect these two static energies to be degenerate in the short-distance limit on symmetry grounds. This behavior can in fact be observed form the lattice data in Fig. 1.
The short-distance degeneracy between (1/2) u and (1/2) u − (3/2) u is a reflection of the degeneracy between the short distance states (3/2) − and (1/2) − , which as far as we know is accidental. The static energies can be computed in the lattice from the large time limit of correlators of appropriate operators that overlap dominantly with the static states. Such operators must have the same quantum numbers as the static states. The static energies are given by where r and R are the relative and the center-of-mass coordinate of the heavy quark pair, respectively, and with ψ the Pauli spinor fields that annihilate a heavy quark, T l ij = lij / √ 2 is a 3 * irreducible tensor [10], C j3 m3 where P is the path-ordering operator. The projectors P κΛ in Eq. (1) act on the light quark spin indices and for spin-1/2 and spin-3/2 take the form with 1 lq n an identity matrix in the light quark spin space of dimension n = 2κ + 1. The spectrum of double heavy baryons corresponds the heavy-quark pair bound states on the static energies defined by Eq. (1), thus the lowest lying states correspond to the static energies in Fig. 1. The binding energies are expected to be smaller than Λ QCD and therefore these bound states can be described in a BO-inspired approach [46]. That is, incorporating an adiabatic expansion in the energy scales of typical binding energies over the one of light quark and gluon degrees of freedom, Λ QCD .
The EFT describing heavy exotic hadrons and double heavy baryons up to 1/m Q for any spin of the light-quark and gluonic degrees of freedom has been presented in Ref. [14]. In the case of double heavy baryons corresponding to the static states of Fig. 1 the Lagrangian is as follows: with the Ψ fields understood as depending on t, r, R, where r = x 1 − x 2 and R = (x 1 + x 2 )/2 are the relative and center-of-mass coordinates of a heavy quark pair. The Ψ fields live both in the light quark and heavy quark pair spin spaces. In the Lagrangian on Eq. (9) we have chosen to leave the spin indices implicit.
The Hamiltonian densities h κ p have the following expansion up to 1/m Q with p = −i∇ r and P = −i∇ R . The kinetic terms in Eq. (10) are diagonal in spin space while the potentials are not. The static potentials, V (0) , are diagonal in the heavy quark spin space, due to heavy quark spin symmetry, while the light quark spin structure is determined by the representations of D ∞h that the κ p quantum numbers can be associated with: with P κΛ the projectors into representations of D ∞h in the spin-κ space. These fulfill the usual projector properties: they are idempotent P 2 κΛ = P κΛ , orthogonal to each other P κΛ P κΛ = δ ΛΛ , and add up to the identity in the spin-κ space Λ P κΛ = 1 2κ+1 .
The subleading potentials V (1) can be split into terms that depend on S QQ or L QQ and terms that do not. The former read with the total heavy quark spin defined as 2S QQ = σ QQ = σ Q1 1 2 Q2 + 1 2 Q1 σ Q2 where the 1 2 are identity matrices in the heavy quark spin space for the heavy-quark labeled in the subindex, L QQ = r × p and the spin-2 irreducible tensor is defined as The heavy quark spin component of the Ψ fields is given by χ Q1 s χ Q2 r with χ s the usual spin-1/2 two-component spinors. The light quark spin component, χ lq α , is a 2-or 4-component spinor for the Ψ (1/2) ± and Ψ (3/2) − fields, respectively.
The matching of the potentials in Eqs. (12), (13), and (14) in terms of static Wilson loops has been presented in Ref. [14]. These Wilson loops are nonperturbative quantities and should be computed on the lattice. The only ones available are the ones corresponding to the static energies, from Refs. [15,16] shown in Fig. 1, which match to the static potentials Furthermore the form of the static potentials can be constrained in the short-and long-distance limits from general grounds. In the short-distance regime, r 1/Λ QCD , one can integrate out the relative momentum scale perturbatively and build weakly coupled pNRQCD [11,12] for double heavy baryons as was done in Ref. [10]. In this regime the potential at leading order in the multipole expansion is the sum of the Coulomb-like potential for two heavy quarks in a triplet state plus a nonperturbative constant with where Q κ p is the light-quark piece of the interpolating operator O κ p . Using quark-diquark symmetry the values of Λ κ p can be obtained from analysis of the D and B meson masses [47,48]. On the other hand, in the long-distance regime, r 1/Λ QCD , we expect the formation of a flux tube that behaves as a quantum string. The formation of such flux tubes has been observed from lattice QCD [49][50][51] in standard and hybrid quarkonium, but is not yet confirmed for double heavy baryons. Nevertheless, the data from Refs. [15,16,52], although it does not strictly reach the long-distance regime, does show a remarkable linear behavior for the larger heavy quark pair distances. Therefore in the long-distance regime we expect the potential to be linear in r plus a possible additive constant depending on the light quark mass

III. DOUBLE HEAVY BARYON SPECTRUM AT LEADING ORDER
The spectrum of double heavy baryons is obtained by solving the Schrödinger equations resulting from the LO Lagrangian The first two terms in Eq. (20) define standard Schrödinger equations, while the last term corresponds to two sets of coupled Schrödinger equations corresponding to the two possible parities of the double heavy baryon states. The coupled Schrödinger equations for the latter case are derived in Appendix A. The main point to keep in mind is that while the states associated with κ p = (1/2) ± are eigenstates of L 2 QQ with eigenvalue l(l + 1), the states associated with κ p = (3/2) − are eigenstates of L 2 = (L QQ + S 3/2 ) 2 with eigenvalue ( + 1). Additionally, in the latter case there are states with positive and negative parity for each with different masses due to the mixing that leads into the coupled Schrödinger equations. This is the so-called Λ-doubling effect known from molecular physics.
Using Eq. (16), the static potentials in Eq. (20) can be obtained by fitting the lattice data from Refs. [15,16] on the corresponding static energies. To fit them we use the following parametrizations of the potentials which interpolate between the expected short-and long-distance behavior: with the ν scale taken as the inverse of the lattice spacing ν lat = 1/a = 2.16 GeV (a = 0.084fm). We constrain the The energy offset c 1 is an additive constant that affects the masses of all the double heavy baryon states; therefore, it is important to determine it accurately. To do so, we fit the short-distance data (r = a, 2a, 3a) to the one-loop expression of the heavy-quark-heavy-quark 3 * potential added to c 1 The value we obtain is and therefore The long-distance behavior is given by the linear term whose constant is fixed at σ = 0.21 GeV 2 [53]. One could obtain similar values by fitting the longer distance lattice points to a straight line plus a constant term; however, there is some correlation between the two parameters that makes this procedure undesirable. The rest of the parameters are obtained by fitting the potentials in Eqs. (21)- (24) to the lattice data. The fits yield the following values: c 2 = 15.782 GeV 2 , c 3 = 9.580 GeV, c 5 = 13.265 GeV 2 , c 6 = 6.560 GeV , In Fig. 1 we show the fitted potentials together with the lattice data. We solve numerically the radial Schrödinger equations in Eqs. (21)- (24) with the short-distance constants of each potential (c 1 , c 4 or b 1 ) subtracted. The effect of this is just to set a common origin of energies so we can compare the binding energies obtained from the Schrödinger equations. To obtain the total masses of the baryons we add to the binding energy two times the heavy quark mass and Λ ( In practice what we are doing is removing the ambiguity on the origin of energies of the lattice data on the static energies by rescaling them by a factor Λ (1/2) + − c 1 so that the short-distance behavior of V (1/2) + matches exactly Eq. (17) The same masses have been used in the kinetic terms of the Schrödinger equations. In the Tables III and IV we present the results for the spectrum of cc and bb double heavy baryons. In Table II we show the full quantum numbers of the double baryon states including mixings and degenerate spin multiplets. We also represent the spectra in terms of j η P states graphically in Figs. 2 and 3 for double charm and double bottom baryons, respectively. The BO approach, corresponding to our LO Lagrangian in Eq. (20), may only be completely consistent for states below the first heavy-meson-heavy-baryon threshold. Above threshold, the states predicted from the Schrödinger equations are expected to acquire widths, corresponding to decays into the particles that form the threshold, and the double heavy baryon masses can also vary due to coupled channel effects. Therefore, above heavy-baryon-meson thresholds our results should be taken with care. In the charm sector the first heavy-baryon-meson threshold is formed by the Λ 0 c and D 0 at ∼ 4.151 GeV. There are only four spin multiplets below this threshold with quantum numbers Λ η = (1/2) g with l = 0, 1, n = 1, Λ η = (1/2) u with l = 0 n = 1 , and Λ η = (3/2) u \(1/2) u with η P = (3/2) + , n = 1. In the bottom sector the first heavy-baryon-meson threshold is formed by Λ 0 b and B 0 at ∼ 10.899 GeV. In this case sixteen spin multiplets are found below threshold: eight multiplets for Λ η = (1/2) g , two for Λ η = (1/2) u , four for Λ η = (3/2) u \(1/2) u and two for (1/2) u .
There are two possible different choices for the origin of energies of the potentials that would affect the number of multiplets below heavy-baryon-meson thresholds. First, instead of using Λ (1/2) + one could adjust the short-distance constant c 1 to reproduce the physical mass of Ξ ++ cc from Ref. [1], which in practice is equivalent to shift down the masses of all the multiplets by 91 MeV. Second, the energy gaps E latt (1a) ( MeV correspond to the mass gaps of the light quark states (3/2) − and (1/2) − with respect to (1/2) + which corresponds the mass difference of the ground and first excited heavy-light mesons. From the PDG [54], these mass gaps read (1/2) + (3/2) + (5/2) + (7/2) + (1/2) -(3/2) -(5/2) -(7/2) -(9/2) -3.   Table III and the corresponding j η P multiplets from Table II. Each line represents a state; the lines with a dot indicate two degenerate states. The color indicates the static energies that generate the state.
The B meson mass gap, which should give the most accurate value, is compatible with the lattice data and any shifting of the states would be small.
To bbq , respectively. However, this uncertainty drops in the mass differences. Additionally, there is the uncertainty on the lattice data for the static energies of Refs. [15,16] and the model dependence on the long-distance parametrization of the static potentials. We estimate the effect of these two sources of uncertainty to be no larger than 10 MeV. Finally, one should keep in mind that M (0) is only the LO contribution to the full double heavy baryon masses. The NLO corrections can be estimated to be of parametrical size Λ 2 QCD /m Q , which amounts to ∼ 100 MeV and ∼ 30 MeV for M ccq and M bbq , respectively. In the next section we study in more detail the heavy quark spin and angular momentum NLO contributions to the double heavy baryon masses. However, one should keep in mind that NLO contributions independent of the spin and angular momentum do also exist.
IV. SPIN SPLITTINGS FOR κ = 1/2 STATES At LO the potential is independent of the spin of the heavy quarks, hence the final j η P states appear in degenerate multiplets. This degeneracy is broken by the spin-and angular-momentum-dependent operators in the Lagrangian in Eqs. (13) and (14). At the moment the potentials corresponding to these spin-dependent operators are unknown   Table IV and the corresponding j η P multiplets from and therefore the full computation of this spin-splitting corrections is not possible. However, in the case of the κ = 1/2 states, due to the relatively small number of operators in the Lagrangian in Eq. (13) it is possible to obtain relations between the masses of states belonging to multiplets of a given l independent of the shape of the potentials. We summarize the states in the multiplets for l = 0, 1, 2 and the corresponding angular expected values of the spin-dependent operators in Table V njl the 1/m Q spin-dependent contributions. Notice that the Pauli exclusion principle constrains s QQ = 0 for odd l and s QQ = 1 for even l; therefore, s QQ is not needed to label the states.
In the case with l = 1 the Pauli principle fixes the heavy quark spin in a singlet state, and both the expected values of S 1/2 · S QQ and S 1/2 · (T 2 · S QQ ) vanish. Only the contribution of S 1/2 · L QQ remains which is found to be that is, the same pattern as the l = 0 states. We can also express the splitting as a mass relation between the two states that form the l = 1 multiplet 2M n 3 l sQQ j S 1/2 · SQQ S 1/2 · (T 2 · SQQ) S 1/2 Notice, that both Eqs. (39) and (41) are equivalent to the statement that the spin average of the l = 0 and l = 1 multiplets is equal to our LO masses M For j = 3/2, 5/2 we have the mixing matrices for = 3/2 and = 5/2 states (1/2) ± n2 + 1 10 V s2 We diagonalize to obtain the physical states Let us consider the following hyperfine splittings among l = 2 which are linear in the expectation values of the potentials These formulas fix V s1 (1/2) ± n2 , V s2 (1/2) ± n2 and V l (1/2) ± n2 in terms of physical masses. Then, we have the following model-independent predictions

V. COMPARISON WITH LATTICE QCD
Let us first discuss double charm baryons. For the ground state spin multiplet (1/2, 3/2) + we get M Ξcc = 3712(63) MeV. If we assign to the (1/2) + state the mass of the Ξ ++ cc found by LHCb [1], M Ξcc = 3621.2 ± 0.7 MeV, and use it in Eq. (39) we obtain a prediction for the mass of the (3/2) + state, M Ξ * cc = 3757(68) MeV. This implies a hyperfine splitting δ hf = 136(44) MeV. We estimate the uncertainty as corrections of O(Λ 3 QCD /m 2 c ) added quadratically to the uncertainties in the masses, for which we take into account that the uncertainties in the Λ (1/2) + and m c vanish for mass differences. This is about 30 MeV larger than what is obtained by assuming heavy quark-diquark symmetry, namely a very compact diquark [9,10,45,55], using the most up to date PDG values. In Table VI we compare this mass splitting to the ones obtained from lattice QCD calculations with fully relativistic charm quarks 3 . We observe that our value is about 50 MeV higher than that of Refs. [31-33, 36, 56] except for Refs. [27,34] for which the discrepancy is higher. In most cases, including the uncertainties, the results are compatible. According to Eq. (39) the spin average of the ground state multiplet should coincide with our LO mass. We can check this by comparing the spin averages of the ground state multiplets from the lattice, summarized in Table VI For the first excitation, we have three spin multiplets very close in mass: two (1/2, 3/2) − doublets, the first one corresponding to the P-wave excitation of the heavy quarks (4062 MeV) and the second one to the S-wave ground state of the (1/2) u BO potential (4095 MeV ), and one (1/2, 3/2, 5/2) − triplet corresponding to the S-wave ground state of the (3/2) u \(1/2) u BO potentials (4066 MeV). This structure qualitatively agrees with the results of Ref. [36] (see Fig. 11 in that reference). It also agrees quantitatively for the first doublet, but the second doublet and the triplet are about 100 MeV and 150 MeV lower, respectively, than in Ref. [36]. In Refs. [33,35], only a doublet is reported. The first reference is compatible with our results whereas the ones of the second reference lie about 100 MeV below ours.
The higher excitations are all above the Λ c −D threshold. The only lattice calculation reporting on them is Ref. [36]. For the positive parity states, our lowest l = 2 spin multiplet (1/2, 3/2, 3/2, 5/2, 5/2, 7/2) + lies very close to the (1/2) u singlet (1/2) + and both are above the n = 2, l = 0 spin doublet (1/2, 3/2) + (see Fig. 2) whereas in Ref. [36] all these states lay in the same energy range. The spin average of the l = 2 multiplet is about 40 MeV below ours, the spin singlet (1/2) + lays 13 MeV above ours and n = 2, l = 0 spin doublet has a spin average 47 MeV above ours. In any case, we can check our model-independent formulas in Eqs. (51) and (52) for the l = 2 spin multiplet against this data. Both formulas are compatible with the lattice masses of Ref. [36]. Recall however that the uncertainties in the masses of the states involved are large. For positive parity, we can qualitatively accommodate the rest of our states below 4450 MeV to those of Ref. [36], though precise identifications are difficult due to the large errors in that reference. Beyond that point we have much less positive parity states than Ref. [36]. This is easy to understand because the first BO potential that we left out starts contributing about this energy. The negative parity states below 4600 MeV (two spin doublets and a spin triplet) agree reasonably well with the spin averages of those of Ref. [36], within 72 MeV. Beyond that point we have more states of higher spin. This is so even if we still miss some states that would arise from the first BO potential left out, which are expected to contribute in this region.
Let us next discuss double bottom baryons, for which no experimental evidence exists yet. The fact that the bottom quark mass is more than three times larger than the charm mass poses further difficulties to direct lattice QCD calculations. A way out is to factor out the bottom quark mass from the problem and do the calculations in lattice NRQCD. So far only results from the ground state spin doublet (1/2, 3/2) + are available. Following Eq. (41) we compare the lattice spin average (s.a.) for the ground state to our value M

VI. CONCLUSIONS
We have put forward an EFT for double heavy baryons that goes beyond the compact diquark approximation. This EFT is built upon two expansions on the small ratios between the characteristic energy scales of double heavy baryon systems. These are: the heavy quark mass expansion, m Q Λ QCD , 1/r, and an adiabatic expansion between the heavy quark and light degree of freedom dynamics, Λ QCD , 1/r E bin . At LO the EFT reproduces the Born-Oppenheimer approximation with a set of static energies. In that sense, it bridges smoothly between QCD and potential model calculations. The spectrum of static energies is obtained from lattice QCD data from Refs. [15,16] and general constraints of the shape of these in the short-and long-distance regimes. We take into account the four lowest-lying static energies below the first heavy baryon-meson threshold, which correspond to the representations (1/2) g , (1/2) u , (3/2) u and (1/2) u of D ∞h . These are plotted together with the lattice data in Fig. 1. In the short distance limit, r → 0, the compact diquark approximation becomes exact and the spectrum of static energies matches the heavy meson spectrum (see Table I for the short-distance quantum numbers).
In Sec. III, we obtained the spectrum of ccq and bbq baryons by numerically solving the Schrödinger equations with the aforementioned static energies including the mixing between (1/2) u \(3/2) u produced by the leading nonadiabatic corrections. A clear advantage of our method in comparison to lattice QCD is that obtaining higher excitations is very easy, even though the systematic errors are less under control. Indeed, the condition E bin 1/r is not fulfilled for some multiplets close to or beyond threshold. For ccq (bbq) only three (ten) of the four (sixteen) multiplets below threshold fulfill it. In the bottom case, there is also a multiplet slightly above threshold that fulfills it. However, E bin ∼ 500 MeV for this multiplet, and for two more below threshold, and hence they fail to meet the E bin Λ QCD condition. In any case, we provide here for the first time the full spectrum below threshold for double bottom baryons based on lattice QCD data.
Our results do not support the heavy quark-diquark approximation, neither for double charm baryons nor for double bottom ones. For ccq, the first angular excitation of the diquark core (negative parity lower blue bands in Fig. 2) is almost degenerate with the two lower light quark excitations (green and red lines in Fig. 2). For bottomonium, the situation is even worse as the first angular excitation of the diquark core (negative parity lower blue bands in Fig. 3) is clearly below the two lower light quark excitations (green and red lines in Fig. 3) and the first principal quantum number excitation has about the same energy as the latter. This is due to the fact that both for ccq and bbq the ground state already lies on the non-Coulombic part of the potential, which is clearly reflected in the positive binding energies displayed in Tables III and IV. Nevertheless, it is still remarkable that the gap between the ground state and the first light quark excitation (383 MeV for ccq and 387 MeV for bbq) is reasonably close to the gap for the spin average of the heavy-light meson systems (427 MeV forcq and 406 MeV forbq).
In Sec. IV, we discussed the hyperfine splittings generated by the NLO heavy quark spin-and angular-momentumdependent operators in the Lagrangian of Eq. (13). Unlike quarkonium, where these corrections start at 1/m 2 Q suppression, in double heavy baryons these start at only 1/m Q suppression. This is a feature shared with exotic heavy hadrons [14], such as heavy hybrids [57][58][59]. These operators are accompanied with, so far, unknown potentials. Despite this, it is possible to obtain relations between the hyperfine splittings of different states that are independent of the specific shape of the potentials. These are presented in Eqs. (39), (41) , (51), and (52).
We have compared our results for the LO masses with the available lattice calculations for the ground state spin multiplet making use of Eq. (39). Both for double charm and bottom baryons we find results compatible with most of the lattice references when uncertainties are taken into account. In the case of double bottom baryons we find a remarkably close agreement. For double charm baryons, the large array of states computed on the lattice in Ref. [36] allows us to compare the general features of the spectrum finding a good qualitative agreement.
The light flavor dependence of the double heavy baryon spectrum enters through the lattice data from which the potentials are obtained. The data of Refs. [15,16] that we use is for degenerate u and d quarks corresponding to a pion mass of 783 MeV. Data at more realistic pion masses would be highly desirable. The EFT used in this paper may be taken in the chiral limit and therefore our results would be independent of the flavor of the light quark forming the double heavy baryon. The implementation of chiral symmetry breaking to this EFT is straightforward and can be done using the standard techniques from Heavy Baryon Chiral Perturbation Theory [60]. However, the couplings with the pseudo-Goldstone bosons would introduce a number of unknown potentials that limit the predictability of the EFT, unless data for different low enough values of the pion mass become available for the potentials. Interestingly, upon implementing chiral symmetry to this EFT and the expansion of the fields in the eigenstates of the Schrördinger equations, our EFT will match the standard chiral hadronic theories. These EFTs have been used abundantly to discuss threshold effects and decays [43,44,61].
The main obstacle in the development of the EFT presented in this paper is the dependence of the potentials in nonperturbative dynamics. In an accompanying paper [14] we have presented the matching of the potentials of the EFT of Sec. II in terms of insertions of NRQCD operators in the Wilson loop. These Wilson loops are suitable to be computed on the lattice with light quarks and gluons only, sidestepping the issue of having to deal with widely separated scales on the lattice. Therefore, combining lattice QCD and EFT the major obstacles of both approaches can be avoided.
The projectors used in the Lagrangians in Eqs. (12) and (14) can be constructed from the projection vectors defined in Eq. (A1).
We can diagonalize the potential matrix for the spin-3/2 field in Eq. (20) by projecting the field into a basis of eigenstates ofr · S 3/2 Note that the projection vectors act onto the light quark spin inside of the Ψ (3/2) − field. In this basis the LO Lagrangian for the Ψ (3/2) − field can be written as Now we want to obtain the matrix P α † 3 2 λ L 2 QQ P α 3 2 λ using the techniques of Refs. [59,62]. Let us start with the following auxiliary formulas. The angular momentum operator in spherical coordinates is The unit vectors in spherical coordinates read r = (sin(θ) cos(φ), sin(θ) sin(φ) , cos(θ)) , θ = (cos(θ) cos(φ), cos(θ) sin(φ) , − sin(θ)) , Let us compute the commutators of the angular momentum operator and the unit vectors in spherical coordinates from which one can obtain Therefore, From an explicit computation we obtain and P α † and we note that If we define L = L QQ + S 3/2 , from Eqs. (A17) and (A16) we can write We are in position now to compute the matrix elements of the square of the heavy-quark angular momentum in between the projector vectors The square of the right-hand side of Eq. (A19) is which is the operator whose eigenfunctions are our angular wave functions, with |m| < and |λ| < . The operators K ± act as the λ-raising and -lowering operators for the angular wave functions v λ m , The mixing terms in Eq. (A20) can be written as a matrix in the λ-λ indices as Introducing Eq. (A27) into Eq. (A20) we arrive at Now, let us look at the expected value of the square angular momentum operator between the angular wave functions The system decouples into two two-state coupled equations in the basis given by the following transformation matrix: Therefore, we arrive at two sets of coupled Schrödinger equations where ψ η P = ±(−1) −1/2 (A44)

Boundary conditions
A system of two linearly coupled differential equations of second order has in general four linearly independent solutions of which two will be singular at the origin. For the case of only one second order differential equation ( = 1/2) there exists two linearly independent solutions, one of which will be singular at the origin. These solutions can be distinguished by their behavior at the origin which will also be important as initial conditions for the numerical solution of the Schödinger equations. To obtain the behavior at the origin of the solutions of the differential equations we note that in the short distance the kinetic term dominates over the potential. In this limit the system can be diagonalized and the two resulting decoupled equations solved by a guess function r b . Two solutions will be obtained for each equation, one of which diverges at the origin and should be discarded. The remaining one for each equation gives us one of the solutions for the system with the weights given by the components of the eigenvectors that diagonalized the system. These are the behaviors at the origin of the good solutions for ≥ 3/2 ψ 3/2+ ψ 1/2+ ∝ − 12 ( + 1) − 9r +1/2 (2l + 3)r +1/2 , (A45)