A Useful one-loop scalar integrals

# $O(α_s v^2)$ correction to pseudoscalar quarkonium decay to two photons

## Abstract

We investigate the correction to the process of pseudoscalar quarkonium decay to two photons in nonrelativistic QCD (NRQCD) factorization framework. The short-distance coefficient associated with the relative-order NRQCD matrix element is determined to next-to-leading order in through the perturbative matching procedure. Some technical subtleties encountered in calculating the QCD amplitude are thoroughly addressed.

###### pacs:
12.38.-t, 12.38.Bx, 13.20.Gd, 13.40.Hq, 14.40.Gx

## I Introduction

Quarkonium inclusive annihilation decays are historically among the earliest applications of perturbative quantum chromodynamics (QCD) (1); (2); (3); (4). At present, it is widely accepted that these quarkonium inclusive decay processes can be systematically described by the nonrelativistic QCD (NRQCD) factorization approach (5), which is based on the effective-field-theory formalism and directly linked with the first principles of QCD.

In NRQCD factorization approach, the inclusive decay rate can be systematically expressed as the sum of product of short-distance coefficients and the NRQCD operator matrix elements. The short-distance coefficients encode the hard effects of quark and antiquark annihilation at length scale of order ( denotes the mass of heavy quark), which therefore can be accessed by perturbation theory owing to asymptotic freedom of QCD. In contrast, the NRQCD matrix elements are sensitive to the nonperturbative dynamics that occurs at distance of or longer ( denotes the typical (anti-)quark velocity inside a quarkonium). An important feature of these nonperturbative matrix elements is that they are universal, and satisfy a definite power-counting in . This attractive feature endows the NRQCD factorization approach with a controlled predictive power.

Among the quarkonium annihilation decay processes, the simplest and cleanest are the quarkonium electromagnetic decays, exemplified by vector quarkonium decay to a lepton pair and pseudoscalar quarkonium decay to two photons. Both of these two processes have been comprehensively studied in theory and experiment for decades. In this work, we will address the correction to the latter process in NRQCD factorization context. We note that investigations on this process from other approaches are also available (e.g., see (6); (7); (8)).

We will generically label a pseudoscalar quarkonium by , where can stand for the or the quark. Through relative order-, the NRQCD factorization formula for the decay rate of reads (9); (10):

 Γ[ηQ→γγ] = F(1S0)m2∣∣⟨0|χ†ψ|ηQ⟩∣∣2+G(1S0)m4Re{⟨ηQ|ψ†χ|0⟩⟨0|χ†(−i2\tensorD)2ψ|ηQ⟩} (1) + H1(1S0)m6⟨ηQ|ψ†(−i2\tensorD)2χ|0⟩⟨0|χ†(−i2\tensorD)2ψ|ηQ⟩ + H2(1S0)m6Re{⟨ηQ|ψ†χ|0⟩⟨0|χ†(−i2\tensorD)4ψ|ηQ⟩}+O(v6Γ),

where and represent Pauli spinor fields that annihilate the heavy quark and heavy antiquark , respectively.

The short-distance coefficients , , in (1) are

 F(1S0)=2πe4Qα2[1+CFαsπ(π24−5)+O(α2s)], (2a) G(1S0)=−8πe4Qα23[1+O(αs)], (2b) H1(1S0)+H2(1S0)=136πe4Qα245[1+O(αs)], (2c)

where denotes the electric charge of the heavy quark . The correction to was first computed in (11); (12); (13). An incomplete correction to this coefficient has also been available about a decade ago (14), which indicates an uncomfortably large negative correction. The coefficient is associated with the order- matrix element, whose tree-level value has been known long ago (15). Recently, the tree-level coefficients associated with the order- matrix elements are also available for the first time (9); (10). Initially only the combination of two coefficients was given (9), as shown in (2c). Later Ref. (10) was able to determine these two coefficients separately: and .

In recent years there have also been many phenomenological investigations on the nonperturbative NRQCD matrix elements for the pseudoscalar quarkonium state (16); (17); (18). Most of them focus on the matrix elements at lowest order (LO) and next-to-leading order (NLO) in . Comparing equation (1) truncated at order with the measured decay rate of , some authors are able to fit the first two NRQCD matrix elements of and find they roughly obey the velocity counting rules (16); (17) 5.

