pNRQCD determination of E1 radiative transitions

This contribution contains the first numerical computation of the complete set of relativistic corrections of relative order $v^{2}$ for electric dipole (E1) transitions in heavy quarkonium; in particular, for the processes $\chi_{bJ}(1P) \to \Upsilon(1S) + \gamma$ with $J=0,\,1,\,2$. We assume that the momentum transfer of the heavy mesons involved in the reactions lies in the weak-coupling regime of the low-energy effective field theory potential non-relativistic QCD (pNRQCD) and thus a full perturbative calculation can be performed.


Introduction
Electromagnetic transitions are often significant decay modes for bottomonium states below BB threshold (10.56 GeV), making them a suitable experimental tool to access the lowest spectra of bottomonia. For instance, the first bb states not directly produced in e + e − collisions were the six triplet-P states, χ b (2P J ) and χ b (1P J ) with J = 0, 1, 2, discovered in radiative decays of the Υ(3S ) and Υ(2S ) in 1982 [1,2] and 1983 [3,4], respectively.
Electric dipole (E1) transitions are defined through the property that they change the orbital angular momentum of the state by one unit, but not the spin. Therefore, the final state has different parity and C-parity than the initial one. Typical E1 quarkonium decays are the ones mentioned above: 2 3 P J → 1 3 S 1 + γ. Here and in the following we denote the states as n 2s+1 J , where n = n r + + 1 corresponds to the principal quantum number with n r = 0, 1, . . . the radial quantum number and the orbital angular momentum. The spin is denoted by s and J is the total angular momentum.
The E1 (and M1) electromagnetic transitions have been treated for a long time by means of potential models that use non-relativistic reductions of QCD-based quark-antiquark interactions (see, e.g., Ref. [6] for a recent application to the bottomonium system). However, the progress made in effective field theories (EFTs) for studying heavy quarkonia [7] and the new large set of accurate experimental data taken in the heavy quark sector by B-factories (BaBar, Belle and CLEO), τ-charm facilities (CLEO-c, BESIII) and even proton-proton colliders (CDF, D0, LHCb, ATLAS, CMS) ask for a systematic and model-independent analysis (see, e.g., Refs. [8,9] for reviews).
Formulae and numerical treatment of M1 transitions within the effective field theory potential NRQCD (pNRQCD) can be found in Refs. [10,11]. Therein, the relativistic corrections to the leading order (LO) expression (which counts as k 3 γ /m 2 where k γ is the photon energy) were computed in two different expansion schemes: (i) strict weak-coupling regime and (ii) including exactly the static potential in the LO Hamiltonian. Within the same theoretical framework, the corresponding formulae for E1 transitions have been presented in Ref. [12]. In this case, the relativistic corrections to the LO decay width (that counts as k 3 γ /(mv) 2 ) are much more involved covering not only higher order terms in the E1 transition operator but also corrections to the initial and final wave function due to higher order potentials and higher order Fock states. These facts have hindered numerical computations of the E1 radiative decays within pNRQCD (for partial calculations see [13]). This contribution aims to close this gap and calculate the decay rate of the reaction 2 3 P J → 1 3 S 1 + γ with J = 0, 1, 2. As a first step, we shall assume that the soft scale lies in the strict weak-coupling regime of pNRQCD and thus a full perturbative calculation can be performed. These proceedings are based on the forthcoming publication [14].

