Decay and electromagnetic production of strongly coupled quarkonia in pNRQCD

We improve the pNRQCD approach to annihilation processes of heavy quarkonium and make first pNRQCD predictions for exclusive electromagnetic production of heavy quarkonium. We consider strongly coupled quarkonia, i.e., quarkonia that are not Coulombic bound states. Possible strongly coupled quarkonia include excited charmonium and bottomonium states. For these, pNRQCD provides expressions for the decay and exclusive electromagnetic production NRQCD matrix elements that depend on the wavefunctions at the origin and few universal gluon field correlators. We compute electromagnetic decay widths and exclusive production cross sections, and inclusive decay widths into light hadrons for P -wave quarkonia at relative order v2 and leading order, respectively. We also compute the decay widths of 2S and 3S bottomonium states into lepton pairs and their ratios with the inclusive widths into light hadrons at relative order v2.


Introduction
Exclusive electromagnetic processes involving heavy quarkonia are good probes of quarkonium production and decay mechanisms: the clean final states enable accurate measurements [1,2] and nonrelativistic effective field theories [3,4] allow to express physical observables through systematic expansions and factorization formulas. Systematic expansions guarantee a control over the accuracy of the theoretical expressions, and factorization casts non perturbative contributions into few long-distance matrix elements (LDMEs). LDMEs are eventually determined from data. Heavy quarkonium annihilations into light hadrons JHEP04(2020)095 are similarly well under control from the theoretical side, but more difficult to determine experimentally as contributions from decay channels into leptons, photons or heavy quarks have to be subtracted from the total width. Heavy quarkonium production in hadron collisions is challenging both theoretically and experimentally. Heavy quarkonia, like charmonia and bottomonia, are nonrelativistic bound states of a heavy quark and a heavy antiquark. Nonrelativistic effective field theories exploit the typical hierarchy of energy scales characterizing such systems. The energy scales are the heavy quark mass, m, the typical momentum and momentum transfer, mv, and the typical kinetic energy, mv 2 , where v is the relative velocity of the heavy quark and antiquark. Because v 1, the above energy scales are hierarchically ordered. Reference values for v are v 2 ≈ 0.3 for charmonia and v 2 ≈ 0.1 for bottomonia.
Nonrelativistic QCD (NRQCD) is the effective field theory, suited to describe states made of a heavy quark and a heavy antiquark, that follows from QCD by integrating out modes of energy and momentum of order m [3]. In NRQCD, heavy quarkonium annihilation rates into photons, leptons or inclusive annihilation rates into light hadrons, Γ, and exclusive electromagnetic production cross sections, σ, are expressed by sums of products of NRQCD LDMEs, O n , with perturbative short distance coefficients, c n : where d n is the mass dimension of the operator O n . The short distance coefficients are process dependent. The LDMEs depend on the quarkonium state, but not on the process. Whereas the short distance coefficients can be computed as a series expansion in the strong coupling constant α s , this is guaranteed by m Λ QCD , the LDMEs are counted in powers of v. In practice, the factorization formula is truncated at a desired order in α s and v. In the case of decay widths, the short distance coefficients c n are dimensionless, while in the case of exclusive electromagnetic production cross sections, they have mass dimension −3 and depend on m and on the center of mass energy of the collision, √ s. The center of mass energy is, besides the heavy quark mass, the other large scale in production processes.
The NRQCD factorization formulas for quarkonium inclusive annihilation widths into light hadrons, electromagnetic annihilation widths and exclusive electromagnetic production cross sections have been proved to all orders in the expansion parameters. Early determinations of several of the short distance coefficients can be found in [3,5,6] and in the review [7]. These have been constantly improved over the last years (see appendix B and references therein).
The NRQCD LDMEs entering quarkonium annihilation and exclusive electromagnetic production are expectation values of four-fermion operators on the quarkonium state. A list of four-fermion operators relevant for the present work is in appendix A. One important feature of NRQCD is that the quarkonium state can contain contributions not only from the leading color-singlet heavy quark-antiquark Fock state, but also from the subleading Fock states that include effects of dynamical gluons. Because gluons carry color, the heavy quark-antiquark pair in the subleading Fock states can be in a color octet state. Hence, four-fermion operators projecting on color octet states contribute to the JHEP04(2020)095 observables. Determinations of these contributions provide important verifications of the NRQCD factorization formalism.
In order to make quantitative statements based on the NRQCD factorization formulas, it is important to be able to determine the LDMEs. This is a difficult task in the standard NRQCD approach, especially for the LDMEs of higher orders in v, which, as a result, are poorly known. Also the power counting of the LDMEs is not unique as they depend on several, still dynamical, energy scales: mv, mv 2 , the typical hadronic scale Λ QCD , . . . Disentangling these energy scales in a suitable nonrelativistic effective field theory of lower energy than NRQCD provides a way to simplify and in some cases compute the NRQCD LDMEs.
Potential NRQCD (pNRQCD) follows from NRQCD by integrating out modes associated with energy scales larger than mv 2 , regardless of these energy scales being perturbative or non perturbative [8]. If all relevant energy scales are perturbative, the LDMEs can be expressed in pNRQCD as series in α s [9]. If mv 2 Λ QCD , then the LDMEs are non perturbative. In the non perturbative, confining, regime a heavy quark-antiquark pair may bind into a quarkonium, i.e., a quark-antiquark pair in a color singlet configuration, a hybrid, i.e., a quark-antiquark pair in a color octet configuration bound to gluons, a quarkonium in the presence of glueballs, a tetraquark, i.e., a heavy quark-antiquark pair bound in different combinations with a light quark-antiquark pair, and so on. We will consider quarkonia that are well below the open flavor threshold, i.e., separated from it by an energy gap of order Λ QCD or larger. Moreover, lattice computations suggest that the quarkonium spectrum may be separated by an energy gap of order Λ QCD or larger from the spectrum of hybrids and quarkonia plus glueballs [10][11][12]. The distribution of energy levels would then be schematically the one shown in figure 1. If the kinematical condition mv 2 Λ QCD is realized, the higher gluonic excitations can be integrated out, and this leaves the quark-antiquark pair in a color singlet configuration as the only dynamical degree of freedom. In this situation, the LDMEs can be factorized in a wavefunction contribution, which encodes information from the quarkonium state, and some universal correlators of gluon fields, which encode contributions coming from higher excitations of the heavy quarkantiquark pair, those induced by gluons or light quarks and separated by an energy gap of order Λ QCD from the quarkonium spectrum [13,14]. The quarkonium wavefunction is obtained by solving the Schrödinger equation that is the equation of motion of pNRQCD. The quarkonium potential may be expressed in terms of Wilson loops and gluon field insertions on it [15,16]. It includes contributions coming from quarkonium modes of order mv and from higher excitations of the heavy quark-antiquark pair induced by gluons or light quarks and separated by an energy gap of order Λ QCD from the quarkonium spectrum. Under the kinematical condition mv 2 Λ QCD , the potential is sensitive to distances of order 1/Λ QCD , where it is non perturbative. Hence, quarkonium satisfying the condition mv 2 Λ QCD is not a Coulombic bound state. To distinguish it from a Coulombic bound state, which is weakly coupled, a non Coulombic quarkonium state is referred to as strongly coupled. Moreover, pNRQCD in the regime mv 2 Λ QCD , is called strongly coupled pNRQCD. In this work, we will focus on decay and exclusive electromagnetic production of quarkonium states for which the condition mv 2 Λ QCD is fulfilled. Under this condition, we will use pNRQCD to express the LDMEs in terms of quarkonium wavefunctions, JHEP04(2020)095 binding energies and gluonic correlators. Ideally the quarkonium wavefunctions and binding energies should be determined from the solution of the pNRQCD Schrödinger equation, which requires knowing the quarkonium potential from lattice QCD. Also the gluonic correlators should be computed on the lattice. In practice, however, due to the incomplete knowledge of these quantities, wavefunctions and binding energies are determined using potential models and the gluonic correlators from the data. Since the gluonic correlators are non perturbative but universal quantities that do not depend on the heavy quark flavor, the non perturbative parameters needed in pNRQCD are in general fewer than the LDMEs needed in NRQCD. When applicable, the pNRQCD factorization formulas have, therefore, more predictive power than the NRQCD ones.
The condition mv 2 Λ QCD is fulfilled by non Coulombic, strongly coupled, quarkonia. 1 Charmonium and bottomonium states that are possibly non Coulombic bound states are higher excited states whose principal quantum number is greater than one. These states will be the subject of our phenomenological investigations in section 4.
The paper is organized in the following way. We start in section 2 by briefly reviewing strongly coupled pNRQCD. In the following section 3, we compute in strongly coupled pNRQCD the relevant LDMEs and give their explicit expressions in terms of quarkonium wavefunctions, binding energies and gluonic correlators. Four-fermion operators are listed in appendix A and NRQCD factorization formulas in appendix B. Details of the computation are in appendix C. Some of the results presented in section 3 correct previous findings of ref. [14]. Using these results and after having determined the wavefunctions and binding energies in several potential models, we fit the unknown gluonic correlators and compute decay widths and exclusive electromagnetic production cross sections for charmonium 1P states and bottomonium 2S, 3S, 1P , 2P and 3P states in section 4. We conclude in section 5.