By far, the contribution (14) and contribution (9); (10) to the decay rate of are known. Assuming , one may naturally wonder what is the actual size of the correction. In order to answer this question, one needs first know the correction to the short-distance coefficient . It is the purpose of this work to compute this correction through the perturbative matching method.

The remainder of this paper is organized as follows. In section II, we outline the perturbative matching strategy which can be utilized to deduce the NRQCD short-distance coefficients for our process. In section III, we elaborate on some technical issues encountered in calculating the correction to with the covariant projection approach. In particular, we specify the prescription of in dimensional regularization adopted in this work. We also mention some technical ambiguities about extracting the -wave amplitude. In section IV, we then employ the covariant projection technique to compute the amplitude to NLO in . The corresponding calculation in the NRQCD side is presented in section V. In section VI, by comparing the NLO QCD amplitude and the respective NRQCD amplitude, we determine the first two short-distance coefficients in (1) through order . Finally in section VII, we present a brief summary. The appendix is devoted to enumerating the analytic expressions of the one-loop scalar integrals encountered in the calculation of the NLO QCD correction to .

## Ii The method of deducing the short-distance coefficients

In many effective field theories, it is a standard procedure to determine the short-distance coefficient through the matching procedure. We will follow this orthodox method in this work. Our calculation is in close analogy with Ref. (21), where the correction for has been deduced also using the matching approach.

It is worth mentioning that the alternative, even more efficient way of deducing the short-distance coefficients exists, exemplified by the method of region developed by Beneke and Smirnov (22). This approach allows one to dissect the loop integral for quarkonium annihilation process into four distinct regions: hard, soft, potential and ultrasoft. The short-distance coefficient only receives the contribution from the hard region, while the NRQCD effective theory characterizes the dynamics from the three low-energy regions. To determine the short-distance coefficient, one can either directly compute the hard region contribution order by order in expansion, or by first calculating the full QCD diagram, then subtracting the contributions from three low energy regions. The second strategy turns out to have some great technical advantage, which enables one to deduce a class of relativistic corrections to arbitrary orders in with ease. Following this spirit, there have recently been progresses in inferring the all-order-in- corrections to  (23) and (24).

### ii.1 NRQCD factorization formula for the decay amplitude and width

Let us label the momenta of and two photons by , and , and the polarization vectors of two photons by and , respectively. The Lorentz and parity invariance dictates the amplitude of uniquely of the structure . For the problem at hand, it is most natural to work in the rest frame. In this frame, the amplitude is then proportional to the kinematic invariant .

When considering quarkonium electromagnetic annihilation decay, one can directly invoke NRQCD factorization at the amplitude level. We need only retain those NRQCD color-singlet operator matrix elements that connect the vacuum to the state. To the order of desired accuracy, one expects that the following factorization formula holds:

 A[ηQ→γγ] = ^k1⋅\boldmathε∗1×\boldmathε∗2[c0⟨0|χ†ψ|ηQ⟩+c2m2⟨0|χ†(−i2\tensorD)2ψ|ηQ⟩+O(v4)], (3)

where () signify the short-distance coefficients associated with the NRQCD matrix elements at LO and NLO in . It is especially convenient in this work to adopt the nonrelativistic normalization for state in both sides of (3). In addition, we have explicitly factored out the kinematic invariant, which is represented by a dimensionless Lorentz pseudoscalar,

 ^k1⋅\boldmathε∗1×\boldmathε∗2 = −2M2ηQϵαβμνPαkβ1ε∗μ1ε∗ν2, (4)

in the right-hand side of (3), where is a unit vector. The separation of this kinematic factor renders () to bear a very simple form.

Squaring the amplitude in (3), summing over the polarizations of photons, and integrating over the phase space and accounting for the indistinguishability of two photons, we can express the decay rate of as

 Γ[ηQ→γγ] = 12!∫d3k1(2π)32k01d3k2(2π)32k02(2π)4δ(4)(P−k1−k2)∑∣∣A[ηQ→γγ]∣∣2 (5) = 18π∣∣∣c0⟨0|χ†ψ|ηQ⟩+c2m2⟨0|χ†(−i2\tensorD)2ψ|ηQ⟩+⋯∣∣∣2,

where the formula has been used.

