Electric dipole transitions of 1P bottomonia

We compute the electric dipole transitions $\chi_{bJ}(1P)\to \gamma\Upsilon(1S)$, with $J=0,1,2$, and $h_{b}(1P)\to \gamma\eta_{b}(1S)$ in a model-independent way. We use potential non-relativistic QCD (pNRQCD) at weak coupling with either the Coulomb potential or the complete static potential incorporated in the leading order Hamiltonian. In the last case, the perturbative series shows very mild scale dependence and a good convergence pattern, allowing predictions for all the transition widths. Assuming $\Lambda_{\text{QCD}} \ll mv^2$, the precision that we reach is $k_{\gamma}^{3}/(mv)^{2} \times \mathcal{O}(v^{2})$, where $k_{\gamma}$ is the photon energy, $m$ is the mass of the heavy quark and $v$ its relative velocity. Our results are: $\Gamma(\chi_{b0}(1P)\to \gamma\Upsilon(1S)) = 28^{+2}_{-2}~\text{keV}$, $\Gamma(\chi_{b1}(1P)\to \gamma\Upsilon(1S)) = 37^{+2}_{-2}~\text{keV}$, $\Gamma(\chi_{b2}(1P)\to \gamma\Upsilon(1S)) = 45^{+3}_{-3}~\text{keV}$ and $\Gamma(h_b(1P)\to \gamma\eta_b(1S)) = 63^{+6}_{-6}~\text{keV}$.


I. INTRODUCTION
Electromagnetic transitions are often a significant decay mode for bottomonium states below the BB threshold (10.56 GeV), making them a suitable experimental tool to access lower states. For instance, the first bb states not directly produced in e þ e − collisions were the six triplet-P states, χ bJ ð2PÞ and χ bJ ð1PÞ, with J ¼ 0, 1, 2, discovered in radiative decays of the ϒð3SÞ and ϒð2SÞ in 1982 [1,2] and 1983 [3,4], respectively.
Electromagnetic transitions can be classified in terms of electric and magnetic multipoles. The most important ones are the E1 (electric dipole) and the M1 (magnetic dipole) transitions; higher order multipole modes E2, M2, E3, … appear in the spectrum, but are suppressed. The width of allowed (hindered) M1 transitions is of order k 3 γ =m 2 (k 3 γ v 2 =m 2 ) where k γ is the photon energy and m is the mass of the heavy quark, whereas the width of E1 transitions is of order k 3 γ =ðmvÞ 2 , where v, which is much smaller than 1, is the relative velocity of the heavy quarks in the quarkonium [5]. Electric dipole transitions happen much more frequently than magnetic dipole transitions. The branching fraction for E1 transitions is indeed significant for the bottomonium states that we shall study in this work [6]: Bðχ b0 ð1PÞ → γϒð1SÞÞ ¼ ð1.94 AE 0.27Þ%, Bðχ b1 ð1PÞ → γϒð1SÞÞ ¼ ð35.0 AE 2.1Þ%, Bðχ b2 ð1PÞ → γϒð1SÞÞ ¼ ð18.8 AE 1.1Þ%, and Bðh b ð1PÞ → γη b ð1SÞÞ ¼ ð52 þ6 −5 Þ%. Even in the χ b0 case this is the largest observed exclusive branching fraction.
Electric dipole transitions are characterized by the fact 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 examples of E1 quarkonium decays are the ones mentioned above: 2 3 P J → 1 3 S 1 þ γ and 2 1 P 1 → 1 1 S 0 þ γ. Here and in the following we denote the states as n 2sþ1 l J , where n ¼ n r þ l þ 1 is the principal quantum number, with n r ¼ 0; 1; … the radial quantum number and l the orbital angular momentum usually represented by a letter: S for l ¼ 0, P for l ¼ 1 and so on. The spin is denoted by s and J is the total angular momentum. We use also the PDG notation, where χ bJ ð1PÞ identifies the state 2 3 P J , and h b ð1PÞ the state 2 1 P 1 . This is to say, in the PDG notation, 1P bottomonia are states with quantum numbers n ¼ 2 and l ¼ 1.
E1 (and M1) electromagnetic transitions between heavy quarkonia have been treated for a long time by means of potential models that use nonrelativistic reductions of QCD-based quark-(anti)quark interactions (see, for instance, Ref. [7] for a recent application to the bottomonium system). However, the release in the last decade of a new large set of accurate experimental data, concerning electromagnetic reactions in the heavy quark sector, by B-factories (BABAR, Belle and CLEO), τ-charm facilities (CLEO-c, BESIII) and even proton-(anti)proton colliders (CDF, D0, LHCb, ATLAS, CMS) [8,9] demands for systematic and model-independent treatments. The aim of this paper is to compute the E1 transitions χ bJ ð1PÞ → γϒð1SÞ, with J ¼ 0, 1, 2, and h b ð1PÞ → γη b ð1SÞ using potential nonrelativistic QCD (pNRQCD). Quarkonium is characterized by the hierarchy of energy scales: where p is the relative momentum of the heavy quarks, proportional to the inverse of the size of the quarkonium, and E is the binding energy. The relative heavy quark velocity, v, is assumed to be v ≪ 1, which qualifies quarkonium as a nonrelativistic bound state. pNRQCD is a nonrelativistic effective field theory that takes advantage of this hierarchy of scales by systematically computing quarkonium observables as expansions in v [10,11] (see Refs. [12,13] for reviews). In the case of radiative transitions another relevant scale is the photon energy, k γ . The photon energy is about the energy gap between the initial and final quarkonium states: for allowed (hindered) M1 transitions it is of the order of mv 4 (mv 2 ), for E1 transitions it is of the order of mv 2 . The theory for M1 transitions in pNRQCD has been developed in [14] and extended to E1 transitions in [15]. Reference [15] provides the theoretical basis for the present study, which aims at computing E1 transitions from 1P bottomonium states at relative order v 2 , i.e., at order k 3 γ =m 2 in the transition width. The specific details of the construction of pNRQCD depend on the relative size of the scale mv 2 with respect to Λ QCD . In this paper, we assume that mv 2 ≫ Λ QCD . 1 The propagation of a color singlet heavy quark-antiquark field, S, is described at relative order v 2 by the Lagrangian density: where r is the quark-antiquark distance parametrizing the color singlet field S and V is the quark-antiquark potential. The operator −i ⃗ ∇ ∼ mv 2 is the center of mass momentum (the derivative acts on the center of mass coordinate), while −i ⃗ ∇ r ∼ mv is the relative momentum (the derivative acts on the distance r). If mv 2 ≫ Λ QCD , the potential V may be computed order by order in perturbation theory and v ∼ α s , where α s is the strong coupling evaluated at the typical momentum transfer scale. At leading order in α s , V is given by the Coulomb potential between static color triplet and color antitriplet sources: V ð0Þ s ¼ −4α s =ð3rÞ. According to the pNRQCD counting V ð0Þ s ∼ mv 2 . E1 transitions are encoded in the part of the pNRQCD Lagrangian, L γpNRQCD , that describes the interaction of the quarkantiquark field S with the electromagnetic field: The displayed term is the leading order electric dipole interaction term (ee Q stands for the electric charge of the heavy quark Q and ⃗ E em for the electric field), whereas the dots stand for higher order operators contributing to the E1 transition at relative order v 2 (or smaller), whose explicit expressions can be read off from Ref. [15].
There seems to be a growing consensus in the literature that the weak-coupling regime mv 2 ≫ Λ QCD may indeed be applied to many physical observables in the bottomonium sector including n ¼ 2 bottomonium states (for early work see [16][17][18], for reviews see [8,9,13], for recent work see [19,20]). In order to reach this conclusion, it is crucial, however, to have a proper treatment for the large terms appearing in the perturbative expansion. As long as α s remains a perturbative coupling, large terms can be due to factorially growing coefficients, which may require renormalon subtraction, or large logarithms in the renormalization scale.
In this work, we adopt methods to deal with both large corrections, eventually achieving a convergent expansion with mild dependence on the renormalization scale. Concerning the renormalon subtraction scheme, we adopt the one of Ref. [21]. Concerning the resummation of large logarithms, we rearrange the perturbative expansion of pNRQCD in such a way that the static potential is exactly included in the leading order (LO) Hamiltonian. This expansion scheme has been applied to the computation of the heavy quarkonium electromagnetic decay ratios in Ref. [22] and to the determination of M1 transitions between low-lying heavy quarkonium states in Ref. [23]. The authors obtain agreement between theory and experiment for the case of the charmonium and bottomonium ground states and for the n ¼ 2 excitations of the bottomonium. Very recently, the same scheme has been applied to the spectrum of n ¼ 2, l ¼ 1 quarkonium states in [20]. Hence, another motivation for the present study is to probe weakly coupled pNRQCD in the context of electric dipole transitions from the spin-triplet and spin-singlet lowest bottomonium P-wave states.
In Ref. [15], the complete set of relativistic corrections of relative order v 2 with respect to the leading order E1 decay width has been derived. In the E1 case, differently from M1 transitions [14,23], the computation of relativistic corrections at relative order v 2 is technically complicated: In addition to the effects due to higher order operators contributing to the E1 transition [the dots in Eq. (3)], one needs to calculate order v and v 2 corrections to the initial and final state wave functions due to higher order 1 The following computations are valid also for mv 2 ∼ Λ QCD . What changes in this case is, however, the parametrical size of the nonperturbative corrections, see Sec. II B 2 and comments in the conclusion. potentials. 2 This complication has hindered so far complete numerical computations of the E1 transitions between lowlying heavy quarkonium states within pNRQCD (for partial calculations see Refs. [24,25]). The present paper aims to close this gap.
The paper is structured in the following way. In Sec. II we discuss the theoretical background of the computation and display the formulas that we use for the decays. In this section, we present also results for the electric dipole transitions when only the LO static potential is incorporated in the Schrödinger equation. Section III is devoted to present the same results but incorporating the complete static potential in the LO Hamiltonian. Renormalon effects and resummation of large logarithms are also taken into account in this part. All of this leads to a good convergence pattern for the studied decay rates and thus to firm predictions for all of them. We summarize our results and conclude in Sec. IV.

