The Born-Oppenheimer approximation in an effective field theory language

The Born--Oppenheimer approximation is the standard tool for the study of molecular systems. It is founded on the observation that the energy scale of the electron dynamics in a molecule is larger than that of the nuclei. A very similar physical picture can be used to describe QCD states containing heavy quarks as well as light-quarks or gluonic excitations. In this work, we derive the Born--Oppenheimer approximation for QED molecular systems in an effective field theory framework by sequentially integrating out degrees of freedom living at energies above the typical energy scale where the dynamics of the heavy degrees of freedom occurs. In particular, we compute the matching coefficients of the effective field theory for the case of the $H^+_2$ diatomic molecule that are relevant to compute its spectrum up to ${\cal O}(m\alpha^5)$. Ultrasoft photon loops contribute at this order, being ultimately responsible for the molecular Lamb shift. In the effective field theory the scaling of all the operators is homogeneous, which facilitates the determination of all the relevant contributions, an observation that may become useful for high-precision calculations. Using the above case as a guidance, we construct under some conditions an effective field theory for QCD states formed by a color-octet heavy quark-antiquark pair bound with a color-octet light-quark pair or excited gluonic state, highlighting the similarities and differences between the QED and QCD systems. Assuming that the multipole expansion is applicable, we construct the heavy-quark potential up to next-to-leading order in the multipole expansion in terms of nonperturbative matching coefficients to be obtained from lattice QCD.