JHEP04(2020)095
2 Effective field theories for strongly coupled quarkonia We compute the LDMEs of NRQCD assuming that the quarkonium states that we consider are well below the open flavor threshold and satisfy the condition mv 2 Λ QCD . We further assume that higher gluonic excitations of the heavy quark-antiquark pair are separated by an energy gap of order Λ QCD or larger from the ground state; this assumption is supported by lattice calculations that show the excitation spectrum of the gluon field around a static quark-antiquark pair separated by a large energy gap from the ground state [10][11][12]. It follows that we can picture the distribution of the energy levels as illustrated in figure 1. Such picture allows us to describe the quarkonium spectrum in an effective field theory where all modes associated to the excitations of the heavy quark-antiquark pair induced by gluons or light quarks and separated by an energy gap of order Λ QCD from the quarkonium spectrum have been integrated out. This effective field theory is strongly coupled pNRQCD [4].
In strongly coupled pNRQCD, the quarkonium potential and the LDMEs are computed by quantum mechanical perturbation theory, order by order in 1/m, around the NRQCD static solution. Each power of 1/m is suppressed by v or Λ QCD /m. The computation of the quarkonium potential in strongly coupled pNRQCD has been first performed in refs. [15,16], and the NRQCD LDMEs have been computed in refs. [13,14,17]. In this section, we briefly review the formalism and compute the potentials relevant for the present work. We compute the LDMEs in section 3.

NRQCD
The degrees of freedom of NRQCD are heavy quark and antiquark fields, ψ and χ, describing modes of energy and momentum smaller than m, gluons and light quarks. The NRQCD Hamiltonian, H NRQCD , can be organized as an expansion in 1/m, so that H NRQCD = H (0) Boldfaced characters indicate three-dimensional vectors. The fields ψ and χ respectively annihilate a heavy quark and create a heavy antiquark, the fields q k are n f massless quark fields, are the chromoelectric and chromomagnetic fields, respectively, G µν a T a = G µν is the gluon field strength tensor, σ i are the Pauli matrices and c F is a short distance coefficient, which is known up to three loops [18]. The fields ψ and χ satisfy the canonical equal time anticommutation relations: The mass m is the heavy quark pole mass.

JHEP04(2020)095
We restrict ourselves to the one-quark-one-antiquark sector of the NRQCD Fock space, where quarkonium states live. In this sector, we denote an energy eigenstate of the NRQCD Hamiltonian by |n; x 1 , x 2 , where n represents a generic set of conserved quantum numbers, and x 1 and x 2 are the positions of the quark and antiquark, respectively. The heavy quark and antiquark positions are conserved quantum numbers in the static limit. We normalize the states as n; . The eigenstates satisfy the Schrödinger equation: The ground state, n = 0, is the Fock state made of a heavy quark-antiquark pair without excitations induced by gluons or light quarks and separated by an energy gap of order Λ QCD from the quarkonium spectrum; other values of n identify heavy quark-antiquark pairs in the presence of such excitations. The functions E n are the corresponding energies. In the static limit the above equation becomes identifies the ground state static energy between a quark and antiquark in a color singlet configuration, which is well approximated by a Cornell-like potential (see section 4.1), E (0) n for n = 0 are the static energies of a quark-antiquark pair in the presence of the excitations described above. The ground state energy and the first gluonic excitation are schematically shown in figure 1. The eigenstates |n; x 1 , x 2 (0) are normalized as (0) n; Having assumed that the gap between the lowest-lying energy, E 0 , and the higher ones is of order Λ QCD mv 2 , we can compute |0; x 1 , x 2 by expanding in 1/m around the static solution |0; x 1 , x 2 (0) : Corrections are obtained from quantum mechanical perturbation theory applied to the NRQCD Hamiltonian. In particular, |0; x 1 , x 2 (1) , which is the correction to the eigenstate at order 1/m, reads where H NRQCD is given in eq. (2.2). The explicit expression for the correction term |0; x 1 , x 2 (1) has been obtained in refs. [15,16]. Here, we briefly list the main ingredients to obtain it, as we will use them to compute the LDMEs. In order to obtain |0; x 1 , x 2 (1) , we need to evaluate the matrix JHEP04(2020)095 . The first step is to make the quark content of the eigenstates explicit: where |n; x 1 , x (0) encodes the light degrees of freedom content of |n; x 1 , x 2 (0) . The states : The states |n; x 1 , x 2 (0) do not contain heavy (anti)quarks, hence they are annihilated by ψ and χ † . This implies the normalization (0) n; x 1 , x 2 |m; x 1 , x 2 (0) = δ nm . Then, the matrix elements (0) n; z 1 , z 2 |H (1) NRQCD |0; x 1 , x 2 (0) can be computed by using Wick's theorem, which removes the quark and antiquark fields leaving delta functions that constrain x 1 = x 1 and x 2 = x 2 . After the quark and antiquark fields have been removed in this way, we can use the following shorthands without ambiguity: Finally, the following identities can be used to simplify the matrix elements and, for n = k, (2.11) Equations (2.9) follows from symmetry considerations, and eqs. (2.10) and (2.11) may be derived from canonical commutation relations [4,15]. The parentheses on the right-hand sides of eqs. (2.10) imply that the derivatives act only on E (0) n . These ingredients are sufficient to derive the explicit expression of |0; x 1 , x 2 (1) .
In this work, the state |0; x 1 , x 2 (1) will turn out to be relevant only for P -wave states.
For P -wave states only terms containing derivatives acting on the wavefunctions give nonvanishing contributions. For our purposes it is sufficient, therefore, to isolate from |0; x 1 , x 2 (1) only this part, which we denote with |0; P -wave . Its explicit expression reads where O 1(2) (t)W r×T means that the fields appearing on the left of W r×T are inserted at a time t on the quark (antiquark) line of the Wilson loop. The Wilson loop and field insertions on it are traced in color space.

Strongly coupled pNRQCD
The degrees of freedom of strongly coupled pNRQCD, i.e., the degrees of freedom that are resolved at an energy scale of order mv 2 are only color singlet heavy quark-antiquark pairs, if we neglect the interaction with light hadrons of energy and momentum of order mv 2 or smaller (these are the Goldstone bosons of the chiral symmetry, for a discussion see ref. [14]). The reason is that, having assumed mv 2 Λ QCD , the scale mv 2 is below the confinement scale, Λ QCD . Hence, the pNRQCD Hamiltonian has the very simple form where S annihilates a color singlet heavy quark-antiquark field. It satisfies the canonical equal time commutation relation: Because of the assumed energy gap of order Λ QCD between the ground state and the higher excitations of the heavy quark-antiquark pair, these have been integrated out when matching to strongly coupled pNRQCD, whereas the NRQCD ground state, |0; x 1 , x 2 , matches the pNRQCD state made of one color singlet heavy quark-antiquark pair, S † (x 1 , x 2 )|vac . If we neglect light hadrons of energy and momentum of order mv 2 or smaller, then |vac is the pNRQCD vacuum.

JHEP04(2020)095
The function h can be organized as an expansion in 1/m: h = h (0) + h (1) /m + . . . . The terms h (0) , h (1) , . . . can be determined by matching h with the NRQCD ground state energy order by order in 1/m: (2.19) and so on. Also in this case, like in the NRQCD case, the mass m should be understood as the heavy quark pole mass. The term h (0) is the static quark-antiquark potential V (0) , which, according to (2.14), can be determined from the large time behavior of the static Wilson loop: The term h (1) contains the quark and antiquark kinetic energies, and the 1/m potential V (1) , which, according to (2.17), may be computed from a Wilson loop with two chromoelectric field insertions: where R = (x 1 + x 2 )/2 is the center of mass coordinate, r = x 1 − x 2 is the quarkantiquark distance, P is the center of mass momentum and ε nJLS is the binding energy. The states, |P |nJLS , are classified according to the center of mass momentum, the principal quantum number, n, the total angular momentum, J, the orbital angular momentum, L, and the spin, S. At leading order, h contains the kinetic energy, −∇ 2 1 /(2m) − ∇ 2 2 /(2m), and the static potential, V (0) , which both count like mv 2 , as a consequence of the scale hierarchy in a nonrelativistic bound state and the virial theorem. 2 Also ε nJLS is of order mv 2 . For a strongly coupled bound state, the potential V (1) /m of eq. (2.21) can be of the same order as the static potential if mv ∼ Λ QCD [15,16]. Under this condition it should be included in the leading order h. Whatever the specific regime that we are describing is, the leading order potential is a central potential, i.e., it depends only on r = |x 1 − x 2 |. Hence, the leading order binding energy may be classified in terms of n and L only. The corresponding Schrödinger equation reads The momenta −i∇1 and −i∇2 may be decomposed in a center of mass momentum and a relative momentum: −i∇1 = −i∇r − i∇ R /2 and −i∇2 = i∇r − i∇ R /2. The relative momentum, −i∇r, scales like mv, while the center of mass momentum, −i∇ R , scales like the momentum of the dynamical low-energy degrees of freedom of the effective theory. Hence, the center of mass momentum scales at most like mv 2 .
In the effective field theory of eq. (2.17) that does not contain dynamical low-energy degrees of freedom besides the quark-antiquark color singlet, the reference frame may be always chosen so that the center of mass momentum is set to zero. The kinetic energy entering the leading order pNRQCD Hamiltonian is, therefore, only the kinetic energy associated with the relative momentum: −∇ 2 r /m.
In summary, a generic strongly coupled quarkonium state |H with quantum numbers n, J, L and S and center of mass momentum P is described in pNRQCD by a state The wavefunction R|P is equal to 1 in the center of mass frame, P = 0. The factor 1/ P = 0|P = 0 normalizes the state. The wavefunction r|nJLS is the solution of the Schrödinger equation (2.22), whose static potential is given by (2.20), 1/m potential by (2.21) and so on [16]. The field S † (x 1 , x 2 ) creates a heavy quark-antiquark pair in a color singlet configuration.