Theoretical set-up Potential non-relativistic QCD (pNRQCD)
Heavy quarkonium systems are characterized by their non-relativistic nature, i.e., the heavy quark bound-state velocity, v, satisfies v 1. This is reasonably fulfilled in bottomonium (v 2 ∼ 0.1) and to a certain extent in charmonium (v 2 ∼ 0.3). Moreover, at least, three widely separated scales appear: the heavy quark mass m (hard scale), the relative momentum of the bound state p ∼ mv (soft scale) and the binding energy E ∼ mv 2 (ultrasoft scale). With v 1, the following hierarchy of scales is satisfied and this allows for a description in terms of EFTs for physical processes taking place at one of the lower scales. The integration out of modes associated with high-energy scales is performed as part of a matching procedure that enforces the equivalence between QCD and the EFT at a given order of the expansion in v. The final result is a factorization at the Lagrangian level between the high-energy modes, which are encoded in the matching coefficients, and the low-energy contributions carried by the dynamical degrees of freedom. The suitable EFT to describe processes that take place at the scale mv such as the E1 radiative transitions between the lowest heavy quarkonium states is potential NRQCD (pNRQCD) [15,16]. It follows by integrating out the modes of order p ∼ 1/r ∼ mv from NRQCD [17,18] which in turn comes from QCD by integrating out the high energy modes of order m. Therefore, pNRQCD takes full advantage of the hierarchy of scales that appear in Eq. (1), and makes a systematic and natural connection between quantum field theory and the Schrödinger equation. Schematically, the pNRQCD equation of motion takes the form  φ( r, t, R ) = 0 + corrections to the potential + interactions with other low-energy degrees of freedom s (r) is the static potential and φ( r ) is the QQ field. Note here that the interactions with other low-energy degrees of freedom produce, among others, non-potential terms that account for singlet Figure 1. Kinematics of the radiative transition H → H γ in the rest frame of the initial-state quarkonium H, taken from [10].
to octet transitions via ultrasoft gluons and provide loop corrections to the leading potential picture. Being induced by low-energy degrees of freedom they encode also non-perturbative effects.
The matching of pNRQCD depends on the relative size between the soft and the Λ QCD scale. Two main situations can be distinguished, namely, the weak-coupling [15,16] (mv Λ QCD ) and the strong-coupling [19] (mv ∼ Λ QCD ) versions of pNRQCD. One major difference between them is that in the former the potential can be computed in perturbation theory unlike in the latter.
It is obvious that the weak-coupling version of pNRQCD is amenable for a theoretically much cleaner analysis. The observables can be computed as an expansion in α s with increasing accuracy. Non-perturbative effects are suppressed by powers of Λ QCD /(mv). Therefore, observables that could be computed with the weak-coupling version of pNRQCD are of the greatest interest.
Decay width of the n 3 P J → n 3 S 1 γ reaction The complete decay rate n 3 P J → n 3 S 1 γ reads up to order k 3 γ /m 2 [12] Γ n 3 P J →n 3 S 1 γ = Γ (0) where R S =1 (J) includes the initial and final state corrections due to higher order potentials and higher order Fock states (see below). The remaining corrections within the brackets are the result of taking into account additional electromagnetic interaction terms in the Lagrangian suppressed by O(v 2 ) [12].
We have displayed terms proportional to the anomalous magnetic moment, κ em Q , however these terms are at least suppressed by α s (m)v 2 and thus go beyond our accuracy and are therefore not considered in the numerical analysis. The LO decay width (∼ k 3 γ /(mv) 2 ) is with α em the electromagnetic fine structure constant, e Q the charge of the heavy quarks in units of the electron charge, and k γ the photon energy determined by the kinematics shown in Fig. 1: The function is a matrix element that involves the radial wave functions of the initial and final states. We shall assume that these states are solutions of the Schrödinger equation with the leading order Hamiltonian in weakly-coupled pNRQCD given by where C F = 4/3. Therefore, ψ (0) n m ( r ) and E (0) n can be written in the hydrogen-like form where m r = m/2 is the reduced mass of the QQ system, ρ n = 2r/na is a dimensionless variable with a = 1/m r C F α s the Bohr radius. Finally, the normalization reads

Relativistic wave-function corrections
Due to higher order potentials and transitions between singlets and octets, the state in Eq. (8) is not an eigenstate of the complete Hamiltonian. Therefore, one has to consider corrections to the wave function, which can contribute to the decay rate at the required order of precision (∼ k 3 γ /m 2 ). To compute these corrections one applies the standard formalism of perturbation theory, either in the language of quantum mechanics or via Feynman diagrams.