II. NUMERICAL ANALYSIS IN PNRQCD
AT WEAK COUPLING: FIXED ORDER CALCULATION

A. Decay width
We aim at computing electric dipole (E1) transitions from 1P bottomonium states at order k 3 γ =m 2 under the condition mv 2 ≫ Λ QCD . The formulas for the decay widths have been derived in Ref. [15]. They read where R S¼1 nn 0 ðJÞ and R S¼0 nn 0 include the initial and final state corrections due to higher order potentials (see Sec. II B 1) and possibly higher order Fock states (see Sec. II B 2). The remaining corrections within the curly brackets are the result of taking into account additional electromagnetic interaction terms in the Lagrangian suppressed by Oðv 2 Þ [the dots in Eq. (3)]. For completeness, we have displayed in the formulas terms proportional to the anomalous magnetic moment, κ em Q . These terms will, however, not be considered in the numerical analyses because they are at least of order α s k 3 γ =m 2 and thus beyond our accuracy. The LO decay width, which scales like k 3 γ =ðmvÞ 2 , is Γ ð0Þ with α em the electromagnetic fine structure constant, e Q the charge of the heavy quark Q in units of the electron charge, and k γ the photon energy determined by the kinematics shown in Fig. 1: The LO decay width follows from the LO electric dipole interaction in the pNRQCD Lagrangian shown in Eq. (3).
All other terms in Eqs. (4) and (5) are of relative order v 2 with respect to the LO decay width. In particular, the function is a matrix element that involves the radial wave functions of the initial and final states. From r ∼ 1=p ∼ 1=ðmvÞ it follows that it scales like ðmvÞ 2þk−N . Under the assumption mv 2 ≫ Λ QCD we can compute the quarkonium potential in perturbation theory, i.e., as an expansion in α s . The wave functions are then the solutions of the Schrödinger equation where H ð0Þ contains the (perturbative) quark-antiquark static potential. More specifically, in this section we take the leading order Hamiltonian as where − ⃗ ∇ 2 r =m is the (nonrelativistic) kinetic energy in the center of mass frame and This means that we include in the static potential only the LO potential in α s , which is the Coulomb potential times the Casimir of the fundamental representation in SU (3), i.e., 4=3. A different choice will be analyzed in Sec. III. With the choice (11), ψ nlm ð⃗ rÞ and E n can be taken from the hydrogen-atom and read where ρ n ¼ 2r=ðnaÞ is a dimensionless variable and a ¼ 3=ð2mα s Þ is the Bohr-like radius. The functions L 2lþ1 n−l−1 and Y lm are the associated Laguerre polynomials and the spherical harmonics, respectively. The normalization reads Finally, if not differently specified, here and in the rest of the paper, α s is understood evaluated at the renormalization scale ν∶ α s ≡ α s ðνÞ. Hence the potential, the Bohr-like radius and, through it, the wave functions depend on ν.

B. Relativistic wave function corrections
The LO wave function (12) gets corrections due to higher order potentials and possibly to higher order Fock states. Corrections due to higher order potentials contribute at relative order v 2 , and therefore have to be included in the analysis to reach a precision of order k 3 γ =m 2 . These corrections will be outlined in the next Sec. II B 1. Corrections due to higher order Fock states will be discussed in Sec. II B 2.