LDMEs
Four-fermion operators show up in NRQCD at order 1/m 2 or higher. Some of them are listed in appendix A. They match into contact terms of pNRQCD. The matching condition reads where O n is a four-fermion operator in the NRQCD Lagrangian, d n is its dimension and V (dn−4) On is a dimension d n − 3 contact term, i.e., a function of δ 3 (r) or its derivatives. 3 From eqs. (2.24) and (3.1) it follows that the LDME of a generic four-fermion operator of the type listed in appendix A, which includes decay and exclusive electromagnetic production LDMEs but excludes hadronic production LDMEs, for a strongly coupled quarkonium H of quantum numbers n, J, L and S, at rest (P = 0), can be expressed in strongly coupled pNRQCD by means of the master formula [14]:

JHEP04(2020)095
Because V (dn−4) On is a contact term, the wavefunctions r|nJLS and their derivatives contribute to eq. (3.2) only at the origin r = r = 0. In the particular case of P -wave states, since their wavefunctions vanish at the origin, the contact term V (dn−4) On must contain a sufficient number of derivatives acting on the wavefunctions in order to give a nonvanishing contribution. In appendix C we compute some relevant LDMEs in strongly coupled pNRQCD from the master formula (3.2). The results are listed and discussed in the following section.

LDMEs in strongly coupled pNRQCD
The LDMEs appearing in the NRQCD factorization formulas for the quarkonium decay widths and electromagnetic production cross sections, see appendix B, involve four-fermion operators of the type listed in appendix A. These LDMEs can be evaluated from the master formula (3.2), which holds when strongly coupled pNRQCD is valid, i.e., for quarkonium states satisfying the condition mv 2 Λ QCD . Equation The chromoelectric fields are evaluated at the same location x 1 = x 2 . The quantity E n is a gluonic matrix element that can be conveniently expressed in terms of a correlator of two chromoelectric fields [14]: where Φ ab (t, 0) is a straight Wilson line in the adjoint representation connecting the points (t, 0) and (0, 0) and T F = 1/2 is the normalization of the color matrices. The Wilson line ensures the gauge invariance of E n . We have used that |0 (0) | x 1 =x 2 = 1 c |vac / √ N c , 1 c being the SU(N c ) identity matrix and N c = 3 the number of colors. Note that the correlator vac|gE i,a (t, 0)Φ ab (t, 0)gE i,b (0, 0)|vac may be understood as the r → 0 limit of a Wilson loop with two chromoelectric field insertions. As we have seen in section 2.2, Wilson loops with field insertions are related to the quarkonium potential. We will exploit this observation in the following.
The power counting of gluon field correlators and their time integrals is obvious, for they depend on only one scale: Λ QCD . Hence they scale like Λ QCD to their dimension. For instance, we have that E n ∼ Λ 3−n QCD . Once the relevant contact terms have been computed and expressed in terms of delta functions at r = 0 and field strength correlators, eq. (3.2) allows to compute the LDMEs. Because of the delta functions at r = 0 the LDMEs will depend on the wavefunction r|nJLS or its derivatives computed at the origin. If the four-fermion operator and/or JHEP04(2020)095 corrections to the state (see eq. (2.12)) contain derivatives, one can generate Laplacian operators acting on the wavefunction at the origin. Such terms can be rewritten in terms of the binding energy of the state by using the Schrödinger equation (2.22): where the dots stand for higher order terms in the 1/m expansion of h. In dimensional regularization, V (0) (r = 0) vanishes as the static potential is purely perturbative at short distances [19] (for a more recent analysis with the same conclusion see [20]) and, therefore, its Fourier transform is scaleless. The situation is different for V (1) (r), which at r = 0 reduces to V (1) (r = 0) = −E 1 (compare eq. (2.21) with eq. (3.4) in the r → 0 limit, taking into account that W r×T | r=0 = 1). Therefore we have in dimensional regularization The first neglected correction in the right-hand side is suppressed by a factor of order Λ QCD /m with respect to E 1 , which is of order Λ 2 QCD . Equation (3.6) corrects an analogous expression obtained and used in ref. [14], where the contribution from V (1) (r = 0) was set to zero. We will show in section 3.1.2 how this modifies some of the LDMEs obtained in [14]. Note that, since mε nJLS is of order (mv) 2 , E 1 is smaller than mε nJLS if mv Λ QCD . Therefore, one can neglect at leading order the term E 1 in the right-hand side of (3.6), if the examined quarkonium state fulfills the kinematical condition mv Λ QCD mv 2 . After matching the contact terms, evaluating the LDMEs with the master formula (3.2) and rewriting the Laplacian of the wavefunction by means of (3.6), the LDMEs are expressed in terms of the quarkonium wavefunctions at the origin, correlators of field strength tensors and the quarkonium binding energies. Because the angular dependence of the wavefunctions is know, we will use the radial parts of the wavefunctions, rather than the wavefunctions. We will denote them with R nJLS (r) (R (0) nL (r) at leading order). We note that the correlators of field strength tensors are universal non perturbative parameters, since they do not depend neither on the heavy quark nor on the quarkonium state. Hence they may be fixed on some set of observables and used in some other one, even involving heavy quarks of different flavor. As it has been noted in ref. [14], this leads eventually to a reduction in the number of non perturbative parameters needed to describe quarkonium decay widths and electromagnetic production cross sections in strongly coupled pNRQCD in comparison to the number of LDMEs required in NRQCD.

P -wave LDMEs
Here we list some relevant P -wave LDMEs computed in strongly coupled pNRQCD. We consider a generic spin one quarkonium state that is a P -wave state with principal quantum number n and total angular momentum J made of a heavy quark-antiquark pair of flavor Q: χ QJ (nP ). 4 Details can be found in appendix C.

JHEP04(2020)095
The hadronic LDME χ QJ (nP )|O 1 ( 3 P J )|χ QJ (nP ) reads in pNRQCD at leading order in the v and Λ QCD /m expansion: where the hadronic operators O 1 ( 3 P J ) have been defined in eqs. (A.18)-(A.20) and R nJLS stands for the derivative of R nJLS . We have computed the corresponding electromagnetic LDME χ QJ (nP )|O em 1 ( 3 P J )|χ QJ (nP ) at next-to-leading order: where the electromagnetic operators O em 1 ( 3 P J ) have been defined in eqs. (A.3)-(A.5). The expressions (3.8) and (3.7) agree at leading order. The leading order expressions are known since ref. [3]. Instead, the correction proportional to iE 2 /m in eq. (3.8) is new. This is the dominant correction to the pure wavefunction contribution. It is of order Λ QCD /m. As we detail in appendix C it originates from the 1/m correction to the quarkonium Fock state given in eq. (2.12).
The electromagnetic LDME χ QJ (nP )|T em 8 ( 3 P J )|χ QJ (nP ) reads in pNRQCD at leading order in the v and Λ QCD /m expansion: where the electromagnetic operators T em 8 ( 3 P J ) have been defined in eqs. (A.11)-(A.13). The result (3.9) agrees with the result of ref. [14] for the J = 0 case. Note that the operator T em 8 ( 3 P J ) has no overlap with the color singlet component of the heavy quarkantiquark pair. Hence, the expression of its matrix element in terms of the quarkonium wavefunction is a specific feature of strongly coupled pNRQCD: the above expression has no equivalent in ref. [3].
The electromagnetic LDME χ QJ (nP )|P em 1 ( 3 P J )|χ QJ (nP ) reads in pNRQCD in dimensional regularization: where the electromagnetic operators P em QCD . The expression follows from having used eq. (3.6). The result is new. Since ), the right-hand sides of eqs. (3.7)-(3.10) are independent of J up to the computed corrections. This verifies the heavy-quark spin symmetry. In particular, the symmetry is realized also for the octet matrix element (3.9).

