Relativistic corrections to exclusive χ cJ + γ production from e + e − annihilation

from e+ e− annihilation Nora Brambilla, Wen Chen, Yu Jia, Vladyslav Shtabovenko, and Antonio Vairo1,∥ Physik-Department, Technische Universität München, James-Franck-Strasse 1, 85748 Garching, Germany Institute for Advanced Study, Technische Universität München, Lichtenbergstrasse 2a, 85748 Garching, Germany Institute of High Energy Physics and Theoretical Physics Center for Science Facilities, Chinese Academy of Sciences, Beijing 100049, China School of Physics, University of Chinese Academy of Science, Beijing 100049, China Center for High Energy Physics, Peking University, Beijing 100871, China Zhejiang Institute of Modern Physics, Department of Physics, Zhejiang University, Hangzhou 310027, China

I.

INTRODUCTION
More than four decades have passed since the experimental groups of Samuel Ting and Burton Richter [1, 2] discovered the J/ψ, the first observed bound state formed by a heavy (charm) quark and a heavy antiquark. This event was of crucial importance for the establishment of Quantum Chromodynamics (QCD) as the theory of the strong interactions. Even though heavy quarkonia are sometimes regarded as the "hydrogen atoms" of QCD, our present understanding of these hadronic systems is not as successful as that of the hydrogen atom in Quantum Electrodynamics (QED). The nonperturbative nature of QCD at low energy makes the theoretical description of heavy quarkonium on the one hand more challenging -and on the other hand also more interesting [3][4][5].
In particular, the theoretical description of the production mechanism for charmonia and bottomonia is complicated. While the creation of a heavy quark pair in a high energy collision can be described in perturbative QCD, this is not the case for the evolution of the pair into a heavy quarkonium, which is governed by nonperturbative long-distance effects.
Effective field theories (EFTs) provide an elegant way to treat this problem by exploiting the nonrelativistic nature of the system and the separation of the relevant scales [6]. The EFT resulting from integrating out modes associated with the energy scale of the heavy-quark mass or larger is known as nonrelativistic QCD (NRQCD) [7]. It conjectures a factorization theorem allowing one to write the quarkonium production cross section as an expansion where each term is the product of a short-distance coefficient and a long-distance matrix element (LDME). The former is computed from matching to perturbative QCD, it is a series in the strong coupling α s , and it incorporates effects from the hard scale m (heavy quark mass) and above. The latter is of nonperturbative nature and can be obtained from fits to experimental data. The LDMEs are assumed to be universal and to obey power counting rules in the relative heavy quark velocity v, which in a nonrelativistic system is much smaller than one. Hence, the importance of each term in the NRQCD-factorized cross section can be estimated by powers of α s and v. This counting gives NRQCD predictive power and allows for a systematic inclusion of radiative and relativistic corrections.
NRQCD also predicts that quarkonium formation is not limited to heavy quark pairs in a color singlet configuration, as it was assumed by the color singlet model [8][9][10]. Instead, the Fock state of a heavy quarkonium can be schematically understood as a sum over different contributions. For example, for the P -wave charmonium χ cJ we have where higher Fock states that involve gluons are suppressed by powers of v. Nevertheless, in higher order calculations such contributions become relevant and must be included to obtain a consistent result. Electromagnetic quarkonium production can be regarded as a relatively simple and clean process that makes it a perfect testing ground for verifying predictions made by NRQCD. Electromagnetic quarkonium production is an important subject of study at the BES-III experiment at the tau-charm factory BEPC-II in China, and at the Belle experiment at the KEKB asymmetric B-factory in Japan. For the latter, one particularly interesting electromagnetic quarkonium production process is the exclusive production of a heavy quarkonium and a hard (|k| m) photon in electron-positron annihilations.
The scope of this work is to improve on the NRQCD calculation of the process e + e − → γ * → χ cJ + γ by including higher Fock state contributions of O(α 0 s v 2 ). The leading cross section of O(α 0 s v 0 ) was obtained in [11]. Subsequently, corrections of O(α s v 0 ) [12,13], partial corrections of O(α 0 s v 2 ) [14,15] and finally partial corrections of O(α s v 2 ) [16] were computed as well. One of the reasons why this process has attracted so much interest from the theory side in the past years is the anticipated connection between C-even quarkonia and some exotic XY Z particles [4,[17][18][19], although such identifications are still controversial [20,21]. In any case, there is clear experimental evidence [22] that the C-even state X(3872) can be measured in the same channel [e + e − → X(3872) + γ] as our process of interest. Therefore, irrespective of the true nature of X(3872), in order to study this state in more detail, it is very useful to have precise and unambiguous predictions for e + e − → γ * → χ cJ + γ. This is all the more important in view of the fact that no experimental measurements for the electromagnetic production of χ cJ and a hard photon are yet available, while good perspectives for this measurement exist at Belle II.
We emphasize that the above-mentioned NRQCD studies considered only operators that contribute through the dominant Fock state |QQ . However, operators that project on the subleading Fock state |QQg already show up at O(α 0 s v 2 ). The importance of this kind of operators was already noticed in the past. In particular, they were explicitly incorporated in studies of the decay χ cJ → γγ [23,24], which is a process very similar to γ * → χ cJ + γ.
In the present work, we will investigate these missing contributions to the production cross section by determining at O(α 0 s v 2 ) in NRQCD the matching coefficients of octet operators containing chromoelectric fields.
The paper is organized in the following way. In Sec. II we discuss the relevant operators and LDMEs, and write the NRQCD-factorized cross sections for e + e − → χ cJ + γ at the desired precision. Different strategies for the determination of the NRQCD matching coefficients are discussed in Sec. III. In Sec. IV we outline the calculation of the NRQCD amplitudes, while in Sec. V we describe the QCD side of the matching and present the short-distance coefficients obtained with two independent methods. A long-standing discrepancy between the results of [23] and [24] regarding some matching coefficients that enter the NRQCD-factorized decay rates for χ cJ → γγ at O(α 0 s v 2 ) is resolved in Sec. VI. Sections VII and VIII contain the final cross sections and numerical predictions. Finally, a summary of the obtained results and their impact is presented in Sec. IX. In Appendix A we present additional short-distance coefficients that arise from the matching. They are not relevant for our calculation but they are new and appear at O(v 4 ) contributing to the dominant Fock state. Appendix B contains the details of the calculation made with the threshold expansion method of Braaten and Chen [25]. Appendix C contains useful formulas about polarization sums of hadrons. Appendix D provides a derivation of the generalized Gremm-Kapustin relations for 3 P J quarkonia.
A generic NRQCD operator of naive scaling dimension d n (which is given by the dimension of the gluonic and fermionic fields only) that describes the production of a heavy quarkonium H can be written as [7] O n = χ † K n ψ X m J |H + X H + X| ψ † K n χ, where K ( ) n consists of a spin matrix (1 or σ), a color matrix (1 or T a ) and a polynomial in B i ≡ ijk G kj /2, E i ≡ G i0 , and D = ∇ − igA with A µ being the gluon field, G µν the field-strength tensor of QCD and g the gauge coupling 1 . Here and throughout the paper, bold fonts specify Cartesian 3-vectors. Furthermore, ψ (χ) denotes a Pauli field that annihilates (creates) a heavy quark (antiquark), while |H + X is a Fock state that contains a quarkonium H and all other particles X that appear in the final state, and have energy and momenta that are typically much smaller than the hard scale m. The summations run over all allowed X states and the magnetic quantum number m J = −J, . . . , J, where J denotes the spin of the quarkonium state. The operatorsÕ n enter the NRQCD factorization formula as vacuum expectation values where |0 denotes the QCD vacuum and the quarkonium states have nonrelativistic normalization, i. e., with P ( ) being the total 3-momentum of the heavy quarkonium. The sum over X in Eq. (3) is one of the reasons why production matrix elements 0|Õ n |0 cannot easily be calculated in lattice simulations [26]. However, in exclusive electromagnetic production this sum is absent, such that the corresponding LDME reads Equation (5) can be further simplified using the rotational invariance of the matrix elements [7,25]. Since all matrix elements that differ only in the quantum number m J are identical, the sum over m J trivially reduces to (2J + 1) and the LDME factorizes into a product of two quarkonium-to-vacuum matrix elements The matrix elements 0| χ † K n ψ |H and H| ψ † K n χ |0 are the same that enter NRQCD factorization formulas for exclusive electromagnetic decays [27]. Hence, up to the prefactor (2J + 1), the LDME for exclusive electromagnetic production 0|O n |0 is identical to the corresponding LDME for electromagnetic decay H|O n em |H with O n em = ψ † K n χ |0 0| χ † K n ψ: 0|O n |0 = (2J + 1) H|O n em |H .
For this reason from now on we will absorb the factor (2J + 1) into the matching coefficients F n : 1 In the case of polarized production K ( ) n may also depend on polarization vectors [7]. Explicit NRQCD calculations of such processes can be found, e. g., in [25]. In this work, we will restrict ourselves only to unpolarized production.

B. Color singlet and color octet production operators
In NRQCD it is customary to classify production operators according to the color configuration of the QQ pair in the leading Fock state through which these operators contribute to the process. In the dominant Fock state |QQ the heavy quark pair is in a color singlet configuration; hence operators that contribute dominantly through this state are denoted as singlet operators. An octet operator is an operator that gives a nonvanishing contribution only when the QQ pair is in a color octet configuration, such as in the subleading Fock state |QQg . An example for a color octet operator is which contributes to the inclusive production of a spin-triplet P -wave quarkonium at leading order in v. In exclusive electromagnetic production there are no color octet operators similar to Eq. (9). This is because a matrix element like 0| χ † σT a ψ |H vanishes due to color conservation. However, if we replace T a with the chromoelectric field multiplied by the gauge coupling, we obtain 0| χ † (gE · σ)ψ |H , which gives a nonvanishing contribution to the exclusive electromagnetic production of a P -wave quarkonium at subleading order in v. The corresponding operator reads where H.c. stands for Hermitian conjugate. According to the usual NRQCD classification, the operator in Eq. (11) is a color octet operator in the sense that this operator necessarily projects on the subleading Fock state |QQg . 2 To our knowledge, the number of studies that considered contributions from operators with chromoelectric or chromomagnetic fields to NRQCD-factorized cross sections and decay rates is very small. One of such studies was concerned with the electromagnetic and hadronic decays of the S-wave heavy quarkonia [28], where LDMEs containing one chromoelectric field were of relative order v 4 as compared to the leading order LDMEs. The authors, however, managed to eliminate such LDMEs from their basis by using the Gremm-Kapustin relations [29]. An explicit determination of the matching coefficients multiplying E-field dependent LDMEs was carried out in [23], where the authors considered electromagnetic decays of the P -wave spin triplet quarkonia. There, the relevant LDMEs contributed already at relative order v 2 . Matching coefficients of various decay LDMEs with chromoelectric and chromomagnetic fields were determined in [24,30], where the authors systematically studied relativistic corrections to electromagnetic and hadronic decays of heavy quarkonia up to relative order v 7 .
As far as the production of heavy quarkonia is concerned, one should mention the study [31] of higher order relativistic corrections to the gluon fragmentation into a spin-triplet S-wave quarkonium. One of the LDMEs contributing at relative order v 4 involves a combination of a covariant derivative and a chromoelectric field. Also in this case the authors 2 Many NRQCD practitioners prefer to use the words "color octet" to denote only operators similar to the one in Eq. (9). To avoid misinterpretations we will speak in the following of operators that contribute through the leading or subleading quarkonium Fock state |QQ or |QQg respectively. could eliminate it from their basis via a field redefinition. We believe that the present work is the first study, in which matching coefficients multiplying LDMEs with chromoelectric fields are explicitly computed for a heavy quarkonium production process and cannot be traded with other operators. 3

C. Production cross section and power counting rules
After having clarified the NRQCD framework and the terminology of this work, let us proceed with applying NRQCD factorization to the exclusive process e + e − → γ * → χ cJ + γ. To assess the relative importance of the NRQCD matrix elements we adopt the power counting rules from [7]. To have a homogeneous counting, we count α s (m) ∼ v 2 , which will be important when combining our results with the existing radiative corrections. Each of our formulas depends on three LDMEs, where two of them are v 2 suppressed as compared to the leading LDME. When we refer to relative O(v 2 ) corrections, we always do so with respect to the absolute size of the leading order matrix element. Hence, at relative O(v 2 ) the NRQCD-factorized cross sections for e + e − → γ * → χ cJ + γ can be written as χ|0 + H.c.
3 As shown in [24], using field redefinitions the relevant LDMEs with chromoelectric fields could, in principle, be traded for LDMEs that involve a time derivative acting on fermion and gluon fields. We have chosen not to include such operators in our basis.
The subscript 1 labels operators contributing through the dominant Fock state |QQ , while the subscript 8 labels operators that project on the subleading Fock state |QQg only. According to the NRQCD power counting rules, the scaling of 0|O 1 ( 3 P J )|0 is v 5 , while 0|P 1 ( 3 P J )|0 scales as v 7 due to the presence of two additional covariant derivatives. In comparison to 0|O 1 ( 3 P J )|0 the matrix element 0|T 8 ( 3 P J )|0 contains one covariant derivative less, but involves a chromoelectric field, which scales as v 3 so that the absolute size of this LDME is also v 7 [7].
It is worth noting that, in general, expectation values of operators like Eq. (9), which contribute through the subleading Fock state |QQg , contain an additional suppression in v, if the QQ pair in the color octet configuration is produced with a different angular momentum than the QQ pair in the dominant Fock state. The reason for this is that in the corresponding short-distance process we create a color octet QQ pair that emits a soft gluon during its nonperturbative evolution into a heavy quarkonium. The emission of the gluon that changes the orbital angular momentum of the pair by one unit and converts it into a color singlet corresponds to an electric dipole (E1) transition, which accounts for the extra suppression. However, according to the power counting rules from [7], for LDMEs that contain chromoelectric fields, the additional O(v) suppression is already accounted for in the power counting rule of E. 4 This is why 0|T 8 ( 3 P J )|0 is only O(v 2 ) suppressed as compared to 0|O 1 ( 3 P J )|0 and hence parametrically as important as 0|P 1 ( 3 P J )|0 .
It is well known that the power counting of [7] is not the most conservative choice in the framework of NRQCD. In particular, the underlying assumption that the typical hadronic scale Λ QCD is of order mv 2 may not be satisfied in all quarkonium production and decay processes. A more cautious choice would be to assume that Λ QCD ∼ mv as was done e. g. in [24] (cf. also [6] and references therein for a detailed discussion of different countings in NRQCD). In the resulting conservative power counting, each chromoelectric field scales only as v 2 (not as v 3 ) but there is also an additional power of v coming from the subleading Fock state |QQg . It is therefore interesting to observe that 0|T 8 ( 3 P J )|0 scales as v 7 both in the perturbative counting of [7] and in the conservative counting of [24].
Equations (12) depend on nine matching coefficients. The matching coefficients F 1 ( 3 P J ) and G 1 ( 3 P J ) at O(α 0 s ) were computed in [11] and [15,32] respectively. In [15,32] (as well as in the order α s analysis of [16]) the authors, although working at O(v 2 ), did not include the matrix element 0|T 8 ( 3 P J )|0 in their cross section formula. Therefore the corresponding matching coefficients T 8 ( 3 P J ) remained unknown. In this work, we compute for the first time the matching coefficients T 8 ( 3 P J ) at O(α 0 s ). This is the last missing piece to have the complete O(v 2 ) corrections for exclusive electromagnetic production of χ cJ and a hard photon. The calculation is essentially a tree level calculation, but nontrivial, as we need to work with a quarkonium composed of two heavy quarks and a gluon. How to determine T 8 ( 3 P J ) from matching QCD to NRQCD will be discussed in the next section. III.

A. Matching conditions
There exist different approaches to calculate the NRQCD matching coefficients in a quarkonium production process. In all of them the matching is done between perturbative QCD and perturbative NRQCD, relying on the same behavior of both theories at low energies. The NRQCD matrix elements are evaluated in such a way that the quarkonium Fock states are replaced by the perturbative Fock states containing on-shell quarkonium constituents, e. g.
The explicit calculation of matrix elements on the right-hand side of Eq. (13) in perturbative NRQCD will be explained in Sec. IV. To derive the values of the coefficients multiplying such matrix elements, we need to impose a relation between suitable quantities (e. g., Green's functions, cross sections, on-shell amplitudes) in QCD and NRQCD. In EFTs such relations are called matching conditions. In [7] the matching condition for quarkonium production was given as which imposes an equality between the partonic cross section to produce an on-shell heavy quark pair in perturbative QCD and the sum of suitable LDMEs multiplied by matching coefficients F n and inverse powers of the heavy quark mass in perturbative NRQCD. Assuming that the relative momentum of the heavy quarks in the rest frame, q, is small, one can expand both sides of Eq. (14) in |q|/m. Within this nonrelativistic expansion one can read off the values of the matching coefficients F n and then substitute them into the NRQCDfactorized production cross section. A more technical explanation of this approach can be found in [33]. One of the difficulties related to practical applications of Eq. (14) is the necessity to perform a nonrelativistic expansion of the phase-space measure on the QCD side of the matching. In processes, where heavy quarkonium is produced together with other particles, such an expansion tends to become complicated and requires great care.
In case of exclusive reactions (such as our process of interest) it is also possible to employ the NRQCD factorization at the amplitude level [34,35], where c i 1 ...i k n is a short-distance coefficient. Equation (15) is an equality between the on-shell amplitude to produce a heavy quark pair in perturbative QCD and a sum of quarkonium-tovacuum matrix elements multiplied or contracted with short-distance coefficients in perturbative NRQCD. Once the short-distance coefficients c n are known, they can be substituted into the NRQCD-factorized production amplitude Squaring Eq. (16) and integrating over the phase space of the physical quarkonium we obtain the NRQCD production cross sections as in Eqs. (12). Notice that in general not all matrix elements in A pert NRQCD give a nonvanishing contribution to A NRQCD . Since QCD amplitudes do not have definite angular momentum (cf. Sec. V), in the matching we will determine more short-distance coefficients than required in Eqs. (12). The advantage of this method is that the matching can be done in a much simpler way, since it does not require one to compute the QCD matrix element squared and expand the phase space measure in |q|/m. In this work we calculate our matching coefficients by employing the NRQCD factorization at the amplitude level.

B. Kinematics
We now make explicit the kinematics that we will use throughout the matching calculation. We distinguish between two frames of reference. The laboratory frame is the center of mass (CM) frame of the colliding leptons, where the heavy quarkonium and the photon fly apart from the interaction point in opposite directions. The rest frame is the frame in which the heavy quarkonium is at rest.
We denote the 4-momenta of the heavy quarks and the gluon in the laboratory frame with p 1 , p 2 and p g respectively; k stands for the 4-momentum of the photon, while the momenta of the colliding leptons are labeled by l 1 and l 2 . We will also make use of the polarization vectors for the photon ε γ , for the gluon ε g and for the dilepton system L µ ≡ ev(l 2 )γ µ u(l 1 ). The direction of the photon 3-momentum is referred to ask ≡ k/|k|. In the matching all the external momenta are put on-shell and the masses of the leptons are neglected as compared to the CM energy √ s. Hence, Finally, the 4-momentum of the heavy quarkonium is denoted by P . In the perturbative matching P is expressed as a sum of the heavy quark and gluon momenta, i. e. P = p 1 + p 2 for a QQ system and P = p 1 + p 2 + p g for a QQg system. In the physical process, however, P is the 4-momentum of the quarkonium with P 2 = M 2 χ cJ , where M χ cJ denotes the mass of χ cJ measured in experiments. We will explicitly state the meaning that we assign to P at the different stages of this work.
To distinguish a laboratory frame vector from a rest frame vector, the latter will be assigned an additional subscript R. For example, when considering the QQ system, it is convenient to introduce the relative momentum of the heavy quarks in the rest frame, defined as q ≡ (p 1,R − p 2,R )/2. In the case of the QQg system, two relative momenta are needed, given by More details regarding the kinematics of the QQ and QQg systems can be found in Appendix B.

C. Automatized expansions
As all matching prescriptions involve nonrelativistic expansions of the perturbative QCD amplitudes, let us also briefly discuss different strategies to realize this task.
In [7], matching calculations for decays were done in such a way that the expanded QCD amplitudes were explicitly rewritten in terms of energies, Cartesian 3-vectors, Pauli matrices and Pauli spinors. Since such objects naturally arise on the NRQCD side of the matching, having both sides of the matching equation expressed as products of 3-vectors and Pauli structures facilitates the extraction of the matching coefficients. This approach to NRQCD matching calculations was developed further and generalized to be applicable both to production and decays in [25] and [36]. Following [36] we will denote it as the threshold expansion method. Applying the technique developed by Braaten and Chen is conceptually simple, yet the calculations tend to become rather cumbersome when one goes beyond tree level or leading order in v. In particular, the necessity to work with noncovariant objects on the QCD side of the matching makes such calculations quite tedious not only when done by hand but also when automatized with FORM [37] or similar symbolic manipulation systems.
The manifest Lorentz covariance of the QCD amplitudes can be preserved if one uses the covariant projector technique [28]. This property is especially useful when the evaluation is automatized using existing software for symbolic calculations in relativistic Quantum Field Theories (QFTs( (e. g. FeynCalc [38,39], Redberry [40] or FDC [41] to name those that are often used in NRQCD calculations). This is one of the main reasons why almost all modern NRQCD matching calculations are done using projectors.
For the present calculation, however, we will not use this approach. The reason is that there is no literature on how to apply the projector technique to extract matching coefficients multiplying matrix elements of the type 0|T 8 ( 3 P J )|0 . On the other hand, the calculation of the matching coefficients T 8 ( 3 P J ) using the threshold expansion method and matching at the amplitude level is straightforward, albeit tedious. Similar calculations have already been done in the investigation of the electromagnetic decay χ cJ → γγ [23,24], such that one can benefit from the existing knowledge on the subject. To automatize our matching calculations we employ FeynArts [42], FeynCalc, FORM and FeynCalcFormLink [43] as well as self-written Mathematica codes. The software framework for automatizing nonrelativistic EFT calculations using FeynCalc will be presented elsewhere [44].

A. NRQCD-factorized production amplitudes and cross sections
In order to apply the threshold expansion method at the amplitude level, we introduce NRQCD production amplitudes at O(v 2 ) relative to the scaling of the leading order quarkonium-to-vacuum matrix element The scaling of these matrix elements in v is estimated according to the power counting rules discussed in Sec. II C. For example, the matrix element As explained in Sec. III, to arrive at production cross sections in the laboratory frame, we need to square the amplitudes in Eqs. (21), average over the polarizations of the leptons, sum over the polarizations of the photon and the heavy quarkonium and integrate over the phase space for producing a heavy quarkonium and a photon, where pols refers to the photon polarizations and the summation over the polarizations of the heavy quarkonium is implicit in the amplitude squared. The prefactor (2J + 1) appears because of our redefinition of the matching coefficients in Eq. (8) and will be precisely canceled by 1/(2J + 1) in the formula for the sum over the heavy quarkonium polarizations (cf. Appendix C). The prefactor 2M χ cJ accounts for the fact that our matrix elements have nonrelativistic normalization. 5 Finally, the phase space measure is given by with the integration done over the solid angle of the photon 3-momentum k. This yields where we used that c 1 c * 3 = c 3 c * 1 and c * 1 d 1 = −c 1 d * 1 . 6 Comparing Eq. (24) with Eq. (12) we can express the matching coefficients F 1 ( 3 P J ), G 1 ( 3 P J ) and T 8 ( 3 P J ) through the short-distance coefficients c 1 , c 3 and d 1 Our task is to calculate all matching coefficients and in particular T 8 ( 3 P J ) at O(α 0 s ), which can be inferred from the knowledge of c 1 and d 1 . To determine the short-distance coefficients c J 1 , c J 3 and d J 1 we need to carry out the matching between perturbative QCD and perturbative NRQCD at the amplitude level.

B. Computation of the NRQCD amplitudes
The perturbative NRQCD matrix elements are most conveniently evaluated in the operator formalism by rewriting heavy quark fields in terms of Pauli spinors and single fermion creation and annihilation operators, as explained in [45]. The evaluation of matrix elements that enter perturbative NRQCD amplitudes is straightforward. For example, replacing |H with |Q(p 1,R )Q(p 2,R ) and Notice that, although matrix elements made only of covariant derivatives contribute through the leading Fock state |QQ , because of the gluon field in the covariant derivative they also give a nonvanishing contribution to the subleading Fock state |QQg . As explained in Sec. III, to match the QCD amplitudes we need, however, more operators and matching coefficients on the NRQCD sides than needed to compute production cross sections alone. For example, apart from P -wave spin triplet quarkonium (χ cJ ) we may also produce S-wave spin singlet quarkonium (η c ). Even though we are actually not interested in the latter, we still have to determine the short-distance coefficients of the corresponding operators. In fact, this property of the threshold expansion method provides an important consistency check of the whole matching calculation: if the short-distance coefficients were determined correctly, the difference A pert QCD − A pert NRQCD would necessarily vanish order by order in the expansion parameters. Notice that the perturbative matching between QCD and NRQCD does not rely on any specific power counting. Instead, we compare A pert QCD to A pert NRQCD order by order in q 1 and q 2 . This follows the approach adopted in [24], where the matching for electromagnetic decays was done order by order in 1/m.
The NRQCD amplitudes relevant for the perturbative matching (including additional matrix elements and Lagrangian insertions allowed by the symmetries) read Some explanations are in order. The matrix elements with the subscript L 2−f contain Lagrangian insertions, where the relevant part of the NRQCD Lagrangian (2-fermion sector only) is In the diagrammatic representation of the NRQCD quarkonium-to-vacuum matrix elements, Lagrangian insertions correspond to the emission of a gluon from one of the external heavy quark lines. For example, for H|ψ † χ|0 and H|ψ † χ|0 L 2−f we have where the quark-gluon vertex on the heavy quark (antiquark) line is understood as a sum of all six vertices from quark-gluon (antiquark-gluon) interaction terms in Eq. (33). Hence, to calculate the diagrams on the left-hand side of Eq. (34b) we need to apply Feynman rules for the ψ † χ operator and for the interaction terms in Eq. (33).
The derivation of such Feynman rules in the operator approach is a simple QFT exercise: we sandwich the operator between the vacuum and a Fock state that contains our fields (so that all momenta are ingoing), calculate the matrix element, amputate external states (Pauli spinors and polarization vectors) and multiply the result by i. For example, for the operator ψ † {D 2 , σ · gB}/(8m 3 )ψ we obtain a, i Once we have the full expressions for the amplitudes in Eq. (32) with |H replaced by the perturbative |QQ and |QQg Fock states, we can expand them in the relative momenta of the quarkonium constituents. For the QQ case we need to expand up to the third order in |q|/m, which allows us to determine the values of c 0 , c 1 , c 2 and c 3 . For the QQg case it is sufficient to expand the amplitude up to the first order in |q 1 |/m and |q 2 |/m, so that we can extract d 0 and d 1 .

V. QCD AMPLITUDES AND MATCHING
In this section, we compute the nonrelativistic expansions for the relevant QCD amplitudes that describe exclusive electromagnetic production of QQ, i. e.
The process in Eq. (36) is described by two QCD Feynman diagrams (cf. left panel of Fig. 1) so that while the reaction in Eq. (37) involves six diagrams (cf. right panel of Fig. 1) and the corresponding amplitude reads where L µ has been defined in Sec. III B and e Q denotes the heavy quark electric charge in units of e. The amplitudes A QQ pert QCD and A QQg pert QCD can be evaluated in the laboratory frame or in the rest frame. Both frames have their advantages and disadvantages. In the course of this study we carry out the calculations in both frames and verify that we arrive at the same results for the final matching coefficients. Figure 1: Representative QCD diagrams for the processes e + e − → QQ + γ and e + e − → QQg+ γ that are relevant for the matching to NRQCD. For the production of QQ there are only two diagrams in total: the one displayed here and a second one where the photon is emitted from the heavy antiquark line. The production of QQg is described by six diagrams that differ in the quark (antiquark) line that emits the photon and/or the gluon.
It is convenient to decompose all Dirac structures involving heavy quarks into scalar, pseudoscalar, vector, axial-vector and tensor components. Since this is a tree level calculation in four dimensions, such decomposition can be done in an unambiguous way. Then the amplitudes can be written as where the coefficients j µ i do not depend on the heavy quark spinors. The vector and axialvector structuresū(p 1 )γ µ v(p 2 ) andū(p 1 )γ µ γ 5 v(p 2 ) can be rewritten in terms of Pauli matrices and spinors using formulas given in Appendix B. Furthermore, all the propagators and scalar products in j µ i have to be expressed in terms of q (for QQ) or of q 1 and q 2 (for QQg). As in Sec. IV, we expand up to the third order in |q|/m and up to the first order in |q 1 |/m and |q 2 |/m.
The NRQCD amplitudes in Eq. (32) are given for specific values of the total angular momentum J, but A QQ pert QCD and A QQg pert QCD should be understood as sums of contributions from different values of J. Therefore, on the QCD side one has to identify Cartesian tensors that correspond to the matrix elements on the NRQCD side and decompose them into irreducible spherical tensors along the lines of [46]. For example, a term that contains the rank 2 tensor σ i q j and no further occurrences of q contributes to from A J=0 pert NRQCD but also to from A J=2 pert NRQCD . We can disentangle these contributions by explicitly projecting out the J = 0, 1, 2 components of the tensor, i. e., Applying such projections to the whole QCD amplitudes allows us to do the matching between QCD and NRQCD for each J-value separately and successfully extract the shortdistance coefficients order by order in the relative momenta. Before presenting our results for the short-distance coefficients, we would like to mention that on the QCD side of the matching we encounter terms that are singular in the limit |p g | → 0, i. e., proportional to 1/|p g |. Such singularities arise when the gluon is emitted from one of the external heavy quark lines and should cancel in the matching. Indeed, matrix elements with the Lagrangian insertions given by Eq. (33) on the NRQCD side of the matching generate terms that precisely reproduce the singularities of the QCD amplitude. The identical infrared behavior of the full QCD and NRQCD amplitudes follows from general principles underlaying the construction of low-energy effective field theories. The cancellation of the 1/|p g | singularities on both sides of the matching has been explicitly checked in both the laboratory frame and the rest frame. A more detailed discussion on this subject can be found in [24] and [30].

A. Matching in the rest frame
Rest frame short-distance coefficients carry an additional subscript R to distinguish them from the coefficients obtained in the laboratory frame. They read with λ ≡ e 2 e 2 Q /s and r ± ≡ (1 ± r) −1 , where r ≡ 4m 2 /s is the kinematic suppression factor. L R is the spatial part of the leptonic current in the rest frame defined in Sec. III B.
The connection to the laboratory frame can be established by means of the following Lorentz boost transformations where the appearance of the physical quarkonium mass stems from the fact that we have employed the kinematics of the physical quarkonium, i. e. used P 2 = M 2 χ cJ . In Sec. V C we will show how the dependence on M χ cJ can be eliminated up to the desired order in v.
Matching in the rest frame of the heavy quarkonium is very similar to the calculation of the O(α 0 s v 2 ) corrections to the decay process χ cJ → γγ. The production of χ cJ at rest can be regarded as the decay of a heavy vector boson γ * into a photon and a quarkonium. As has already been observed in [11], in the limit s → 0, the production short-distance coefficients from Eqs. (46) should reduce into the short distance coefficients for the decay of χ cJ into two photons. In our calculation this is indeed the case.

B. Matching in the laboratory frame
In the laboratory frame it is straightforward to carry out the phase space integrations in Eq. (25) and thus arrive at the final matching coefficients. The more complicated part of the calculation is the manipulation of the QCD amplitudes. Obviously, potentially large laboratory frame momenta of the heavy quarks and gluon are not suitable parameters for doing a nonrelativistic expansion. The correct expansion can be done only after these momenta have been rewritten in terms of small rest frame momenta, which greatly increases the number of intermediate terms.
Lorentz boost transformations that relate potentially large laboratory frame momenta of a QQ system to the small rest frame momenta in a way that is useful for NRQCD matching calculations were introduced in [25]. In this work we generalize those formulas for a QQg system and summarize them in Appendix B.
The matching procedure yields the following values for the short-distance coefficients:

C. Relations between matrix elements and heavy quarkonium masses
From the NRQCD point of view, the LDMEs should be regarded as independent nonperturbative parameters. Nevertheless, symmetries and approximations can be used to establish relations between different LDMEs that are valid up to a certain order in v.
Particularly useful relations can be derived from the equations of motion of NRQCD and are known as Gremm-Kapustin relations [29]. Since our LDMEs directly follow from those that appear in the electromagnetic decay χ cJ → γγ, we can use the Gremm-Kapustin relations 7 that were first obtained in [23], where E χ cJ = M χ cJ −2m is the binding energy of χ cJ . Equation (49) is valid up to corrections of O(v 4 ) and it is useful for two reasons. First, we can solve it for M χ cJ and use Eq. (50) to eliminate the dependence of the matching coefficients on the heavy quarkonium mass in Eqs. (25). This is a crucial ingredient to check that the matching coefficients that we will list explicitly in Sec. VII are the same if obtained from the expansion of the QCD amplitudes in the laboratory frame or in the rest frame. Second, Eq. (49) allows us to reduce the number of unknown matrix elements in the final NRQCD-factorized production cross sections at the cost of making the numerical predictions depend on the heavy quarkonium binding energy. We will see this in the following Sec. VIII.

VI. ELECTROMAGNETIC DECAYS OF χ cJ INTO TWO PHOTONS
In this section, we show that our matching calculations in the rest frame also contribute to the resolution of the known discrepancy between the results of [23] and [24] on the values of some matching coefficients entering the NRQCD-factorized decay rates for χ cJ → γγ at O(α 0 s v 2 ). The formulas for these decay rates read (using a remark made in Sec. II A we express eventually the widths in terms of production LDMEs 8 ) The factor 1/2 in front of m 0|T 8 ( 3 P J )|0 accounts for the fact that 0|T 8 ( 3 P J )|0 = 2 χ cJ |F EM ( 3 P J )|χ cJ , with χ cJ |F EM ( 3 P J )|χ cJ being the chromoelectric LDME from [23]. 8 With the definitions of the decay and production LDMEs from [24] and Eqs. (12) respectively, we have 0|{O, P} 1 ( 3 P J )|0 = χ cJ |{O, P} 1 em ( 3 P J )|χ cJ and 0|T 8 ( 3 P J )|0 = − χ cJ |T 8 em ( 3 P J )|χ cJ .
While both collaborations obtained the same expressions for Im g em ( 3 P 0 ) = −7α 2 e 4 Q π, (52c) they disagree on the values of the other coefficients; cf. Table I. In 2013, it was reported [47,48] that an independent investigation of the χ cJ → γγ decays confirmed the results of [24]. In the course of this study we could repeat the calculation of the O(α 0 s v 2 ) corrections to these decay processes. Matching at the amplitude level and at the level of the total decay rates we also found agreement with the values of the matching coefficients obtained in [24]. These are Note that in [23] the decay rate of χ c2 contains the additional matrix element Since in [24] it was shown that the contribution of this operator to the total decay rate is already accounted for by the inclusion of 0|P 1 ( 3 P 2 )|0 , we do not include it in the expression for the decay width.
Plugging the laboratory frame short-distance coefficients from Eqs. (48) into Eqs. (25) we obtain the values of the matching coefficients that enter the factorized cross sections given in Eqs. (12). If we define where θ denotes the angle between the beam line and the outgoing photon in the laboratory frame, the results read and In order to write the production cross sections, for consistency we combine our findings with the known O(α s v 0 ) corrections to F 1 ( 3 P J ) computed in [12,17]. The explicit values of the coefficients C i j (r) can be found in Appendix B of [12]. The O(α s v 2 ) corrections from [16] are, however, not included. The first reason is that the corrections reported in [16] are incomplete, as they do not include contributions from 0|T 8 ( 3 P J )|0 . The second one is that, to be consistent with the NRQCD power counting rules at that order in α s and v, one would also need to include O(α 0 s v 4 ) and O(α 2 s v 0 ) corrections, which are not available yet. The differential and total cross sections read where we recall that r ≡ 4m 2 /s. The dependence on the physical quarkonium mass arises due to the presence of M χ cJ in Eq. (22) and Eq. (23). This dependence can easily be eliminated via the Gremm-Kapustin relations. In particular, using Eq. (50) we arrive at cross sections that do not explicitly depend on the mass of the χ cJ . These are and σ(e + e − → χ c2 + γ) = (4πα) 3 e 4 Q (1 + 3r + 6r 2 ) 9πm 3 s 2 (1 − r) Comparing these results with the literature, we observe that the matching coefficients F 1 ( 3 P J ) agree with those reported in [11] and in all subsequent studies. Furthermore, our values for G 1 ( 3 P J ) are consistent with the K-factors given in [14]. The values of the matching coefficientsT 8 ( 3 P J ) and T 8 ( 3 P J ) are original of this work.

VIII. NUMERICAL RESULTS
The main limitation in the determination of the e + e − → γ * → χ cJ + γ production cross sections comes from the uncertainties in the LDMEs. Nevertheless we will see that we can express the LDMEs in such a way to reduce those uncertainties in the final results.
The error corresponds to an O(v 2 ) uncertainty that we estimate to be 30% of the central value. 9 This uncertainty accounts also for the breaking of the heavy-quark spin symmetry. Further, we notice that within the error bounds the choice in Eq. (64) is consistent with the most accurate values of 0|O 1 ( 3 P 0,2 )|0 determined in Table I of the recent study [50] devoted to O(α 2 s v 0 ) corrections of the χ cJ → 2γ decay. The number of the remaining independent LDMEs can be reduced by applying the heavyquark spin symmetry [7], which relates matrix elements with the same orbital momenta but different spin values and is valid up to corrections of O(v 2 ). Since 0|P 1 ( 3 P J )|0 and 0|T 8 ( 3 P J )|0 are already v 2 suppressed as compared to 0|O 1 ( 3 P J )|0 , the corrections from the breaking of the heavy-quark spin symmetry to the total cross section are of O(v 4 ) and therefore irrelevant for the accuracy that we are aiming at: This is clearly different from the case of 0|O 1 ( 3 P J )|0 where spin symmetry breaking effects, being relevant at O(v 2 ), cannot be neglected at our accuracy and have been included in the uncertainty of Eq. (64).
The values of 0|P 1 ( 3 P 0 )|0 and 0|T 8 ( 3 P 0 )|0 can be determined from the electromagnetic decays χ c0 → γγ and χ c2 → γγ that have recently been measured by the BES-III experiment [51], where the first error is statistical, and the second one consists of systematic errors and the error in the branching fraction combined in quadrature. Combining the NRQCD results for the electromagnetic decay rates in Eqs. (51) with the Gremm-Kapustin relation in Eq. (49) we obtain and The decay rate formulas in Eqs. (69a) and (70a) contain not only corrections of O(α 0 s v 0 ) and O(α 0 s v 2 ), but also the O(α s ) correction to Im f em ( 3 P 0,2 ) [52,53]. This is required for consistency with the NRQCD power counting, where corrections of O(α s v 0 ) and O(α 0 s v 2 ) are assumed to be equally important.
The mass m appearing in the above equations is the charm pole mass. The pole mass is a poorly known quantity that, however, can be traded for the physical masses of the χ cJ , which are experimentally known, and the binding energy, which is one of our unknowns, according to In order to be consistent with the NRQCD power counting, we replace m with (M χ cJ − E χ cJ )/2 in contributions that are of order v 0 and use m = M χ cJ /2 in contributions that are of order v 2 (with α s ∼ v 2 ). Solving the system for the unknown LDMEs and expanding in E χ cJ ∼ v 2 up to relative order v 2 yield for the determination using Eqs. (69) and for the determination using Eqs. (70). Furthermore, we choose where for M χ cJ we use the most recent PDG values [54], while the values of the strong coupling constant (at 1-loop accuracy) were obtained with RunDec [55][56][57]. The errors of α s from varying M χ cJ are of relative order 10 −5 and hence can safely be neglected in comparison to other error sources in the present analysis. The matrix elements 0|P 1 ( 3 P 0 )|0 and 0|T 8 ( 3 P 0 )|0 depend on the binding energy, E χ cJ . Since of the binding energy we only know the nonrelativistic scaling, E χ cJ ∼ mv 2 , where v 2 ≈ 0.3 for charmonia [7], it may be useful to parametrize the binding energy in terms of a dimensionless parameter x J of O(1) and write it as The LDMEs 0|P 1 ( 3 P 0 )|0 and 0|T 8 ( 3 P 0 )|0 as a function of x J are shown in Fig. 2. The uncertainties coming from M χ cJ , Γ(χ cJ → γγ) and 0|O 1 ( 3 P 0 )|0 have been estimated using Gaussian error propagation and are shown by the band. Because the Buchmüller-Tye potential or the Cornell potential or a purely confining potential would all give a positive binding energy, 10 we will restrict x J to positive values centered around the natural value x J = 1.0. Allowing for a 100% uncertainty, we take With this choice we obtain 0|P 1 ( 3 P 0 )|0 = (0.103 ± 0.092 ± 0.070 ± 0.026) GeV 7 , (81a) from the determination using Eqs. (69) and from the determination using Eqs. (70), where the uncertainties were estimated using the Gaussian error propagation formula. The first error originates from the uncertainty in 0|O 1 ( 3 P 0 )|0 , the second one originates from the uncertainty in x J , and the third one is the combined (in quadrature) uncertainty in M χ cJ and Γ(χ cJ → γγ).
The values of the LDMEs given in Eqs. (81) and Eqs. (82) are consistent with each other within the error bounds, as expected from the heavy-quark spin symmetry. Under the assumption that the heavy-quark spin symmetry holds at leading order in the velocity expansion, we can average the two determinations, which yields our best determination for these LDMEs 0|P 1 ( 3 P 0 )|0 = (0.058 ± 0.074 ± 0.086 ± 0.033) GeV 7 , (83a) where the uncertainties in 0|O 1 ( 3 P 0 )|0 and x J were added linearly [first and second uncertainties in Eqs. (83) respectively], while those in M χ cJ and Γ(χ cJ → γγ) (third uncertainty) are regarded as independent and were therefore added quadratically. The values given in Eqs. (83) represent our best determination of 0|P 1 ( 3 P 0 )|0 and 0|T 8 ( 3 P 0 )|0 using the heavy-quark spin symmetry.  Figure 3: Cross sections for the exclusive production of a χ cJ and a hard photon in the energy region of Belle II. The matrix elements 0|P 1 ( 3 P 0 )|0 and 0|T 8 ( 3 P 0 )|0 were determined from Eqs. (69). The dash-dotted curve shows only the contribution from the 0|O 1 ( 3 P 0 )|0 matrix element multiplied with the tree level matching coefficient. The dotted curve also includes the loop correction to the matching coefficient of 0|O 1 ( 3 P 0 )|0 . The dashed curve incorporates the contribution from 0|P 1 ( 3 P 0 )|0 but not 0|T 8 ( 3 P 0 )|0 . Finally, the solid curve displays our final result that contains all the relevant contributions.
In order to reduce the uncertainties in the total cross sections σ(e + e − → χ cJ + γ), it is useful to compute them by replacing the matrix elements 0|P 1 ( 3 P 0 )|0 and 0|T 8 ( 3 P 0 )|0 with their expressions from Eqs. (72) or Eqs. (73). In this way each cross section will have five independent sources of uncertainty, which are M χ cJ , Γ(χ cJ → γγ), 0|O 1 ( 3 P 0 )|0 ), x J and the renormalization scale µ that enters through α s .
The process e + e − → χ cJ + γ can be measured at B-factories with a sufficiently high CM energy to produce a χ cJ and a hard photon, e. g., at Belle II. Hence, we select our input parameters as √ s = 10.6 GeV, In principle, one power of α that originates from the photon emission of the heavy quark (cf. Fig. 1) could run as low as α(M χ cJ /2). However, for simplicity we will ignore this potential uncertainty and evaluate all α at the same value. To estimate the uncertainty in the choice of the renormalization scale of the strong coupling, we use four possible values of µ   Here the only relevant error source is the uncertainty in 0|O 1 ( 3 P 0 )|0 . The loop corrections are added in the second column. In the third column we include contributions from 0|P 1 ( 3 P 0 )|0 but not from 0|T 8 ( 3 P 0 )|0 . In the fourth column all three LDMEs are included. The last column displays the ratio of the two final cross sections, when 0|T 8 ( 3 P 0 )|0 is included (σ 8 ) or omitted (σ 1 ). frame). In Fig. 7 we also show the ratios of the total cross sections, where the LDME 0|T 8 ( 3 P 0 )|0 is included (denoted as σ 8 ) or omitted (denoted as σ 1 ) as a function of √ s.
As for the influence of the additional LDME 0|T 8 ( 3 P J )|0 on the total cross section, we refer to the last columns of Table II and III and to Fig. 7. We see that the contribution of 0|T 8 ( 3 P J )|0 is small, but, at its central value, not negligible and thus could be potentially significant for all χ cJ .  Tables II and III. While our analytic results for the matching coefficients are unambiguous, the numerical results strongly depend on the estimated size of the quarkonium binding energy and the assumed value of 0|O 1 ( 3 P J )|0 . Moreover, the experimental errors in the knowledge of Γ(χ cJ → γγ) are also not negligible. A reliable determination of 0|O 1 ( 3 P J )|0 , 0|P 1 ( 3 P 0 )|0 and 0|T 8 ( 3 P 0 )|0 from experimental data or lattice calculations is therefore highly desirable and would provide a major improvement of the present cross section determinations.

IX. SUMMARY
NRQCD factorization provides a systematic prescription to compute production cross sections of heavy quarkonia order by order in α s and v. Higher order corrections in α s arise from loop effects, while higher order corrections in v are generated by the inclusion of higher-dimensional operators. A consistent computation of higher order corrections requires both loop corrections to matching coefficients of lower dimensional operators and relativistic corrections from higher dimensional operators.
In this work we have computed the O(α 0 s v 2 ) contribution induced by the higher Fock state |QQg to the e + e − → γ * → χ cJ + γ cross section, in this way completing the order v 2 corrections to this P -wave quarkonium production process. The new pieces of information contained in this work are the leading order expressions of the matching coefficients T 8 ( 3 P J ) multiplying LDMEs that depend on the chromoelectric field.
We match QCD to NRQCD perturbatively at the amplitude level, where we treat the final state gluon as a part of the quarkonium system, together with the heavy quark pair. In the course of the calculation we encountered infrared singularities induced by the emission of the gluon from a heavy quark line. We explicitly verified that such terms cancel exactly in the matching. As a further consistency check we performed the matching calculation in two different reference frames: the CM frame of the colliding leptons and the rest frame of the heavy quarkonium. Although these two calculations are technically very different, they lead to the same matching coefficients.
The process e + e − → γ * → χ cJ + γ should be observable at B-factories with sufficiently high CM energy, e. g. at Belle II. As no published experimental results are available so far, the phenomenological relevance of contributions from higher Fock states remains to be tested.
for J = 2 states. We provide their values in the laboratory frame and in the rest frame as a reference for future studies. For the notation we refer to Sec. V.
Appendix B: Threshold expansion method

2-body system
We consider, first, the case of a quarkonium state described by a QQ pair whose quark and antiquark carry momenta p 1 and p 2 . We follow [25]. These momenta can be written in terms of the Jacobi momenta of a 2-body system as where P = p 1 + p 2 denotes the total momentum of the quarkonium, while Q = (p 1 − p 2 )/2 stands for the relative momentum between the heavy quark and the heavy antiquark.
Here q is indeed soft and can be used as an expansion parameter. The two frames are related by a Lorentz transformation, which means that we can express Q and P in terms of q. According to [25] the Lorentz transformation reads where P 0 = 4E 2 q + P 2 .
The Dirac spinors in the rest frame are given by where N is the normalization factor. In case of relativistic normalization, it is given by E q + m, while in the nonrelativistic case one has N = (E q + m)/(2E q ). When boosted to the laboratory frame they become u(p 1 ) = 1 4E q (P 0 + 2E q ) (2E q + / P γ 0 )N ξ For practical purposes, it is more convenient to work with boosted Dirac bilinears. Only vector and axial vector bilinears appear in our process of interest; they can be written as

3-body system
We now turn to the case of a quarkonium state composed by a heavy quark, a heavy antiquark and a gluon. The momenta of the three constituents can be written in terms of the Jacobi momenta of a 3-body system as p 1 = 1 3 P + Q 1 − Q 2 , p 2 = 1 3 P − Q 1 − Q 2 , p g = 1 3 P + 2Q 2 .