Comparing (5) with (1), one can express the short-distance coefficients and that appear in (1), the standard NRQCD factorization formula for the decay rate, in terms of and :

 F(1S0) = m28π|c0|2, (6a) G(1S0) = m24πRe[c0c∗2]. (6b)

The goal of this work is to determine through the order . Therefore, we first need determine both of the coefficients and to order .

### ii.2 Matching at the amplitude level

To determine the values of and , we follow the moral that these short-distance coefficients, encapsulating the hard quantum fluctuations that occur at the length scale , are insensitive to the long-distance hadronic dynamics. As a convenient calculational device, one can replace the physical meson by a free pair of the quantum number , so that both the full amplitude, , and the NRQCD operator matrix elements can be accessed by perturbation theory. The short-distance coefficients then can be solved by equating the QCD amplitude and the corresponding NRQCD amplitude , order by order in . This procedure is commonly referred to as perturbative matching. Analogous to (3), one can write down the pertubative matching formula for the process:

 A[Q¯¯¯¯Q(1S[1]0)→γγ]=^k1⋅\boldmathε∗1×% \boldmathε∗2ANRQCD, (7a) ANRQCD=c0⟨0|χ†ψ|Q¯¯¯¯Q(1S[1]0)⟩+c2m2⟨0|χ†(−i2\tensorD)2ψ|Q¯¯¯¯Q(1S[1]0)⟩+⋯. (7b)

In Eq. (7), we again adopt the nonrelativistic normalization for the and states in the computations of the full QCD amplitude and the NRQCD matrix elements.

One can organize the full amplitude in powers of the relative momentum between and , denoted by :

 A[Q¯¯¯¯Q(1S[1]0)→γγ]=^k1⋅\boldmathε∗1×% \boldmathε∗2[A0+q2m2A2+O(q4)]. (8)

For convenience, we have again factored out the kinematic invariant in the amplitude. To the desired accuracy, one can truncate the series at , with the first two Taylor coefficients denoted by and . To our purpose, we need compute both through NLO in . This will be conducted in Section IV.

The NRQCD amplitude in (7), or equivalently, the NRQCD vacuum-to- matrix elements, need also be worked out in perturbation theory through . The encountered NRQCD matrix elements at LO in are particularly simple 6:

 ⟨0|χ†ψ|Q¯¯¯¯Q(1S[1]0)⟩(0)=√2Nc, (9a) ⟨0|χ†(−i2\tensorD)2ψ|Q¯¯¯¯Q(1S[1]0)⟩(0)=√2Ncq2, (9b)

where the factor is due to the spin and color factors of the normalized state. The computation of these matrix elements to will be addressed in Section V.

## Iii Techniques about computing full QCD amplitude

In this section, we outline some necessary techniques about calculating the amplitude for .

### iii.1 Kinematic setup

Let and represent the momenta carried by and , respectively. It is customary to decompose them in the following form:

 p = 12P+q, (10a) ¯p = 12P−q, (10b)

where is the total momentum of the pair with invariant mass , is the relative momentum. Enforcing and to stay on their mass shells, one requires that and .

In the rest frame of the pair, which is our default choice, the explicit forms of all the momenta are given by

 Pμ = (2E,0), (11a) qμ = (0,q), (11b) pμ = (E,q), (11c) ¯pμ = (E,−q). (11d)

Hence the total momentum becomes purely timelike, while the relative momentum is purely spacelike, and one has and . In this frame, the momenta carried by both photons have a magnitude of .

For latter use, we define two velocity variables:

 β ≡ |q|E, (12a) v ≡ |q|m. (12b)

We will distinguish these two variables even in the nonrelativistic limit.

### iii.2 Covariant projection approach

#### Projection of spin-singlet Q¯¯¯¯Q state

We start with the quark amplitude , with the momenta of and defined in (10):

 ¯u(p)Tv(¯p)=Tr[v(¯p)¯u(p)T]. (13)

Here denotes a matrix in Dirac-color space.

To proceed, we need first to project the amplitude (13) onto the spin-singlet color-singlet state, by replacing the with a suitable projection matrix. The projector that is valid to all orders in for the spin-singlet color-singlet channel, denoted by , was first derived in (9):

 Π(1)1(p,¯p) = ∑s1,s2u(p,s1)¯v(¯p,s2)⟨12,s1;12,s2|00⟩⊗1c√Nc (14) = 18√2E2(E+m)(/p+m)(/P+2E)γ5(/¯p−m)⊗1c√Nc,