Corrections due to higher order potentials
To account for Oðv 2 Þ corrections to the decay width due to higher order potentials, we need to consider the Hamiltonian The quark-antiquark static potential up to next-to-next-toleading order (NNLO) is given by where the coefficients of the Oðα s Þ and Oðα 2 s Þ radiative corrections to the LO static potential are 3 . The coefficients β i are the coefficients of the β-function with β 0 ¼ 11 − 2n f =3 and β 1 ¼ 102 − 38n f =3; n f is the number of massless flavors. The Oðα s Þ correction was computed in Ref. [26] and the Oðα 2 s Þ one in Ref. [27]. In this section, we consider higher order corrections to the static potential as perturbations around the leading order solution of Sec. II A. Hence, the order α s correction contributes to the transition width at relative order v in first order quantum mechanical perturbation theory and at relative order v 2 in second order quantum mechanical perturbation theory, whereas the order α 2 s correction contributes at relative order v 2 in first order quantum mechanical perturbation theory. On the other hand, the Oðα 3 s Þ correction, which is also known from Refs. [28][29][30], would give a contribution to the E1 decay rate of relative order v 3 , which is beyond our precision. Therefore, we will not include Oðα 3 s Þ corrections in this part of our analysis.
The term δH contains relativistic corrections to the potential and to the kinetic energy. They can be organized as an expansion in the inverse of the heavy quark mass, m. At the order we are interested in, such an expansion includes 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 have distinguished between spinindependent (SI) and spin-dependent (SD) terms: where f; g stands for the anticommutator. The above potentials read at leading (nonvanishing) order in perturbation theory (see, e.g., Ref. [12]): All these potentials contribute through first order quantum mechanical perturbation theory at relative order v 2 to the E1 width.
Using quantum-mechanical perturbation theory, we compute the first and, for one term, the second order correction, induced by δV ¼ ðV s − V ð0Þ s Þ þ δH, to the wave function ψ nlm ð⃗ rÞ ≡ h⃗ rjnlmi of energy E n . The second order correction to the wave function is only needed when the perturbation is given by the next-to-leading order (NLO) term in the static potential, i.e., the one proportional to a 1 ðν; rÞ. The (normalized) corrections to the wave function are at first order and at second order The operator P n 0 ≠n;l 0 ;m 0 jn 0 l 0 m 0 ihn 0 l 0 m 0 j E n −E n 0 appearing in the Eqs. (23) and (24) can be rewritten as and it can thus be identified with the pole-subtracted Coulomb Green function. In coordinate space, it reads G 0 n ð⃗ r 1 ; ⃗ r 2 Þ ≡ h⃗ r 1 j 1 ðE n − HÞ 0 j⃗ r 2 i ¼ lim where Gð⃗ r 1 ; ⃗ r 2 Þ is the Coulomb Green function [31,32]: E ≡ −4mα 2 s =ð9λ 2 Þ and ρ λ;i ¼ 2r i =ðλaÞ. In calculations it may be useful to set λ ¼ n= , since in this way we have E ¼ E n ð1 − ϵÞ and E → E n for ϵ → 0. Therefore, the first order and second order corrections to the expectation values of an arbitrary operator O may be written as (note that, for the sake of simplicity, only initial state corrections are shown herein, but the same corrections affect also the final state): where δE ð1Þ V is the first order correction to the energy: δE ð1Þ V ≡ R d 3 r ψ Ã nlm ð⃗ rÞδVð⃗ rÞψ nlm ð⃗ rÞ. As a final remark we note that, although in (19) we have included the center of mass kinetic energy, − ⃗ ∇ 2 =4m, this term does not contribute at our accuracy. 3 The reason is that, even if the center of mass kinetic energy scales like a term of relative order v 2 , nevertheless, its contribution vanishes at first order in quantum mechanical perturbation theory, Eq. (23), as the states are eigenstates (in fact simple plane waves) of the center of mass momentum.

Corrections due to higher order Fock states
The LO correction to E1 transitions due to higher order Fock states comes from diagrams in which a heavy quark-antiquark color singlet state is coupled to a heavy quark-antiquark color octet state via emission and reabsorption of gluons whose energy and momentum are of order mv 2 or Λ QCD . The coupling of the color singlet field S with the color octet field O is encoded in the pNRQCD Lagrangian in a chromoelectric dipole interaction term: The relevant Feynman diagrams in pNRQCD are shown in Fig. 8 of Ref. [15]: they are diagrams corresponding to the normalization of the initial and final state wave functions, diagrams accounting for the corrections to the initial and final state wave functions due to the presence of octet states, and a diagram representing an electric dipole transition mediated by an intermediate octet state. According to the power counting of pNRQCD, those diagrams contribute to relative order α s v 2 if the gluons carry an energy and a momentum of order mv 2 . They contribute to relative order Λ 2 QCD =ðmvÞ 2 or Λ 3 QCD =ðm 3 v 4 Þ if the gluons are nonperturbative and carry an energy and a momentum of order Λ QCD . In the first case, their contribution is smaller than v 2 by a factor α s and hence beyond our accuracy. In the second case, it is also smaller than v 2 if mv 2 ≫ Λ QCD , which is what we have assumed. It should be remarked, however, that it suffices mv 2 ∼ Λ QCD for the nonperturbative contributions to be of the same relative order, v 2 , as the ones coming from higher order potentials.