S-wave LDMEs
The inclusion of the term E 1 in eq. (3.6) modifies the expression of some of the S-wave color singlet LDMEs computed in ref. [14]. The modification is relevant at relative order JHEP04(2020)095 where V Q (nS) (P Q (nS)) is an S-wave vector (pseudoscalar) quarkonium state made of a heavy quark and antiquark of flavor Q. The operators are four-chromoelectric field correlators of mass dimension two, whose definition can be found in ref. [14] but is irrelevant for the present work.
The (leading order) binding energy ε (0) n0 scales like mv 2 , the correlator E 3 is a scaleless constant, whereas all other correlators in eqs. (3.11)-(3.15) scale like Λ 2 QCD . Hence the term proportional to the binding energy is the dominant correction in eqs. (3.11)-(3.14) if the quarkonium state satisfies the condition mv Λ QCD mv 2 . All corrections are of the same order if mv ∼ Λ QCD .
The inclusion of the term E 1 in eq. (3.6) has modified eqs. (3.11)-(3.14) with respect to ref. [14] by adding the term proportional to −2E 1 E 3 /(9m 2 ). This term is of order (Λ QCD /m) 2 . It has also modified eq. (3.15). Differently from the version in ref. [14], eq. (3.15) does not contain the term −N c |R We note that eqs. (3.11)-(3.14) are accurate up to relative order v 2 , hence also the wavefunctions at the origin include corrections of relative order v 2 . These corrections distinguish the vector from the pseudoscalar radial wavefunctions. On the contrary, eq. (3.15) is accurate only at leading order, hence the wavefunction appearing there needs not to be more accurate than that. At leading order the radial parts of the vector and pseudoscalar wavefunctions are equal.

JHEP04(2020)095
For completeness we reproduce here also the S-wave color octet LDMEs of the op- These relations are valid at leading order in the velocity and Λ QCD /m expansion.

Gremm-Kapustin relations
The NRQCD equations of motion imply that some LDMEs are related at leading order in the velocity expansion. These relations are often referred to as Gremm-Kapustin relations [22]. Over the years, following the same method, more relations have been derived, see, for instance, refs. [23][24][25].
The Gremm-Kapustin relations are automatically satisfied by the expressions of the LDMEs derived in strongly coupled pNRQCD in sections 3.1.1 and 3.1.2. The reason is that at the level of pNRQCD the information encoded in the equations of motion has been implemented through the Schrödinger equation, or, more specifically through eq. (3.6).
In particular, for electromagnetic P -wave LDMEs, from eqs. (3.8), (3.9) and (3.10) it follows that This relation was first derived in ref. [23]. According to the expressions of the LDMEs in strongly coupled pNRQCD, it holds at the orders |R nJ11 (0)| 2 × (mv) 2 and |R nJ11 (0)| 2 × Λ 2 QCD . This is the leading order if mv ∼ Λ QCD , but goes beyond it if mv Λ QCD mv 2 . For S-wave LDMEs, from eqs. (3.11)-(3.15) it follows that at leading order in both regimes, mv ∼ Λ QCD and mv Λ QCD mv 2 , This relation was first derived in ref. [22]. We note that ref. [14] could reproduce this relation only in the regime mv Λ QCD mv 2 .

JHEP04(2020)095 4 Fits, analyses and results
We apply now the pNRQCD factorization of the LDMEs computed in the previous section to the analysis of some quarkonium decay and production observables, in particular electromagnetic ones. The strategy that we will pursue is the following: first we determine the quarkonium wavefunctions and binding energies by means of models, then we fit the relevant chromoelectric field correlators on charmonium data and finally we compute the observables. Because the correlators are universal we are in the position to predict observables in the bottomonium sector. We focus on two sets of observables that depend on two distinct sets of correlators. In the first part of the section, we compute quarkonium P -wave electromagnetic decay widths and production cross sections. In the second part, we analyze quarkonium P -wave widths for inclusive decays into light hadrons, and bottomonium S-wave widths for decays into lepton pairs. We use the strongly coupled pNRQCD factorization formulas in their regime of validity, i.e., for non Coulombic quarkonium states. For this reason we limit ourselves to states with principal quantum number greater than one.
The section is organized as follows. In section 4.1 we establish some reasonable values for the quarkonium wavefunctions at the origin and the binding energies by comparing several potential models. In section 4.2 we fit the correlators E 1 and iE 2 on the P -wave charmonium electromagnetic decay widths and the recently measured cross section σ(e + e − → χ c1 (1P ) + γ). We also compute these quantities within our framework. In section 4.3, we compute the P -wave bottomonium electromagnetic decay widths and electromagnetic cross sections from the determined correlators. In section 4.4, we compute the P -wave charmonium widths for inclusive decays into light hadrons and fit the correlator E 3 . We also compute with this information P -wave bottomonium widths for inclusive decays into light hadrons. Finally, in section 4.5, we use the determination of the correlator E 3 to compute, under some assumptions, the leptonic decay widths of the bottomonium S-wave states Υ(2S) and Υ(3S).

Potential models
The first ingredients entering the LDMEs are the quarkonium wavefunctions at the origin and the binding energies. In particular the wavefunctions at the origin are very important as they affect all LDMEs in the pNRQCD formulation and contribute to widths and cross sections at leading order. The uncertainty of the wavefunction at the origin is typically the major source of uncertainty for these observables.
Ideally, quarkonium wavefunctions and binding energies should follow from the solution of the Schrödinger equation (2.22) with the potentials computed within lattice QCD from the corresponding Wilson loops, like the one in eq. (2.20) for the static potential or the one in eq. (2.21) for the 1/m potential. The knowledge of the lattice potentials beyond the static one is, however, incomplete and sometimes poor [26][27][28]. In practice, therefore, one uses potential models, with the idea that tuning the potential model parameters on some observables may provide enough input to mimic the full real dynamics, the one that lattice computations are not yet in the position to provide. Clearly, the use of potential models introduces JHEP04(2020)095  Table 1. For the four potential models described in the text, we list the squared derivative of the radial wavefunction at the origin, |R possibly large and, to some extent, uncontrolled uncertainties. Nevertheless, it has also proved to be successful in many cases, besides being the only available solution at present.
We employ four different potential models to compute wavefunctions at the origin and binding energies for charmonium states. We label these potential models A, B, C, and D. A common feature of them is that they reduce to a liner rising potential at long distances.
Model A is the Cornell potential model of refs. [29,30], where V (0) (r) = −κ/r + σr, with κ = 0.52 and σ = 0.1826 GeV 2 ; σ may be identified with a string tension. In this model, the charm and bottom quark mass parameters are taken to be m c = 1.84 GeV and m b = 5.18 GeV, respectively.
Model B is the frozen α s model of refs. [31,32]. It is similar to model A, except that now κ depends on r. The r-dependent values of κ are tabulated in ref. [31].
Model C is the Buchmüller-Tye potential model [33]. Here the charm and bottom quark mass parameters are taken to be m c = 1.48 GeV and m b = 4.88 GeV, respectively.
Model D is the Cornell potential model in the version of ref. [34]. In this version, the parameters are set to be κ = 0.538, σ = 0.1682 GeV 2 and m c = 1.44 GeV in order to reproduce the mass difference of the J/ψ and ψ(2S), and the leptonic width of the J/ψ. Similarly, the bottom quark mass parameter is taken to be m b = 3.98 GeV in order to reproduce the mass difference of the Υ(1S) and Υ(2S), and the leptonic width of the Υ(3S) [35].
We have determined from these potential models the binding energy, ε 21 , and squared derivative of the radial wavefunction at the origin, |R  [30,32].
Furthermore, we have used models A, B, C, and D to compute the wavefunctions at the origin and binding energies of some nP bottomonium states. The results are listed in table 2. The values of the squared derivatives of the radial wavefunctions at the origin for the models A, B and C are taken from refs. [30,32]. For bottomonium we employ also a potential model E.
Model E is similar to model D. The only difference is that we take κ = 0.508 and m b = 4.68 GeV so to reproduce the mass difference of the Υ(2S) and Υ(3S), and the leptonic width of the Υ(3S) in the calculation of section 4.5.
In order to compute electromagnetic decay widths and production cross sections of P -wave quarkonia at relative order v 2 , we would need to include also order v 2 corrections to the wavefunctions at the origin. Since the model dependence of the values of |R that we employ exceeds v 2 , their inclusion would not improve, however, the accuracy of the model determinations. Hence, we account for the order v 2 corrections to the wavefunctions at the origin only in our final error budget. The sole corrections that we add explicitly are those that depend on the total angular momentum, since they contribute to distinguish between decay and production of χ QJ (nP ) states with different J. Nevertheless, also in this case the corrections are smaller than the systematic uncertainty due to the models. The total angular momentum corrections originate from one single 1/m 2 spin-orbit potential, V LS /m 2 . In the short range, this potential generates corrections to the wavefunction at the origin of relative order α 2 s that are divergent. The renormalization of these divergences requires introducing order α 2 s short distance coefficients for P -wave quarkonium production and decay processes. Since these are unknown at present, in this work we will not include corrections due to the short distance part of the spin-orbit potential. In the long range, the behavior of the spin-orbit potential is entirely described by the string tension σ and fixed by Lorentz symmetry [36,37]. This behavior is confirmed by lattice calculations [27]. The spin-orbit potential in the long range reads: where L and S are the total orbital angular momentum and spin, respectively. For P -wave spin-triplet states, the corrections to |R LS /m 2 have the following form where ∆ LS (nP ) depends on the radial excitation. We have listed the values of ∆ LS (nP ) in the different models for the charmonium 1P state in table 1, and for the bottomonium 1P , 2P , and 3P states in table 2. We see explicitly that the correction induced by the spin-orbit potential is smaller than the intrinsic potential model uncertainty, which we infer from the spread of the different wavefunction determinations. Finally, we observe JHEP04(2020)095  Table 3. For the five potential models described in the text, we list the squared radial wavefunctions at the origin, |R n0 (0)| 2 , and the binding energies, ε n0 , of the bottomonium 2S and 3S states.
that uncalculated corrections of relative order v 2 coming from the quantum-mechanical 1/m expansion of the quarkonium Fock state, in particular for the LDME of eq. (3.8), can be spin and angular momentum dependent as well. These uncalculated corrections are included in the error budget of the LDME, although the tuning of the potential model parameters may effectively reduce their size.
With the same five potential models described above we have also determined at leading order in v the squared radial wavefunctions at the origin and the binding energies of the 2S and 3S bottomonium states. The results are listed in table 3. The values of the wavefunctions at the origin for the models A, B, and C are taken from refs. [30,32].