where is the unit matrix in the fundamental representation of the color group. The above spin-singlet projector is derived by assuming the nonrelativistic normalization convention for Dirac spinor. Applying this projector to (13), one obtains the amplitude for a spin/color-singlet pair annihilation decay into two photons:

 Asing[Q¯¯¯¯Q→γγ]=Tr{Π(1)1(p,¯p)T}, (15)

where the trace is understood to act on both Dirac and color spaces.

#### Projection of S-wave amplitude

In the amplitude in (15), the pair is warranted to be in the spin-singlet state, but not necessarily in the -wave orbital-angular-momentum state. To project out the -wave amplitude, one needs average the amplitude over all the directions of the relative momentum in the rest frame. Afterwards one expands the resulting -wave amplitude in powers of and truncate the series at the desired order.

A question may concern us immediately– should the relative momentum be taken as a 4-dimensional or a -dimensional vector upon -wave angular averaging? Since for our process the matching can be conducted solely at amplitude level, there is no a priori criterion to tell which treatment is superior. At first glancing, our incapability to pick up a “unique” and “correct” scheme may look troublesome, since different treatments may, conceivably, lead to different answers for the full QCD amplitude once beyond LO in expansion. Nevertheless, it is possible that both schemes are equally acceptable, provided that ultimately they yield identical short-distance coefficients. In this respect, the ambiguity about the dimensionality of may turn into a virtue, in that it can serve as a useful consistency check against our calculation. In Section IV, we will separately treat to be 4-dimensional (dubbed scheme), and -dimensional (labeled scheme) 7.

Since we are only concerned with the first-order relativistic correction, we will follow a standard shortcut to extract the -wave amplitude rather than literally perform the angular integration. In scheme, we first expand the spin-singlet amplitude in through the quadratic order, then make the following replacement:

 qμqν → q23Πμν(P)(q4scheme), (16)

where

 Πμν(P) ≡ −gμν+PμPνP2.

Subsequently we would be able to identify () as indicated in (8).

If is assumed to be a -dimensional Lorentz vector, one can follow the same step as described above, except one make the following substitution:

 qμqν → q2D−1Πμν(P)(qDscheme). (17)

We stress that in both schemes, the momenta , , and polarization vectors , are always assumed to reside only in physical spacetime dimensions.

### iii.3 γ5-prescription in Dimensional Regularization

When applying the covariant spin-singlet projector as given in Sec. III.2.1, one needs to deal with the trace of with a string of Dirac -matrices. This will pose one notorious problem, that a definite prescription of must be specified if the spacetime dimension is deformed from four.

We first point out that, for our process, in principle there exists no any technical subtlety about in the scheme. One can always choose to first calculate the quark amplitude through NLO in . As usual, dimensional regularization (DR) can be chosen to regularize both ultraviolet (UV) and infrared (IR) divergences. After the loop integration is done, one will end up with the -matrix in (13) that only depends on the external kinematic variables, e.g. the momenta of quarks and photons, as well as the polarization vectors of photons, which are all 4-dimensional objects. Upon projecting out the spin-singlet amplitude, it is obviously legitimate to use the standard 4-dimensional trace formula involving in (15).

In the scheme, the problem cannot be circumvented even if one first carries out the loop integration for the quark amplitude, because the quark momenta and appearing in the spin-singlet projector (14) can now penetrate into unphysical dimensions. In this case, the rule about -dimensional trace operation involving in (15) must be specified.

In practical computation, it is simpler to apply (15) prior to carrying out the loop integration. Since the internal fermion propagators and vertices are all -dimensional objects, and all of them will enter into the trace, it becomes compulsory to specify the prescription of in DR, for both and schemes.

In literature, there are two popular prescriptions about in DR, the naive dimensional regularization (NDR) (25); (26) and ’t Hooft-Veltman dimensional regularization (HVDR) (27); (28). In the former prescription, one assumes for all . In the latter, one explicitly constructs , which anticommutes with for but commutes with for .