C. Numerical analysis
We specify, first, the parameters that enter in the determination of the bottomonium E1 transition widths. We have 4 where n f is the number of massless flavors, 5 e b is the electric charge of the bottom quark in units of the electron charge e and α em is the electromagnetic fine structure constant.
The masses of the initial and final quarkonium states are chosen to be the ones reported by the PDG [6], listed in Table I. The photon energies are determined by the kinematics of the two body decay, Eq. (7), and are given by Our reference value for the strong coupling constant is α ðn f ¼3Þ s ð1 GeVÞ ¼ 0.480. We obtain this value by using the RUNDEC package [34] to run down from α ðn f ¼5Þ s ðM Z ¼ 91. 19 GeVÞ ¼ 0.118 at four-loop accuracy. We then run α s to the typical scales of the bound state.
We fix the bottom quark pole mass using the experimental mass of the ϒð1SÞ state and the leading order binding energy. This means that if the bottom mass is which is the expression that goes into the wave function. Higher order terms are beyond our accuracy. Indeed, even the Oðα 2 s Þ term given above is beyond our accuracy if used for higher order corrections in the 1=m expansion. For those corrections we set the bottom quark mass to be Notation η b ð1SÞ ϒð1SÞ h b ð1PÞ χ b0 ð1PÞ χ b1 ð1PÞ χ b2 ð1PÞ n 2sþ1 l J 1 1 S 0 1 3 S 1 2 1 P 1 2 3 P 0 2 3 P 1 2 3 P 2 Mass [6] 9.399 9.460 9.899 9.859 9.893 9.912 3 Potentials depending on the center of mass momentum contribute, instead, at relative order v 2 to M1 transitions [14]. 4 The value of α em ð400 MeVÞ, where 400 MeV is a reference value for the photon energy in a typical n ¼ 2, l ¼ 1 bottomonium E1 transition, has been computed using the ALPHAQED package [33]. 5 At the typical momentum transfer inside the bb system the charm quark decouples [17].
1. χ bJ ð1PÞ → γϒð1SÞ with J = 0, 1, 2 We begin the numerical analysis of the electric dipole transitions χ bJ ð1PÞ → γϒð1SÞ, with J ¼ 0, 1, 2, focusing on the contributions that appear in Eq. (4) and come from higher order electromagnetic operators in the pNRQCD Lagrangian. As one can see in Fig. 2, the leading order decay width depends strongly on the renormalization scale ν. This is due to the scale dependence of the Bohr-like radius that enters the wave functions. The effects from higher order electromagnetic operators are small. The correction to the LO decay width is at most ≈1%, ≈2% and ≈5% when the initial state is a χ b0 , χ b1 and χ b2 , respectively. This can be understood analyzing each contribution separately: The contributions almost cancel for J ¼ 0 but this is not the case for J ¼ 1 and J ¼ 2.
The radiative corrections to the LO static potential [the terms in the sum of Eq. (16)] lead to first order and second order quantum-mechanical corrections to the decay widths. These terms are proportional to (soft) logs [like ln n ðνrÞ] and thus one expects a significant scale dependence of the resulting matrix elements. This is indeed the case as shown in Fig. 3. 6 The plotted matrix elements, M, stand for the first order and second order corrections to the matrix elements of the specified potentials, according to Eqs. (28) and (29). 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. Among the features shown by the panels, the following are of particular interest: (i) The matrix elements clearly exceed the value of the LO one. To some extent, this is due to the factors stemming from the β-function in Eqs. (17) and (18) that are large. (ii) The matrix elements depend strongly on the scale ν, especially for small ν. A similar behavior shows up in some matrix elements contributing to the M1 transitions [23]. (iii) The zero crossing in some of the matrix elements comes from the logarithms in the Eqs. (17) and (18). The scale where this effect occurs is ν ≈ 1.2 GeV. (iv) Initial and final state corrections partially cancel each other, order by order.
The corrections to the matrix element of the LO electric dipole operator (3), due to the relativistic corrections to the bottomonium wave functions discussed in Sec. II B 1, are shown in Fig. 4. These corrections contribute to the term R S¼1 21 ðJÞ in Eq. (4). As one can see, most of the contributions are small, except for the final state correction induced by V ð1Þ and the correction due to V S 2 . The overall dependence on the scale ν is weak in all cases but a slight trend toward larger values by decreasing scale can be observed.
We sum the matrix elements that include radiative corrections to the static potential (see Fig. 3) and higher order relativistic corrections to the potential and kinetic energy (see Fig. 4) at each order and the result is shown in Fig. 5, first row. The corresponding decay widths are displayed in the second row. From both plot sequences we can see that each LO, NLO and NNLO contribution depends strongly on the renormalization scale and also that subleading contributions may be of similar size to the leading one. Moreover, the overall impact of the corrections decreases with increasing total angular momentum: For J ¼ 0 the NLO þ NNLO curves exceed for some ν the LO curve, for J ¼ 1 they touch each other and for J ¼ 2 they stay slightly below. The kink, visible in the NNLO matrix element (and subsequently also in the NLO þ NNLO matrix element) at about 1.2 GeV, can be traced back either to the zero crossing or to the maximum in the matrix elements of Fig. 3. Also the NLO and NNLO matrix element sums show a zero crossing, and the combined NLO þ NNLO matrix element has a clear maximum. The zero crossings yield vanishing contribution to the respective decay widths, as visible in the second row.
The results that follow from summing up all previous corrections, i.e., those that contribute to the term R S¼1 21 ðJÞ in (4) (we recall that these are radiative corrections to the static LS potential (dotted green for initial state correction), due to the V ð2Þ S 2 potential (dot-dashed red for final state correction), and due to the V ð2Þ S 12 potential (dashed violet and dotted brown for initial and final state corrections, respectively). Third row: Matrix element at LO (solid blue), and at NNLO: Including relativistic corrections due to the V ð2Þ p 2 potential (dashed orange and green for initial and final state corrections, respectively), and due to the kinetic energy term − ⃗ ∇ 4 =4m 3 (dotted red and violet for initial and final state corrections, respectively).
potential, due to the one and two loop corrections in (16), and higher order relativistic corrections to the potential and the kinetic energy, due to (19); their combined effect to the E1 transition width is shown in Fig. 5) and those that contribute to the terms other than R S¼1 21 ðJÞ in (4) (these are due to higher order electromagnetic operators in the pNRQCD Lagrangian (3); their effect to the E1 transition width is shown in Fig. 2), are shown in Fig. 6. The renormalization scale dependence of the decay widths is reduced as the NLO and NNLO corrections are included. For instance, varying the renormalization scale from 1 GeV to 3 GeV for the J ¼ 1 case, the LO spans over the range ð17-74Þ keV, incorporating the NLO contribution shrinks the range to ð35-75Þ keV, and adding the NNLO corrections results in a range of ð32-79Þ keV. Although a slight shift toward higher upper bounds is noticeable, the whole range and thus the overall scale dependence somewhat decreases from the LO.
Another feature of the panels in Fig. 6 is that by setting the terms proportional to a 1 ðν; rÞ and a 2 ðν; rÞ to zero, the decay width exhibits a different ν-dependence in the low ν region (dotted green curve). This suggests that the terms proportional to the logs in Eqs. (17) and (18) give rise to non-negligible contributions, whose dependence on the renormalization scale needs to be treated carefully, as we shall see in the next section. [keV] FIG. 6. Decay widths of the electric dipole transitions χ bJ ð1PÞ → γϒð1SÞ. The three panels refer to the three cases J ¼ 0, 1, 2, respectively. In each panel, the dashed blue curve is the LO decay width, the dot-dashed orange one incorporates LO þ NLO corrections and the solid black curve incorporates LO þ NLO þ NNLO corrections coming from higher order electromagnetic operators, radiative corrections to the static potential and higher order relativistic corrections to the potential and the kinetic energy. The dotted green curve is similar to the black one but it omits all corrections to the decay width due to radiative corrections of the static potential (one and two loops). We take our central value at ν ¼ 1.25 GeV, whereas the gray band indicates the associated uncertainty. The scale setting and the uncertainty estimate are explained in the text.
The convergence of the perturbative series is poor. This can be seen by looking at the difference between the LO and NLO, and between the NLO and NNLO results. Also the strong scale dependence in the range 1 GeV ≤ ν ≤ 3 GeV is a consequence of having large higher order corrections. As a consequence, it is difficult (if not impossible) to get a reliable result using fixed order perturbation theory. Nevertheless, in the following we will produce a first rough determination of the 1P bottomonium dipole electric transitions, with a large error reflecting the large uncertainty. We will overcome this difficulty and provide a reliable determination with a small uncertainty in the next section.
We choose to set the central value of the decay widths at the renormalization scale that self-consistently solves the Bohr-like radius equation: This scale is ν ¼ 1=a ¼ 1.25 GeV. 7 We estimate the uncertainty associated to the central value in a twofold way: (i) First, we vary the renormalization scale from 1 GeV to 3 GeV, which is a conservative interval including the lowest scale where perturbation theory may be still applicable and more than twice the inverse of the Bohr radius. (ii) Second, we estimate the uncertainty associated with truncating the perturbative series at NNLO and the fact that the series is poorly converging by taking one half of the maximum difference between the LO and the NNLO decay width. For the final error we choose the largest of these two values, which is indicated in the plots by a gray band. Further sources of uncertainties are given by the input parameters, these being the masses of initial and final states, and the value of α s . If we assume that these quantities are accurate within ≲ð1-3Þ%, their uncertainty is largely inside the final error.
Hence, a fixed order determination at NNLO gives for the E1 transition widths of the χ bJ ð1PÞ: As anticipated, the errors are large, reflecting the poor convergence of the perturbative series. In Sec. III, we will see how resumming the known terms of the perturbative expansion of the static potential into the wave functions will enormously improve the above determinations providing convergent expansions with tiny theoretical uncertainties.
We apply now the former analysis to the electric dipole transition h b ð1PÞ → γη b ð1SÞ. Figure 7 shows the LO decay rate and its correction due to higher order electromagnetic operators in the pNRQCD Lagrangian. In comparison with the χ bJ radiative transitions, the LO transition width and the correction induced by higher order operators are larger in this case. This is because the photon energy, k γ , is larger for increasing J in the χ bJ ð1PÞÞ → γϒð1SÞ transitions and even larger in the h b ð1PÞ → γη b ð1SÞ transition, see Eq. (31). The fact that the photon energy enters with the third power in the expression of the decay width explains then the overall increasing effect. The correction due to higher order operators is about 10%.
The corrections to the decay width due to the radiative corrections to the static potential (16) are the same as the ones shown in Fig. 3. As already mentioned, this is so because none of these potentials depends on either the spin s or the total angular momentum J. Figure 8 shows the corrections to the matrix element of the LO electric dipole operator (3), due to the relativistic corrections to the bottomonium wave functions discussed in Sec. II B 1. These corrections contribute to the term R S¼0 21 in Eq. (5). Since the initial and final states in the transition are now spin-singlet states, corrections to the wave functions due to the spin-orbit, spin-spin, and tensor potentials are absent. This has a major impact on the total NNLO matrix element because the correction induced by the V ð2Þ S 2 potential is zero now, whereas in the χ bJ case it is large (and negative) especially in the low ν region.
The left (middle) panel of Fig. 9 shows for each order the sum of all matrix elements (decay widths) including radiative corrections to the static potential (see Fig. 3) and higher order relativistic corrections to the potential and kinetic energy (see Fig. 8). The kink, visible in the NNLO and thus in the NLO þ NNLO matrix elements at about 1.2 GeV, can be traced back to the zero crossing or maximum in the matrix elements that account for the radiative corrections to the static potential. The absence of several negative contributions at NNLO yields a stronger dependence on the scale ν for values ν ≲ 1.5 GeV than in the χ bJ case. In this region of ν, the NLO þ NNLO matrix element and the subsequent decay width clearly exceed the leading order ones. The result that follows from summing up all previous corrections, i.e., those that contribute to the term R S¼0 21 in (5), shown in the first two panels of Fig. 9, and those that contribute to the terms other than R S¼0 21 in (5), shown in Fig. 7, is shown in the right panel of Fig. 9. The renormalization scale dependence of the decay width is reduced when the NLO and NNLO corrections are included: by varying the renormalization scale from 1 GeV to 3 GeV the LO decay width spans over the range ð27-114Þ keV, incorporating the NLO contribution shrinks the range to ð64-115Þ keV, and incorporating the NNLO correction shrinks further the range to ð97-127Þ keV. This comes at the cost of an even worse convergence pattern of the perturbative series than in the χ bJ case. We observe again a slight shift towards higher upper bounds, but the whole range and thus the overall scale dependence decreases.
Omitting the corrections to the decay width induced by the radiative corrections to the static potential results in a curve (green-dotted curve in the right panel of Fig. 9) that is quite close to the LO one at large values of ν and whose ν-scale dependence is weaker than the complete result at low values of ν. This is in contrast with the effect observed for the χ bJ states, but understandable since several additional contributions appear in the χ bJ case that are not present here.  We choose to set the central value of the decay width at ν ¼ 1.25 GeV, following the same prescription discussed for the χ bJ case. The main differences in comparison to the χ bJ transition width curves in Fig. 6 are the overall weaker scale dependence, the worse convergence of the perturbative series, and the shape of the curve for large values of the renormalization scale ν. Assigning the error to the transition width as in the χ bJ → γϒð1SÞ case discussed above, a fixed order determination at NNLO gives for the E1 transition width of the h b : In the following Sec. III, we will see how to improve also this determination by resumming the known terms of the perturbative expansion of the static potential into the wave function.