P -wave charmonium electromagnetic decay and production
In this section, we compute the charmonium decay widths Γ(χ cJ (1P ) → γγ) and the cross sections σ(e + e − → χ cJ (1P )+γ) using the NRQCD factorization formulas (B.4) and (B.19), which are valid up to order v 2 , and rewriting the LDMEs according to the strongly coupled pNRQCD factorization formulas (3.8)-(3.10). We determine the gluonic correlators E 1 and iE 2 by fitting the available data.
The experimental inputs that we use are the χ c0 (1P ) and χ c2 (1P ) two photon decay widths and the cross section σ(e + e − → χ c1 (1P ) + γ). The BESIII measurements for the former give [38] Γ(χ c0 (1P ) → γγ) BESIII = 2.33 ± 0.20 ± 0.22 keV , For the latter, very recently Belle has observed the process e + e − → χ c1 (1P ) + γ and measured at √ s = 10.6 GeV [39] σ(e + e − → χ c1 (1P ) + γ) Belle = 17.3 +4.2 −3.9 ± 1.7 fb . From the theoretical side, rather than using the NRQCD factorization formulas for electromagnetic processes in their original form (see appendix B) we prefer using NRQCD factorization formulas at the amplitude level. So that the matching, the velocity expansion and the power counting are done for the amplitudes rather than for the decay widths or cross sections. In practice, one moves from the original factorization formulas to the ones JHEP04(2020)095 at the amplitude level through replacements of the type (4.6) The advantage is that in this way, without losing any systematicity, one is effectively including some potentially large contributions of order v 4 , α s v 2 and α 2 s in the expressions of the observables. Moreover, we replace the uncertain heavy quark pole mass with the spin average of the masses of the P states. This is defined in the case of charmonium 1P states as We take the 1P charmonium masses from ref. [21]. At our accuracy it is sufficient to use the following relation between M 1Pc and the charm pole mass m c : which is valid up to order v 2 . Afterwards we expand in powers of the binding energy up to relative order v 2 accuracy in the amplitude. Eventually, the theoretical expressions for the two photon decay widths that we use in the numerical analyses are where e c = 2/3. Note that Γ(χ c2 (1P ) → γγ) does not depend explicitly on the binding energy ε 21 . Similarly, the expression for the cross sections σ(e + e − → χ cJ (1P ) + γ) at the center of mass energy √ s that we use is  Table 4. Results of the fit of the gluonic correlators E 1 and iE 2 , when the wavefunctions and binding energies are computed within the potential models described in section 4.1.
In eqs. (4.9)-(4.11) the corrections proportional to α s come from the short distance coefficients. Hence α s should be understood as evaluated at a high energy scale and counting parametrically like a v 2 correction in the velocity expansion. For the decay widths, we use α s = 0.282, which is evaluated at the scale M 1Pc /2, and α = 1/137 reflecting the fact that the photons in the final state are on-shell. For the cross section at √ s = 10.6 GeV, we evaluate α s at the scale √ s/2 and take α s = 0.200. Of the fine structure constants appearing in eqs. (B.20)-(B.22), two originate from the virtual production mechanism, and are evaluated at virtuality √ s, and one originates from the real photon emission, and is evaluated at zero virtuality. For α at the scale √ s we take α = 1/131, while for α at zero virtuality we take α = 1/137. In the fit, we take the uncertainties in the decay rates and in the cross section to be 0.3 times the central values, for the order v 2 corrections that we have not included, and α 2 s times the central values, for the uncalculated corrections of higher orders in α s . We also add the experimental errors.
We determine E 1 and iE 2 by a least squares fit. The results for E 1 and iE 2 , when wavefunctions and binding energies are computed by means of the potential models described in section 4.1, are shown in table 4. The errors include the theoretical errors due to higher order corrections in v 2 and α s as described above, and the experimental errors in the data. Taking the averages over the different models, we get E 1 = −0.20 +0.14 −0.14 ± 0.90 GeV 2 , (4.12) where the first uncertainty comes from the potential model dependence and the second one is the average of the uncertainties in each potential model determination. The correlators E 1 and iE 2 have a size that is consistent, within uncertainties, with their naive scaling in powers of Λ QCD . The uncertainties are, however, large, reflecting the large uncertainty carried by the potential models. Vanishing small correlators are also consistent with our determinations.
One may wonder if it would not be possible to fit also the quarkonium wavefunctions eliminating in this way a major source of uncertainty. For the considered observables, see eqs. (4.9)-(4.11), it is not possible to disentangle the contribution of the wavefunction from the one of the correlator iE 2 . Hence, a fit would be able to determine a combination of the wavefunction and iE 2 , but not each of the two. Since in section 4.3 we aim at making some predictions for P -wave bottomonium electromagnetic decay and production, we have chosen to add the information on the wavefunction coming from potential models and gain some insight in the universal correlator iE 2 .
The correlators undergo renormalization and therefore depend on a subtraction scheme and a renormalization scale. The considered observables, at their present accuracy, are however insensitive to the renormalization of E 1 and iE 2 , and, in particular, they are insensitive to the renormalization scale of the correlators. We may reasonably expect that the obtained values refer to a renormalization scale that is of the order of the typical hadronic scale, but at this point further specifications are not possible. In section 4.4, we will see instead a case where the observable is sensitive to the renormalization of the involved correlator, so that a proper renormalization scale can be fixed, at least at leading logarithmic accuracy.
Combining the values of E 1 and iE 2 with the potential model calculations of the wavefunctions and binding energies in table 1, we obtain the χ cJ (1P ) LDMEs listed in table 5, whose errors are due to the uncertainties in the correlators E 1 and iE 2 . Here and in the following we take into account that the uncertainties in E 1 and iE 2 are correlated. The averages read where the first uncertainties come from the potential model dependence and the second ones from the uncertainties in E 1 and iE 2 . In χ cJ (1P )| T em 8 ( 3 P J ) |χ cJ (1P ) and χ cJ (1P )| P em 1 ( 3 P J ) |χ cJ (1P ) , we have ignored the total angular momentum dependent corrections to the wavefunction at the origin, because they are of higher order in v. The obtained value for the octet LDME is compatible with zero.