In an arbitrary dimension, the definition and are incompatible (28), therefore in NDR is an ambiguous object. In order to obtain consistent predictions in this scheme, one must impose some additional rules, e.g., to give up the cyclicity property of trace and to place in a fixed position called “reading point” (26). Nevertheless, for its technical simplicity, the NDR scheme has been widely utilized in computing the NLO QCD corrections to quarkonium decay and production processes (29), though most of which are at the LO accuracy in only 8.

In contrast to NDR, the HVDR scheme turns to be a mathematically consistent scheme, in which is a well-defined and unique object. For example, the HVDR scheme can automatically guarantee to recover the celebrated Adler-Bell-Jackiw chiral anomaly, while in the NDR scheme, some ad hoc prescription has to be imposed to achieve this.

In HVDR, the -matrices in dimension obey the following anticommutation algebra:

 {γμ,γν} = 2gμν,{γ5,γμ}=2^gμνγ5γν. (18)

denotes the projection of the metric tensor onto the unphysical dimensions, which equals for , and equals otherwise. Some useful relations about this tensor are , . A nuisance of the HVDR scheme is that, since the first four dimensions are singled out as special, Lorentz covariance has been sacrificed in the intermediate stage. Moreover, one may get the impression that the messy anticommutation rule for in dimension would render the practical calculation a formidable task.

In this work, in favor of its internal consistency, we choose to work with the HVDR scheme. It is necessary to spell out the recipe of -dimensional trace operation involving in this scheme. In arbitrary spacetime dimension, the trace of a with odd number of -matrices always vanishes. Starting from (18), West has derived a recursive formula for the trace of with an even number of -matrices (32):

 Tr[γ5γμ1γμ2γμ3γμ4]=−4iϵμ1μ2μ3μ4, (19a) Tr[γ5γμ1γμ2⋯γμn]=2n−4n∑i=2i−1∑j=1(−1)i+j+1gμiμjTr⎡⎣γ5n∏k=1(≠i,j)γμk⎤⎦(forn≥6),

where the Levi-Civita tensor is a 4-dimensional object in HVDR, i.e., .

One attractive point of West’s trace formula is that, it involves only the 4-dimensional antisymmetric tensor together with the -dimensional metric tensor (not the evanescent metric tensor !). As a consequence, this recursive algorithm can be readily implemented in the Mathematica packages specialized to high energy physics, such as FeynCalc (33).

Equation (19) is valid in any dimension, of course also in , however it is superficially much more involved than the familiar 4-dimensional trace formula 9.

To make use of West’s formula when projecting out the spin-singlet amplitude, we can employ the cyclicity of trace to move in (15) to the leftmost, since this property persists to be a valid operation in the HVDR scheme:

 Asing(Q¯¯¯¯Q→γγ)=18√2NcE2(E+m)Tr{γ5(/¯p−m)T(/p+m)(/P+2E)}, (20)

and the color trace has been implicit.

For obvious reason, we would not bother to use the anti-commutation relation of to realize this goal.

## Iv QCD amplitude of Q¯¯¯¯Q(1S[1]0)→γγ through NLO in \boldmathαs

In this section we employ the covariant projection technique described in the preceding section to compute through NLO in . The HVDR scheme will be used throughout.

### iv.1 Tree-level amplitude and matching coefficients

There are two diagrams for , one of which is illustrated in Fig. 1). The corresponding LO -matrix reads:

 T(0) = −ie2e2Q[⧸ε∗2⧸p−⧸k1+m−2p⋅k1⧸ε∗1+⧸ε∗1−⧸¯p+⧸k1+m−2¯p⋅k1⧸ε∗2]⊗1c. (21)

If lives only in physical dimensions, one can directly substitute (21) into (15) to project out the spin-singlet amplitude, and use the 4-dimensional trace formula to obtain (9):

 Asing(0)(Q¯¯¯¯Q→γγ) = e2e2Q√2Nc^k1⋅% \boldmathε∗1×\boldmathε∗2mE2E4−(k1⋅q)2. (22)

If is instead allowed to leak into the unphysical dimensions, one needs substitute (21) into (20), and utilize (19) to carry out the -dimensional trace:

 Asing(0)(Q¯¯¯¯Q→γγ)=e2e2Q√Nc21(E+m)(E4−(k1⋅q)2) (23) × [E(E+m)ϵμναβk1α+ϵμαβγk1αqγqν−ϵναβγk1αqγqμ+ϵμναβqαk1⋅q]Pβε∗1με∗2ν.