A. Log resummation and renormalon subtraction
We have seen in the previous section that the electric dipole transitions from the lowest-lying P-wave bottomonium states are not reliably described by fixed order calculations. The reason is that, even if these states are weakly coupled and the potential computable in perturbation theory, considering them at LO as Coulombic bound states is inadequate. Indeed, expanding around the Coulomb potential, V ð0Þ s , has led to a poor convergence of the perturbative series, resulting in a strong dependence on the renormalization scale and large theoretical uncertainties.
We deal with this problem by rearranging the perturbative expansion in pNRQCD in such a way that the static potential is exactly included in the LO Hamiltonian. One motivation for this reorganization of the perturbative series is the observation, originating from Ref. [35] (for more recent studies see, for instance, [36]), that, when comparing the static potential with lattice perturbation theory at short distances, the inclusion of higher order corrections is necessary to get a good agreement. An accurate treatment of the potential is particularly important for those observables, like the electric dipole transition widths, that are sensitive to the precise form of the wave function.
The new expansion scheme was applied in Ref. [22] to study electromagnetic decays of heavy quarkonium, and in Ref. [23] to compute magnetic dipole transitions between low-lying heavy quarkonia. The effect of the new rearrangement was found to be large. In particular, the exact treatment of the soft logarithms of the static potential made the renormalization scale dependence much weaker. We proceed herein to apply the same scheme to the electric dipole transitions under study. Like in the magnetic dipole transition computation performed in Ref. [23], an improvement in the convergence of the perturbative expansion is expected. The perturbative expansion will consist of just two terms: A leading order term, incorporating exactly the static potential, and a term incorporating the remaining corrections coming from higher order electromagnetic operators, and higher order relativistic corrections to the wave functions.
We follow the same setup of Ref. [23]. The leading order Hamiltonian reads now: where the static potential is ideally summed to all orders in perturbation theory. In practice, it is only known up to order α 4 s , hence we take 8 The analytical expressions of a 1 ðν; rÞ and a 2 ðν; rÞ have been given in Eqs. (17) and (18), respectively. The term a 3 ðν; ν us ; rÞ is known from Refs. [29,30]: where β 2 ¼ 2857=2 − 5033n f =18 þ 325n 2 f =54, a 3 may be read from the original literature or, for instance, from [32], and δa us 3 ðν; ν us Þ, encoding the subtraction of ultrasoft corrections from the static potential, is taken [28] δa us Ultrasoft corrections to the static potential are due to gluons carrying energy and momentum of order α s =r; the scale ν us is the factorization scale separating the ultrasoft energy and momentum region from higher ones. We will not resum here ultrasoft logs, like the one appearing in (43), although the result is known at leading [37] and next-to-leading accuracy [38]. The reason is that their numerical effect is small with respect to other sources of error. 8 To keep the notation simple, we will not explicitly write the dependence on the scale for quantities where this is due only to the truncation of the perturbative expansion.
The perturbative expansion (41) does not converge due to factorially growing terms that, once Borel resummed, give rise to singularities in the Borel plane, known as renormalons. The leading order renormalon affecting the static potential, V s , cancels against twice the pole mass m [39][40][41]. To make this cancellation explicit one adds/ subtracts the same renormalon contribution from twice the pole mass/the static potential ensuring that both are expressed in series of α s to the same power and at the same scale, e.g., ν: where ν Þα kþ1 s ðνÞ encodes the pole mass renormalon contribution, ν f is the renormalon factorization scale and X stands for the chosen renormalon subtraction scheme. For the renormalon subtraction scheme we use here the RS 0 scheme [21], 9 which amounts at choosing where Sðn; bÞ ¼ The mass that we are using for the bottom quark is m b;RS 0 ðν f ¼ 1.0 GeVÞ ¼ 4.859 GeV. It can be translated into the MS-mass: m b ðm b Þ ¼ 4.19 GeV [42]. Our reference value for N m is N m ¼ 0.574974 (for three light flavors) from Ref. [21]. 10 As in the previous section, our reference value for α s is α ðn f ¼3Þ s ð1 GeVÞ ¼ 0.480, and, like there, the running is implemented with four-loop accuracy. We set ν us ¼ ν f . This choice is motivated by the fact that ν us has to be smaller than the typical momentum transfer scale, i.e., ν us < p ∼ 1=a ¼ 1.25 GeV on the one hand, and ν us has to be larger than the scale where perturbation theory breaks down, say 0.7 GeV. Varying ν us from 0.7 GeV to 1.25 GeV induces a change from þ4% to −2% in the coefficient δa us 3 ðν; ν us Þ. The numerical impact of this change in the three loop coefficient of the static potential is negligible with respect to the dependence on the scale ν. This is not surprising as ultrasoft corrections are beyond the accuracy of the present study.
In the short range, it is possible to further improve the static potential by resumming potentially large logs of the type lnðνrÞ by setting the scale ν ¼ 1=r and yet achieve renormalon cancellation order by order in α s ð1=rÞ (see [35]). Following [23], we finally define our renormalon subtracted static potential in the RS 0 scheme as 9 We have checked against the RS [21] and the potential subtracted (PS) [41] schemes that the LO matrix element depends only mildly on the adopted renormalon subtraction scheme. 10 In the literature, there is an updated value, N m ¼ 0.563126, from Ref. [43], as well as other recent determinations, like N m ¼ 0.535 AE 0.010 from Ref. [44]. Since we have verified that these different determinations vary our results well inside the final errors, we will neglect in the following the uncertainty of N m .
The scale ν r separates short distances, where logs are resummed in the coupling (r < 1=ν r ), from long distances, where the coupling is evaluated at the fixed scale ν (r > 1=ν r ). If ν r ¼ ∞, this is equivalent to compute with a fixed scale over all distances; if ν r ¼ 0, this is equivalent to compute the coupling at 1=r over the full distance range.
The renormalon factorization scale, ν f , must be chosen low enough that the subtracted mass, δm RS 0 , does not jeopardize the power counting, i.e., δm RS 0 must be of order mv 2 or smaller, but also large enough that δm RS 0 encompasses the renormalon, i.e., the renormalon subtracted series converges, and perturbation theory holds. In our analysis, we observe that we can use the rather low value ν f ¼ 1.0 GeV and yet achieve renormalon cancellation.
Other choices of ν f are possible, but, given the above constraints, the allowed range of variation for ν f is even more restricted than for ν us . In Refs. [20,23] the effect of taking ν f ¼ 0.7 GeV has been considered. The impact on the bottomonium mass is at most 1%. We consider this to be a reasonable upper limit also for the transition widths. The uncertainty coming from the scale ν f (as well as the one from the scale ν us considered before) is, therefore, negligible with respect to the one coming from the scale ν, which is, on the overall, the largest theoretical uncertainty in our computation.
We can look at the effects on the leading order transition width, i.e., the matrix element of the leading order E1 operator (3), when incorporating the static potential (50) at different perturbative orders into the exact solution of the Schrödinger equation. Differently from the previous section, now the Schrödinger equation with the potential (50) can be solved, beyond LO, only numerically. We provide some details on the numerical solution of the Schrödinger equation in Appendix.
Let us consider, as an example, the transition χ b1 ð1PÞ → γϒð1SÞ; the other transitions at leading order follow from this one just by rescaling all the curves by the constant factor ðk γ =423 MeVÞ 3 , which corrects for the photon energy. The left panel of Fig. 10 shows the leading order transition rate when the coupling in the static potential is computed at the fixed scale ν, corresponding to the case ν r ¼ ∞ in Eq. (50). Solving the Schrödinger equation with only the Coulomb-like term in the static potential gives back the same LO result as in Sec. II. This decay rate (solid blue curve) depends strongly on the renormalization scale: It ranges from 18 keV to 72 keV when running ν from 1 GeV to 3 GeV. However, the ν-scale dependence becomes mild as NLO (dashed orange curve), NNLO (dot-dashed green curve) and NNNLO (dotted red curve) radiative corrections to the static potential are added to the Schrödinger equation. Indeed, the decay rate changes only of about 4 keV over the considered ν-range, when the three loop static potential is considered. Moreover, the convergence of the perturbative series has improved with respect to the fixed order case. Convergence tends to worsen only for low ν.
The right panel of Fig. 10 shows the same quantity when the coupling in the static potential is computed at the scale 1=r for r < 1.0 GeV −1 and at the scale ν for r > 1.0 GeV −1 , corresponding to the case ν r ¼ 1.0 GeV in Eq. (50). The perturbative series appears to converge over the whole range 1 GeV ≤ ν ≤ 3 GeV, and in particular for low ν. As for the curves in the left panel, also for the curves shown in the right panel the dependence on the renormalization scale becomes mild with increasing order: At NNNLO the decay rate changes by less than 8 keV when ν goes from 1 to 3 GeV, which is slightly more than for the corresponding NNNLO decay rate in the left panel.