Corrections due to higher order potentials
In order to account for corrections to the decay width of relative order v 2 , we need to consider the complete Hamiltonian The static potential is given by where, as mentioned above, V (0) s (r) = −C F α s /r, is the leading order potential or Coulomb-like potential that goes directly in the Schrödinger equation. The O(α s ) and O(α 2 s ) radiative corrections to the LO static potential are (the constants shown herein can be found, e.g., in Appendix C of Ref. [20]): a 2 (ν, r) = a 2 + π 2 3 β 2 0 + (4a 1 β 0 + 2β 1 )ln(νe γ E r) + 4β 2 0 ln 2 (νe γ E r) .
The O(α s ) term was computed in Ref. [21] and the O(α 2 s ) in Ref. [22]. The static potential is known up to order O(α 4 s ) with the O(α 3 s ) radiative correction computed in Refs. [23][24][25][26]. However, already O(α 3 s ) terms would give a contribution to the E1 decay rate that goes beyond present precision. The term δH encodes the relativistic corrections which can be organized as an expansion in the inverse of the heavy quark mass, m. At the order we are interested in, such expansion covers all the 1/m and 1/m 2 potentials and, at order 1/m 3 , the first relativistic correction to the kinetic energy: At order 1/m 2 , we can split the contributions into spin-independent (SI) and spin-dependent (SD) terms [7] where S = S 1 + S 2 = ( σ 1 + σ 2 )/2, L = r× p and S 12 = 3(r· σ 1 )(r· σ 2 )− σ 1 · σ 2 are, respectively, the total spin, total orbital angular momentum and tensor operators acting on the system. In the weak-coupling case, the above potentials read at leading (non-vanishing) order in perturbation theory We now make use of standard quantum mechanical perturbation theory and compute the first and second order correction, induced by a potential V, to a state |n (0) ≡ |n . The second order correction to the wave function is only needed when the perturbation is given by the static potential proportional to the a 1 (ν, r) term. The normalised corrected wave-function is for the first order, and for the second one.
As one can see in Eq. (20), a particular re-arrangement of the terms allows us to have a key expression that can be re-written as n n , |n n | for the first order, and for the second order. The term δE (1) V in Eq. (24) is the first order correction to the energy induced by a potential V: δE (1) V ≡ d 3 r ψ * n m ( r ) V( r ) ψ n m ( r ); and G n ( r 1 , r 2 ) has the following expression where G( r 1 , r 2 , E) is the Coulomb Green function

Corrections due to higher order Fock states
The weakly coupled quarkonia may also get corrections from the coupling of the heavy quarkantiquark pair to other low-energy degrees of freedom. In particular, the leading order electromagnetic dipole transition may get a correction from diagrams (see Fig. 8 in [12]) in which a singlet state is 1 In order to perform the computation it is specially useful to use such expression, because for λ = n coupled to an octet state due to the emission and re-absorption of an ultrasoft gluon. These diagrams come from terms of the pNRQCD Lagrangian like [7] where S = S 1 c / √ N c is a quark-antiquark field that transforms as a singlet under S U(3) c and U(1) em , O = √ 2O a T a is a quark-antiquark field which transforms as an octet under S U(3) c and as a singlet under U(1) em , and E is the chromo-electric field.
The first two diagrams in Fig. 8 of [12] correspond to the renormalization of the initial and final wave function. The diagrams 2, 3a and 3b account for the correction of the initial and final wave functions due to the presence of octet states. The diagram 4 represents an electric dipole transition mediated by the intermediate octet state. According to the power counting, the first two diagrams contribute to relative order Λ 2 QCD /(mv) 2 whereas the remaining ones scales as Λ 3 QCD /(mv 2 )/(mv) 2 . We shall not consider these contributions herein because in the strict weak-coupling regime, E ∼ mv 2 Λ QCD , one can argue that they should be negligible.
It is noteworthy that, in contrast to the E1 transitions, the colour-octet contributions for allowed M1 transitions cancel [10]. This is a consequence of the fact that the magnetic dipole operator behaves as an identity operator in position space.