We have dropped those terms linear in , which trivially vanish due to the constraints and . Note this expression is much more complicated than (22), though they should be exactly identical when is a 4-dimensional vector.

The Levi-Civita tensor, , and are all 4-dimensional quantities. For any vector in (23) that contracts with them, only its first four components can contribute, that is, one can make the replacement , without affecting the answer. The unphysical components of , , contribute to (23) only implicitly through the factor .

Expand (22) and (23) in , then apply (16) and (17), one can pick up the corresponding tree-level -wave amplitudes in the and schemes. It is then straightforward to identify as introduced in (8):

 A(0)0=√2Nc4πe2Qαm, (24a) A(0)2∣∣∣q4scheme=−√2Nc8πe2Qα3m, (24b) A(0)2∣∣∣qDscheme=−√2Nc2Dπe2Qα(D−1)m. (24c)

Not surprisingly, the -wave amplitudes differ in and schemes.

The NRQCD matarix elements at LO in have been given in (9). According to (7b), one can write down the Born-order perturbative NRQCD amplitude:

 A(0)NRQCD=√2Nc[c(0)0+c(0)2v2+⋯]. (25)

Combining (7), (8), (24) and (25), one easily recognizes the short-distance coefficients ():

 c(0)0=4πe2Qαm, (26a) c(0)2∣∣∣q4scheme=−8πe2Qα3m, (26b) c(0)2∣∣∣qDscheme=−2−ϵ3−2ϵ4πe2Qαm, (26c)

where as usual, . The coefficient in scheme differs from that in scheme by an constant. At first sight, it seems no need to retain this extra piece since no any divergence emerges at this order. However, as will become clear later, this piece plays a key role for ultimately obtaining the scheme-independent short-distance coefficients.

### iv.2 QCD amplitude at NLO in \boldmathαs

We proceed to compute NLO QCD correction to the process. At , there are eight one-loop diagrams, including two self-energy diagrams, four triangle diagrams and two box diagrams. Half of these diagrams have been shown in Fig. 1) through 1).

Each individual NLO diagram may contain UV or IR divergences. For example, the self-energy and triangle diagrams contain UV divergences, and the box diagrams possess IR divergence. We will choose DR as a convenient regulator to regularize both types of divergences. Since our concern is to calculate the gauge-invariant on-shell amplitude, for simplicity we will work with Feynman gauge in this Section.

In accordance with the LSZ reduction formula, we need multiply the tree-level amplitude in Fig. 1) by the residue of the heavy quark propagator at its pole, . This contribution is represented by Fig. 1). In Feynman gauge, the residue is given by

 ZQ = 1−CFαs4π(1ϵUV+2ϵIR−3γE+3ln4πμ2m2+4)+O(α2s), (27)

where is the Euler constant, and is the Casmir for the fundamental representation of the group.

In addition, we also need replace the bare quark mass in the quark propagator in Fig. 1) by

 mbare = m[1−CFαs4π(3ϵUV−3γE+3ln4πμ2m2+4)+O(α2s)]. (28)

It is straightforward to write down the respective -matrix for each NLO diagram in Fig. 1. We then substitute them into (20), and for simplicity, use (19) to carry out the -dimensional trace prior to performing loop integration 10.

Practically, we resort to the Mathematica package FeynCalc (33) to accomplish the abovementioned trace operation, because West’s formula (19) is its built-in algorithm for calculating the trace involving . We continue to use FeynCalc to reduce all the encountered one-loop tensor integrals to the one-loop scalar integrals. It turns out that only a couple of two-point, three-point scalar integrals and one four-point scalar integral are required. For reader’s convenience, the closed-form expressions for these scalar integrals have been collected in Appendix A.

By far we have obtained the analytic expression for , the spin-singlet amplitude for . We proceed to expand it to second order in , then apply (16) and (17) to extract the corresponding -wave amplitudes defined in (8), in both and schemes, respectively. The intermediate steps are straightforward but cumbersome, and some special care has to be paid when dealing with the box diagrams. Fortunately, almost all these manipulations can be handled by computer.