B. Numerical analysis
We are now in the position to discuss the final determinations of the electric dipole transitions χ bJ ð1PÞ → γϒð1SÞ with J ¼ 0, 1, 2 and h b ð1PÞ → γη b ð1SÞ. We use wave functions obtained from the solution of the Schrödinger equation with the full static potential (50). The static potential is taken at three loops, Eq. (41), including ultrasoft effects. The leading order renormalon is subtracted according to the RS 0 scheme defined in Eqs. (45)- (49). The relevant factorization scales are set to be ν us ¼ ν r ¼ ν f ¼ 1.0 GeV. Higher order corrections of relative order v 2 come from higher order electromagnetic operators in the pNRQCD Lagrangian, terms in (4) and (5) other than R S¼1 21 ðJÞ and R S¼0 21 , respectively, and from higher order relativistic corrections affecting initial and final states, terms contributing to R S¼1 21 ðJÞ and R S¼0 21 in (4) and (5), respectively, and stemming from Eq. (19).
The decay width for the χ b0 ð1PÞ → γϒð1SÞ transition is shown in Fig. 11. The leading order (full V s ) nonrelativistic decay rate is the dashed blue curve, the dot-dashed orange curve includes relativistic contributions stemming from higher order electromagnetic operators and the solid black one includes both contributions from higher order electromagnetic operators and relativistic corrections to the wave functions of the initial and final states.
The leading order decay width depends weakly on the renormalization scale: It varies from Γ ≈ 26 keV at ν ¼ 1 GeV to Γ ≈ 31 keV at ν ¼ 3 GeV. This feature is preserved when higher order electromagnetic operators are included and also in the final result. In fact, the ν-dependence of the final result, which is about 3 keV, is weaker than that of the leading order result and also weaker than that obtained from including only higher order electromagnetic operators. A variation of 3 keV over a central value of about 28 keV represents an uncertainty of about 11% in our determination of the decay rate. Moreover, higher order electromagnetic operators and relativistic corrections to the initial and final states provide relatively small changes to the LO transition width.
The gray error band accounts for the uncertainty due to the unknown higher order terms in the perturbative expansion. This is computed, here and in the following plots, by taking the largest between the variation of the result with the scale and one half of the maximum difference between the leading order and the final result, as described in the previous section after Eq. (35).
An interesting feature of Fig. 11 is that the corrections induced by higher order electromagnetic operators diminish the LO decay rate, whereas relativistic corrections to the initial and final states increase it. As a result, at the renormalization scale ν ¼ 1.25 GeV, the value of the decay width Γðχ b0 ð1PÞ → γϒð1SÞÞ turns out to be very similar to the LO result. This will not be the case for the other transitions.
We have performed the same analysis for the electric dipole transitions χ b1 ð1PÞ → γϒð1SÞ and χ b2 ð1PÞ → γϒð1SÞ in Figs. 12 and 13, respectively. Similar features, as the one observed in the χ b0 ð1PÞ → γϒð1SÞ case, are seen here, too. However, we notice that the effect due to relativistic corrections to the initial and final states is a factor 2-3 larger in these cases. We also observe that the final decay rates for χ b1 ð1PÞ → γϒð1SÞ and χ b2 ð1PÞ → γϒð1SÞ show a weaker dependence on the renormalization scale than for χ b0 ð1PÞ → γϒð1SÞ. The scale variation for χ b1 ð1PÞ → γϒð1SÞ is ≲8%, and the scale variation for χ b2 ð1PÞ → γϒð1SÞ is ≲5%. Finally, we remark that for the decay width Γðχ b2 ð1PÞ → γϒð1SÞÞ, the LO result at the renormalization scale ν ¼ 1.25 GeV is outside the final result error band.
2. h b ð1PÞ → γη b ð1SÞ Figure 14 shows the results for the h b ð1PÞ → γη b ð1SÞ transition. The corrections to the decay width induced by higher order electromagnetic operators are very similar to the ones obtained in the previous cases. Their effect is to reduce the LO decay rate by about 2-3 keV. However, the effect due to relativistic corrections to the initial and final state wave functions is larger for the h b ð1PÞ → γη b ð1SÞ transition than for the three transitions considered before. In particular, the decay width changes from about 52 keV (dot-dashed orange curve) to about 63 keV (solid black curve) at ν ¼ 1.25 GeV. This is because the initial and final state bottomonia in the transition are spin-singlet states and thus many corrections to the wave functions, like those induced by the spin-orbit, spin-spin and tensor potentials, are absent. In the case of the χ bJ states, since they are spintriplets, these corrections appear and tend to compensate other relativistic corrections due to different relative signs.
Similarly to the case of the χ b1 ð1PÞ and χ b2 ð1PÞ electric dipole transitions, also the decay width of the h b ð1PÞ → γη b ð1SÞ transition displays a very weak dependence on ν. The rate varies by a mere 1 keV along the whole range of the renormalization scale studied herein. For this reason and for the one given in the paragraph above, the h b ð1PÞ → γη b ð1SÞ decay width appears to be a well suited observable for studying relativistic corrections to the heavy quarkonium wave function. However, the uncertainty due to possible higher order corrections in the perturbative expansion, estimated by looking at one half of the maximum difference between the leading order and the final result, is about six times larger than the one coming from the scale variation in the transition. It is also larger than in the case of the χ bJ transitions. This reflects in a larger final theoretical uncertainty. A related feature is that for the decay width Γðh b ð1PÞ → γη b ð1SÞÞ, the LO result is outside the final result error band for ν ¼ 1.25 GeV. As in the previous section, we choose to set the central value of the decay widths at the scale that self-consistently solves the Bohr-like radius equation (35). This scale is ν ¼ 1=a ¼ 1.25 GeV.