Results
We discuss in detail our theoretical result for the χ b1 (1P) → Υ(1S )γ reaction, but a similar analysis has been performed for the transitions χ bJ (1P) → Υ(1S )γ with J = 0, 2. The mean value for the decay width and an estimate of its theoretical error will be given at the end of this Section for all transitions. Figure 2 shows the χ b1 (1P) → Υ(1S )γ LO decay rate and its relativistic correction due to higher order electromagnetic interactions that appear in the pNRQCD Lagrangian. In other words, we are analysing Eq. (2) without the factor R S =1 (J = 1). As one can see in the left panel of 1 to 3 GeV. This range encompasses the typical momentum transfer in the bottomonium system, still being consistent with perturbation theory. Let us focus now our attention to the computation of the wave function corrections due to higher order potentials, which are encoded in the factor R S =1 (J = 1) of Eq. (2). The upper panels of Fig. 3 show the matrix elements correcting the E1 decay rate up to O(v 2 ) and coming from the radiative corrections to the static potential, Eq. (12). The left and middle panels refer to the first order initial and final wave function corrections coming from a 1 (ν, r) and a 2 (ν, r), respectively. The right panel refers to the second order correction due to the a 1 (ν, r) term of the static potential. Amongst the features shown by the panels, the following are of particular interest: (i) the matrix elements clearly exceed the value of the LO one. (ii) The matrix elements depend quite dramatically on the scale ν, especially for small ν; in some sense, we expected such behaviour from the numerical analysis of the M1 transitions in Refs. [10,11]. (iii) The zero crossing in some of the matrix elements comes from the logarithms in (13) and (14).
The lower panels of Fig. 3 show the remaining matrix element contributions coming from δH, Eq. (15). One can see that only few of them are relevant corrections to the LO decay rate. Moreover, the ν-dependence of every matrix element is smaller than in the case of the radiative corrections.  Summing all the contributions discussed in the paragraph above, the left panel of Fig. 4 shows the next-to-leading order (NLO), NNLO and NLO+NNLO matrix elements and compares them with the LO term. The most important features have been already mentioned: the subleading matrix elements are of the same order of magnitude than the leading one and the dependence with ν in the logs dominates the picture. In the right panel of Fig. 4, we draw the decay rate associated with the χ b1 (1P) → Υ(1S )γ reaction at LO, NLO and NNLO. It is worth to remark that the NLO contribution is negligible at large-ν but multiplies by a factor of 2 the LO decay width at ν = 1 GeV. A big correction to the decay rate is due to the NNLO contribution. One can see in the Figure that the theoretical result depends slightly on the scale for ν 1.75 GeV, whereas the ν-dependence is dramatic for lower values due to the logarithmic functions. This fact is demonstrated by the additional curve (dotted green) where we omitted the contributions coming from the radiative corrections to the static potential, hence set the a 1 (ν, r) and a 2 (ν, r) terms to zero. Note that the relativistic corrections to the leading order E1 transition operator are included in the NNLO curve.
Finally, our theoretical results for the decay rates of the transitions under consideration are obtained by choosing the value at ν = 1.5 GeV, yielding: where we have chosen a very conservative error estimation that includes the total range of our final result, obtained by varying ν from (1-3) GeV.

Epilogue
We have presented the first numerical determination of the decay rate χ bJ (1P) → Υ(1S )γ with J = 0, 1, 2 within potential NRQCD. We have assumed that the momentum scale of the heavy quarkonium involved lies in the strict weak-coupling regime of pNRQCD and non-perturbative effects are suppressed, such that a full perturbative calculation can be performed. Relativistic corrections of relative order v 2 to the LO decay rate are included. The analysis separates those contributions that account for the higher order electromagnetic interaction terms in the pNRQCD Lagrangian and those that account for quarkonium state corrections due to higher order potentials and transitions between singlets and octets.