Potential model
The results for the two photon decay widths of the charmonium P -wave states χ c0 (1P ) and χ c2 (1P ) for each potential model determination of the wavefunction and binding energy are listed in table 6. The errors are due to the uncertainties in the correlators E 1 and iE 2 . The averages of these determinations read where the first uncertainty comes from the potential model dependence and the second one is the average of the uncertainties from each potential model. The determined values of E 1 and iE 2 allow us to make predictions for the cross sections σ(e + e − → χ cJ (1P ) + γ). In table 7, we list for each potential model the results at √ s = 10.6 GeV. The uncertainties in table 7 are computed from the uncertainties of E 1 and iE 2 , which already account for the uncertainties originating from the missing corrections of relative order v 2 and α 2 s , and from adding in quadrature the uncertainty that comes from varying α s between α s ( where the first uncertainty is from the model dependence and the second one is the average of the uncertainties in table 7. The obtained cross sections are consistent, inside errors, with the results of ref. [25]. It is worthwhile emphasizing that, although the measured two photon decay widths of the χ c0 (1P ) and χ c2 (1P ) states and the cross section σ(e + e − → χ c1 (1P ) + γ) have been used as an experimental input, the theoretical results for these quantities, eqs. ( Table 8. Results for the two photon decay widths of the states χ b0 (nP ) and χ b2 (nP ), indicated with Γ γγ χc0(nP ) and Γ γγ χc2(nP ) respectively, for each of the potential models described in section 4.1.
and (4.22), and their agreement with the data, eqs. (4.3)-(4.5), is nevertheless significant. The reason is that the two correlators E 1 and iE 2 are the result of a least squares fit of three data and not of a fine tuning of some of them.

P -wave bottomonium electromagnetic decay and production
With the values of E 1 and iE 2 determined in the previous section we can make predictions for P -wave electromagnetic decay widths and production cross sections of bottomonium states. Wavefunctions and binding energies are computed according to the potential model results listed in table 2. We consider, first, the two photon decay rates of the states χ b0 (nP ) and χ b2 (nP ) with n = 1, 2 and 3. Following the same procedure discussed in section 4.2 for charmonium states, we replace in our theoretical expressions for the χ bJ (nP ) decay widths the bottom pole mass with the spin average of the nP bottomonium masses. We use the equivalent of eq. (4.7) and eq. (4.8). The masses of the 1P , 2P and 3P bottomonium states are taken from ref. [21]. 5 Moreover, we take α = 1/137 reflecting the fact that the photons in the final state are on shell, and we take α s = 0.200 at the scale of half the spin averaged masses. The results for each choice of potential model used to compute wavefunctions and binding energies are shown in table 8. The uncertainties in table 8 come from the correlated uncertainties in E 1 and iE 2 , as well as from the uncertainties stemming from uncalculated corrections of order v 2 and α 2 s in the bottomonium sector, which we estimate to be 0.1 and α 2 s times the central values, respectively. The uncertainties are added in quadrature. After averaging over the determinations from the different potential models, we obtain the following predictions  Γ(χ b0 (3P ) → γγ) = 46.1 +11.9 −11.8 ± 6.3 eV, (4.28) where the first uncertainty comes from the potential model dependence, and the second one is the average of the uncertainties in table 8. Using the same input as for the two photon decay widths, we can also make predictions for the cross sections σ(e + e − → χ bJ (nP ) + γ). As has been pointed out in this context in ref. [41] and mentioned at the end of appendix B, the perturbative expression of the electromagnetic cross section becomes singular when the center of mass energy approaches the heavy quark-antiquark pair production threshold. In the bottomonium case, this threshold is around 10 GeV. Therefore, in order to make predictions for σ(e + e − → χ bJ (nP ) + γ) using the factorization formulas provided in appendix B, the center of mass energy has to JHEP04(2020)095 be significantly larger than 10 GeV. We look at the energy range 20 GeV< √ s <500 GeV that encompasses the energies of a possible future e + e − collider. We evaluate α s at the scale √ s/2. Furthermore, we fix α = 1/128 and neglect the running, as the running of α only affects the cross section by less than a few percent, which is negligible compared to other uncertainties. The results for the electromagnetic cross sections of the 1P , 2P , and 3P bottomonium states are shown in figures 2, 3 and 4, respectively. The central values are obtained by averaging over the determinations of the wavefunctions and binding energies from the different potential models described in section 4.1. The bands account for the uncertainties, which include potential model dependence, uncertainties in E 1 and iE 2 , uncertainties from uncalculated corrections of order v 2 and α 2 s , which we estimate to be 0.1 and α 2 s times the central values, and uncertainties coming from varying α s between α s ( √ s/4) and α s ( √ s). We add these uncertainties in quadrature. In particular, at √ s = 20 GeV the cross sections are σ(e + e − → χ b0 (1P ) + γ) = (2.47 ± 0.83 ± 0.56) × 10 −3 fb , where the first uncertainties come from the dependence on the potential models, and the second ones account for the other uncertainties that we have mentioned above eq. (4.30).

P -wave charmonium and bottomonium decay into light hadrons
In this section, we analyze inclusive P -wave quarkonium decays into light hadrons (LH) at leading order in the velocity expansion. The NRQCD factorization formula has been first derived in [42] and we reproduce it in eq. (B.14). It depends on two LDMEs. The color singlet LDME has been factorized in strongly coupled pNRQCD at leading order in v in eq. (3.7). The color octet matrix element can be written in strongly coupled pNRQCD at leading order in v as [13] where the gluonic correlator E 3 has been defined in eq. (3.4). The expression of the decay width of a P -wave quarkonium into light hadrons under the conditions of validity of strongly coupled pNRQCD reads therefore at leading order in v: In eq. (4.49) we have emphasized that both Imf 1 ( 3 P J )(µ Λ ) and E 3 (µ Λ ) depend on a cutoff µ Λ . We use the MS scheme for both quantities. The correlator E 3 (µ Λ ) is dimensionless, and depends logarithmically on the scale µ Λ . This dependence cancels in the decay width (4.49) against the µ Λ dependence of the short distance coefficient Imf 1 ( 3 P J )(µ Λ ) [13]. Also the one loop running with respect to the scale µ Λ of Imf 1 ( 3 P J )(µ Λ ) and E 3 (µ Λ ) is known.
theoretical expressions for the decay rates that we use are valid up to next-to-leading order in α s , except for Γ(χ c1 (1P ) → LH), which is known only at leading order in α s . The decay rates Γ(χ cJ (1P ) → LH) have been obtained by subtracting radiative decay rates and transition rates into other charmonia from the total widths of χ cJ (1P ) given in ref. [21]. This result is compatible, within errors, with a previous determination in ref. [13]. Nevertheless, we note that the uncertainties are smaller, despite the determination in [13] did not include theoretical uncertainties. 6 From eq. (4.53), we can compute E 3 (µ Λ ) at different scales by using the one loop renormalization group improved expression [13]: where

JHEP04(2020)095
where the first uncertainty comes from the potential model dependence, and the second one is the average of the uncertainties in table 9. From eq. (4.53) and eq. (4.49), we can compute the decay widths of the χ cJ (1P ) states into light hadrons. Rather than computing the decay rates directly using eq. (4.49), we determine Γ(χ cJ (1P ) → LH) for J = 0 and 2 states by combining the ratios Γ(χ cJ (1P ) → LH)/ Γ(χ cJ (1P ) → γγ) with the two photon decay widths computed in eq. (4.19) and eq. (4.20). Similarly, we compute Γ(χ c1 (1P ) → LH) by combining these determinations of Γ(χ c0 (1P ) → LH) and Γ(χ c2 (1P ) → LH) with the ratios Γ(χ c0 (1P ) → LH)/Γ(χ c1 (1P ) → LH) and Γ(χ c1 (1P ) → LH)/Γ(χ c2 (1P ) → LH). This approach has the advantage that the ratios do not depend on the choice of potential model, a fact that reduces significantly the uncertainties. Our results for the ratios Γ(χ cJ (1P ) → LH)/ Γ(χ cJ (1P ) → γγ) for J = 0 and 2 are where the uncertainties come from uncalculated corrections of relative order v 2 and α 2 s , which are taken to be 0.3 and α 2 s times the central values, respectively. These uncertainties are added in quadrature. Since the uncertainty in E 3 is dominated by uncertainties from uncalculated higher order corrections to the theoretical expressions of the ratios, we do not include the uncertainty in E 3 to avoid double counting. Using eqs.

Υ(2S) and Υ(3S) decay into lepton pairs
The NRQCD factorization formula for the decay width of a vector S-wave quarkonium state into a lepton pair at relative order v 2 is given by eq. (B.1). It depends on two LDMEs: V Q (nS)|O em 1 ( 3 S 1 )|V Q (nS) , whose factorization in strongly coupled pNRQCD at relative order (Λ QCD /m) 2 is in (3.13), and V Q (nS)|P em 1 ( 3 S 1 )|V Q (nS) , whose factorization in strongly coupled pNRQCD at leading order is in (3.15).
At relative order (Λ QCD /m) 2 the matrix element V Q (nS)|O em 1 ( 3 S 1 )|V Q (nS) depends, besides on E 1 and E 3 , also on a correlator involving four chromoelectric fields and a correlator involving two chromomagnetic fields. In this section, also to avoid dealing with correlators about which practically nothing is known, we will explore the leptonic decays of the bottomonium states Υ(2S) and Υ(3S) assuming that these states satisfy the kinematical condition Under this condition the matrix element Υ(nS)|O em 1 ( 3 S 1 )|Υ(nS) for n = 2, 3 can be written in strongly coupled pNRQCD at JHEP04(2020)095  Table 10. Results for the leptonic decay widths of the states Υ(2S) and Υ(3S), indicated with Γ e + e − Υ(2S) and Γ e + e − Υ(3S) for short. Wavefunctions at the origin and binding energies have been computed within the potential models of section 4.1.
relative order v 2 as as all other contributions in eq. (3.13) are of relative order (Λ QCD /m) 2 and, therefore, suppressed.
Under the assumption that the states Υ(2S) and Υ(3S) satisfy the condition m b v Λ QCD m b v 2 , their decay width into a lepton pair can be written up to relative order v 2 in strongly coupled pNRQCD as n0 , which is valid up to order v 2 , and expanded in the leading order binding energy, ε (0) n0 , up to relative order v 2 . Hence, under the above assumptions, the decay widths Γ(Υ(nS) → e + e − ) for n = 2, 3 depend on the Υ(nS) wavefunctions at the origin, the binding energies and the chromoelectric correlator E 3 . The wavefunctions at the origin and the binding energies have been computed in the potential models of section 4.1 and are listed in table 3. The correlator E 3 has been computed in the previous section from the decays of P -wave charmonia and its value at 1 GeV is given in eq. (4.53). 7 The next-to-leading order expression of Im gee( 3 S1) contributes at relative order αsv 2 , which is beyond our accuracy. It is worth noting, however, that the next-to-leading order expression of Im gee( 3 S1) depends on the cutoff µΛ and that this dependence cancels against the µΛ dependence of E3 in the expression of the dilepton decay width of S-wave quarkonia:
We take α = 1/131 and compute α s at the scale of the meson mass, which gives Under the same kinematical conditions considered above, we can also compute the inclusive decay widths of the Υ(2S) and Υ(3S) states into light hadrons. The NRQCD factorization formula valid at relative order v 2 is eq. (B.11). It depends on the LDMEs: Υ(nS)|O 1 ( 3 S 1 )|Υ(nS) and Υ(nS)|P 1 ( 3 S 1 )|Υ(nS) for n = 2, 3. At relative order v 2 and under the condition m b v Λ QCD m b v 2 , from the comparison of eq. (3.11) with eq. (3.13) it follows that Υ(nS)|O 1 ( 3 S 1 )|Υ(nS) is equal to Υ(nS)|O em 1 ( 3 S 1 )|Υ(nS) . Moreover, Υ(nS)|P 1 ( 3 S 1 )|Υ(nS) is equal to Υ(nS)|P em 1 ( 3 S 1 )|Υ(nS) at leading order in the velocity and Λ QCD /m expansion. Hence, by using the same strongly coupled pNRQCD factorization formulas for LDMEs employed above, we can write at relative order v 2 and neglecting corrections of order Λ 2 QCD /m 2 The expressions for the short distance coefficients are in eq. (B.12) and eq. (B.13), and we have again expressed the bottom mass in terms of the Υ(nS) mass and the corresponding binding energy up to relative order v 2 . 8 If we consider the ratios of the leptonic decay widths, eq. (4.79), with the corresponding decay widths into light hadrons, eq. (4.88), the wavefunction dependence drops out and 8 Using the same reasoning of footnote 7, from requiring that the right-hand side of eq. (4.88) is independent of the factorization scale µΛ, it follows that Img1( 3 S1) must develop a µΛ dependence at order α 4 s that exactly cancels the one in E3. The µΛ dependent part of Img1( 3 S1) at order α 4 s then reads JHEP04(2020)095  the ratio depends only on the binding energies, whose potential model dependent values are in table 3, and on the chromoelectric field correlator E 3 , which has been determined in eq. (4.53) and whose running at leading logarithmic accuracy is described by eq. (4.54). Furthermore, if we expand the ratio in powers of v 2 , the dependence on E 3 also drops out, and we obtain an expression that depends only on the binding energy. To relative order v 2 accuracy, it reads we get an expression that is valid up to relative order v 2 and whose order v 2 correction is better under control. The results for this ratio are listed in table 12  the ratios Γ(Υ(nS) → e + e − )/Γ(Υ(nS) → LH). In order to compare the result in eq. (4.91) with measurements, we compute Γ(Υ(2S) → LH) and Γ(Υ(3S) → LH) by subtracting the radiative decay widths and the transition widths into other bottomonia from the total Υ(2S) and Υ(3S) widths given in ref. [21]; also the decay widths into e + e − are taken from ref. [21]. We obtain the following experimental value Compared to eq. (4.91), the theoretical result is compatible with the experimental value within uncertainties. The result (4.91) improves a more qualitative determination that can be found in ref. [14].
Other quarkonium S-wave vector states are the ψ(2S), the Υ(1S) and the Υ(4S). In the present analysis, we have not considered the states ψ(2S) and Υ(4S), because they are very close or above the open flavor threshold, respectively. Effects due to degrees of freedom that are relevant above or close to the open flavor threshold have not been included in the pNRQCD Hamiltonian (2.17). Hence, the effective field theory as formulated in this work is not suited to treat quarkonia like the ψ(2S) and Υ(4S). In this analysis, we have not considered the Υ(1S) too, because the hierarchy mv 2 Λ QCD is unlikely to be realized for this state, which is commonly treated assuming mv 2 Λ QCD (see, for instance, refs. [1,2,46]). Hence, the only S-wave quarkonia that satisfy possibly the condition mv Λ QCD mv 2 within the effective field theory (2.17) are the Υ(2S) and Υ(3S). If they really realize this kinematical condition may be eventually established only at the hand of phenomenological analyses of the kind presented here.

Summary and conclusion
In the paper, we have computed decay widths and exclusive electromagnetic production cross sections of charmonia and bottomonia based on strongly coupled pNRQCD. In strongly coupled pNRQCD, nonperturbative LDMEs are expressed in terms of quarkonium wavefunctions, binding energies and gluonic correlators. Wavefunctions and binding energies are the solutions of the equation of motion of pNRQCD. The gluonic correlators are nonperturbative parameters, which are independent of the quarkonium state and of the flavor of the heavy quark. Owing to the universal nature of the gluonic correlators, the number of nonperturbative unknowns needed in pNRQCD to describe decay widths and exclusive electromagnetic production cross sections of charmonia and bottomonia is smaller than the number of LDMEs needed in NRQCD, as they depend on the quarkonium state and on the flavor of the heavy quark. This enables specific predictions in pNRQCD that are not possible in NRQCD. Since strongly coupled pNRQCD is suited to describe non Coulombic quarkonium, we have restricted our applications to charmonium 1P states and bottomonium 2S, 3S, 1P , 2P and 3P states, which are possibly non Coulombic bound states.
The calculation of the NRQCD LDMEs in strongly coupled pNRQCD was first done in refs. [13,14]. We have computed new corrections to P -wave LDMEs and revised the computation of S-wave LDMEs by adding some new contributions proportional to Λ 2 QCD /m 2 . The newly computed corrections are expressed in terms of gluonic correlators. Our results

JHEP04(2020)095
for P -wave and S-wave LDMEs in the strongly coupled pNRQCD factorization are listed in sections 3.1.1 and 3.1.2, respectively. The LDMEs that we have computed satisfy the Gremm-Kapustin relations, as discussed in section 3.2.
We have applied strongly coupled pNRQCD for the first time to exclusive electromagnetic production cross sections. In particular, we have computed the cross sections e + e − → χ cJ (1P ) + γ in section 4.2 and e + e − → χ bJ (nP ) + γ for n = 1, 2 and 3 in section 4.3. Although straightforward in the case of exclusive electromagnetic production, the application of strongly coupled pNRQCD has never been attempted before for any production process. This has enabled us to make first and so far unique predictions for exclusive electromagnetic production cross sections of bottomonia. Furthermore, in section 4.2 we have computed the decay widths χ c0 (1P ) → γγ and χ c2 (1P ) → γγ, in section 4.3 the decay widths χ b0 (nP ) → γγ and χ b2 (nP ) → γγ for n = 1, 2 and 3 and in section 4.5, under some assumptions, the dilepton decay widths Υ(2S) → e + e − and Υ(3S) → e + e − . We have also considered hadronic annihilations. In From the theoretical side, our expressions are generally accurate up to relative order v 2 in the velocity expansion, with the exception of the P -wave inclusive annihilation widths into light hadrons, where we have truncated our expansions in v at leading order. The results that we obtain are in agreement, within errors, with experimental data, when available. The determinations of the e + e − → χ c0 (1P ) + γ, e + e − → χ c2 (1P ) + γ, e + e − → χ bJ (nP ) + γ cross sections, and of the χ b0 (nP ) → γγ, χ b2 (nP ) → γγ and χ b (nP ) → LH decay widths for n = 1, 2 and 3 are predictions. These predictions were made possible by the universal nature of the potential and gluonic correlators that determine the LDMEs in strongly coupled pNRQCD. The gluonic correlators can be computed, in principle, in lattice QCD. However, since lattice QCD determinations of the gluonic correlators are not available yet, we have determined them from the available data on decay and production of charmonia and used to compute bottomonium observables. This procedure should be contrasted with NRQCD, where one cannot, in general, infer the bottomonium LDMEs from the charmonium ones.
Our results rely on potential models for determining the quarkonium wavefunctions at the origin and the binding energies. We rely on potential models because a first principle determination of the wavefunctions and binding energies from the equation of motion of pNRQCD is hindered by the poor or incomplete knowledge of the corrections to the wavefunction at higher orders in v stemming from 1/m suppressed terms in the potential. Difficulties come from the limited accuracy in our knowledge of the potential beyond the static term, and also from the renormalization of the divergences in the wavefunction at the origin due to the potential at short distances. Concerning the 1/m corrections to the potential, we remark that, although these corrections can be expressed in terms of Wilson loops and gauge field strength insertions on them, not all of them have been computed in lattice JHEP04(2020)095 QCD and with enough precision. Indeed, besides the static potential, no 1/m suppressed correction to the quarkonium potential has been computed in full (unquenched) QCD.
The lack of a reliable determination of the quarkonium wavefunctions at the origin and the binding energies reflects in the wide spread of potential model results presented and discussed in section 4.1. As the quarkonium wavefunction enters most quarkonium observables at leading order, the poor knowledge of the quarkonium wavefunction is the main source of uncertainty for most observables computed in NRQCD, and, in particular, for the observables computed in this work that are not ratios of suitably chosen decay widths. Parametrically, the uncertainty in the wavefunction affects these observables at least at relative order v 2 assuming a perfect knowledge of the leading order potential. It can be argued that this is the case if the leading order potential coincides with the static potential; the relative order v 2 uncertainty stems, then, from the 1/m suppressed terms in the potential that have not been included neither in an accurate nor in a complete form. The uncertainty may be larger, however, if the term V (1) /m contributes at leading order to the potential (see the discussion in section 2.2).
The poor knowledge of the quarkonium wavefunctions at the origin and the binding energies is the main limitation in the phenomenological analyses done in the present (and similar) works. It could be overcome by improving our knowledge of the quarkonium potential, ideally via lattice QCD, but also by using all available short distance and long distance information on the potential in a comprehensive analysis. Wavefunctions at the origin and binding energies could be determined, in principle, also from a global fit of quarkonium observables versus data. The fit would then determine these parameters together with the field strength correlators encoding the universal non perturbative parts of the LDMEs. As we discussed in section 4.2, for the set of observables considered in the present work, it is not possible, even in principle, to disentangle the wavefunction at the origin from all field strength correlators. It remains an open and interesting question to answer if an enlarged set of observables may be able to solve the problem and fix on the data all non perturbative parameters.
Finally, we note that the strategy used in this work to compute the nonperturbative LDMEs could be possibly applied to study inclusive hadroproduction of heavy quarkonia too. This may improve our understanding of the inclusive production mechanism of heavy quarkonium that remains elusive to this day. In particular, the pNRQCD calculations of color-octet LDMEs could lead to a reduction of the number of nonperturbative unknowns.

A Four-fermion operators
In the following, we list the four-fermion operators in the NRQCD Lagrangian relevant for the present analysis. Extensive lists of four-fermion operators can be found in [47][48][49]. The specified quantum numbers identify the state on which the operator projects dominantly. The electromagnetic (em) operators are where the fields are defined as in section 2.1, |vac vac| projects on the QCD vacuum state, T (ij) ≡ (T ij + T ji )/2 − T kk δ ij /3, and H.c. stands for Hermitian conjugate. The subscript "1" labels four-fermion operators that project dominantly on a color singlet heavy quark-antiquark pair, whereas the subscript "8" labels four-fermion operators that project dominantly on a color octet heavy quark-antiquark pair. The relevant four-fermion hadronic operators are JHEP04(2020)095

B NRQCD factorization formulas
In this appendix, we list the NRQCD factorization formulas for quarkonium decay widths and exclusive electromagnetic production cross sections for which we provide the strongly coupled pNRQCD factorization formulas in the main body of the paper. We list first the decay formulas and conclude with the electromagnetic production ones. Leptonic decay widths of S-wave quarkonium vector states are described in NRQCD up to relative order v 2 by two LDMEs [3]: Γ(V Q (nS) → e + e − ) = 2 Im f ee ( 3 S 1 ) m 2 V Q (nS)|O em 1 ( 3 S 1 )|V Q (nS) + 2 Im g ee ( 3 S 1 ) m 4 V Q (nS)|P em 1 ( 3 S 1 )|V Q (nS) . (B.1) The short distance coefficients, which are the imaginary parts of the coefficients of the corresponding four fermion operators in the NRQCD Lagrangian, read at next-to-leading order [50][51][52][53] (the next-to-next-to-leading order correction to Im f ee ( 3 S 1 ) has been computed in refs. [54,55], and the next-to-next-to-next-to leading order correction to Im f ee ( 3 S 1 ) has been computed in ref. [56]) where e Q is the fractional electric charge of a heavy quark of flavor Q, α is the fine structure constant, α s is the strong coupling in the MS scheme and C F = (N 2 c − 1)/(2N c ) (= 4/3

JHEP04(2020)095
in QCD) is the Casimir of the fundamental representation of SU(N c ). The infrared cutoff µ Λ arises from renormalization of the short distance coefficients in NRQCD. In here and in the following, the short distance coefficients are renormalized in the MS scheme. Two photon decay widths of spin one and J = 0, 2 P -wave quarkonia are described in NRQCD up to relative order v 2 by the factorization formulas [23,47]: The coefficients Im f em ( 3 P J ) for J = 0, 2 have been computed at order α s in refs. [57,58], and the coefficients Im g em ( 3 P J ) and Im t 8 em ( 3 P J ) for J = 0, 2 have been computed at leading order in ref. [47]. Inclusive decay widths into light hadrons (LH) of S-wave quarkonium vector states are described in NRQCD up to relative order v 2 by the factorization formula [14]: Γ(V Q (nS) → LH) = 2 Im f 1 ( 3 S 1 ) m 2 V Q (nS)|O 1 ( 3 S 1 )|V Q (nS) (B.11) The approximation in the last line holds when mv Λ QCD mv 2 and when neglecting, at relative order v 2 , all contributions that scale like (Λ QCD /m) 2 (one should consider in the power counting of strongly coupled pNRQCD the expressions of the LDMEs in eqs. (3.11), (3.15), (3.16), (3.17) and the expression of V Q (nS)|O 8 ( 3 S 1 )|V Q (nS) given in JHEP04(2020)095 C Derivation of some LDMEs in strongly coupled pNRQCD We give here few more details on the computation of the electromagnetic P -wave LDMEs appearing in section 3.1.1.

C.1
χ QJ (nP )|O em 1 ( 3 P J )|χ QJ (nP ) We start by computing the matrix element (0) n; where we have used eq. (2.7) for the first equality, Wick's theorem for the second one and |0 For n = k = 0 and from eqs. (2.9), it follows For n = 0, k = 0 and from eqs. (2.11), it follows where we have also used the fact that (0) n|gE 1 |0 (0) and (0) n|gE T 2 |0 (0) are equal at r = 0. Finally, for completeness we give the result for n = 0 and k = 0, which reads (0) n; x 1 , We also need the correction to the matrix element (C.5) stemming from next-to-leading order in the quantum-mechanical expansion in 1/m. In particular, we need the part that projects on P waves. It is given by where we have used eq. (2.12), eq. (C.6) and, in the last equality, the definition (

C.2
χ QJ (nP )|T em 8 ( 3 P J )|χ QJ (nP ) We compute here the matrix element 0; x 1 , x 2 | d 3 x T em 8 ( 3 P J )(x)|0; x 1 , x 2 . Differently from the matrix elements computed in the previous and in the next section, at leading order in the 1/m expansion this matrix element does not contribute to P -wave quarkonium states. Hence, the first non vanishing contribution of the matrix element of the color octet operator T em 8 ( 3 P J ) on a P -wave state appears only at next-to-leading order in 1/m. Indeed, at leading order the matrix element reads where we have used eqs. (2.10). The function δ (3) (r) brings the derivative of the static energy, ∇ r E 0 (x 1 , x 2 ), into the perturbative regime. At r = 0 this vanishes in dimensional JHEP04(2020)095 regularization, because it is scaleless. Moreover, eq. (C.10) defines through (3.1) a contact term V (5) T em 8 ( 3 P J ) with only one derivative left to act on the wavefunctions. This is not enough in eq. (3.2) to produce a non vanishing contribution for P -wave states. At least two derivatives are necessary.
At next-to-leading order in the quantum-mechanical 1/m expansion, the missing derivative comes from the 1/m correction to the state shown in eq. (2.12): Plugging eq. (C.11) into eq. (3.1) one gets a contact term V 2) the expression of the octet LDME, χ QJ (nP )|T em 8 ( 3 P J ) |χ QJ (nP ) , given in eq. (3.9) follows.

C.3
χ QJ (nP )|P em 1 ( 3 P J ) |χ QJ (nP ) Finally, we compute the matrix element 0; x 1 , x 2 | d 3 x P em 1 ( 3 P J )(x)|0; x 1 , x 2 . For our purposes, it is sufficient to consider it at leading order in the quantum-mechanical 1/m expansion: The corresponding contact term from eq. (3.1) reads (C.14) From V (6) P em 1 ( 3 P J ) and eq. (3.2) the expression of the LDME, χ QJ (nP )|P em 1 ( 3 P J ) |χ QJ (nP ) , written in eq. (3.10) follows if we replace the Laplacian acting on the wavefunction at the origin with the expression given in eq. (3.6). Note that the two derivatives in (C.14), ∇ i r and ∇ j r , need to act on the wavefunctions to ensure a non vanishing contribution for P -wave states to the right-hand side of eq. (3.2).
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.