C. Summary and comparisons
Our final results for the electric dipole transitions χ bJ ð1PÞ → γϒð1SÞ, with J ¼ 0, 1, 2, and h b ð1PÞ → γη b ð1SÞ, at relative order v 2 in the counting scheme adopted in this section that consists in treating the whole static potential as a leading order contribution, read Because of the very mild dependence on the renormalization scale, and the good convergence of the perturbative series, the results appear solid and their associated uncertainties are small. The uncertainties correspond to the gray bands shown in Figs. 11-14, and have been computed as described after Eq. (35). In the plots, the errors have not been rounded. We compare our results with those obtained in several other theoretical approaches in Table II. These are a nonrelativistic constituent quark model (CQM) [7], a relativistic quark model (R) [45], a study based on the Godfrey-Isgur model (GI) [46], a study based on the Buchmüller-Tye potential model (BT) [47], a light-front quark model (LFQM) [48], and a screened potential model with zeroth-order wave functions (SNR 0 ) and first-order relativistically corrected wave functions (SNR 1 ) [49]. Reference [47] does not provide a prediction for the h b ð1PÞ → γη b ð1SÞ width, whereas Ref. [48] is restricted to the study of the h b ð1PÞ → γη b ð1SÞ transition only. Our results agree well with those of other approaches for the χ b0 ð1PÞ → γϒð1SÞ and χ b1 ð1PÞ → γϒð1SÞ transitions. In the case of the χ b2 ð1PÞ → γϒð1SÞ transition our result is slightly larger than the bulk of the other predictions, whereas in the case of the h b ð1PÞ → γη b ð1SÞ transition it is significantly larger. The reasons for the differences may be diverse, and follow from the theoretical approaches  Table II being, to various degrees, phenomenological models that neither include QCD corrections in a systematic way, nor derive their parameters from QCD. Hence, they differ from our model independent determination in more than one way. For example, Refs. [7,46] do not include spin-independent 1=m and 1=m 2 potentials, while Refs. [47,49] miss the 1=m potential. Our final results (51)-(54) are predictions, as the bottomonium P-wave E1 transition widths have not been measured so far. In fact, for these electromagnetic transitions only the branching fractions are known, while there are no measurements of any of the total decay widths of the χ bJ , with J ¼ 0, 1, 2, and h b states. Nevertheless, we can use the branching fractions given by the PDG [6] and our results for the decay rates of the electric dipole transitions to predict the total decay widths of the χ bJ ð1PÞ and h b ð1PÞ bottomonia. The results are given in Table III, where the errors are obtained via standard Gaussian uncertainty propagation. The Belle collaboration has reported an upper limit on the total decay width of the χ b0 ð1PÞ at 90% confidence level [50]: Γðχ b0 ð1PÞÞ < 2.4 MeV, which is compatible with our prediction.