I. INTRODUCTION AND MOTIVATION
The discovery in the last decade of the XY Z mesons has brought into QCD challenges enduring since the early days of molecular physics in QED -for a recent overview, see Ref. [1]. A great variety of possible models have been introduced to explain the observed pattern of new mesons. A recent proposal [2,3] (see also [4]) advocates the use of the Born-Oppenheimer (BO) approximation [5][6][7][8], familiar to QED molecular physics, as a starting point for a coherent description of the new QCD structures. The rational for this being that many of the new mesons contain a heavy quark-antiquark pair, and the time scale for the evolution of the gluon and lightquark fields is small compared to that for the motion of the heavy quarks. Although the BO approximation has been used in the past to study heavy hybrids by means of quenched lattice data for gluonic static potentials [9][10][11] 1 , the new aspect of the proposal in Refs. [2,3] is the * nora.brambilla@ph.tum.de † gkrein@ift.unesp.br ‡ jaume.tarrus@tum.de § antonio.vairo@ph.tum.de 1 Models have also been used for the determinations of the gluonic static potentials and heavy hybrids in a BO framework, see for example Refs. [12,13] recognition that the BO approximation can also be applied to mesons with light quark and antiquark flavors when input from lattice simulations becomes available.
In the present paper we go one step further in this proposal and develop an effective field theory (EFT) that allows to calculate in a systematic and controlled manner corrections to the BO approximation for QED and QCD molecular systems. An EFT is built by sequentially integrating out degrees of freedom induced by energy scales higher than the energy scale we are interested in. For QED molecules, such a sequential process proceeds as follows: (A) integrating out hard modes associated with the masses of the charged particles leading to nonrelativistic QED (NRQED) [14,15], (B) integrating out soft modes associated with the relative momenta between electrons and nuclei in NRQED leading to potential NRQED (pN-RQED) [16,17], and (C) exploiting the fact that the nuclei move much slower than the electrons due to their heavier masses, modes associated with the electron and photon dynamics at the electron binding energy scale, the ultrasoft scale, can be integrated out leading to an EFT for the motion of the nuclei only. In QED these steps can be done in perturbation theory.
In the present paper we compute this ultimate EFT in the simple case of a QED molecule formed by two heavy nuclei and one electron, like the H + 2 ion molecule. Because the BO approximation emerges as the leading-order approximation in this EFT, we call it Born-Oppenheimer EFT (BOEFT). Furthermore we show how the EFT allows to systematically improve on the leading-order approximation by calculating corrections in the inverse of the mass of the nuclei as well as electromagnetic corrections. We give explicit analytical expressions, regularized in dimensional regularization when needed, for the different contributions to the binding energy of the two nuclei plus one electron molecule up to O(mα 5 ). It is at this order that the Lamb shift is generated.
The BOEFT that we construct is new, although NRQED has been applied in atomic and molecular physics for nearly two decades [15,18]. In particular, NRQED has been used for computing the leading relativistic, recoil and radiative corrections to the energy levels of the H + 2 molecule in Ref. [19] and for computing higher-order corrections in Refs. [20][21][22][23][24]. The new and distinctive aspect of our approach is that we carry out the full EFT program for the diatomic molecule, integrating out not only the hard scale, as in NRQED, but also the soft and ultrasoft scales. The advantage is that each term in the Lagrangian has a unique size and the scaling of Feynman diagrams is homogeneous. This greatly facilitates the determination of all the relevant contributions to a given observable up to a given precision, a feature that is particularly useful for higher-order calculations.
An analog EFT for QCD states containing a heavy quark-antiquark pair in a color-octet state bound with light quarks or a gluonic color-octet state can be built following a similar path. However, unlike QED molecules, the QCD states are determined by nonperturbative interactions. The hard scale set by the heavy-quark mass can always be integrated out perturbatively, leading to nonrelativistic QCD (NRQCD) [14,25]. At short enough distances the relative momentum of the heavy quarks can also be integrated out perturbatively resulting in potential nonrelativistic QCD (pNRQCD) [16,[26][27][28]. 2 . Similarly to the diatomic molecule case, the heavy quarks move slower than the light degrees of freedom, whose spectrum is assumed to appear at the scale Λ QCD . Thus, one can construct an EFT for these QCD "molecular" states by integrating out the scale Λ QCD . Since this is the scale of nonperturbative physics, the matching coefficients will be nonperturbative quantities to be determined, for instance, by lattice calculations. When light quarks are neglected, one regains in this way the EFT recently constructed for quarkonium hybrids [32].
The paper is organized as follows. In Sec. II we construct the pNRQED Lagrangian for two nuclei and one electron. In Sec. III we proceed with integrating out the ultrasoft scale and constructing the molecular EFT, BOEFT. Section IV is devoted to the power counting of the BOEFT, which we use to assess the importance of the nonadiabatic coupling and other corrections to the molecular energy levels. The EFT for the QCD analog of the diatomic molecule, quarkonium hybrids and tetraquark mesons built out of a heavy quark and antiquark, is developed in Sec. V. Section VI contains the conclusions and an outlook for future developments. The Appendix presents a detailed calculation of the Lamb shift for the H + 2 molecule.

II. pNRQED
We aim at building an EFT for a molecular system containing heavy and light particles: the heavy particles (nuclei) have electric charge +Ze and mass M and the light particles (electrons) have electric charge −e and mass m, with M ≫ m. Both kinds of particles are nonrelativistic. Such a molecular system has several wellseparated energy scales, as we will see more in detail in the following. From the highest to the lowest one the relevant energy scales are the masses of the heavy and light constituents (hard scales), the typical relative momentum p = p ∼ mv between heavy and light particles (soft scale) and the binding energy of the light particles E ∼ mv 2 (ultrasoft scale). For a Coulomb-type interaction it holds that v ∼ α with α = e 2 4π ∼ 1 137 the fine structure constant. Finally, specific of molecules an extra low-energy scale appears: the binding energy of the heavy nuclei.
The EFT suitable for describing QED bound states at the ultrasoft scale is pNRQED. In Ref. [17] it was worked out for the hydrogen atom, in this section we extend pNRQED to describe systems with two nuclei and one electron. In Sec. III we will integrate out the ultrasoft modes and build the EFT suitable to describe the molecular states.
The Lagrangian of pNRQED can be written in terms of the light and heavy fermion fields, ψ(t, x) and N (t, x) respectively, and the ultrasoft-photon field, A µ (t, x). The meaning of A µ (t, x) being ultrasoft is that it must be multipole expanded (e.g., about the position of the center of mass (c.m.) of the constituents). The operators of the pNRQED Lagrangian can be organized in an expansion in α and m M . In order to homogenize the counting in these two expansion parameters, we will use that m M is numerically similar to ∼ α 3 2 . Then, the pNRQED Lagrangian relevant to compute the spectrum up to order O(mα 5 ) reads where F µν = ∂ µ A ν − ∂ ν A µ and all photons are ultrasoft. Moreover we have used where D q is the covariant derivative, with q = −e for the electron and q = +Ze for the nuclei: The electron-nucleus potential V Ze (x, σ) is given by where LO (leading order) and NLO (next-to-leading order) refer to the order mα 2 and mα 4 contributions to the spectrum respectively. The LO potential is the Coulomb potential while the NLO one is the sum of a contact and spin-orbit interaction with where c D , c S and d 2 are matching coefficients that up to order α read The coefficient c D has been renormalized in the MS scheme. The scale µ is the dimensional regularization scale that in the case of c D acts as an infrared factorization scale. Finally, the V ZZ potential in Eq. (1) contains the LO nucleus-nucleus Coulomb potential: Further contributions to (5) and (11), which can be found in Ref. [33], are beyond our accuracy. Next, we project the Lagrangian in Eq. (1) on the subspace of one electron and two nuclei. This is similar to the pNRQED bound state calculations for the hydrogen atom [16,17], but since the projection for one light and two heavy particles with different charges has not been done so far in the literature, we present the procedure with some detail. The subspace of one electron and two nuclei is spanned by Fock-space states of the form where ϕ(t, x, y 1 , y 2 ) is the wave function of the system and US⟩ is the Fock-space state containing no hard particles (electrons or nuclei) and an arbitrary number of ultrasoft ones (photons). The corresponding projected Lagrangian, adequate for calculating the spectrum up to O(mα 5 ), is where we have promoted ϕ(t, x, y 1 , y 2 ) to a tri-local field.
To ensure that the photon fields A µ are ultrasoft one may multipole expand them about the c.m. of the system. The task is facilitated by defining an appropriate c.m. and relative coordinates. The c.m. coordinate R of the system is given by To describe the motion of the electron relative to the positions y 1 and y 2 of the nuclei we use and for the relative coordinate of the nuclei The multipole expansion spoils manifest gauge invariance. It is important, however to recall that we have an EFT for ultrasoft gauge fields, hence gauge transformations must not introduce into the EFT gauge fields with large-momentum components; that is, the allowed gauge transformations are those that produce fields that still are within the EFT. One can recover manifest (ultrasoft) gauge invariance at least for charge neutral systems by introducing the field redefinition: where U q is the Wilson line Under a gauge transformation A 0 (t, R) → A 0 (t, R) − ∂ t θ(t, R) and A(t, R) → A(t, R) + ∇ R θ(t, R), the field S(t, R, r, z) transforms as where e tot is the total charge: For a charge-neutral system, e tot = 0, and the field S(t, R, r, z) is gauge invariant.
The Lagrangian in terms of the field S is given by where with M tot being the total mass is the electric field and e eff is the effective charge: The sizes of the different terms that appear in the Lagrangian (21) are as follows.
1. Relative electron-nuclei momentum −i∇ z and inverse relative distance 1 z have size mα.
2. Photon fields, derivatives acting on photon fields, the time derivative, and c.m. momentum, −i∇ R , acting on S have size mα 2 .
3. As we shall discuss in Sec. IV, the inverse relative nuclei-nuclei distance is 1 r ∼ mα, whereas the radial part of the derivative ∇ r ∼ (M m) 1 4 mα ∼ mα 5 8 when acting on the nuclei, but ∇ r ∼ mα when acting on the electron cloud. This implies that the kinetic energy associated with the relative motion of the nuclei is −∇ 2 r M ∼ mα 2 m M ∼ mα 11 4 .
Using this counting, and disregarding operators that produce emission or absorption of photons that contribute only in loops, the leading-order operators in Eq. (21) are h 0 (r, z) + V LO ZZ (r), which are of O mα 2 . Since the kinetic energy associated with the relative motion of the two nuclei, −∇ 2 r M , is of O mα 11 4 , at leading order the nuclei are static and V LO ZZ (r) is just a constant. Therefore, at leading order, the Euler-Lagrange equation from the Lagrangian (21) is nothing else than a Schrödinger equation for the electronic energy levels with Hamiltonian h 0 (r, z). Corrections to these energy levels can be obtained in perturbation theory. Parametrically, the first of such corrections is given by the recoil term, ∇ 2 z 4M , which is O mα 7 2 , and the second one by ∇ 4 z 8m 3 + V NLO Ze , which starts at O mα 4 . The O mα 5 corrections include the Lamb shift, and originate from ultrasoft photon loops and subleading contributions to the NLO potentials.
To obtain the molecular energy levels we need to solve the dynamics of the r coordinate. In principle we could do this by adding subleading terms to the Hamiltonian, . . , and solving the corresponding Schrödinger equation. However, in this paper, following the logic of EFTs, we will integrate out from pNRQED the ultrasoft degrees of freedom to obtain an EFT at the energy scale of the two-nuclei dynamics. The Euler-Lagrange equation of this EFT provides a Schrödinger equation for the molecular energy levels. We will develop this EFT, which we call BOEFT, in the following section.
Since the c.m. motion does not affect the internal dynamics of the molecule, we can simply work in the c.m. frame and ignore the dependence on R of the field S. We also use the notation A 0 (t, 0) and E(t, 0) to indicate quantities defined at the origin of the coordinate system, i.e., R = 0.

III. BORN-OPPENHEIMER EFT FOR DIATOMIC MOLECULES
Our purpose is to build the BOEFT, an EFT for the diatomic molecule at the energy scale of the two-nuclei dynamics. This EFT is obtained by integrating out the ultrasoft scale, mα 2 , from pNRQED for two nuclei and one electron given in Sec. II. We will include effects that contribute to the binding energy of the molecule up to O(mα 5 ).
Since the electron dynamics occurs at the ultrasoft scale, integrating out this scale entails that all the electronic degrees of freedom are integrated out. Moreover, also ultrasoft photons are integrated out. Therefore, the degrees of freedom of the BOEFT are nuclei and photons with energies of O mα 11 4 or smaller.
The tree-level matching contributions can be easily obtained by expanding the field S(t, r, z) in the pNRQED Lagrangian of Eq. (21) in eigenfunctions of the leadingorder Hamiltonian h 0 (r, z) of Eq. (22). This corresponds in expanding the field S(t, r, z) as where φ κ (r; z) = ⟨z r, κ⟩ satisfy the electronic eigenvalue equation The eigenvalues V light κ (r) are the static energies, with κ representing the set of quantum numbers spec-ifying the electronic state for a fixed separation r of the nuclei. The r in the state vector r, κ⟩ emphasizes that eigenvalues labeled by κ refer to a given nuclei separation r. The eigenfunctions φ κ (r; z) are orthonormal: The static electronic energies V light κ (r) scale like mα 2 .
The set of quantum numbers κ is familiar from molecular physics and corresponds to representations of the symmetry group of a diatomic molecule [34]: the eigenvalue λ = 0, ±1, ⋯ of the projection of the electron angular momentum on the axis joining the two nuclei,r, traditionally denoted by Λ = λ and conventionally labeled by Σ, Π, ∆, . . . for Λ = 0, 1, 2, . . . ; the total electronic spin S, with the number of states (multiplicity) for a given S being 2S + 1, and indicated with an index, like 2S+1 Σ; additionally, for the Σ state, there is a symmetry under reflection in any plane passing through the axisr, the eigenvalues of the corresponding symmetry operator being ±1 and indicated as Σ ± ; and, in the situation of identical heavy nuclei, the eigenvalues ±1 of the parity operator of reflections through the midpoint between the two nuclei, denoted by g = +1 and u = −1. 3 In this way, a possible ground state is denoted by κ = 1 Σ + g . The tree-level matching is sufficient up to terms in the Lagrangian of O(mα 4 ). Ultrasoft photon loops start contributing at O(mα 5 ) and are responsible for the Lamb shift of the diatomic molecule. We detail the calculation of the leading ultrasoft loop in Appendix A.
The BOEFT Lagrangian up to O(mα 5 ) reads The photon fields carry energies and momenta of O mα 11 4 or smaller. The operator H (0) κ is the leading-order nuclei-nuclei Hamiltonian: and δE κ (r) is the sum of the tree-level and second order recoil, Breit-Pauli corrections as well as the one-loop ultrasoft one: The counting of H (0) κ will be justified in the next section, but we have already anticipated that the eigenvalues of H The different contributions to δE κ (r) read which is of order mα 2 m M ∼ mα 7 2 , which is of order which starts at order mα 4 , and where and ρ κ (r) is the electron density at the positions of the nuclei The ultrasoft contribution is of order mα 5 log(α) and mα 5 . Note that the ultrasoft contribution has been renormalized in the MS scheme and its µ dependence cancels against that one of the matching coefficient c D [see Eq. (10)] in the NLO potential of Eq. (33). Finally, C nad κκ ′ (r) is the nonadiabatic coupling [8,35]: The first integral in the second line is the matrix element of the kinetic energy operator of the relative motion of the nuclei, it is of order mα 2 m M ∼ mα 7 2 , and the second integral involves the momentum of their relative motion, it is of order mα 2 (m M ) 3 4 ∼ mα 25 8 . When the φ κ 's are real and κ = κ ′ , the second integral vanishes. We conclude by commenting on some general features of the BOEFT. First, we would like to notice that there is no extra approximation by writing S(t, r, z) as in Eq. (25), since the eigenfunctions φ κ (r; z) form a complete set and the Ψ κ (t, r) play the role of time-dependent expansion coefficients. However, as it is well-known in treatments employing the Born-Oppenheimer approximation, this is useful in practice only when the dynamics of the heavy degrees of freedom (with mass M ) is much slower than the dynamics of the light degrees of freedom (with mass m), a feature that permits to define an adiabatic dynamics for the heavy particles and to treat departure from adiabaticity using perturbation theory in the small parameter m M ≪ 1, as we have done above. Otherwise, when M ≃ m, the concept of adiabatic motion for one of the particles loses sense and an expansion like Eq. (25) would be useless. A way to see this is by noticing that mixing terms in the energy levels of the BOEFT would count like mα 2 , a fact that would prevent the separation of the electron from the nuclei dynamics.
Under the adiabatic assumption the molecular energy levels are distributed as sketched in Fig. 1. Electronic excitations define for each nuclei separation a potential V light κ (r). These potentials are separated by large gaps of order mα 2 . For each electronic excitation, the nuclei motion induces smaller excitations of order mα 2 m M . We can compute these smaller excitations in the BOEFT for each electronic potential V light κ (r). They are at lead- ing order the eigenvalues of H κ . It is astounding that the wave functions of these nuclear vibrational modes can not only be computed but experimentally directly visualized: for the H + 2 ground state potential V light 0 (r) see [36].

IV. POWER COUNTING IN THE BOEFT
In this section we examine in detail the power counting of the BOEFT that we have just developed. The main aim is to substantiate the starting assumption in the construction of the BOEFT, namely that the kinetic term −∇ 2 r M ≪ mα 2 . Also of interest is the size of the nonadiabatic coupling.
The derivative ∇ r can act on the nuclei fields Ψ κ (t, r) as well as on the electronic wave functions φ κ (r; z). The size of the derivative turns out to be different for nuclei and electrons. In the case of ∇ r acting on φ κ (r; z), it scales like ∼ mv. Since the electron is bound to the nuclei through Coulomb interactions, we have that v ∼ α. In the case that the derivative acts on Ψ κ (t, r), it scales like ∼ M w, where w is the relative velocity of the nuclei. Therefore, our goal is to asses the size of w.
Since the system is bound, the nuclei will have a stable equilibrium arrangement and oscillate around an average separation r 0 . Without the electron the two nuclei would not form a bound state, hence r 0 is an emergent scale, whose size needs to be determined. Let us consider the ground-state electron energy (κ = 0) and expand the total potential V (r) = V LO ZZ (r) + V light 0 (r) around the equilibrium position r 0 (we have adjusted the potential so that its minimum is zero): The Hamiltonian of the relative motion is that of a harmonic oscillator. The ground-state energy E 0 is given by The equilibrium position r 0 of the nuclei is determined from Because V light 0 (r 0 ) is the ground state energy of Eq. (26), it is of order mα 2 (O(Z 2 ) ∼ 1). Hence Eq. (40) implies That is, the average size of the nuclei separation is of the same order as the electron-nucleus separation. Clearly, this is a particular feature of the Coulomb interaction between the nuclei; for a different r dependence of the nucleus-nucleus interaction, r 0 may be not of the order of the Bohr radius. From the above result it follows that and that the ground-state vibrational energy is Transitions between low-lying vibrational states are also of order mα 2 m M . We note that the scaling behavior of E 0 implies a large cancellation between V LO ZZ (r) and V light 0 (r) near the equilibrium position, since each of these two potentials scales like mα 2 .
The virial theorem for the harmonic oscillator relates the expectation value of the kinetic energy with the total energy, from where the size of the kinetic-energy operator acting on Ψ follows Our initial assumption was that the kinetic energy associated with the relative motion of the nuclei is small compared to the ultrasoft scale, from there we integrated out the latter and matched pNRQED to the BOEFT. The above analysis shows that the energy scale associated with the relative motion of the nuclei is indeed largely suppressed by a factor m M ∼ α 3 4 ≈ 0.025 with respect to the ultrasoft scale, which justifies the initial assumption. The size of ∇ r acting on Ψ and the relative velocity of the nuclei follows from (45): w ∼ α m M 3 4 .
A more detailed look reveals, however, that the counting of Eq. (46) applies only to the radial component of ∇ r . Indeed, in spherical coordinates we have ∇ r = (∂ r , ∂ θ r, ∂ φ (r sin θ)), and since the angles are dimensionless variables, the size of the last two components is determined by r ∼ r 0 ∼ 1 (mα). This implies also that the counting (45) is appropriate for the radial part of the kinetic energy, whereas −2 (M r) ∂ ∂r ∼ mα 2 (m M ) 3 4 and the angular part L 2 (M r 2 ) scales like mα 2 (m M ).
The size of the kinetic term in Eq. (45) sets the energy scale for the BOEFT. Hence it determines the scaling of photon fields and derivatives acting on them. The last ingredient to complete the counting rules for the BOEFT is the scaling of ∇ z ∼ 1 z ∼ mα, which is inherited from pNRQED of Sec. II. The molecular energy scales are summarized in Fig. 2. We apply now the counting rules to the nonadiabatic coupling C nad (r) defined in (37). The largest contribution comes from the radial piece of the second term, which is of O mα 2 (m M ) 3 4 , while the first term and the angular piece of the second one are O mα 2 (m M ) . Therefore, at leading order the nonadiabatic coupling can be neglected and the equation of motion for the field Ψ κ (t, r) reads which is nothing else than the Schrödinger equation that describes the motion of the heavy particles in the Born-Oppenheimer approximation [5][6][7]. Equation (48) produces the leading-order energy eigenvalues for the diatomic molecule, but it does not describe well the angular wave functions [8]. This is a consequence of the angular piece of the kinetic term being of the same size as the angular parts of C nad κκ . The adiabatic approximation [8,35] corresponds to including in the above Schrödinger equation the diagonal term C nad κκ (r) One can use an iterative procedure to solve the problem: starting from the zeroth-order solution in which the nonadiabatic coupling C nad is neglected, one can treat C nad as a perturbation [37] since its contribution to the energy is suppressed by an amount (m M ) 1 4 ≈ 0.15 with respect to the zeroth-order energy. We emphasize again that this relies on the Coulomb nature of the nucleusnucleus interaction and on the smallness of the ratio m M . Let Ψ The leading-order correction E κn comes from the diagonal nonadiabatic coupling and reads It is of order mα 2 (m M ) 3 4 ∼ mα 25 8 . The nondiagonal nonadiabatic coupling provides mixing with different electronic excitations. The first contribution appears at order mα 2 (m M ) 3 2 ∼ mα 17 4 and reads More important than the mixing with states belonging to different electronic excitations is the mixing with states in the same one. The mixing is in this case suppressed by a mere factor (m M ) 1 4 ∼ α 3 8 . We will not display here explicitly this kind of contributions that follow straightforwardly from time-independent quantum-mechanical perturbation theory. We add that the recoil corrections to the electronic levels (31) and (32) contribute first at order mα 2 (m M ) ∼ mα 7 2 and mα 2 (m M ) 2 ∼ mα 5 respectively. Finally, the NLO corrections to the electronic levels (33) contribute first at order mα 4 , while the ultrasoft corrections (34) contribute first at order mα 5 log(α) and mα 5 . Let us now summarize the steps necessary for a numerical evaluation of the molecular energy levels using the BOEFT. First, the electronic static energies V light κ and wave functions φ κ are obtained by solving the eigenvalue equation (26) (see, for example, Ref. [38]). The BOEFT matching coefficients in Eqs. (31)- (34) and (37) can then be evaluated. The nuclei wave functions Ψ

V. THE BOEFT FOR QCD: HEAVY HYBRIDS AND ADJOINT TETRAQUARK MESONS
In the context of QCD, it exists a system analog to the QED diatomic molecule. It is the system formed by a heavy quark-antiquark pair and some light degrees of freedom that can be either gluonic or light quark in nature. Similarly to the QED bound state, the QCD system develops three well separated energy scales: the heavyquark mass M (hard scale), the relative momentum M w (soft scale), where w is the heavy-quark relative velocity, and the binding energy M w 2 . Furthermore, there is the scale associated with nonperturbative physics, Λ QCD that plays the role of the ultrasoft scale in the hadronic case. Restricting ourselves to the case M w ≫ Λ QCD , we can use weakly-coupled pNRQCD [16,27] to describe the heavy quark-antiquark pair, which is called quarkonium if bound, pretty much in the same way as pNRQED, described in Sec. II, can be used to describe electromagnetic bound states. However, a situation that has no analog in pNRQED, the heavy quark-antiquark fields can appear in pNRQCD either in a color-octet or in a color-singlet configuration.
At energies of the order of Λ QCD , the spectrum of QCD is formed by color-singlet hadronic states that are nonperturbative in nature. An interesting case it that one of exotic hadrons made of a color-octet heavy quarkantiquark pair bound with light degrees of freedom. Such a system can be studied similarly to the QED diatomic molecules. The heavy quarks play the role of the nuclei and the gluons and light quarks play the role of the electrons.
In a diatomic molecule the electrons are nonrelativistic with energies of the order of the ultrasoft scale, mα 2 , whereas, as we have seen, the nuclei have a smaller energy due to their heavier mass. In a hadron made of a color-octet heavy quark-antiquark pair, the light degrees of freedom are relativistic with a typical energy and momentum of order Λ QCD . This implies that the typical size of the hadron is of the order of 1 Λ QCD . If the mass of the heavy quarks is much larger than Λ QCD , there may be cases where also the typical momentum M w of the heavy quarks in the hadron is larger than Λ QCD . The scaling of the typical distance of the heavy quarkantiquark pair depends on the details of the full interquark potential, which has a long-range nonperturbative part and a short-range Coulomb interaction. It may therefore happen that the heavy quark and antiquark are more closely bound than the light degrees of freedom. This situation is interesting because the hadron would present a hierarchy between the distance of the quarkantiquark pair and the typical size of the light degrees of freedom that does not exist in the diatomic molecular case where the electron cloud and the two nuclei have the same size. A consequence of this is that while the molecule is characterized by a cylindrical symmetry, the symmetry group of the hadron would be a much stronger spherical symmetry at leading order in a (multipole) ex-pansion in the distance of the heavy quark-antiquark pair. This modifies significantly the power counting of the hadronic BOEFT with respect to the molecular one leading to new effects. In order to emphasize the difference between the hadronic and molecular case, we will assume in the following that the typical distance between the heavy quark and antiquark is of order 1 (M w).
The kinetic energy associated with the relative motion of the quark-antiquark pair scales like M w 2 . If we look at hadrons that are in the ground state or in the first excited states only, we may require that M w 2 ≪ Λ QCD . As we have seen discussing the diatomic molecule, in order for a Born-Oppenheimer picture to emerge and for the BOEFT to provide a valuable theory it is crucial that the excitations between the heavy particles happen at an energy scale that is smaller than the energy scale of the light degrees of freedom. In summary, we will require the following hierarchy of energy scales to hold true: M w ≫ Λ QCD ≫ M w 2 [27]. The different energy scales are shown in Fig. 3. After integrating out the hard and soft scales from QCD and projecting on quarkonium states, one arrives at the pNRQCD Lagrangian in the weakly-coupled regime, which at leading order in 1 M and at O(r) in the multipole expansion is (we neglect the light-quark masses and higher-order radiative corrections to the dipole operators) where S and O are the heavy quark-antiquark colorsinglet and color-octet fields respectively normalized with respect to color. They depend on t, r, the relative coordinate, and R, the c.m. position of the heavy quarkantiquark pair. All the fields of the light degrees of freedom in Eq. (54) are evaluated at R and t; in particular, G µν a = G µν a (R, t), q i = q i (R, t) and The field E is the chromoelectric field, G µν a the gluonic field strength tensor and q i are light-quark fields appearing in n f flavors. The singlet and octet Hamiltonians read (in the c.m. frame) where V s (r) = −4α s (3r) + . . . and V o (r) = α s (6r) + . . . are the color-singlet and color-octet potentials respectively; α s is the strong coupling. The Lagrangian (54) is the analog of the Lagrangian (21) for diatomic molecules. The difference is that in the Lagrangian (54) the number of gluons and light quarks is not fixed as the number of electrons is in (21). This stems from the fact that the electrons are nonrelativistic, which implies that their number is conserved at the low energy of pNRQED, while gluons and light quarks are massless relativistic particles and thus their creation and annihilation are still allowed in the Lagrangian (54).
The Hamiltonian density corresponding to the light degrees of freedom at leading order in 1 M and in the multipole expansion is It plays the same role as the Hamiltonian density of Eq. (22) does for the diatomic molecule. As anticipated, the symmetry groups of the two Hamiltonians are, nevertheless, different: the Hamiltonian density in Eq. (22) has a cylindrical symmetry, while Eq. (57) has a spherical symmetry. The color-octet G ia κ (R) operators that generate the eigenstates of h 0 (R) form a basis of octet light degrees of freedom operators, labeled by the light-flavor f and J P C quantum numbers, and an extra label i for states belonging to the same J P C representation. Note that the energy eigenvalue Λ κ is in general a complex number, whose imaginary part accounts for the possible decay of the state. If we introduce the states which are eigenstates of the octet sector of the pNRQCD Hamiltonian at leading order in the multipole expansion with eigenvalues h o + Λ κ , we can now project the Lagrangian of (54) onto the Fock subspace spanned by This step is the equivalent for the hadronic system to the projection on the state of Eq. (12) and the expansion (25) for the diatomic molecule. Using Eq. (61) and integrating out light degrees of freedom of energy of order Λ QCD we derive the BOEFT Lagrangian that describes the heavy quark-antiquark pair physics at the scale M w 2 . Since we are interested in bound states we will not consider sectors of the Lagrangian that describe transitions between states with different κ and decays into singlet states. Up to next-toleading order in the multipole expansion the Lagrangian reads where P i κλ are projection operators along the heavyquark axis of the light degrees of freedom operator (an implicit sum is understood over repeated i, j indices). There is one projection operator for each − j ≤ λ ≤ j . These operators select different polarizations of the wave function Ψ iκ . For example, in the case of J = 1 the operators are given by withr = (sin(θ) cos(φ), sin(θ) sin(φ) , cos(θ)) T , θ = (cos(θ) cos(φ), cos(θ) sin(φ) , − sin(θ)) T , For higher J the projection operators can be built by multiplying j powers of (63) and (64) with appropriate symmetrization of the indices (see also [39]). The projection operators are necessary to organize the states in Eq. (60) according to the quantum numbers of the exotic hadron. In particular they project the light degrees of freedom operator onto the heavy quark-antiquark axis. The quantum numbers of the exotic hadron are the same as the ones of the diatomic molecule presented in Sec. III plus charge conjugation: as we discussed, at leading order in the multipole expansion the symmetry of the hadron is spherical, hence the projectors commute with the eigenstates of h 0 (the equivalent statement is not true in the molecular case), but higher-order terms break this symmetry to the original cylindrical one. In Eq. (62), the next-to-leading order term in the multipole expansion is P i κλ b κλ r 2 P j † κλ , whereas the dots stand for higher-order terms.
The specific value of the next-to-leading-order term, P i κλ b κλ r 2 P j † κλ , depends on nonperturbative physics and is unknown, however some of its characteristics can be determined on general grounds. This term has its origin in the chromoelectric dipole interactions of Eq. (54), which couple the light degrees of freedom operator G ia κ to the octet field giving corrections to the (static) energy of the system. That this kind of corrections shows up for the static energy is a specific feature of QCD [26,27], however, for nonstatic nuclei dipole interactions are also responsible for the Lamb shift of the diatomic molecule, as we have seen. The r 2 dependence arises from the necessity of having at least two chromoelectric dipoles in order to conserve the J P C quantum numbers of G ia κ . Cylindrical symmetry and charge conjugation also imply b κλ = b κ−λ = b κΛ . In Fig. 4 we show static potentials for the case of quarkonium hybrids, that is, for the case in which the considered light degrees of freedom are purely gluonic. The potentials correspond to κ = 1 +− and are compared to the static energies computed on the lattice in the quenched approximation. The values of b κλ are fitted to the lattice data for r ≲ 0.5fm.  Figure 4. Comparison of the hybrid quarkonium static energies generated by the lowest mass gluelump (κ = 1 +− ) computed on the lattice in Refs. [40] (red squares) and [41] (green dots) compared to the BOEFT static potential up to next-toleading-order (solid black line), V κλ = Vo(r) + Λκ + b κλ r 2 . The octet potential is taken in the Renormalon Subtracted (RS) scheme and up to α 3 s . The mass of lowest laying gluelump is computed also in the RS scheme Λ RS 1 +− = 0.87 GeV [40]. The b κλ coefficients are fitted to the lattice data for r ≲ 0.5 fm yielding the values b10 = 1.112 GeV/fm 2 and b1±1 = 0.110 GeV/fm 2 . For lattice determinations of higher laying gluelump masses and static energies see Refs. [9,10,[41][42][43][44][45][46][47].
Defining the projected wave function as and using we can rewrite Eq. (62) as The last term can be split into a kinetic operator acting on the heavy quark-antiquark field and a nonadiabatic coupling with being the nonadiabatic coupling analog to Eq. (37) for the diatomic molecule. At this point it is important to review the sizes of the different terms appearing in Eq. (68). All dimensional quantities that arose from integrating out Λ QCD are of order Λ QCD to their dimension. Hence Λ κ is of order Λ QCD and b κλ is of order Λ 3 QCD . The temporal derivative, the kinetic term and the potential up to the constant shift Λ κ are of order M w 2 . Unlike in the diatomic molecule case, ∇ r has the same size for radial and angular pieces, because the momentum of the heavy quark is taken to scale like the inverse of the distance, r, between the quark and the antiquark. For the nonadiabatic coupling C nad κλλ ′ , the radial piece of the derivative ∇ r acting on the projection operators P i κλ ′ vanishes, since they do not depend on r . According to our counting, the size of the angular piece L 2 (M r 2 ), P i κλ ′ is M w 2 , i.e., of the same order as the kinetic operator of the heavy quarks. This is different from the diatomic molecular case.
The equations of motion for the fields Ψ κλ (t, r, R) that follow from the Euler-Lagrange equation at leading order are nothing else than a set of coupled Schrödinger equations By solving them we obtain the eigenvalues E N that give the masses M N of the states as In summary, the spectrum of exotic hadrons that are sufficiently tightly bound that our hierarchy of scales, and in particular the multipole expansion, applies is similar to that one of diatomic molecules illustrated in Fig. 1. The quantum number κ identifies, through different shifts Λ κ , different excitations of the light degrees of freedom.
The gap between different excitations is (at least for the lower states) of order Λ QCD . In the case of the diatomic molecule the different electronic excitations are separated by a gap of order mα 2 . For each BO potential the vibrational modes of the heavy quark-antiquark pair generate a fine structure of levels, E N , separated for fixed κ by small gaps of order M w 2 . Similarly, in the molecular case the vibrational modes of the nuclei induce small splittings of order mα 2 m M . There are, however, also noteworthy differences. In the hadronic case, if the size of the hadron is much larger than the distance between the heavy quark and antiquark, then κ labels spherically symmetric states. Because the symmetry of the hadron is cylindrical, this means that at short distances some excitations of the light degrees of freedom turn out to be degenerate. As a consequence the equations of motion are the coupled Schrödinger equations of Eq. (71) that mix different excitations, labeled by λ, λ ′ , with the same κ. The mixing happens through the nonadiabatic coupling, which under our assumptions counts like the quark-antiquark kinetic energy. A physical consequence of the mixing is the so-called Λ-doubling, i.e., a lifting of degeneracy between states with the same parity [32]. In the molecular case, the size of the molecule and the typical distance between the nuclei is of the same order. Because there is no special hierarchy between these two lengths there is neither a special symmetry at short distance nor a corresponding degeneracy pattern. The equation of motion for the molecular case is the simple Schrödinger equation (48) [or (49) in the adiabatic approximation]. In this case, different electronic excitations do not mix at leading order. Moreover, the nonadiabatic coupling is subleading with respect to the relative kinetic energy of the nuclei.
The masses for heavy hybrid states have been obtained in Ref. [32] following the method just described. There, the light-quark part of h 0 was omitted. In Fig. 5 we reproduce the results of Ref. [32] compared with an updated list of possible experimental candidates. Tetraquarks were discussed in Ref. [3] in the context of the BO approximation (see also [39]). In [3], preliminary estimates for their masses were given assuming that the tetraquark static energies have the same shape as the hybrid ones and using values for Λ κ from Ref. [48]. One major difficulty is the lack of knowledge of the static energies carrying light-quark flavor quantum numbers. One expects that lattice QCD will soon provide results on these and other crucial nonperturbative matrix elements to be used in the BOEFT developed here.

VI. CONCLUSIONS AND PERSPECTIVES
The Born-Oppenheimer approximation is the usual tool for solving the Schrödinger equation of molecules. It relies on the movement of the nuclei being much slower than that of the electrons, a circumstance that allows to study the electronic eigenstates and energy levels for fixed positions of the nuclei, the so-called static energies. The wave functions of the molecule can then be expanded in terms of these electronic eigenfunctions resulting in a Schrödinger equation describing the molecular energy levels. We have used this hierarchy of scales to build an EFT that systematically describes the energy levels of the simplest diatomic molecule, H + 2 . Our starting point has been an EFT of QED for the ultrasoft scale, pNRQED, adapted to the case of two nuclei and one electron. Since pNRQED for two heavy and one light particle has not been presented in the literature before, we have worked out its derivation in some detail. Particular care has been put in including all the relevant operators suppressed in powers of m M , where m and M are the electron and nuclei masses respectively. Counting m M ∼ α 3 2 we have derived the pNRQED Lagrangian relevant to compute the spectrum up to O(mα 5 ).
The assumption that the nuclei move slower than the electrons, which is at the basis of the Born-Oppenheimer approximation, is equivalent to take the kinetic term of the nuclei to be of a smaller size than the energy scale of the electron dynamics, the ultrasoft scale. Being these two scales well separated, it is natural in an EFT framework to integrate out the ultrasoft degrees of freedom in order to obtain an EFT that describes the molecular degrees of freedom only. We have carried out this integration obtaining a molecular EFT that we have named Born-Oppenheimer EFT (BOEFT). Up to O mα 4 it is sufficient to match pNRQED and BOEFT at tree level, or equivalently, to expand the matter field in the pN-RQED Lagrangian in eigenfunctions of the leading-order Hamiltonian for the electron, as it is done in the Born-Oppenheimer approximation of the Schrödinger equation. Loop diagrams involving ultrasoft photons start contributing at O mα 5 , the first of such contributions being responsible for the H + 2 molecular Lamb shift. We have computed the leading ultrasoft loop and obtained the BOEFT Lagrangian relevant to compute the spectrum up to O mα 5 .
The precise size of the nuclei kinetic operator has been obtained using the virial theorem to relate it to the potential acting on the nuclei. At leading order this potential is formed by the repulsive Coulomb potential between the nuclei and the attractive electronic static energies. Since the system is bound, the nuclei do not move over the whole size of the molecule, but oscillate around the minimum of the potential. The size of the kinetic operator of the nuclei is of the order of mα 2 m M , which is smaller than the ultrasoft scale mα 2 . This is consistent with the original statement that the two nuclei dynamics occurs at a lower energy scale than the electronic one. The size of the nonadiabatic coupling could also be assessed resulting in the conclusion that for diatomic molecules its contribution to the energy levels is suppressed by a factor (m M ) 1 4 .
In the present paper we have derived the BOEFT Lagrangian for the H + 2 molecule up to operators relevant for the spectrum up to O mα 5 . This can be system-  atically improved by including higher-order operators in the power counting detailed in Sec. IV, and computing their corresponding matching coefficients. Similarly, all the relevant contributions up to a certain precision to a specific observable can be determined with the help of the power counting, which may be of crucial importance to handle high-precision calculations.
Having set the general framework for constructing the BOEFT in QED, we have analyzed systems in QCD analog to the diatomic molecule. These are systems made of a heavy quark-antiquark pair, which plays the role of the heavy degrees of freedom, bound with light-quarks or excited gluonic states, playing the role of the light degrees of freedom. In particular, we have studied the case in which the quark-antiquark pair appears in a color-octet state. In the short distance regime, r ≪ 1 Λ QCD , the multipole expansion is applicable and the system can be described using weakly-coupled pNRQCD.
The energy scale of the leading-order light degrees of freedom dynamics is Λ QCD , while, as in the molecular case, the heavy degrees of freedom dynamics, in this case that of the heavy quark-antiquark pair, takes place at the lower energy scale M w 2 . We have identified the leading-order Hamiltonian in the multipole and 1 M expansions for the light degrees of freedom, h 0 , and defined a basis of color-octet light degrees of freedom operators, which, together with the heavy quark-antiquark octet field, generate hadronic (color-singlet) eigenstates of the pNRQCD Hamiltonian. The Λ QCD scale has been integrated out and pNRQCD matched into a QCD version of the BOEFT. At LO in the multipole expansion the matching can be done by just projecting the octet sector of the pNRQCD Lagrangian on the basis of eigenstates of h 0 . At NLO the matching requires a full non-perturbative computation, nevertheless, some constraints on the form of the NLO term can be obtained from the multipole expansion itself and the cylindrical symmetry that the system possesses at finite separation between the heavy quarks. As in the diatomic molecular case, a nonadiabatic coupling between the heavy quarks and the light degrees of freedom arises from the matching procedure, however, unlike in the molecular case, this does not need to be suppressed with respect to the kinetic operator. Furthermore, the nonadiabatic coupling mixes states that in the short distance limit have degenerate potentials, therefore the mixing has to be taken into account when solving the set of Schrödinger equations that result from the Euler-Lagrange equations of the BOEFT. As a result the phenomenon known as Λ-doubling in molecular physics [34] is more prominent in the QCD case [32].
The BOEFT has been used to obtain the masses of the quarkonium hybrids in Ref. [32] (see also [49]). Preliminary studies on quarkonium tetraquarks using a similar framework based on the BO approximation were carried out in Ref. [3]. A further analysis is in preparation [39]. The EFT presented here could be straightforwardly extended to describe any system made of two heavy quarks bound adiabatically with some light degrees of freedom. An example are doubly heavy baryons, i.e., states with two heavy quarks and one light-quark. Experimentally, doubly heavy baryons have been first observed at the LHCb [50]. For a study of this system in the framework of pNRQCD, we refer to [51]. Another example are pentaquark states made of two heavy quarks and three lightquarks. Candidates have been observed at the LHCb [52], but a pNRQCD based study of these systems is still to be done.