IV. CONCLUSION
We have computed the electric dipole transitions χ bJ ð1PÞ → γϒð1SÞ, with J ¼ 0, 1, 2, and h b ð1PÞ → γη b ð1SÞ, within potential nonrelativistic QCD, assuming that the typical binding energy scale, mv 2 , is much larger than Λ QCD , where m is the mass of the heavy quark and v its relative velocity. Consequences of this assumption are that n ¼ 2, l ¼ 1 bottomonia are taken as weakly coupled bound states, and that nonperturbative effects are smaller than the accuracy reached in the calculation. This assumption would not be suited for n ¼ 2, l ¼ 1 charmonia.
The precision that we have reached in this paper is k 3 γ =ðmvÞ 2 × Oðv 2 Þ, k γ being the photon energy. At relative order v 2 we have included higher order electromagnetic interactions in the pNRQCD Lagrangian and higher order relativistic corrections to the initial and final state bottomonia, due to 1=m and 1=m 2 potentials, and 1=m 3 relativistic corrections to the kinetic energy. Concerning radiative corrections to the static potential, we have included them in two different counting schemes: in Sec. II, perturbatively, counting higher order corrections as perturbations of the leading order Coulomb-like potential, and, in Sec. III, nonperturbatively, counting all known terms in the perturbative expansion of the static potential as leading order and including them in the numerical solution of the Schrödinger equation for the initial and final state wave functions.
We summarize the main conclusions drawn from the first scheme. (i) The decay widths show a strong dependence on the renormalization scale ν. At leading order, the strong dependence is due to the running of α s ðνÞ, which affects primarily the Bohr-like radius entering the initial and final state wave functions. At higher orders a significant ν-dependence persists, due to the corrections to the initial and final state wave functions induced by the radiative corrections of the static potential. The static potential contains terms proportional to powers of lnðνrÞ that become large at low values of ν. (ii) Most of the corrections to the decay rates induced by the 1=m and 1=m 2 potentials are relatively small and do not change much as a function of the renormalization scale ν. The largest contributions come from the 1=m potential and the spin-spin one, especially for low values of the scale ν. (iii) The convergence of the perturbative series for all the studied electric dipole transitions is poor. This indicates that bottomonium 1P  [11][12][13][14] and the column NNLO our final results (solid black curves in Figs. [11][12][13][14], both taken at ν ¼ 1.25 GeV. We compare them with those reported by a nonrelativistic constituent quark model (CQM) [7], a relativistic quark model (R) [45], a study based on the Godfrey-Isgur model (GI) [46], a study based on the Buchmüller-Tye potential (BT) [47], a light-front quark model (LFQM) [48], and a screened potential model with zeroth-order wave functions (SNR 0 ) and first-order relativistically corrected wave functions (SNR 1 ) [49]. All decay widths are given in units of keV. states are difficult to accommodate in this scheme, an observation that led us to adopt for our final analysis the second scheme.
In the second scheme, the Schrödinger equation is solved at leading order with all known terms of the perturbative static potential included, i.e., up to three loops. Further, we subtract to the static potential the leading order renormalon and resum at short distances potentially large logs of the type lnðνrÞ. The main effects are: (i) The leading order decay rates depend weakly on the renormalization scale and this is also so when higher order electromagnetic operators and relativistic corrections to the initial and final states are taken into account at relative order v 2 . (ii) Both Oðv 2 Þ corrections do not change much as functions of the renormalization scale and produce corrections to the leading order decay widths that are relatively small. (iii) The corrections induced by higher order electromagnetic operators tend to diminish the leading order decay rates, whereas the opposite effect is found for the relativistic corrections to the initial and final state wave functions. These observations support our initial assumptions on the nature of the 1P bottomonia. Because the perturbative series appears convergent and only mildly dependent on the renormalization scale, the final results are affected by small uncertainties.
If the most critical of our assumptions, mv 2 ≫ Λ QCD , is relaxed to mv 2 ∼ Λ QCD , then nonperturbative corrections may become as large as the Oðv 2 Þ corrections considered above. Since the uncertainties on our final results have been chosen to include one half of the Oðv 2 Þ corrections, the effect of assuming mv 2 ∼ Λ QCD would be (at least) to double our final errors. A challenging alternative is to compute the nonperturbative contributions listed in Ref. [15].

ACKNOWLEDGMENTS
We thank Nora Brambilla, Yuichiro Kiyo, Clara Peset, Antonio Pineda and Yukinari Sumino for numerous informative discussions. This work has been supported by the DFG and the NSFC through funds provided to the Sino-German CRC 110 "Symmetries and the Emergence of Structure in QCD", and by the DFG cluster of excellence "Origin and structure of the universe" (www.universecluster.de). J. S. acknowledges the financial support from the Alexander von Humboldt Foundation and thanks the Technische Universität München for hospitality while most of this work was carried out.

APPENDIX: SOLVING THE SCHRÖDINGER EQUATION
In a generic central potential, VðrÞ, the Schrödinger equation for the reduced wave function, u nl ðrÞ ¼ rR nl ðrÞ, has the form − 1 m d 2 u nl ðrÞ dr 2 þ VðrÞ þ lðl þ 1Þ mr 2 u nl ðrÞ ¼ E nl u nl ðrÞ; in the case of two particles of mass m. This is a one dimensional Schrödinger equation, which has significance only for positive values of r, and must be supplemented by a boundary condition at r ¼ 0. We require that the radial function R nl ðrÞ remains finite at the origin, which implies that u nl ð0Þ ¼ 0.
If close to the origin the potential VðrÞ has the form where p is an integer such that p ≥ −1, we can expand the solution u nl ðrÞ in the vicinity of the origin as Since we are interested in finding bound states, we also impose that where k ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffi ffi mjE nl j p is the wave function number. In numerical applications we introduce short-and longdistance cutoffs, denoted by r in and r fi , respectively, for which we require u nl ðr in Þ ¼ r lþ1 in ; ðA6Þ u nl ðr fi Þ ¼ e −kr fi : The dependence of physical observables on the short distance cutoff is quite sensible, and it can hinder the numerical search of the ground state and its excitations due to the dominance of the irregular solutions at very small values of r. In order to improve on this, we can use, for two different energies E nl ≠ E n 0 l , the orthogonality relation between their bound state wave functions: which follows from multiplying Eq. (A1) by u n 0 l ðrÞ and later subtracting the same equation, but with n and n 0 exchanged. The regularity condition at the origin u nl ðr in Þ ¼ 0 for r in → 0 makes the states automatically orthogonal in the r in → 0 limit. We can further enforce orthogonality also for finite r in by requiring for any two states, meaning that the logarithmic derivative at short distances becomes independent of the principal quantum number. This condition has many advantages, such as the possibility of working with singular potentials at r ¼ 0, like r 2 VðrÞ ¼ AE∞ for r → 0. Moreover, in order to avoid pollution from the irregular solutions, we can use Eq. (A9) to match at an intermediate distance r me the  Fig. 16, but for the radial wave function, R 21 ðrÞ, of the lowest P-wave state. solutions of the one dimensional Schrödinger equation obtained when integrating it from r in to r me , with boundary (A6), and from r me to r fi with boundary (A7). Finally, Eq. (A9) is ideal to find excited states because, as we remarked, the logarithmic derivative at short distances becomes independent of n. For further details we refer to Ref. [51].
In order to solve the differential equation (A1) for the potential (50), we use the fourth order Runge-Kutta algorithm with adaptive step size implemented in FORTRAN77. This implementation automatically takes care of convergence and numerical accuracy. The numerical implementation of the Green function, Eqs. (26) and (27)