aainstitutetext: HEP Division, Argonne National Laboratory,
9700 South Cass Avenue, Argonne, IL 60439, U.S.A.bbinstitutetext: Department
of Physics, Korea University,
Seoul 136-713, Korea.

# Higher-order relativistic corrections to gluon fragmentation into spin-
triplet $\bm{S}$-wave quarkonium

Geoffrey T. Bodwin b U-Rae Kim b and Jungil Lee<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We compute the relative-order-$v^{4}$ contribution to gluon fragmentation into
quarkonium in the ${}^{3}S_{1}$ color-singlet channel, using the
nonrelativistic QCD (NRQCD) factorization approach. The QCD fragmentation
process contains infrared divergences that produce single and double poles in
$\epsilon$ in $4-2\epsilon$ dimensions. We devise subtractions that isolate
the pole contributions, which ultimately are absorbed into long-distance NRQCD
matrix elements in the NRQCD matching procedure. The matching procedure
involves two-loop renormalizations of the NRQCD operators. The subtractions
are integrated over the phase space analytically in $4-2\epsilon$ dimensions,
and the remainder is integrated over the phase-space numerically. We find that
the order-$v^{4}$ contribution is enhanced relative to the order-$v^{0}$
contribution. However, the order-$v^{4}$ contribution is not important
numerically at the current level of precision of quarkonium-hadroproduction
phenomenology. We also estimate the contribution to hadroproduction from gluon
fragmentation into quarkonium in the ${}^{3}P_{J}$ color-octet channel and
find that it is significant in comparison to the complete next-to-leading-
order-in-$\alpha_{s}$ contribution in that channel.

###### Keywords:

quarkonium, fragmentation, NRQCD, relativistic corrections

††arxiv: 1208.5301††preprint: ANL-HEP-PR-12-65

##  1 Introduction

Heavy-quarkonium production in hard-scattering collisions has a long and rich
history of experimental measurements and theoretical calculations
Brambilla:2010cs . Intense efforts in this area are expected to continue as
the Large Hadron Collider (LHC) makes available data with unprecedented
momentum transfers and statistics.

In recent years, a great deal of theoretical effort has been focused on the
nonrelativistic QCD (NRQCD) factorization approach BBL  to calculations of
quarkonium production rates. In this approach, it is conjectured that the
inclusive quarkonium production cross section at large transverse momentum
(${{p}}_{T}$) can be written as a sum of products of short-distance
coefficients and long-distance matrix elements (LDMEs):

$\sigma(H)=\sum_{n}F_{n}(\mu_{\Lambda})\langle 0|{\cal
O}_{n}^{H}(\mu_{\Lambda})|0\rangle.$ (1)

Here, $\mu_{\Lambda}$ is the factorization scale, which is the cutoff of the
effective field theory NRQCD. A short-distance coefficient
$F_{n}(\mu_{\Lambda})$ is, in essence, the partonic cross section to produce a
heavy-quark-antiquark ($Q\bar{Q}$) pair with certain quantum numbers,
convolved with parton distribution functions. The short-distance coefficients
can be calculated as perturbation series in the strong-coupling constant
$\alpha_{s}$. A production LDME ${\langle 0|\cal
O}_{n}^{H}(\mu_{\Lambda})|0\rangle$ is the probability for a $Q\bar{Q}$ pair
with certain quantum numbers to evolve into a particular heavy-quarkonium
state. It is expressed as the vacuum expectation value of a four-fermion
operator

${\cal O}_{n}^{H}(\mu_{\Lambda})=\langle 0|\chi^{\dagger}\kappa_{n}\psi{\cal
P}_{H(P)}\psi^{\dagger}\kappa^{\prime}_{n}\chi|0\rangle,$ (2)

where $\psi^{\dagger}$ and $\chi$ are two-component (Pauli) fields that create
a heavy quark and a heavy antiquark, respectively, and $\kappa_{n}$ and
$\kappa_{n}^{\prime}$ are combinations of Pauli and color matrices.111 It was
pointed out by Nayak, Qiu, and Sterman that gauge invariance requires that the
definitions of the NRQCD LDMEs include Wilson lines that run from the quark
and antiquark fields to infinity Nayak:2005rw ; Nayak:2005rt . For simplicity,
we have omitted these Wilson lines here.

${\cal P}_{H(P)}=\sum_{X}|H(P)+X\rangle\langle H(P)+X|$ (3)

is a projection onto a state consisting of a quarkonium $H$, with four-
momentum $P$, plus anything. ${\cal P}_{H(P)}$ contains a sum over any
quarkonium polarization quantum numbers that are not specified explicitly. The
NRQCD LDMEs are evaluated in the rest frame of the quarkonium, in which
$P=(M,\bm{0})$, where $M$ is the quarkonium mass. In the remainder of this
paper, we suppress the momentum argument of ${\cal P}_{H(P)}$ in NRQCD LDMEs.

The production LDMEs for the evolution of color-singlet $Q\bar{Q}$ pairs into
a quarkonium state are related to the color-singlet quarkonium decay LDMEs.
These color-singlet production LDMEs can be determined from comparison of
theory with quarkonium production or decay data or from lattice QCD
calculations. However, the production LDMEs for the evolution of color-octet
$Q\bar{Q}$ pairs into a quarkonium state can be determined, at least at
present, only through comparison of theory with experimental quarkonium-
production data.

Complete calculations of short-distance coefficients in the NRQCD
factorization approach now exist through next-to-leading order (NLO) in
$\alpha_{s}$ for production of the $J/\psi$ and the $\psi(2S)$ in $e^{+}e^{-}$
collisions, in $ep$ collisions, and in $p\bar{p}$ and $pp$ collisions
Brambilla:2010cs ; Butenschoen:2009zy ; Ma:2010yw ; Butenschoen:2010rq ;
Butenschoen:2010px ; Ma:2010jj ; Butenschoen:2011yh ; Butenschoen:2011ks ;
Butenschoen:2012px ; Chao:2012iv ; Butenschoen:2012qh ; Gong:2012ug . These
calculations include the contributions from all of the color-octet channels
through relative order $v^{4}$, as well as the contribution of the color-
singlet channel at leading order (LO) in $v$. Specifically, the calculations
include the contributions of the ${}^{3}S_{1}$, ${}^{1}S_{0}$ and
${}^{3}P_{J}$ color-octet channels and the ${}^{3}S_{1}$ color-singlet channel
at the leading nontrivial order in $v$ in each channel. Here, $v$ is half the
relative velocity of the heavy quark and the heavy antiquark in the quarkonium
rest frame. $v^{2}\approx 0.22$ for the $J/\psi$, and $v^{2}\approx 0.1$ for
the $\Upsilon$. These theoretical results are generally compatible with
experimental measurements of quarkonium production cross sections. However,
significant discrepancies remain between theoretical predictions for
quarkonium polarization and experimental measurements Chao:2012iv ;
Butenschoen:2012qh ; Gong:2012ug . These discrepancies might point to as-yet-
uncalculated theoretical contributions, to experimental difficulties, to a
failure of convergence of the NRQCD series in $\alpha_{s}$ or $v$, or to a
failure of the NRQCD factorization conjecture itself.

The calculations at NLO in $\alpha_{s}$ have revealed very large corrections
in that order to the ${}^{3}P_{J}$ color-octet channel and the ${}^{3}S_{1}$
color-singlet channel. Large corrections have also been found in a calculation
of the real-emission contributions to $\Upsilon$ hadroproduction at next-to-
next-to-leading order in $\alpha_{s}$ Artoisenet:2008fc . The large
corrections are the result of kinematic enhancements of the higher-order cross
sections at high ${{p}}_{T}$ relative to the LO cross sections. The sizes of
these corrections have cast some doubt on the convergence of the perturbation
series.

It has been suggested recently that the large higher-order corrections to
quarkonium production can be brought under control by re-organizing the
perturbation series according to the ${{p}}_{T}$ behavior of the various
contributions Kang:2011zza . In this approach, the cross section can be shown
to factorize into convolutions of hard-scattering cross sections with
fragmentation functions. The factorization holds up to corrections of relative
order $m_{c}^{4}/{{p}}_{T}^{4}$, where $m_{c}$ is the charm-quark mass. The
factorized cross section consists of a leading contribution, which arises from
single-particle fragmentation into a quarkonium and falls as $1/{{p}}_{T}^{4}$
in the partonic cross section and a first subleading contribution, which
arises from two-particle fragmentation into a quarkonium and falls as
$1/{{p}}_{T}^{6}$ in the partonic cross section. If NRQCD factorization holds,
then the various fragmentation functions can be expressed in terms of a sum of
products of short-distance coefficients and NRQCD LDMEs. This picture has been
shown to account for the large corrections at NLO in $\alpha_{s}$ in the
${}^{3}S_{1}$ color-singlet channel Kang:2011mg .

In this paper, we compute the NRQCD short-distance coefficient for gluon
fragmentation into a ${}^{3}S_{1}$ color-singlet $Q\bar{Q}$ pair in relative
order $v^{4}$. The short-distance coefficients for gluon fragmentation in this
channel have been computed in relative order $v^{0}$ Braaten:1993rw ;
Braaten:1995cj  and relative order $v^{2}$ Bodwin:2003wh . In both cases, the
contributions are not important phenomenologically. Nevertheless, it is
worthwhile to consider the order-$v^{4}$ contribution for two reasons. First,
this contribution is interesting theoretically because it is in order $v^{4}$
that the ${}^{3}S_{1}$ color-singlet fragmentation channel first develops soft
divergences in full QCD. As we shall see, these soft divergences in full QCD
correspond to soft divergences in the LDMEs for the ${}^{3}S_{1}$ and
${}^{3}P_{J}$ color-octet channels and cancel in the short-distance
coefficients, as is required by NRQCD factorization. A second motivation for
examining the order-$v^{4}$ contribution is that it is potentially large.
Contributions from gluon fragmentation into the ${}^{3}S_{1}$ and
${}^{3}P_{J}$ color-octet channels are known to be important
phenomenologically. The ${}^{3}S_{1}$ color-singlet channel mixes with these
channels in order $v^{4}$, and the partitioning of the various contributions
is controlled by single and double logarithms of the factorization scale.
Therefore, it is plausible that the order-$v^{4}$ contributions to the color-
singlet channel could be large.

Our method of calculation is based on the Collins-Soper definition
Collins:1981uw  of the fragmentation function for a gluon fragmenting into a
quarkonium. We assume that NRQCD factorization holds, that is, that the
fragmentation function can be decomposed into a sum of products of short-
distance coefficients and NRQCD LDMEs. We then compute the full-QCD
fragmentation functions for a gluon fragmenting into free $Q\bar{Q}$ states
with various quantum numbers. The ultimate aim is to match these full-QCD
fragmentation functions to the corresponding NRQCD fragmentation functions in
order to determine the NRQCD short-distance coefficients. Some additional
details of this approach can be found in Ref. Bodwin:2003wh .

This paper is organized as follows. We give the Collins-Soper definition of
the fragmentation function in Sec. 2. Section 3 contains the NRQCD
factorization formula for the fragmentation function and also contains a
discussion of the NRQCD LDMEs and short-distance coefficients that are
relevant through relative order $v^{4}$. The kinematics and variables that we
use in our calculation are described in Sec. 4. In Sec. 5, we discuss the
calculation of the fragmentation processes in full QCD. As we have mentioned,
an important feature of the present calculation is that soft divergences arise
in the ${}^{3}S_{1}$ color-singlet channel in both full QCD and NRQCD. These
divergences ultimately cancel in the short-distance coefficients when we carry
out the matching between full QCD and NRQCD. In both full QCD and NRQCD, we
regulate the divergences dimensionally. In the case of full QCD, we devise
subtractions that remove the divergent terms from the integrand, and we
compute the subtraction contributions analytically. This computation is
described in Sec. 5.3. After we remove the subtraction terms, we calculate the
remainder of the full-QCD contribution in four dimensions, carrying out the
integration numerically. We compute the relevant NRQCD LDMEs for free
$Q\bar{Q}$ states analytically in dimensional regularization. These
calculations are described in Sec. 6. We also determine the evolution
equations for the LDMEs and find a discrepancy with a result in Ref.
Gremm:1997dq . In Sec. 7, we match the NRQCD and full-QCD fragmentation
functions to obtain the short-distance coefficients, and we present numerical
results for them in Sec. 8. Finally, in Sec. 9, we summarize our results.

##  2 Collins-Soper definition of the fragmentation function

Here, and throughout this paper, we use the following light-cone coordinates
for a four-vector $V$:

$\displaystyle V$ $\displaystyle=$
$\displaystyle(V^{+},V^{-},\bm{V}_{\bot})=(V^{+},V^{-},V^{1},V^{2}),$ (4a)
$\displaystyle V^{+}$ $\displaystyle=$ $\displaystyle(V^{0}+V^{3})/\sqrt{2},$
(4b) $\displaystyle V^{-}$ $\displaystyle=$
$\displaystyle(V^{0}-V^{3})/\sqrt{2}.$ (4c)

The scalar product of two four-vectors $V$ and $W$ is then

$V\cdot W=V^{+}W^{-}+V^{-}W^{+}-\bm{V}_{\bot}\cdot\bm{W}_{\bot}.$ (5)

The Collins-Soper definition for the fragmentation function for a gluon
fragmenting into a hadron (quarkonium) $H$ Collins:1981uw  is

$\displaystyle D[g\to H](z,\mu_{\Lambda})$ $\displaystyle=$
$\displaystyle\frac{-g_{\mu\nu}z^{d-3}}{2\pi
k^{+}(N_{c}^{2}-1)(d-2)}\int_{-\infty}^{+\infty}dx^{-}e^{-ik^{+}x^{-}}$
$\displaystyle\times\langle
0|G^{+\mu}_{c}(0)\mathcal{E}^{\dagger}(0,0,\bm{0}_{\perp})_{cb}\;\mathcal{P}_{H(P)}\;\mathcal{E}(0,x^{-},\bm{0}_{\perp})_{ba}G^{+\nu}_{a}(0,x^{-},\bm{0}_{\perp})|0\rangle\,.$

Here, $z$ is the fraction of the gluon’s $+$ component of momentum that is
carried by the hadron, $G_{\mu\nu}$ is the gluon field-strength operator, $k$
is the momentum of the field-strength operator, $\mu_{\Lambda}$ is the
factorization scale, and $d=4-2\epsilon$ is the number of space-time
dimensions. There is an implicit average over the color and polarization
states of the initial gluon. The projection $\mathcal{P}_{H(P)}$ is given in
Eq. (3). The fragmentation function is evaluated in the frame in which the
hadron has zero transverse momentum:
$P=[zk^{+},M^{2}/(2zk^{+}),\bm{0}_{\perp}]$. The operator
$\mathcal{E}(0,x^{-},\bm{0}_{\perp})$ is a path-ordered exponential of the
gluon field:

$\mathcal{E}(0,x^{-},\bm{0}_{\perp})_{ba}\;=\;\textrm{P}\exp\left[+ig_{s}\int_{x^{-}}^{\infty}dz^{-}A^{+}(0,z^{-},\bm{0}_{\perp})\right]_{ba},$
(7)

where $g_{s}=\sqrt{4\pi\alpha_{s}}$ is the QCD coupling constant and
$A^{\mu}(x)$ is the gluon field. Both $A_{\mu}$ and $G_{\mu\nu}$ are SU(3)
matrices in the adjoint representation. The expression (LABEL:eq:D-def) is
manifestly gauge invariant. We use the Feynman gauge in our calculation.

The Feynman rules for the perturbative expansion of Eq. (LABEL:eq:D-def) are
given in Ref. Collins:1981uw . The quantity
$\mathcal{E}(0,x^{-},\bm{0}_{\perp})$ appears in the Feynman rules as an
eikonal line. Owing to the charge-conjugation properties of the $Q\bar{Q}$
states that we consider and the Landau-Yang theorem landau-thm ; Yang:1950rg ,
gluon attachments to the eikonal lines from
$\mathcal{E}(0,x^{-},\bm{0}_{\perp})$ do not appear in our calculation. Hence,
we need only the standard QCD Feynman rules, an overall factor

$C_{\rm frag}=\frac{z^{d-3}k^{+}}{2\pi(N_{c}^{2}-1)(d-2)}$ (8)

from Eq. (LABEL:eq:D-def), and the special Feynman rule for the vertex that
creates a gluon and an eikonal line. That vertex is shown in Fig. 1. Its
Feynman rule, in momentum space, is a factor

$+i\left(g^{\nu\alpha}-\frac{Q^{\nu}n^{\alpha}}{k^{+}}\right)\delta_{ab},$ (9)

where $k$ is the sum of the momenta of the gluon and the eikonal line, $Q$ is
the momentum of the gluon, $\alpha$ is the polarization index of the gluon,
and $a$ and $b$ are the color indices, respectively, of the gluon and the
eikonal line. In the absence of interactions with the eikonal lines, $k=Q$.
$n$ is a light-like vector whose components are given by
$n=(0,1,\bm{0}_{\perp})$.

Figure 1: Feynman diagram for the vertex that creates a gluon and an eikonal
line. The circle represents the operator $G_{a}^{+\nu}$, which creates a gluon
with momentum $Q$ and polarization and color indices $\alpha$ and $a$,
respectively. $b$ is the color index for the eikonal line. The operator
momentum $k$ is the sum of $Q$ and the momentum of the eikonal line.

The final-state phase space that is implied by Eq. (LABEL:eq:D-def) is

$d\Phi_{n}=\frac{4\pi
M}{S}\;\delta\left(k^{+}-P^{+}-\sum_{i=1}^{n}k_{i}^{+}\right)\theta(k^{+})\prod_{i=1}^{n}\frac{dk_{i}^{+}}{4\pi
k^{+}_{i}}\,\frac{d^{d-2}\bm{k}_{i\perp}}{(2\pi)^{d-2}}\,\theta(k_{i}^{+}),$
(10)

where $S$ is the statistical factor for identical particles in the final
state, $k_{i}$ is the momentum of the $i$th final-state particle, and the
product is over all of the final-state particles except $H$. We use
nonrelativistic normalization for the state $H$, and so a factor $2M$ appears
in the phase space in order to cancel the relativistic normalization of $H$ in
the definition (LABEL:eq:D-def). We use relativistic normalization for all
particles other than $H$.

##  3 NRQCD factorization

We assume that the fragmentation function for a gluon fragmenting into a
quarkonium $H$ satisfies NRQCD factorization. Then, in analogy with Eq. (1),
we have

$D[g\to H](z)=\sum_{n}d_{n}(z)\langle 0|\mathcal{O}_{n}^{H}|0\rangle,$ (11)

where the $\langle 0|\mathcal{O}_{n}^{H}|0\rangle$ are NRQCD LDMEs and the
$d_{n}(z)$ are the fragmentation short-distance coefficients. We have
suppressed the dependences of $d_{n}(z)$ and $\langle
0|\mathcal{O}_{n}^{H}|0\rangle$ on the factorization scale $\mu_{\Lambda}$. In
discussing specific cases, we use the notation $\langle
0|\mathcal{O}_{n_{v},n^{\prime}}^{H}({}^{2s+1}l_{j}^{[c]})|0\rangle$ for the
LDMEs and the notation $d_{n_{v},n^{\prime}}({}^{2s+1}l_{j}^{[c]})(z)$ for the
short-distance coefficients, where ${}^{2s+1}l_{j}$ is the standard
spectroscopic notation for the angular-momentum quantum numbers of the
corresponding NRQCD operator, and $c$ is the color quantum number of the NRQCD
operator ($1$ or $8$). $n_{v}$ is the order in $v$, relative to the leading
order, of the field operators and derivatives in
$\mathcal{O}_{n_{v},n^{\prime}}^{H}$, excluding factors of $v$ from the
projection onto the final state $H$. $n^{\prime}$ is an integer that is used
to distinguish operators that have the same quantum numbers and order in $v$.
We denote the contribution of order $v^{k}$ to $D[g\to H]$ by $D_{k}[g\to H]$.

We can write the fragmentation functions for gluon fragmentation into free
$Q\bar{Q}$ states as

$D[g\to Q\bar{Q}](z)=\sum_{n}d_{n}(z)\langle
0|\mathcal{O}_{n}^{Q\bar{Q}}|0\rangle.$ (12)

Since the short-distance coefficients $d_{n}(z)$ are independent of the
specifics of the hadronic states, the $d_{n}(z)$ in Eq. (12) are identical to
the $d_{n}(z)$ in Eq. (11). We determine the $d_{n}(z)$ by computing the left
side of Eq. (12) in full QCD and comparing it with the right side, in which
the free $Q\bar{Q}$ LDMEs are computed in NRQCD. Since we choose a
factorization scale $\mu_{\Lambda}$ of order the heavy-quark mass $m$, we can
carry out this computation in perturbation theory. We will denote the
contribution of order $v^{k}$ to $D[g\to Q\bar{Q}]$ by $D_{k}[g\to Q\bar{Q}]$.

If $H$ is a ${}^{3}S_{1}$ quarkonium state, such as the $J/\psi$, then, in LO
in $v$, we must consider the LDME

$\langle 0|\mathcal{O}_{0}^{H}({}^{3}S_{1}^{[1]})|0\rangle=\langle
0|\chi^{\dagger}\sigma^{i}\psi\;{\cal
P}_{H}\;\psi^{\dagger}\sigma^{i}\chi|0\rangle.$ (13)

In relative order $v^{2}$, we must consider the LDME

$\langle 0|\mathcal{O}_{2}^{H}({}^{3}S_{1}^{[1]})|0\rangle=\frac{1}{2}\langle
0|\chi^{\dagger}\sigma^{i}(-\tfrac{i}{2}\\!\\!\stackrel{{\scriptstyle\leftrightarrow}}{{\bm{D}}})^{2}\psi\;{\cal
P}_{H}\;\psi^{\dagger}\sigma^{i}\chi+\textrm{H.~{}c.}|0\rangle.$ (14)

In relative order $v^{3}$, we must consider the LDME

$\langle 0|{\cal O}_{0}^{H}({}^{1}S_{0}^{[8]})|0\rangle=\langle
0|\chi^{\dagger}T^{a}\psi\;{\cal P}_{H}\;\psi^{\dagger}T^{a}\chi|0\rangle.$
(15)

In relative order $v^{4}$, we must consider the LDMEs

$\displaystyle\langle 0|{\cal O}_{4,1}^{H}({}^{3}S_{1}^{[1]})|0\rangle$
$\displaystyle=$ $\displaystyle\langle
0|\chi^{\dagger}\sigma^{i}(-\tfrac{i}{2}\\!\\!\stackrel{{\scriptstyle\leftrightarrow}}{{\bm{D}}})^{2}\psi\;{\cal
P}_{H}\;\psi^{\dagger}\sigma^{i}(-\tfrac{i}{2}\\!\\!\stackrel{{\scriptstyle\leftrightarrow}}{{\bm{D}}})^{2}\chi|0\rangle,$
(16a) $\displaystyle\langle
0|\mathcal{O}_{4,2}^{H}({}^{3}S_{1}^{[1]})|0\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{2}\langle
0|\chi^{\dagger}\sigma^{i}(-\tfrac{i}{2}\\!\\!\stackrel{{\scriptstyle\leftrightarrow}}{{\bm{D}}})^{4}\psi\;{\cal
P}_{H}\;\psi^{\dagger}\sigma^{i}\chi+\textrm{H.~{}c.}|0\rangle,$ (16b)
$\displaystyle\langle 0|\mathcal{O}_{4,3}^{H}({}^{3}S_{1}^{[1]})|0\rangle$
$\displaystyle=$ $\displaystyle\frac{1}{2}\langle
0|\chi^{\dagger}\sigma^{i}\psi\;{\cal
P}_{H}\;\psi^{\dagger}\sigma^{i}(\stackrel{{\scriptstyle\leftrightarrow}}{{\bm{D}}}\\!\cdot
g_{s}\bm{E}+g_{s}\bm{E}\cdot\\!\stackrel{{\scriptstyle\leftrightarrow}}{{\bm{D}}})\chi$
(16c)
$\displaystyle\qquad-\chi^{\dagger}\sigma^{i}(\stackrel{{\scriptstyle\leftrightarrow}}{{\bm{D}}}\\!\cdot
g_{s}\bm{E}+g_{s}\bm{E}\cdot\\!\stackrel{{\scriptstyle\leftrightarrow}}{{\bm{D}}})\psi\;{\cal
P}_{H}\;\psi^{\dagger}\sigma^{i}\chi|0\rangle,$ $\displaystyle\langle 0|{\cal
O}_{0}^{H}({}^{3}S_{1}^{[8]})|0\rangle$ $\displaystyle=$ $\displaystyle\langle
0|\chi^{\dagger}\sigma^{i}T^{a}\psi\;{\cal
P}_{H}\;\psi^{\dagger}\sigma^{i}T^{a}\chi|0\rangle,$ (16d)
$\displaystyle\langle 0|\mathcal{O}^{H}_{0}({}^{3}P_{0}^{[1]})|0\rangle$
$\displaystyle=$ $\displaystyle\frac{1}{d-1}\langle
0|\chi^{\dagger}(-\tfrac{i}{2}\\!\\!\stackrel{{\scriptstyle\leftrightarrow}}{{\bm{D}}}\cdot\bm{\sigma})\psi\;{\cal
P}_{H}\;\psi^{\dagger}(-\tfrac{i}{2}\\!\\!\stackrel{{\scriptstyle\leftrightarrow}}{{\bm{D}}}\cdot\bm{\sigma})\chi|0\rangle,\phantom{xxx}$
(16e) $\displaystyle\langle 0|\mathcal{O}^{H}_{0}({}^{3}P_{1}^{[1]})|0\rangle$
$\displaystyle=$ $\displaystyle\langle
0|\chi^{\dagger}(-\tfrac{i}{2}\\!\\!\stackrel{{\scriptstyle\leftrightarrow}}{{D}}\phantom{}\\!\\!^{[i}\sigma^{j]})\psi\;{\cal
P}_{H}\;\psi^{\dagger}(-\tfrac{i}{2}\\!\\!\stackrel{{\scriptstyle\leftrightarrow}}{{D}}\phantom{}\\!\\!^{[i}\sigma^{j]})\chi|0\rangle,\phantom{xxx}$
(16f) $\displaystyle\langle 0|\mathcal{O}^{H}_{0}({}^{3}P_{2}^{[1]})|0\rangle$
$\displaystyle=$ $\displaystyle\langle
0|\chi^{\dagger}(-\tfrac{i}{2}\\!\\!\stackrel{{\scriptstyle\leftrightarrow}}{{D}}\\!\\!\phantom{}^{(i}\sigma^{j)})\psi\;{\cal
P}_{H}\;\psi^{\dagger}(-\tfrac{i}{2}\\!\\!\stackrel{{\scriptstyle\leftrightarrow}}{{D}}\\!\\!\phantom{}^{(i}\sigma^{j)})\chi|0\rangle.$
(16g)

Here, the symmetric traceless product is defined by

$A^{(i}B^{j)}=\frac{1}{2}(A^{i}B^{j}+A^{j}B^{i})-\frac{1}{d-1}\,\delta^{ij}A^{k}B^{k},$
(17)

and the antisymmetric product is defined by

$A^{[i}B^{j]}=\frac{1}{2}(A^{i}B^{j}-A^{j}B^{i}).$ (18)

For purposes of our calculation, it is also useful to define

$\displaystyle\langle 0|{\cal O}_{0}^{H}({}^{3}P^{[8]})|0\rangle$
$\displaystyle=$ $\displaystyle\sum_{J=0,\,1,\,2}\langle 0|{\cal
O}_{0}^{H}({}^{3}P_{J}^{[8]})|0\rangle$ (19) $\displaystyle=$
$\displaystyle\langle
0|\chi^{\dagger}(-\tfrac{i}{2}\\!\\!\stackrel{{\scriptstyle\leftrightarrow}}{{\bm{D}}})^{r}\sigma^{n}T^{a}\psi\;{\cal
P}_{H}\;\psi^{\dagger}(-\tfrac{i}{2}\\!\\!\stackrel{{\scriptstyle\leftrightarrow}}{{\bm{D}}})^{r}\sigma^{n}T^{a}\chi|0\rangle.$

It was shown in Ref. Bodwin:2002hg  that, by making use of the NRQCD equations
of motion, one can express the LDME $\langle
0|\mathcal{O}_{4,3}^{H}({}^{3}S_{1}^{[1]})|0\rangle$ in terms of the LDMEs
$\langle 0|\mathcal{O}_{4,1}^{H}({}^{3}S_{1}^{[1]})|0\rangle$ and $\langle
0|\mathcal{O}_{4,2}^{H}({}^{3}S_{1}^{[1]})|0\rangle$. (Equivalently, one can
eliminate the LDME $\langle
0|\mathcal{O}_{4,3}^{H}({}^{3}S_{1}^{[1]})|0\rangle$ by making use of a field
redefinition Brambilla:2008zg .) Hence, we need not consider $\langle
0|\mathcal{O}_{4,3}^{H}({}^{3}S_{1}^{[1]})|0\rangle$ in our analysis.

In the LDMEs $\langle 0|\mathcal{O}_{4,1}^{H}({}^{3}S_{1}^{[1]})|0\rangle$ and
$\langle 0|\mathcal{O}_{4,2}^{H}({}^{3}S_{1}^{[1]})|0\rangle$, one can replace
$|H+X\rangle\langle H+X|$ in the projector ${\cal P}_{H(M,\bm{0})}$ with
$|H\rangle\langle H|$ (vacuum-saturation approximation), making an error of
relative order $v^{4}$. If one takes this approximation and evaluates the
LDMEs in dimensional regularization in a potential model, then they are equal
Bodwin:2006dn . Since the static potential model is valid up to corrections of
order $v^{2}$, we have

$\langle 0|\mathcal{O}_{4,1}^{H}({}^{3}S_{1}^{[1]})|0\rangle=\langle
0|\mathcal{O}_{4,2}^{H}({}^{3}S_{1}^{[1]})|0\rangle+{\cal O}(v^{2}).$ (20)

Hence, up to corrections of relative order $v^{2}$, only the sum of short-
distance coefficients
$d_{4,1}({}^{3}S_{1}^{[1]})(z)+d_{4,2}({}^{3}S_{1}^{[1]})(z)$ appears in the
fragmentation function.

###  3.1 NRQCD factorization formulas for ${g\to J/\psi}$ through order
$v^{4}$

In summary, we have the following NRQCD factorization formulas for gluon
fragmentation into $J/\psi$ through relative order $v^{4}$.

In relative order $v^{0}$ we have

$D_{0}[g\to J/\psi]=d_{0}[g\to Q\bar{Q}({}^{3}S_{1}^{[1]})]\langle
0|\mathcal{O}_{0}^{J/\psi}({}^{3}S_{1}^{[1]})|0\rangle.$ (21)

The short-distance coefficient $d_{0}[g\to Q\bar{Q}({}^{3}S_{1}^{[1]})]$ was
calculated in Refs. Braaten:1993rw ; Braaten:1995cj .

In relative order $v^{2}$ we have

$D_{2}[g\to J/\psi]=d_{2}[g\to Q\bar{Q}({}^{3}S_{1}^{[1]})]\langle
0|\mathcal{O}_{2}^{J/\psi}({}^{3}S^{[1]})|0\rangle.$ (22)

The short-distance coefficient $d_{2}[g\to Q\bar{Q}({}^{3}S_{1}^{[1]})]$ was
calculated in Ref. Bodwin:2003wh .

In relative order $v^{3}$ we have

$D_{3}[g\to J/\psi]=d_{0}[g\to Q\bar{Q}({}^{1}S_{0}^{[8]})]\langle
0|\mathcal{O}_{0}^{J/\psi}({}^{1}S_{0}^{[8]})|0\rangle.$ (23)

The short-distance coefficient $d_{0}[g\to Q\bar{Q}({}^{1}S_{0}^{[1]})]$ was
calculated in Ref. Braaten:1996rp  and differs from the short-distance
coefficient $d_{0}[g\to Q\bar{Q}({}^{1}S_{0}^{[8]})]$ in Eq. (23) only by a
color factor, which we provide in Sec. 7.

In relative order $v^{4}$ we have

$\displaystyle D_{4}[g\to J/\psi]$ $\displaystyle=$
$\displaystyle\big{\\{}\,d_{4,1}[g\to
Q\bar{Q}({}^{3}S_{1}^{[1]})]+d_{4,2}[g\to
Q\bar{Q}({}^{3}S_{1}^{[1]})]\,\big{\\}}\,\langle
0|\mathcal{O}_{4}^{J/\psi}({}^{3}S_{1}^{[1]})|0\rangle$ (24) $\displaystyle+$
$\displaystyle d_{0}[g\to Q\bar{Q}({}^{3}P^{[8]})]\langle
0|\mathcal{O}_{0}^{J/\psi}({}^{3}P^{[8]})|0\rangle$ $\displaystyle+$
$\displaystyle d_{0}[g\to Q\bar{Q}({}^{3}S_{1}^{[8]})]\langle
0|\mathcal{O}_{0}^{J/\psi}({}^{3}S_{1}^{[8]})|0\rangle.$

The short-distance coefficient $d_{0}[g\to Q\bar{Q}({}^{3}S_{1}^{[8]})]$ was
calculated at LO in $\alpha_{s}$ in Refs. Braaten:1996rp ; Bodwin:2003wh  and
at NLO in $\alpha_{s}$ in Refs. BL:gfrag-NLO ; Lee:2005jw . We verify the LO
calculations in Refs. Braaten:1996rp ; Bodwin:2003wh  in the present paper,
giving our result in Sec. 7. We compute $d_{0}[g\to Q\bar{Q}({}^{3}P^{[8]})]$
in this paper, giving the result in Sec. 7. The short-distance coefficient
$d_{0}[g\to Q\bar{Q}({}^{3}P^{[1]})]$ was calculated in Refs. Braaten:1996rp
and differs from the short-distance coefficient $d_{0}[g\to
Q\bar{Q}({}^{3}P^{[8]})]$ in Eq. (24) only by a color factor. The computation
of the combination of short-distance coefficients $d_{4,1}[g\to
Q\bar{Q}({}^{3}S_{1}^{[1]})]+d_{4,2}[g\to Q\bar{Q}({}^{3}S_{1}^{[1]})]$ is the
main goal of this paper. The result of that computation is given in Sec. 7.

##  4 Kinematics

In the calculations to follow, in both full QCD and NRQCD, we employ the
following kinematics.

We take the $Q$ and the $\bar{Q}$ to be free (on-shell) states with momenta

$\displaystyle p$ $\displaystyle=$ $\displaystyle\tfrac{1}{2}P+q,$ (25a)
$\displaystyle\bar{p}$ $\displaystyle=$ $\displaystyle\tfrac{1}{2}P-q,$ (25b)

respectively. The heavy quark has three-momentum $\bm{q}$ in the $Q\bar{Q}$
rest frame, and, so, the invariant mass of the $Q\bar{Q}$ state is

$P^{2}=M^{2}=4E^{2},$ (26)

where

$E=\sqrt{m^{2}+\bm{q}^{2}}.$ (27)

We work in the frame in which the transverse momentum of the $Q\bar{Q}$ pair
vanishes. In this frame, the initial-state gluon, the final-state $Q\bar{Q}$
pair and the final-state gluons, respectively, have the momenta

$\displaystyle k$ $\displaystyle=$
$\displaystyle\left(k^{+},k^{-}=\frac{k^{2}+(P_{\perp}/z)^{2}}{2k^{+}},-\frac{\bm{P}_{\perp}}{z}\right),$
(28a) $\displaystyle P$ $\displaystyle=$
$\displaystyle\left(zk^{+},\frac{M^{2}}{2zk^{+}},\bm{0}_{\perp}\right),$ (28b)
$\displaystyle k_{1}$ $\displaystyle=$
$\displaystyle\left(z_{1}k^{+},\frac{k_{1\perp}^{2}}{2z_{1}k^{+}},\bm{k}_{1\perp}\right),$
(28c) $\displaystyle k_{2}$ $\displaystyle=$
$\displaystyle\left(z_{2}k^{+},\frac{k_{2\perp}^{2}}{2z_{2}k^{+}},\bm{k}_{2\perp}\right),$
(28d)

where we have introduced the longitudinal momentum fractions

$\displaystyle z$ $\displaystyle=$ $\displaystyle\frac{P^{+}}{k^{+}},$ (29a)
$\displaystyle z_{1}$ $\displaystyle=$ $\displaystyle\frac{k_{1}^{+}}{k^{+}},$
(29b) $\displaystyle z_{2}$ $\displaystyle=$
$\displaystyle\frac{k_{2}^{+}}{k^{+}}.$ (29c)

Because of the conservation of four-momentum, $k=P+k_{1}+k_{2}$, the momenta
$k$, $k_{1}$ and $k_{2}$ depend implicitly on $P$, and, therefore, on $q$. We
can make the dependence on $q$ explicit by writing quantities in terms of
dimensionless momenta

$\displaystyle\bar{P}$ $\displaystyle=$ $\displaystyle\frac{P}{\sqrt{P^{2}}},$
(30a) $\displaystyle\bar{k}$ $\displaystyle=$
$\displaystyle\frac{k}{\sqrt{P^{2}}},$ (30b) $\displaystyle\bar{k}_{1}$
$\displaystyle=$ $\displaystyle\frac{k_{1}}{\sqrt{P^{2}}},$ (30c)
$\displaystyle\bar{k}_{2}$ $\displaystyle=$
$\displaystyle\frac{k_{2}}{\sqrt{P^{2}}}.$ (30d)

It is also useful to express the Lorentz invariants in terms of the following
dimensionless variables:

$\displaystyle e_{1}$ $\displaystyle=$ $\displaystyle\bar{k}_{1}\cdot\bar{P},$
(31a) $\displaystyle e_{2}$ $\displaystyle=$
$\displaystyle\bar{k}_{2}\cdot\bar{P},$ (31b) $\displaystyle x$
$\displaystyle=$
$\displaystyle\bar{k}_{1}\cdot\bar{k}_{2}=e_{1}e_{2}(1-\hat{\bm{k}}_{1}\cdot\hat{\bm{k}}_{2}),$
(31c)

where $\hat{\bm{k}}_{i}$ is the unit vector that is parallel to the three-
vector $\bm{k}_{i}$ in the $Q\bar{Q}$ rest frame.

The phase space in Eq. (10) can be expressed in terms of the dimensionless
variables as

$\displaystyle d\Phi_{0}$ $\displaystyle=$ $\displaystyle\frac{4\pi
M}{k^{+}}\,\delta(1-z),$ (32a) $\displaystyle d\Phi_{1}$ $\displaystyle=$
$\displaystyle\frac{4\pi
M^{d-1}}{k^{+}}\,\theta(z_{1})\,\delta(1-z-z_{1})\,\frac{dz_{1}}{4\pi
z_{1}}\frac{d^{d-2}\bar{\bm{k}}_{1\perp}}{(2\pi)^{d-2}},$ (32b) $\displaystyle
d\Phi_{2}$ $\displaystyle=$ $\displaystyle\frac{4\pi
M^{2d-3}}{Sk^{+}}\,\theta(z_{1})\,\theta(z_{2})\,\delta(1-z-z_{1}-z_{2})\,\frac{dz_{1}}{4\pi
z_{1}}\frac{dz_{2}}{4\pi
z_{2}}\frac{d^{d-2}\bar{\bm{k}}_{1\perp}}{(2\pi)^{d-2}}\frac{d^{d-2}\bar{\bm{k}}_{2\perp}}{(2\pi)^{d-2}}.\phantom{xxx}$
(32c)

We have not replaced the overall factor $1/k^{+}$ in Eq. (32) with
$1/(\bar{k}^{+}\sqrt{P^{2}})$, because it ultimately will be cancelled by the
factor $k^{+}$ in $C_{\rm frag}$ in Eq. (8). The ranges of the variables $z$,
$z_{1}$ and $z_{2}$ are completely determined by the $\delta$ and $\theta$
functions. When we expand the fragmentation function in powers of $q$, it is
convenient to make use of the phase space at LO in $q$,

$d\tilde{\Phi}_{n}=d\Phi_{n}\big{|}_{\bm{q}\to\bm{0}},$ (33)

where $\bm{q}\to\bm{0}$ means that, in the phase space in Eq. (32), we replace
$M$ with $2m$. Then $d\Phi_{n}$ can be expressed in terms of
$d\tilde{\Phi}_{n}$ as follows:

$\displaystyle d\Phi_{0}$ $\displaystyle=$
$\displaystyle\frac{E}{m}d\tilde{\Phi}_{0},$ (34a) $\displaystyle d\Phi_{1}$
$\displaystyle=$
$\displaystyle\left(\frac{E}{m}\right)^{3-2\epsilon}d\tilde{\Phi}_{1},$ (34b)
$\displaystyle d\Phi_{2}$ $\displaystyle=$
$\displaystyle\left(\frac{E}{m}\right)^{5-4\epsilon}d\tilde{\Phi}_{2}.$ (34c)

We express our results for the fragmentation contributions in terms of
integrals over the phase spaces $d\tilde{\Phi}_{n}$:

$\tilde{D}[g\to Q\bar{Q}]=D[g\to Q\bar{Q}]\big{|}_{d\Phi_{n}\to
d\tilde{\Phi}_{n}}.$ (35)

The factors of $E/m$ in Eq. (34) are then an additional source of relativistic
corrections.

##  5 Full-QCD calculations

In this section, we compute the relevant fragmentation functions for free
$Q\bar{Q}$ states in full QCD. We have carried out the calculations by writing
independent codes using reduce REDUCE  and using the feyncalc package
Mertig:an  in mathematica MATH . At each stage of the calculations we have
checked that the independent codes give identical results.

The computations are carried out in $d=4-2\epsilon$ dimensions with
dimensional-regulariza-tion scale $\mu$. We use the modified-minimal-
subtraction ($\overline{\rm MS}$) scheme throughout. Then, in $d=4-2\epsilon$
dimensions, there is a factor $[\mu^{2}\exp({{\gamma}}_{\rm
E})/(4\pi)]^{\epsilon}$ that is associated with each factor of the strong
coupling $\alpha_{s}$, where ${{\gamma}}_{\rm E}$ is the Euler-Mascheroni
constant.

In computing the $Q\bar{Q}$ fragmentation functions, it is convenient to make
use of projection operators for the spin and color states of the $Q\bar{Q}$
pair. The projection operators for a $Q\bar{Q}$ pair in the color-singlet and
color-octet configurations are

$\displaystyle\Lambda_{1}$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{N_{c}}}\mathbbm{1},$ (36a)
$\displaystyle\Lambda^{a}_{8}$ $\displaystyle=$ $\displaystyle\sqrt{2}T^{a},$
(36b)

where $\mathbbm{1}$ and $T^{a}$ are the identity matrix and the generator of
the fundamental (triplet) representation of SU(3), $a$ is the adjoint-
representation color index ($a=1,\ldots,N_{c}^{2}-1$), and $N_{c}=3$. Spin-
projection operators at LO in $v$ were first given in Refs. Barbieri:1975am ;
Barbieri:1976fp ; Chang:1979nn ; Guberina:1980dc ; Berger:1980ni . Projectors
accurate to all orders in $v$ were given in Ref. Bodwin:2002hg . For the spin-
triplet state, the projection operator, correct to all orders in $v$, is

$\Lambda(P,q,\epsilon_{S}^{*})=N(\not\\!{\overline{p}}-m)\not\\!{\epsilon}^{*}_{S}\,\frac{\not\\!{P}+2E}{4E}(\not\\!{p}+m),$
(37)

where $\epsilon_{S}$ is the spin polarization of the $Q\bar{Q}$ pair, and
$N=[2\sqrt{2}E(E+m)]^{-1}$. Note that we use nonrelativistic normalization for
the heavy-quark spinors.

The use of the spin projection (37) in $d$ dimensions requires some
justification. It accounts for only the $d-1$ vector polarization states,
which, in the $Q\bar{Q}$ rest frame, correspond to the $d-1$ Pauli matrices
$\sigma_{i}$. In general, in $d$ dimensions, one must consider states that
correspond to products of the $\sigma_{i}$ that are linearly independent of
the $\sigma_{i}$ Braaten:1996rp . These additional states vanish as $\epsilon$
goes to zero. Hence, they can contribute only in conjunction with a pole in
$\epsilon$. The poles in $\epsilon$ in our calculation correspond to soft
divergences. The divergent parts of soft interactions arise from the
convection current on fermions lines, and, hence, do not change the fermion
spin. Therefore, the additional states that correspond to products of the
$\sigma_{i}$ never mix in our calculation with the vector states that
correspond to the $\sigma_{i}$. Consequently, we need consider only the $d-1$
vector states in our calculation.222Some elements of this argument were
presented in Ref. Petrelli:1997ge .

The spin-triplet, color-singlet part of an amplitude $\mathcal{C}$ is

$\mathcal{M}=\textrm{Tr}[\,\mathcal{C}(\Lambda\otimes\Lambda_{1})\,],$ (38)

and the spin-triplet, color-octet part of an amplitude $\mathcal{C}$ is

$\mathcal{M}^{a}=\textrm{Tr}[\,\mathcal{C}(\Lambda\otimes\Lambda_{8}^{a})\,],$
(39)

where the traces are over the Dirac and color indices. The amplitude
$\mathcal{C}$ includes the propagator of the initial gluon, as well as the
associated polarization factor in Eq. (9). In our calculation, the amplitudes
$\mathcal{C}$, $\mathcal{M}$, and $\mathcal{M}^{a}$ are all expressed in terms
of the dimensionless variables in Eq. (30) or invariants that are formed from
them, and so the dependence on $q$ is explicit.

The $S$-wave part of $\mathcal{M}$ (with color index suppressed in the color-
octet case) can be written as an expansion in powers of
$v^{2}=\bm{q}^{2}/m^{2}$:

$\mathcal{M}_{S}=\mathcal{M}_{S0}+\mathcal{M}_{S2}+\mathcal{M}_{S4}+O(\bm{q}^{6}/m^{6}),$
(40)

where

$\displaystyle\mathcal{M}_{S0}$ $\displaystyle=$
$\displaystyle\left(\mathcal{M}\right)_{\bm{q}\to\bm{0}}\;,$ (41a)
$\displaystyle\mathcal{M}_{S2}$ $\displaystyle=$
$\displaystyle\frac{\bm{q}^{2}}{2!(d-1)}\,I^{\alpha\beta}\left(\frac{\partial^{2}\mathcal{M}}{\partial
q^{\alpha}\partial q^{\beta}}\right)_{\bm{q}\to\bm{0}}\;,$ (41b)
$\displaystyle\mathcal{M}_{S4}$ $\displaystyle=$
$\displaystyle\frac{\bm{q}^{4}}{4!(d-1)(d+1)}\,I^{\alpha\beta\gamma\delta}\left(\frac{\partial^{4}\mathcal{M}}{\partial
q^{\alpha}\partial q^{\beta}\partial q^{\gamma}\partial
q^{\delta}}\right)_{\bm{q}\to\bm{0}}\;,$ (41c)

and

$\displaystyle I^{\alpha\beta}$ $\displaystyle=$
$\displaystyle-g^{\alpha\beta}+P^{\alpha}P^{\beta}/(4E^{2}),$ (42a)
$\displaystyle I^{\alpha\beta\gamma\delta}$ $\displaystyle=$ $\displaystyle
I^{\alpha\beta}I^{\gamma\delta}+I^{\alpha\gamma}I^{\beta\delta}+I^{\alpha\delta}I^{\beta\gamma}.$
(42b)

In order to project out the $P$-wave part of the amplitude ${\cal M}$, we
multiply ${\cal M}$ by the $P$-wave orbital-angular-momentum state
$-\sqrt{d-1}\,\epsilon_{L}^{*}\cdot\hat{q}$ and average over the direction of
$\bm{q}$.333 In some calculations in NRQCD, the $P$-wave orbital-angular-
momentum state is normalized as $-\epsilon_{L}^{*}\cdot\hat{q}$. Here,
$\epsilon_{L}$ is the polarization vector for the orbital-angular-momentum
state, and $\hat{q}=(0,\hat{\bm{q}})$ in the rest frame of the $Q\bar{Q}$
pair. Then, the $P$-wave part of the amplitude is

$\mathcal{M}_{P}=\mathcal{M}_{P1}+O(\bm{q}^{3}/m^{3}),$ (43)

where

$\mathcal{M}_{P1}=-\frac{|\bm{q}|}{\sqrt{d-1}}\,\epsilon^{*}_{L\alpha}\,I^{\alpha\beta}\left(\frac{\partial\mathcal{M}}{\partial
q^{\beta}}\right)_{\bm{q}\to\bm{0}}\;.$ (44)

We define squared amplitudes for the color-singlet and color-octet states as

$\displaystyle\mathcal{A}({}^{2s+1}l_{j}^{[1]})$ $\displaystyle=$
$\displaystyle C_{\rm frag}|{\cal M}({}^{2s+1}l_{j}^{[1]})|^{2},$ (45a)
$\displaystyle\mathcal{A}({}^{2s+1}l_{j}^{[8]})$ $\displaystyle=$
$\displaystyle C_{\rm frag}\sum_{a}|{\cal M}^{a}({}^{2s+1}l_{j}^{[8]})|^{2},$
(45b)

where $C_{\rm frag}$ is given in Eq. (8), and it is implicit that there are
sums over the spin and orbital-angular-momentum polarizations of the
$Q\bar{Q}$ states and sums over the polarizations of the initial and final
gluons. Note that

$\sum_{\lambda}\epsilon_{S}^{\alpha*}(\lambda)\epsilon_{S}^{\beta}(\lambda)=\sum_{\lambda}\epsilon_{L}^{\alpha*}(\lambda)\epsilon_{L}^{\beta}(\lambda)=I^{\alpha\beta}.$
(46)

We denote the order-$v^{k}$ contribution to $\mathcal{A}$ by
$\mathcal{A}_{k}$.

###  5.1 $D_{0}[g\to Q\bar{Q}({}^{3}S_{1}^{[8]})]$

Figure 2: Feynman diagram for the fragmentation process $g\to
Q\bar{Q}({}^{3}S_{1}^{[8]})$. The dashed line represents the final-state cut.
The momenta for $Q$ and $\bar{Q}$ on the left side of the cut are
$p=\tfrac{1}{2}P+q$ and $\bar{p}=\tfrac{1}{2}P-q$, respectively. The momenta
on the right side of the cut are $\tfrac{1}{2}P+q^{\prime}$ and
$\tfrac{1}{2}P-q^{\prime}$, respectively. Here, $|\bm{q}^{\prime}|=|\bm{q}|$
in the rest frame of the $Q\bar{Q}$ pair, but we distinguish the directions of
$\bm{q}$ and $\bm{q}^{\prime}$ in order to be able to project out orbital-
angular-momentum states in the amplitude and its complex conjugate.

The diagram for gluon fragmentation into a ${}^{3}S_{1}$ color-octet
$Q\bar{Q}$ pair at order $v^{0}$ and $\alpha_{s}^{1}$ is shown in Fig. 2. A
straightforward computation yields

$\mathcal{A}_{0}({}^{3}S_{1}^{[8]})=\frac{\alpha_{s}k^{+}}{8m^{4}}\left(\frac{\mu^{2}}{4\pi}e^{{{\gamma}}_{\rm
E}}\right)^{\epsilon}.$ (47)

Carrying out the trivial integration over the phase space $d\tilde{\Phi}_{0}$
in Eqs. (32a) and (34a), we obtain

$D_{0}[g\to
Q\bar{Q}({}^{3}S_{1}^{[8]})]=\frac{\pi\alpha_{s}}{m^{3}}\left(\frac{\mu^{2}}{4\pi}e^{{{\gamma}}_{\rm
E}}\right)^{\epsilon}\delta(1-z).$ (48)

###  5.2 $D_{2}[g\to Q\bar{Q}({}^{3}P^{[8]})]$

Figure 3:  One of the four Feynman diagrams for the fragmentation process
$g\to Q\bar{Q}({}^{3}P^{[8]})$. Three additional diagrams can be obtained by
permuting the gluon-fermion vertices on the left and right sides of the cut.

The diagrams for gluon fragmentation into a ${}^{3}P^{[8]}$ color-octet
$Q\bar{Q}$ pair at LO in $\alpha_{s}$ and $v$ are shown in Fig. 3. We obtain

$\mathcal{A}_{1}({}^{3}P^{[8]})=\frac{k^{+}\pi\alpha_{s}^{2}(N_{c}^{2}-4)\bm{q}^{2}z^{1-2\epsilon}}{8N_{c}(d-1)(d-2)m^{8}}\left(\frac{\mu^{2}}{4\pi}e^{{{\gamma}}_{\rm
E}}\right)^{2\epsilon}\sum_{n=0}^{3}\frac{\rho_{n}(z)}{e_{1}^{n}(1+2e_{1})^{2}},$
(49)

where $\rho_{n}(z)$ are given by

$\displaystyle\rho_{0}(z)$ $\displaystyle=$
$\displaystyle(3-4\epsilon)(2-2\epsilon-4z+4z^{2}),$ (50a)
$\displaystyle\rho_{1}(z)$ $\displaystyle=$ $\displaystyle
2\big{[}5-10\epsilon+4\epsilon^{2}-z(5-12\epsilon)+2z^{2}(1-4\epsilon)\big{]},$
(50b) $\displaystyle\rho_{2}(z)$ $\displaystyle=$ $\displaystyle
3-12\epsilon+4\epsilon^{2}+2z(3+4\epsilon)-z^{2}(5+4\epsilon),$ (50c)
$\displaystyle\rho_{3}(z)$ $\displaystyle=$ $\displaystyle-2(1-z)^{2}.$ (50d)

We carry out the integration over the phase space $d\tilde{\Phi}_{1}$ in Eqs.
(32b) and (34b) by making use of the methods that are described in Appendix A.
Then, we obtain

$\displaystyle D_{2}[g\to Q\bar{Q}({}^{3}P^{[8]})]$ $\displaystyle=$
$\displaystyle\frac{8\alpha_{s}^{2}\,\bm{q}^{2}}{(d-1)m^{5}}\,\frac{N_{c}^{2}-4}{4N_{c}}(1-\epsilon)\Gamma(1+\epsilon)\left(\frac{\mu^{2}}{4\pi}e^{{{\gamma}}_{\rm
E}}\right)^{\epsilon}\left(\frac{\mu^{2}}{4m^{2}}e^{{{\gamma}}_{\rm
E}}\right)^{\epsilon}$ (51a)
$\displaystyle\times\left[-\frac{1}{2\epsilon_{\rm
IR}}\,\delta(1-z)+f(z)\right],$ where the finite function $f(z)$ is defined by
$\displaystyle f(z)$ $\displaystyle=$
$\displaystyle\left[\frac{1}{(1-z)^{1+2\epsilon}}\right]_{+}+\frac{1}{4(1-\epsilon)^{2}}\Bigg{\\{}\frac{(1-z)^{-\epsilon}-(1-z)^{-2\epsilon}}{\epsilon}\,(13-7z)$
(51b)
$\displaystyle+(1-z)^{-2\epsilon}[10+3z-5z^{2}+2\epsilon(2-2z+z^{2})-4\epsilon^{2}]$
$\displaystyle+\frac{1}{2}(1-z)^{-\epsilon}[-28+15z+\epsilon(8-11z)+4\epsilon^{2}z]\Bigg{\\}}.$

Here, the distribution $[g(z)]_{+}$ is defined by

$\int_{0}^{1}dz\,h(z)[g(z)]_{+}\equiv\int_{0}^{1}dz\,[h(z)-h(1)]g(z).$ (52)

In extracting the pole in Eq. (51), we have made use of the identity,

$\frac{1}{(1-z)^{1+n\epsilon}}=-\frac{1}{n\epsilon}\,\delta(1-z)+\left[\frac{1}{(1-z)^{1+n\epsilon}}\right]_{+},$
(53)

which applies when the domain of integration is $0\leq z\leq 1$. The
expression in Eq. (51b) gives the exact $\epsilon$ dependence. We can expand
the plus function $[1/(1-z)^{1+n\epsilon}]_{+}$ as

$\left[\frac{1}{(1-z)^{1+n\epsilon}}\right]_{+}=\sum_{k=0}^{\infty}\frac{(-n\epsilon)^{k}}{k!}\left[\frac{\log^{k}(1-z)}{1-z}\right]_{+}.$
(54)

In the analysis of $f(z)$, we need to keep only terms through order
$\epsilon^{1}$. Then, we can simplify $f(z)$ as follows:

$\displaystyle f(z)\\!$ $\displaystyle=$
$\displaystyle\\!\\!\left(\frac{1}{1-z}\right)_{+}\\!\\!-2\epsilon\left[\frac{\log(1-z)}{1-z}\right]_{+}\\!\\!+\frac{1}{8}\big{[}\\!-8+21z-10z^{2}+2(13-7z)\log(1-z)\big{]}$
$\displaystyle+$
$\displaystyle\frac{\epsilon}{8}\big{[}(23-16z)z+5(8-11z+4z^{2})\log(1-z)-3(13-7z)\log^{2}(1-z)\big{]}+O(\epsilon^{2}).$

###  5.3 $D_{4}[g\to Q\bar{Q}({}^{3}S_{1}^{[1]})]$

Figure 4:  One of the 36 Feynman diagrams for the fragmentation process $g\to
Q\bar{Q}({}^{3}S_{1}^{[1]})$. Thirty-five additional diagrams can be obtained
by permuting the gluon-fermion vertices on the left and right sides of the
cut.

The diagrams for gluon fragmentation into a ${}^{3}S_{1}$ color-singlet
$Q\bar{Q}$ pair at LO in $\alpha_{s}$ are shown in Fig. 4. Through relative
order $v^{4}$, the relevant squared amplitudes are

$\displaystyle\mathcal{A}_{0}({}^{3}S_{1}^{[1]})$ $\displaystyle=$
$\displaystyle C_{\rm frag}\,|{\cal M}_{0}({}^{3}S_{1}^{[1]})|^{2},$ (56a)
$\displaystyle\mathcal{A}_{2}({}^{3}S_{1}^{[1]})$ $\displaystyle=$
$\displaystyle 2\,C_{\rm frag}\,{\rm Re}[{\cal M}_{2}({}^{3}S_{1}^{[1]}){\cal
M}_{0}^{*}({}^{3}S_{1}^{[1]})],$ (56b)
$\displaystyle\mathcal{A}_{4}({}^{3}S_{1}^{[1]})$ $\displaystyle=$
$\displaystyle C_{\rm frag}\left\\{|{\cal M}_{2}({}^{3}S_{1}^{[1]})|^{2}+2{\rm
Re}[{\cal M}_{4}({}^{3}S_{1}^{[1]}){\cal
M}_{0}^{*}({}^{3}S_{1}^{[1]})]\right\\}.$ (56c)

The order-$v^{0}$ contribution $\mathcal{A}_{0}({}^{3}S_{1}^{[1]})$ and the
order-$v^{2}$ contribution $\mathcal{A}_{2}({}^{3}S_{1}^{[1]})$ have been
computed previously in Ref. Bodwin:2003wh . Here, we wish to compute the
order-$v^{4}$ contribution $\mathcal{A}_{4}({}^{3}S_{1}^{[1]})$.

The integration of $\mathcal{A}_{4}({}^{3}S_{1}^{[1]})$ over the phase space
contains soft divergences. These can arise when one or both of the final-state
gluons become soft. We identify the divergent part of the integrand in Eq.
(56c) that arises when both gluons become soft by making the substitutions
$k_{1}\to k_{1}\lambda$ and $k_{2}\to k_{2}\lambda$, multiplying by
$\lambda^{4}$, and taking the limit $\lambda\to 0$. The result is

$\mathcal{S}_{12}=\frac{C_{{}_{\\!\mathcal{S}}}(d-2)}{e_{1}^{4}e_{2}^{4}}\left[(d-2)e_{1}^{2}e_{2}^{2}-2e_{1}e_{2}\,x+x^{2}\right],$
(57)

where

$C_{{}_{\\!\mathcal{S}}}=\frac{k^{+}\pi^{2}\alpha_{s}^{3}z^{d-3}}{(d-1)^{2}(d-2)m^{8}}\,\frac{N_{c}^{2}-4}{4N_{c}^{2}}\left(\frac{\mu^{2}}{4\pi}e^{{{\gamma}}_{\rm
E}}\right)^{3\epsilon}\frac{\bm{q}^{4}}{m^{4}}.$ (58)

We identify the divergent part of the integrand in Eq. (56c) that arises when
only $k_{1}$ ($k_{2}$) becomes soft by subtracting $\mathcal{S}_{12}$, making
the substitution $k_{1}\to k_{1}\lambda$ ($k_{2}\to k_{2}\lambda$),
multiplying by $\lambda^{2}$, and taking the limit $\lambda\to 0$. The result
is

$\displaystyle\mathcal{S}_{1}$ $\displaystyle=$
$\displaystyle\frac{C_{{}_{\\!\mathcal{S}}}}{e_{1}^{4}e_{2}^{4}(1+2e_{2})^{2}}\sum_{n=0}^{2}x^{n}\,h_{n}(e_{1},e_{2},z_{1},z_{2}),$
(59a) $\displaystyle\mathcal{S}_{2}$ $\displaystyle=$
$\displaystyle\frac{C_{{}_{\\!\mathcal{S}}}}{e_{2}^{4}e_{1}^{4}(1+2e_{1})^{2}}\sum_{n=0}^{2}x^{n}\,h_{n}(e_{2},e_{1},z_{2},z_{1}),$
(59b)

where the $h_{n}(e_{1},e_{2},z_{1},z_{2})$ are given by

$\displaystyle h_{0}(e_{1},e_{2},z_{1},z_{2})$ $\displaystyle=$ $\displaystyle
2e_{2}(1+2e_{2})\big{\\{}-2(3-2\epsilon)(1-\epsilon)e_{1}^{2}e_{2}^{2}+z_{1}e_{1}e_{2}^{2}+z_{2}e_{1}^{2}e_{2}$
$\displaystyle\times[1-2(1-2\epsilon)e_{2}]-z_{1}^{2}e_{2}^{2}(1+\epsilon
e_{2})+z_{1}z_{2}e_{1}e_{2}[1-2(1-\epsilon)e_{2}]$
$\displaystyle+z_{2}^{2}e_{1}^{2}[-1+(1-3\epsilon)e_{2}+2(1-2\epsilon)e_{2}^{2}\,]\,\big{\\}}+4(3-2\epsilon)(1-\epsilon)e_{1}^{2}e_{2}^{4},$
$\displaystyle h_{1}(e_{1},e_{2},z_{1},z_{2})$ $\displaystyle=$ $\displaystyle
2e_{2}(1+2e_{2})\big{\\{}z_{1}e_{2}(1+2\epsilon
e_{2})+z_{2}e_{1}[1+(3-2\epsilon)e_{2}]-z_{1}z_{2}(1+2\epsilon e_{2}$
$\displaystyle+2\epsilon
e_{2}^{2})-z_{2}^{2}e_{1}[3-2\epsilon+2(1-\epsilon)e_{2}]\big{\\}}+2e_{1}e_{2}^{2}[5-4\epsilon+6(1-\epsilon)e_{2}],$
$\displaystyle h_{2}(e_{1},e_{2},z_{1},z_{2})$ $\displaystyle=$
$\displaystyle-(1+2e_{2})\big{\\{}e_{2}(4-3\epsilon+2\epsilon
e_{2})+2z_{2}e_{2}(1-2\epsilon-2\epsilon e_{2})$ (60c)
$\displaystyle-2z_{2}^{2}(1+e_{2})(1-\epsilon-\epsilon
e_{2})\big{\\}}-(2-\epsilon)e_{2}.$

We carry out the integrations of $\mathcal{S}_{12}$, $\mathcal{S}_{1}$ and
$\mathcal{S}_{2}$ over the phase space $d\tilde{\Phi}_{2}$ in Eq. (34c) [see
also Eq. (32c)] by making use of the methods that are described in Appendix A.
Then, we obtain

$\displaystyle I[\mathcal{S}_{12}]$ $\displaystyle=$
$\displaystyle\left\\{\frac{1}{8\epsilon_{\rm
IR}^{2}}\,\delta(1-z)-\frac{1}{2\epsilon_{\rm
IR}}\left[\frac{1}{(1-z)^{1+4\epsilon}}\right]_{+}+\frac{1-z^{1+2\epsilon}}{2\epsilon_{\rm
IR}(1-z)^{1+4\epsilon}}\right\\}$ (61)
$\displaystyle\times\left(\frac{8\alpha_{s}}{3\pi
m^{2}}\right)^{2}\frac{N_{c}^{2}-4}{16N_{c}^{2}}\,\frac{\pi\alpha_{s}}{(d-1)m^{3}}\left(\frac{\mu^{2}}{4\pi}e^{{{\gamma}}_{\rm
E}}\right)^{\epsilon}\frac{\bm{q}^{4}}{d-1}$
$\displaystyle\times\left(\frac{\mu^{2}}{4m^{2}}e^{{{\gamma}}_{\rm
E}}\right)^{2\epsilon}\frac{\Gamma^{2}(1+\epsilon)\Gamma^{2}(1-2\epsilon)}{\Gamma(1-4\epsilon)}(1-\epsilon)(6-2\epsilon-\epsilon^{2}-2\epsilon^{3}),$

where we have used the identity in Eq. (53). If we expand Eq. (61) in powers
of $\epsilon$, then we find that

$\displaystyle I[\mathcal{S}_{12}]$ $\displaystyle=$
$\displaystyle\bigg{\\{}\frac{1}{8\epsilon_{\rm
IR}^{2}}\,\delta(1-z)-\frac{1}{2\epsilon_{\rm
IR}}\left[\delta(1-z)\left(\frac{1}{3}-\log\frac{\mu}{2m}\right)+\left(\frac{1}{1-z}\right)_{+}\\!\\!-1\right]$
(62)
$\displaystyle\;+\delta(1-z)\bigg{(}\frac{1-3\pi^{2}}{48}-\frac{2}{3}\log\frac{\mu}{2m}+\log^{2}\frac{\mu}{2m}\bigg{)}+\left(\frac{1}{1-z}\right)_{+}\\!\\!\left(\frac{2}{3}-2\log\frac{\mu}{2m}\right)$
$\displaystyle\;+2\left[\frac{\log(1-z)}{1-z}\right]_{+}\\!\\!-2\left(\frac{1}{3}-\log\frac{\mu}{2m}\right)-\frac{z\log
z}{1-z}-2\log(1-z)\bigg{\\}}$
$\displaystyle\times\left(\frac{8\alpha_{s}}{3\pi
m^{2}}\right)^{2}\frac{N_{c}^{2}-4}{16N_{c}^{2}}\,\frac{\pi\alpha_{s}}{(d-1)m^{3}}\left(\frac{\mu^{2}}{4\pi}e^{{{\gamma}}_{\rm
E}}\right)^{\epsilon}\frac{6\bm{q}^{4}}{d-1}+O(\epsilon).$

We also find that

$\displaystyle I[\mathcal{S}_{1}]$ $\displaystyle=$ $\displaystyle
I[\mathcal{S}_{2}]=\left(-\frac{\tau_{1}}{2\epsilon_{\rm
IR}}+\tau_{0}\right)\left(\frac{\mu^{2}}{4m^{2}}e^{{{\gamma}}_{\rm
E}}\right)^{2\epsilon}\frac{z^{-2+2\epsilon}(1-z)^{-4\epsilon}\Gamma^{2}(1+\epsilon)}{48(1-\epsilon)}$
(63) $\displaystyle\times$ $\displaystyle\left(\frac{8\alpha_{s}}{3\pi
m^{2}}\right)^{2}\frac{N_{c}^{2}-4}{16N_{c}^{2}}\left(\frac{\mu^{2}}{4\pi}e^{{{\gamma}}_{\rm
E}}\right)^{\epsilon}\frac{\pi\alpha_{s}}{(d-1)m^{3}}\frac{\bm{q}^{4}}{d-1}+O(\epsilon),$

where

$\displaystyle\tau_{0}$ $\displaystyle=$ $\displaystyle
3\bigg{\\{}12(13-7z)z^{2}\,\text{Li}_{2}(1-z)+3(7z-13)z^{2}\log^{2}(1-z)$
(64a) $\displaystyle-2z\left[z^{3}-(43+7\pi^{2})z^{2}+6(2z+9)z^{2}\log
z+(48+13\pi^{2})z-42\right]$
$\displaystyle+3(8z^{4}-7z^{3}+34z^{2}-44z+28)\log(1-z)\bigg{\\}},$
$\displaystyle{}\tau_{1}$ $\displaystyle=$ $\displaystyle
18z^{2}\left[z(21-10z)+2(13-7z)\log(1-z)\right].$ (64b)

Here, $\text{Li}_{2}(z)$ is the Spence function, which is defined by

$\text{Li}_{2}(x)=-\int_{0}^{x}\,dt\,\frac{\log(1-t)}{t}=\sum_{k=1}^{\infty}\frac{x^{k}}{k^{2}}.$
(65)

The fragmentation-function contribution of
$\mathcal{A}_{4}({}^{3}S_{1}^{[1]})$ is then

$\tilde{D}_{4}[g\to Q\bar{Q}({}^{3}S_{1}^{[1]})]=\tilde{D}_{4}[g\to
Q\bar{Q}({}^{3}S_{1}^{[1]})]^{\rm
finite}+I[\mathcal{S}_{12}]+I[\mathcal{S}_{1}]+I[\mathcal{S}_{2}],$ (66)

where

$\tilde{D}_{4}[g\to Q\bar{Q}({}^{3}S_{1}^{[1]})]^{\rm finite}=\int
d\tilde{\Phi}_{2}[\mathcal{A}_{4}({}^{3}S_{1}^{[1]})-\mathcal{S}_{12}-\mathcal{S}_{1}-\mathcal{S}_{2}].$
(67)

Since $\tilde{D}_{4}[g\to Q\bar{Q}({}^{3}S_{1}^{[1]})]^{\rm finite}$ contains
no soft divergences, we compute it in $d=4$ dimensions, using numerical
integration over the phase space. This computation is described in Sec. 8.

As was mentioned previously, $D_{4}[g\to Q\bar{Q}({}^{3}S_{1}^{[1]})]^{\rm
finite}$ also contains contributions that arise from the difference between
the fully relativistic phase space $d\Phi_{2}$ and the order-$v^{0}$ phase
space $d\tilde{\Phi}_{2}$ [Eq. (34c)]. Hence, we write

$D_{4}[g\to Q\bar{Q}({}^{3}S_{1}^{[1]})]=D_{4}[g\to
Q\bar{Q}({}^{3}S_{1}^{[1]})]^{\rm
finite}+I[\mathcal{S}_{12}]+I[\mathcal{S}_{1}]+I[\mathcal{S}_{2}],$ (68)

where

$\displaystyle D_{4}[g\to Q\bar{Q}({}^{3}S_{1}^{[1]})]^{\rm finite}$
$\displaystyle=$ $\displaystyle\tilde{D}_{4}[g\to
Q\bar{Q}({}^{3}S_{1}^{[1]})]^{\rm
finite}+\frac{5\bm{q}^{2}}{2m^{2}}\tilde{D}_{2}[g\to
Q\bar{Q}({}^{3}S_{1}^{[1]})]$ (69)
$\displaystyle+\frac{15\bm{q}^{4}}{8m^{4}}\tilde{D}_{0}[g\to
Q\bar{Q}({}^{3}S_{1}^{[1]})],$ $\displaystyle\tilde{D}_{0}[g\to
Q\bar{Q}({}^{3}S_{1}^{[1]})]$ $\displaystyle=$
$\displaystyle\int\\!d\tilde{\Phi}_{2}\,\mathcal{A}_{0}({}^{3}S_{1}^{[1]}),$
(70a) $\displaystyle\tilde{D}_{2}[g\to Q\bar{Q}({}^{3}S_{1}^{[1]})]$
$\displaystyle=$
$\displaystyle\int\\!d\tilde{\Phi}_{2}\,\mathcal{A}_{2}({}^{3}S_{1}^{[1]}),$
(70b)

and we have used the expansion

$\left(\frac{E}{m}\right)^{5}=1+\frac{5\bm{q}^{2}}{2m^{2}}+\frac{15\bm{q}^{4}}{8m^{4}}+O(\bm{q}^{6}/m^{6}).$
(71)

$\tilde{D}_{0}[g\to Q\bar{Q}({}^{3}S_{1}^{[1]})]$ and $\tilde{D}_{2}[g\to
Q\bar{Q}({}^{3}S_{1}^{[1]})]$ are finite, and they have been evaluated in
$d=4$ dimensions in Ref. Bodwin:2003wh . We have checked those calculations,
computing the phase-space integrations in $\tilde{D}_{0}[g\to
Q\bar{Q}({}^{3}S_{1}^{[1]})]$ and $\tilde{D}_{2}[g\to
Q\bar{Q}({}^{3}S_{1}^{[1]})]$ numerically.

##  6 NRQCD LDMEs

In this section we compute the NRQCD LDMEs for free
$Q\bar{Q}({}^{3}S_{1}^{[1]})$ states that are relevant through relative order
$v^{4}$. These computations are carried out in each case at the leading
nontrivial order in $\alpha_{s}$ in $d=4-2\epsilon$ dimensions, with
dimensional-regularization scale $\mu$. We remind the reader that, because we
use the $\overline{\rm MS}$ scheme in computing the QCD corrections to NRQCD
LDMEs, there is a factor $[\mu^{2}\exp({{\gamma}}_{\rm E})/(4\pi)]^{\epsilon}$
that is associated with each factor of the strong coupling $\alpha_{s}$ in $d$
dimensions.

###  6.1 Order $\alpha_{s}^{0}$

The matrix elements of the $Q\bar{Q}$ NRQCD operators at order
$\alpha_{s}^{0}$ are normalized as

$\displaystyle\langle
0|\mathcal{O}_{0}^{Q\bar{Q}({}^{1}S_{0}^{[8]})}({}^{1}S_{0}^{[8]})|0\rangle^{(0)}$
$\displaystyle=$ $\displaystyle(N_{c}^{2}-1),$ (72a) $\displaystyle\langle
0|\mathcal{O}_{0}^{Q\bar{Q}({}^{3}S_{1}^{[1]})}({}^{3}S_{1}^{[1]})|0\rangle^{(0)}$
$\displaystyle=$ $\displaystyle 2(d-1)N_{c},$ (72b) $\displaystyle\langle
0|\mathcal{O}_{0}^{Q\bar{Q}({}^{3}S_{1}^{[8]})}({}^{3}S_{1}^{[8]})|0\rangle^{(0)}$
$\displaystyle=$ $\displaystyle(d-1)(N_{c}^{2}-1),$ (72c)
$\displaystyle\langle
0|\mathcal{O}_{2}^{Q\bar{Q}({}^{3}S_{1}^{[n]})}({}^{3}S_{1}^{[n]})|0\rangle^{(0)}$
$\displaystyle=$ $\displaystyle\bm{q}^{2}\langle
0|\mathcal{O}_{0}^{Q\bar{Q}({}^{3}S_{1}^{[n]})}({}^{3}S_{1}^{[n]})|0\rangle^{(0)},$
(72d) $\displaystyle\langle
0|\mathcal{O}_{4}^{Q\bar{Q}({}^{3}S_{1}^{[n]})}({}^{3}S_{1}^{[n]})|0\rangle^{(0)}$
$\displaystyle=$ $\displaystyle\bm{q}^{4}\langle
0|\mathcal{O}_{0}^{Q\bar{Q}({}^{3}S_{1}^{[n]})}({}^{3}S_{1}^{[n]})|0\rangle^{(0)},$
(72e) $\displaystyle\langle
0|\mathcal{O}_{0}^{Q\bar{Q}({}^{3}P^{[8]})}({}^{3}P^{[8]})|0\rangle^{(0)}$
$\displaystyle=$ $\displaystyle\bm{q}^{2}(d-1)(N_{c}^{2}-1),$ (72f)

where a sum over the final-state polarizations is implied. The superscript
$(k)$ indicates the order in $\alpha_{s}$.

###  6.2 Order $\alpha_{s}$

####  6.2.1 ${}^{3}P^{[8]}\to{}^{3}S_{1}^{[1]}$

Figure 5:  One of the four Feynman diagrams for the computation of the LDME
$\langle
0|\mathcal{O}_{0}^{Q\bar{Q}({}^{3}S_{1}^{[1]})}[Q\bar{Q}({}^{3}P^{[8]})]|0\rangle^{(1)}$.
The solid circles represent the $Q\bar{Q}$ operators in the LDME. As in the
full-QCD calculation, we take the free $Q$ and $\bar{Q}$ states to have
momenta $\tfrac{1}{2}P+q$ and $\tfrac{1}{2}P-q$ on the left side of the cut
and $\tfrac{1}{2}P+q^{\prime}$ and $\tfrac{1}{2}P-q^{\prime}$ on the right
side of the cut, where, $|\bm{q}^{\prime}|=|\bm{q}|$ in the rest frame of the
$Q\bar{Q}$ pair, but we distinguish the directions of $\bm{q}$ and
$\bm{q}^{\prime}$ in order to be able to project out orbital-angular-momentum
states in the amplitude and its complex conjugate. Three additional diagrams
can be obtained by permuting the gluon-fermion and operator vertices on the
left and right sides of the cut.

In order $\alpha_{s}$, the ${}^{3}P^{[8]}$ $Q\bar{Q}$ operator can couple to
the ${}^{3}S_{1}^{[1]}$ state through the diagrams that are shown in Fig.
5.444We suppress Wilson lines in diagrams involving NRQCD production operators
Nayak:2005rw ; Nayak:2005rt  because the diagrams involving interactions with
Wilson lines vanish for the orders in $\alpha_{s}$ and the operator and final-
state quantum numbers that we consider.

The set of diagrams in Fig. 5 is gauge invariant. We find it convenient to
work in the $Q\bar{Q}$ center-of-momentum frame and to compute the diagrams in
the Coulomb gauge. Then, the real gluons in the final state must be
transverse. Through relative order $v^{2}$, a straightforward computation in
dimensional regularization gives

$\langle 0|{\cal
O}_{0}^{Q\bar{Q}({}^{3}S_{1}^{[1]})}({}^{3}P^{[8]})|0\rangle^{(1)}=\left(\mathcal{M}_{a}^{(1)}+\mathcal{M}_{b}^{(1)}+\mathcal{M}_{c}^{(1)}+\mathcal{M}_{d}^{(1)}\right)_{{}^{3}S_{1}^{[1]}},$
(73)

where

$\displaystyle\mathcal{M}_{a}^{(1)}$ $\displaystyle=$
$\displaystyle\pi\alpha_{s}\left(\frac{\mu^{2}}{4\pi}e^{{{\gamma}}_{\rm
E}}\right)^{\epsilon}\frac{N_{c}^{2}-1}{4N_{c}^{2}}\,\xi^{\dagger}\sigma^{k}\eta\,\eta^{\dagger}\sigma^{k}\xi$
(74a)
$\displaystyle\times\int\frac{d^{d-1}\bm{k}}{2|\bm{k}|(2\pi)^{d-1}}\frac{(k+2q)^{l}(k+2q)^{i}(\delta^{ij}-\hat{\bm{k}}^{i}\hat{\bm{k}}^{j})(k+2q^{\prime})^{j}(k+2q^{\prime})^{l}}{(\bm{k}^{2}+2\bm{k}\cdot\bm{q}+2m|\bm{k}|)(\bm{k}^{2}+2\bm{k}\cdot\bm{q}^{\prime}+2m|\bm{k}|)},\phantom{xxx}$
$\displaystyle\mathcal{M}_{b}^{(1)}$ $\displaystyle=$
$\displaystyle\pi\alpha_{s}\left(\frac{\mu^{2}}{4\pi}e^{{{\gamma}}_{\rm
E}}\right)^{\epsilon}\frac{N_{c}^{2}-1}{4N_{c}^{2}}\,\xi^{\dagger}\sigma^{k}\eta\,\eta^{\dagger}\sigma^{k}\xi$
(74b)
$\displaystyle\times\int\frac{d^{d-1}\bm{k}}{2|\bm{k}|(2\pi)^{d-1}}\frac{(-k+2q)^{l}(-k+2q)^{i}(\delta^{ij}-\hat{\bm{k}}^{i}\hat{\bm{k}}^{j})(k+2q^{\prime})^{j}(k+2q^{\prime})^{l}}{(\bm{k}^{2}-2\bm{k}\cdot\bm{q}+2m|\bm{k}|)(\bm{k}^{2}+2\bm{k}\cdot\bm{q}^{\prime}+2m|\bm{k}|)},\phantom{xxx}$
$\displaystyle\mathcal{M}_{c}^{(1)}$ $\displaystyle=$
$\displaystyle\pi\alpha_{s}\left(\frac{\mu^{2}}{4\pi}e^{{{\gamma}}_{\rm
E}}\right)^{\epsilon}\frac{N_{c}^{2}-1}{4N_{c}^{2}}\,\xi^{\dagger}\sigma^{k}\eta\,\eta^{\dagger}\sigma^{k}\xi$
(74c)
$\displaystyle\times\int\frac{d^{d-1}\bm{k}}{2|\bm{k}|(2\pi)^{d-1}}\frac{(-k+2q)^{l}(-k+2q)^{i}(\delta^{ij}-\hat{\bm{k}}^{i}\hat{\bm{k}}^{j})(-k+2q^{\prime})^{j}(-k+2q^{\prime})^{l}}{(\bm{k}^{2}-2\bm{k}\cdot\bm{q}+2m|\bm{k}|)(\bm{k}^{2}-2\bm{k}\cdot\bm{q}^{\prime}+2m|\bm{k}|)},\phantom{xxx}$
$\displaystyle\mathcal{M}_{d}^{(1)}$ $\displaystyle=$
$\displaystyle\pi\alpha_{s}\left(\frac{\mu^{2}}{4\pi}e^{{{\gamma}}_{\rm
E}}\right)^{\epsilon}\frac{N_{c}^{2}-1}{4N_{c}^{2}}\,\xi^{\dagger}\sigma^{k}\eta\,\eta^{\dagger}\sigma^{k}\xi$
(74d)
$\displaystyle\times\int\frac{d^{d-1}\bm{k}}{2|\bm{k}|(2\pi)^{d-1}}\frac{(k+2q)^{l}(k+2q)^{i}(\delta^{ij}-\hat{\bm{k}}^{i}\hat{\bm{k}}^{j})(-k+2q^{\prime})^{j}(-k+2q^{\prime})^{l}}{(\bm{k}^{2}+2\bm{k}\cdot\bm{q}+2m|\bm{k}|)(\bm{k}^{2}-2\bm{k}\cdot\bm{q}^{\prime}+2m|\bm{k}|)}.\phantom{xxx}$

In Eq. (74), $\xi$ and $\eta^{\dagger}$ are the Pauli spinors for the free $Q$
and $\bar{Q}$ states, respectively. The subscript ${}^{3}S_{1}^{[1]}$ in Eq.
(73) indicates that the bispinors $\xi^{\dagger}\eta$ and $\eta^{\dagger}\xi$
are in color-singlet, spin-triplet states and that we project onto $S$-wave
states by averaging over the directions of $\bm{q}$ and $\bm{q}^{\prime}$. A
sum over the polarizations of the spin-triplet $Q\bar{Q}$ pair is implicit.

We expand the integrands in Eq. (74) in powers of $1/m$. In dimensional
regularization, only the leading power contributes because the expressions for
higher powers in $1/m$ produce power-divergent, homogeneous integrals.555This
approach was first used in Appendix B of Ref. BBL . It has been discussed
subsequently in Refs. Braaten:1996rp ; Manohar:1997qy ; Beneke:1997av . In
Ref. Bodwin:1998mn , it was pointed out that this approach allocates
contributions that are infrared finite to the short-distance coefficients,
rather than to the LDMEs, and is, therefore, the NRQCD analogue of the
standard methods for computing dimensionally regulated short-distance
coefficients for hard-scattering processes in collinear factorization in QCD.
The result is

$\displaystyle\langle 0|{\cal
O}_{0}^{Q\bar{Q}({}^{3}S_{1}^{[1]})}({}^{3}P^{[8]})|0\rangle^{(1)}$
$\displaystyle=$
$\displaystyle\frac{8\pi\alpha_{s}}{m^{2}}\,\frac{d-2}{d-1}\left(\frac{\mu^{2}}{4\pi}e^{{{\gamma}}_{\rm
E}}\right)^{\epsilon}\frac{N_{c}^{2}-1}{4N_{c}^{2}}\left(\xi^{\dagger}q^{\prime
i}q^{\prime
l}\sigma^{k}\eta\,\eta^{\dagger}q^{i}q^{l}\sigma^{k}\xi\right)_{{}^{3}S_{1}^{[1]}}$
(75) $\displaystyle\times\int\frac{d^{d-1}\bm{k}}{|\bm{k}|^{3}(2\pi)^{d-1}}.$

The remaining integration over the spatial components of $k$ is
straightforward and yields

$\displaystyle\langle 0|{\cal
O}_{0}^{Q\bar{Q}({}^{3}S_{1}^{[1]})}({}^{3}P^{[8]})|0\rangle^{(1)}=\frac{8\alpha_{s}c(\epsilon)}{3\pi
m^{2}}\left(\frac{1}{2\epsilon_{\rm UV}}-\frac{1}{2\epsilon_{\rm
IR}}\right)\frac{N_{c}^{2}-1}{4N_{c}^{2}}\left(\xi^{\dagger}q^{\prime
i}q^{\prime
l}\sigma^{k}\eta\,\eta^{\dagger}q^{i}q^{l}\sigma^{k}\xi\right)_{{}^{3}S_{1}^{[1]}}.$

Here, we have separated the ultraviolet (UV) and infrared (IR) divergent
contributions of the scaleless integral. The quantity $c(\epsilon)$ is given
by

$c(\epsilon)=\frac{(\mu^{2}e^{{{\gamma}}_{\rm
E}})^{\epsilon}(1-\epsilon)\Gamma(\tfrac{1}{2})}{\left(1-\frac{2}{3}\epsilon\right)(1-2\epsilon)\Gamma(\tfrac{1}{2}-\epsilon)}.$
(77)

Note that

$c(0)=1.$ (78)

Now,

$\displaystyle\frac{c(\epsilon)}{\epsilon_{\rm UV}}$ $\displaystyle=$
$\displaystyle\frac{1}{\epsilon_{\rm UV}}+\frac{c(\epsilon)-1}{\epsilon},$
(79a) $\displaystyle\frac{c(\epsilon)}{\epsilon_{\rm IR}}$ $\displaystyle=$
$\displaystyle\frac{1}{\epsilon_{\rm IR}}+\frac{c(\epsilon)-1}{\epsilon},$
(79b)

where we have dropped the subscripts “UV” and “IR” in the second terms of the
above equations because those terms are finite. Hence, we have

$c(\epsilon)\left(\frac{1}{2\epsilon_{\rm UV}}-\frac{1}{2\epsilon_{\rm
IR}}\right)=\frac{1}{2\epsilon_{\rm UV}}-\frac{1}{2\epsilon_{\rm IR}}.$ (80)

Therefore,

$\displaystyle\langle 0|{\cal
O}_{0}^{Q\bar{Q}({}^{3}S_{1}^{[1]})}({}^{3}P^{[8]})|0\rangle^{(1)}=\frac{8\alpha_{s}}{3\pi
m^{2}}\left(\frac{1}{2\epsilon_{\rm UV}}-\frac{1}{2\epsilon_{\rm
IR}}\right)\frac{N_{c}^{2}-1}{4N_{c}^{2}}\left(\xi^{\dagger}q^{\prime
i}q^{\prime
l}\sigma^{k}\eta\,\eta^{\dagger}q^{i}q^{l}\sigma^{k}\xi\right)_{{}^{3}S_{1}^{[1]}}.$

We renormalize $\langle 0|{\cal
O}_{0}^{Q\bar{Q}({}^{3}S_{1}^{[1]})}({}^{3}P^{[8]})|0\rangle^{(1)}$ in the
$\overline{\rm MS}$ scheme. In the $\overline{\rm MS}$ scheme, one constructs
the counterterm for a UV-divergent subdiagram by subtracting the poles in
$\epsilon_{\rm UV}$ that appear in that subdiagram.666In some versions of the
$\overline{\rm MS}$ scheme, one subtracts constants, as well as poles, in
constructing the counterterms. These constants are accounted for in method
that we use in this paper by the factors $(\mu^{2}e^{\gamma_{\rm
E}})^{2\epsilon}$ that are associated with $g_{s}^{2}$. In order to insure
that the counterterm contribution to an LDME removes precisely the
contribution that is proportional to the UV divergence in the divergent
subdiagram, one must compute all factors that are external to the divergent
subdiagram, such as the projections of external momenta onto particular
angular-momentum states, in $d=4-2\epsilon$ dimensions. Hence, we find that
the $\overline{\rm MS}$-counterterm contribution to $\langle 0|{\cal
O}_{0}^{Q\bar{Q}({}^{3}S_{1}^{[1]})}({}^{3}P^{[8]})|0\rangle^{(1)}$ is

$\displaystyle\delta\langle 0|{\cal
O}_{0}^{Q\bar{Q}({}^{3}S_{1}^{[1]})}({}^{3}P^{[8]})|0\rangle^{(1)}$
$\displaystyle=$ $\displaystyle\frac{8\alpha_{s}}{3\pi
m^{2}}\left(\frac{-1}{2\epsilon_{\rm
UV}}\right)\frac{N_{c}^{2}-1}{4N_{c}^{2}}\left(\xi^{\dagger}q^{\prime
i}q^{\prime
l}\sigma^{k}\eta\,\eta^{\dagger}q^{i}q^{l}\sigma^{k}\xi\right)_{{}^{3}S_{1}^{[1]}}$
(82) $\displaystyle=$
$\displaystyle\left[Z^{(1)}({}^{3}P^{[8]}\to{}^{3}S_{1}^{[1]})-1\right]\left(\xi^{\dagger}q^{\prime
i}q^{\prime
l}\sigma^{k}\eta\,\eta^{\dagger}q^{i}q^{l}\sigma^{k}\xi\right)_{{}^{3}S_{1}^{[1]}}.\phantom{xxxxx}$

Then, we have

$\displaystyle\langle 0|{\cal
O}_{0}^{Q\bar{Q}({}^{3}S_{1}^{[1]})}({}^{3}P^{[8]})|0\rangle_{\overline{\rm
MS}}^{(1)}$ $\displaystyle=$ $\displaystyle\langle 0|{\cal
O}_{0}^{Q\bar{Q}({}^{3}S_{1}^{[1]})}({}^{3}P^{[8]})|0\rangle^{(1)}+\delta\langle
0|{\cal O}_{0}^{Q\bar{Q}({}^{3}S_{1}^{[1]})}({}^{3}P^{[8]})|0\rangle^{(1)}$
(83) $\displaystyle=$ $\displaystyle\frac{8\alpha_{s}}{3\pi
m^{2}}\left(\frac{-1}{2\epsilon_{\rm
IR}}\right)\frac{N_{c}^{2}-1}{4N_{c}^{2}}\left(\xi^{\dagger}q^{\prime
i}q^{\prime
l}\sigma^{k}\eta\,\eta^{\dagger}q^{i}q^{l}\sigma^{k}\xi\right)_{{}^{3}S_{1}^{[1]}}.\phantom{xxxxx}$

Now we extract the $S$-wave part by averaging over the angles of $\bm{q}$ and
$\bm{q}^{\prime}$:

$\displaystyle\langle 0|{\cal
O}_{0}^{Q\bar{Q}({}^{3}S_{1}^{[1]})}({}^{3}P^{[8]})|0\rangle_{\overline{\rm
MS}}^{(1)}$ $\displaystyle=$ $\displaystyle\frac{8\alpha_{s}}{3\pi
m^{2}}\left(\frac{-1}{2\epsilon_{\rm
IR}}\right)\frac{N_{c}^{2}-1}{4N_{c}^{2}}\,\frac{1}{d-1}\langle 0|{\cal
O}_{4,1}^{Q\bar{Q}({}^{3}S_{1}^{[1]})}({}^{3}S_{1}^{[1]})|0\rangle^{(0)}.$

####  6.2.2 ${}^{3}S_{1}^{[8]}\to{}^{3}P^{[8]}$

In order $\alpha_{s}$, the ${}^{3}S_{1}^{[8]}$ $Q\bar{Q}$ operator can couple
to the $Q\bar{Q}({}^{3}P^{[8]})$ state. By carrying out a calculation that is
very similar to the one in the preceding section, we find that

$Z^{(1)}({}^{3}S_{1}^{[8]}\to{}^{3}P^{[8]})-1=\frac{8\alpha_{s}}{3\pi
m^{2}}\left(\frac{-1}{2\epsilon_{\rm UV}}\right)\frac{N_{c}^{2}-4}{4N_{c}}$
(85)

and that

$\langle 0|{\cal
O}_{0}^{Q\bar{Q}({}^{3}P^{[8]})}({}^{3}S_{1}^{[8]})|0\rangle_{\overline{\rm
MS}}^{(1)}=\frac{8\alpha_{s}}{3\pi m^{2}}\left(\frac{-1}{2\epsilon_{\rm
IR}}\right)\frac{N_{c}^{2}-4}{4N_{c}}\langle 0|{\cal
O}_{0}^{Q\bar{Q}({}^{3}P^{[8]})}({}^{3}P^{[8]})|0\rangle^{(0)}.\phantom{xxx}$
(86)

The equivalent result for NRQCD decay LDMEs was obtained in Ref.
Petrelli:1997ge .

###  6.3 Order $\alpha_{s}^{2}$

Figure 6:  One of the 36 Feynman diagrams for computing the LDME $\langle
0|\mathcal{O}_{n}^{Q\bar{Q}({}^{3}S_{1}^{[1]})}[Q\bar{Q}({}^{3}S_{1}^{[8]})]|0\rangle^{(2)}$.
The solid circles represent the $Q\bar{Q}$ operators in the LDME. Thirty-five
additional diagrams can be obtained by permuting the gluon-fermion and
operator vertices on the left and right sides of the cut.

In order $\alpha_{s}^{2}$, the ${}^{3}S_{1}^{[8]}$ $Q\bar{Q}$ operator can
couple to the $Q\bar{Q}({}^{3}S_{1}^{[1]})$ state through the diagrams that
are shown in Fig. 6. As in the order-$\alpha_{s}$ case, we expand the
integrands for these diagrams in powers of $1/m$ and, again, only the leading
power contributes in dimensional regularization. Now we combine contributions
that differ in the order of the gluon vertices on a quark or an antiquark line
by making use of the fact that the color factor for the color-singlet
contribution is symmetric under the interchange of the gluon color indices and
by making use of the identity

$\frac{1}{|\bm{k}|}\frac{1}{|\bm{k}|+|\bm{\ell}|}+\frac{1}{|\bm{\ell}|}\frac{1}{|\bm{k}|+|\bm{\ell}|}=\frac{1}{|\bm{k}|}\frac{1}{|\bm{\ell}|}.$
(87)

After we combine the contributions of all of the diagrams in this way, the
loop integrations decouple, and we have

$\displaystyle\langle 0|{\cal
O}_{0}^{Q\bar{Q}({}^{3}S_{1}^{[1]})}({}^{3}S_{1}^{[8]})|0\rangle^{(2)}$
$\displaystyle=$
$\displaystyle\frac{(4\pi\alpha_{s})^{2}}{m^{4}}\left(\frac{\mu^{2}}{4\pi}e^{{{\gamma}}_{\rm
E}}\right)^{2\epsilon}\frac{(N_{c}^{2}-4)(N_{c}^{2}-1)}{8N_{c}^{3}}$ (88)
$\displaystyle\times\left(\xi^{\dagger}q^{\prime i}q^{\prime
r}\sigma^{n}\eta\,\eta^{\dagger}q^{j}q^{s}\sigma^{n}\xi\right)_{{}^{3}S_{1}^{[1]}}$
$\displaystyle\times\int\frac{d^{d-1}k}{(2\pi)^{d-1}}\frac{\delta^{ij}-\hat{\bm{k}}^{i}\hat{\bm{k}}^{j}}{|\bm{k}|^{3}}\int\frac{d^{d-1}\ell}{(2\pi)^{d-1}}\frac{\delta^{rs}-\hat{\bm{\ell}}^{r}\hat{\bm{\ell}}^{s}}{|\bm{\ell}|^{3}}$
$\displaystyle=$ $\displaystyle\frac{1}{2}\left[\frac{8\alpha_{s}}{3\pi
m^{2}}\left(\frac{1}{2\epsilon_{\rm UV}}-\frac{1}{2\epsilon_{\rm
IR}}\right)\right]^{2}\frac{(N_{c}^{2}-1)(N_{c}^{2}-4)}{16N_{c}^{3}}$
$\displaystyle\times\left(\xi^{\dagger}q^{\prime i}q^{\prime
r}\sigma^{n}\eta\,\eta^{\dagger}q^{i}q^{r}\sigma^{n}\xi\right)_{{}^{3}S_{1}^{[1]}},$

where we have made use of Eq. (80).

Figure 7: One of the four Feynman diagrams for computing the one-loop
correction to the counterterm $Z({}^{3}S_{1}^{[8]}\to{}^{3}P^{[8]})-1$. The
symbols $\otimes$ represent the counterterm
$Z({}^{3}S_{1}^{[8]}\to{}^{3}P^{[8]})-1$. Three additional diagrams can be
obtained by permuting the gluon-fermion and counterterm vertices on the left
and right sides of the cut.

We carry out the renormalization of $\langle 0|{\cal
O}_{0}^{Q\bar{Q}({}^{3}S_{1}^{[1]})}({}^{3}S_{1}^{[8]})|0\rangle^{(2)}$ in the
$\overline{\rm MS}$ scheme. First, we add the contribution of the one-loop
diagrams involving the counterterm
$Z^{(1)}({}^{3}P^{[8]}\to{}^{3}S_{1}^{[1]})-1$, which are shown in Fig. 7. The
contribution of these diagrams is

$\delta_{1}\langle 0|{\cal
O}_{0}^{Q\bar{Q}({}^{3}S_{1}^{[1]})}({}^{3}S_{1}^{[8]})|0\rangle^{(2)}=\left[Z^{(1)}({}^{3}S_{1}^{[8]}\to{}^{3}P^{[8]})-1\right]\langle
0|{\cal O}_{0}^{Q\bar{Q}({}^{3}S_{1}^{[1]})}({}^{3}P^{[8]})|0\rangle^{(1)},$
(89)

where $\langle 0|{\cal
O}_{0}^{Q\bar{Q}({}^{3}S_{1}^{[1]})}({}^{3}P^{[8]})|0\rangle^{(1)}$ is given
in Eq. (6.2.1) and $Z^{(1)}({}^{3}S_{1}^{[8]}\to{}^{3}P^{[8]})-1$ is given in
Eq. (85). Thus,

$\displaystyle\delta_{1}\langle 0|{\cal
O}_{0}^{Q\bar{Q}({}^{3}S_{1}^{[1]})}({}^{3}S_{1}^{[8]})|0\rangle^{(2)}$
$\displaystyle=$ $\displaystyle-\frac{8\alpha_{s}}{3\pi
m^{2}}\left(\frac{1}{2\epsilon_{\rm UV}}-\frac{1}{2\epsilon_{\rm
IR}}\right)\frac{(N_{c}^{2}-1)(N_{c}^{2}-4)}{16N_{c}^{3}}$ (90)
$\displaystyle\times\frac{8\alpha_{s}}{3\pi m^{2}}\,\frac{1}{2\epsilon_{\rm
UV}}\left(\xi^{\dagger}q^{\prime i}q^{\prime
r}\sigma^{n}\eta\,\eta^{\dagger}q^{i}q^{r}\sigma^{n}\xi\right)_{{}^{3}S_{1}^{[1]}}.$

Adding this counterterm contribution to Eq. (88), we obtain

$\displaystyle\langle 0|{\cal
O}_{0}^{Q\bar{Q}({}^{3}S_{1}^{[1]})}({}^{3}S_{1}^{[8]})|0\rangle^{(2)}+\delta_{1}\langle
0|{\cal
O}_{0}^{Q\bar{Q}({}^{3}S_{1}^{[1]})}({}^{3}S_{1}^{[8]})|0\rangle^{(2)}$
$\displaystyle=\;\frac{1}{8}\left(\frac{8\alpha_{s}}{3\pi
m^{2}}\right)^{2}\biggl{(}\frac{-1}{\epsilon_{\rm
UV}^{2}}+\frac{1}{\epsilon_{\rm
IR}^{2}}\biggr{)}\frac{(N_{c}^{2}-1)(N_{c}^{2}-4)}{16N_{c}^{3}}\left(\xi^{\dagger}q^{\prime
i}q^{\prime
r}\sigma^{n}\eta\,\eta^{\dagger}q^{i}q^{r}\sigma^{n}\xi\right)_{{}^{3}S_{1}^{[1]}}.\phantom{xxx}$
(91)

We see that the counterterm contribution removes the cross term between the
pole in $\epsilon_{\rm UV}$ and the pole in $\epsilon_{\rm IR}$. Hence, the
remaining overall UV divergence is not coupled to infrared contributions, as
is required for the consistency of the renormalization program. There are no
single poles in $\epsilon_{\rm UV}$. The remaining double pole in
$\epsilon_{\rm UV}$ is removed by adding the $\overline{\rm MS}$-counterterm
contribution

$\displaystyle\delta_{2}\langle 0|{\cal
O}_{0}^{Q\bar{Q}({}^{3}S_{1}^{[1]})}({}^{3}S_{1}^{[8]})|0\rangle^{(2)}$
$\displaystyle=$ $\displaystyle\left(\frac{8\alpha_{s}}{3\pi
m^{2}}\right)^{2}\frac{1}{8\epsilon_{\rm
UV}^{2}}\,\frac{(N_{c}^{2}-1)(N_{c}^{2}-4)}{16N_{c}^{3}}$
$\displaystyle\times\left(\xi^{\dagger}q^{\prime i}q^{\prime
r}\sigma^{n}\eta\,\eta^{\dagger}q^{i}q^{r}\sigma^{n}\xi\right)_{{}^{3}S_{1}^{[1]}}$
$\displaystyle=$
$\displaystyle\left[Z^{(2)}({}^{3}S_{1}^{[8]}\to{}^{3}S_{1}^{[1]})-1\right]\left(\xi^{\dagger}q^{\prime
i}q^{\prime
r}\sigma^{n}\eta\,\eta^{\dagger}q^{i}q^{r}\sigma^{n}\xi\right)_{{}^{3}S_{1}^{[1]}}.$

Adding this counterterm contribution to Eq. (91), we obtain

$\displaystyle\langle 0|{\cal
O}_{0}^{Q\bar{Q}({}^{3}S_{1}^{[1]})}({}^{3}S_{1}^{[8]})|0\rangle_{\overline{\rm
MS}}^{(2)}$ $\displaystyle=$ $\displaystyle\langle 0|{\cal
O}_{0}^{Q\bar{Q}({}^{3}S_{1}^{[1]})}({}^{3}S_{1}^{[8]})|0\rangle^{(2)}+\delta_{1}\langle
0|{\cal
O}_{0}^{Q\bar{Q}({}^{3}S_{1}^{[1]})}({}^{3}S_{1}^{[8]})|0\rangle^{(2)}$ (93)
$\displaystyle+\delta_{2}\langle 0|{\cal
O}_{0}^{Q\bar{Q}({}^{3}S_{1}^{[1]})}({}^{3}S_{1}^{[8]})|0\rangle^{(2)}$
$\displaystyle=$ $\displaystyle\left(\frac{8\alpha_{s}}{3\pi
m^{2}}\right)^{2}\frac{1}{8\epsilon_{\rm
IR}^{2}}\,\frac{(N_{c}^{2}-1)(N_{c}^{2}-4)}{16N_{c}^{3}}$
$\displaystyle\times\left(\xi^{\dagger}q^{\prime i}q^{\prime
r}\sigma^{n}\eta\,\eta^{\dagger}q^{i}q^{r}\sigma^{n}\xi\right)_{{}^{3}S_{1}^{[1]}}.$

We evaluate the $S$-wave part by carrying out the average over the directions
of $\bm{q}$ and $\bm{q}^{\prime}$, with the result

$\displaystyle\langle 0|{\cal
O}_{0}^{Q\bar{Q}({}^{3}S_{1}^{[1]})}({}^{3}S_{1}^{[8]})|0\rangle_{\overline{\rm
MS}}^{(2)}$ $\displaystyle=$ $\displaystyle\frac{1}{8\epsilon_{\rm
IR}^{2}}\left(\frac{8\alpha_{s}}{3\pi
m^{2}}\right)^{2}\frac{(N_{c}^{2}-1)(N_{c}^{2}-4)}{16N_{c}^{3}}$ (94)
$\displaystyle\times\frac{1}{d-1}\langle 0|{\cal
O}_{4,1}^{Q\bar{Q}({}^{3}S_{1}^{[1]})}({}^{3}S_{1}^{[1]})|0\rangle^{(0)}.$

###  6.4 Renormalization-group evolution

In the renormalized LDMEs, we can identify the dimensional-regularization
scale $\mu$ with the NRQCD factorization scale $\mu_{\Lambda}$. We now work
out the renormalization-group evolution with respect to $\mu_{\Lambda}$ of the
LDMEs $\langle 0|{\cal O}_{0}^{Q\bar{Q}({}^{3}S_{1}^{[1]})}$
$({}^{3}P^{[8]})|0\rangle_{\overline{\rm MS}}^{(1)}$, $\langle 0|{\cal
O}_{0}^{Q\bar{Q}({}^{3}P^{[8]})}({}^{3}S_{1}^{[8]})|0\rangle_{\overline{\rm
MS}}^{(1)}$, and $\langle 0|{\cal
O}_{0}^{Q\bar{Q}({}^{3}S_{1}^{[1]})}({}^{3}S_{1}^{[8]})|0\rangle_{\overline{\rm
MS}}^{(2)}$. First, we take $d/d\log{\mu_{\Lambda}}$ of Eq.
(LABEL:3pj3s1msbar):

$\displaystyle\frac{d}{d\log\mu_{\Lambda}}\langle 0|{\cal
O}_{0}^{Q\bar{Q}({}^{3}S_{1}^{[1]})}({}^{3}P^{[8]})|0\rangle_{\overline{\rm
MS}}^{(1)}$ $\displaystyle=$
$\displaystyle\frac{d\alpha_{s}}{d\log\mu_{\Lambda}}\frac{d}{d\alpha_{s}}\langle
0|{\cal
O}_{0}^{Q\bar{Q}({}^{3}S_{1}^{[1]})}({}^{3}P^{[8]})|0\rangle_{\overline{\rm
MS}}^{(1)}$ $\displaystyle=$ $\displaystyle\frac{8\alpha_{s}}{3\pi
m^{2}}\,\frac{N_{c}^{2}-1}{4N_{c}^{2}}\,\frac{1}{d-1}\langle 0|{\cal
O}_{4,1}^{Q\bar{Q}({}^{3}S_{1}^{[1]})}({}^{3}S_{1}^{[1]})|0\rangle^{(0)},$

where we have used

$\frac{d\alpha_{s}}{d\log\mu_{\Lambda}}=-2\epsilon\alpha_{s}+O(\alpha_{s}^{2}).$
(96)

The result in Eq. (LABEL:3pj3s1-evo) agrees with the corresponding results of
Refs. Gremm:1997dq ; z-g-he .

Similarly, taking $d/d\log{\mu_{\Lambda}}$ of Eq. (86), we obtain

$\displaystyle\frac{d}{d\log\mu_{\Lambda}}\langle 0|{\cal
O}_{0}^{Q\bar{Q}({}^{3}P^{[8]})}({}^{3}S_{1}^{[8]})|0\rangle_{\overline{\rm
MS}}^{(1)}$ $\displaystyle=$ $\displaystyle\frac{8\alpha_{s}}{3\pi
m^{2}}\frac{N_{c}^{2}-4}{4N_{c}}\langle 0|{\cal
O}_{0}^{Q\bar{Q}({}^{3}P^{[8]})}({}^{3}P^{[8]})|0\rangle^{(0)}.\phantom{xxx}$
(97)

Taking $d/d\log{\mu_{\Lambda}}$ of Eq. (94), we also obtain

$\displaystyle\frac{d}{d\log\mu_{\Lambda}}\langle 0|{\cal
O}_{0}^{Q\bar{Q}({}^{3}S_{1}^{[1]})}({}^{3}S_{1}^{[8]})|0\rangle_{\overline{\rm
MS}}^{(2)}$ $\displaystyle=$ $\displaystyle\left(\frac{8\alpha_{s}}{3\pi
m^{2}}\right)^{2}\bigg{(}\frac{-1}{2\epsilon_{\rm
IR}}\bigg{)}\frac{(N_{c}^{2}-1)(N_{c}^{2}-4)}{16N_{c}^{3}}$ (98)
$\displaystyle\times\frac{1}{d-1}\langle 0|{\cal
O}_{4,1}^{Q\bar{Q}({}^{3}S_{1}^{[1]})}({}^{3}S_{1}^{[1]})|0\rangle^{(0)}.$

Substituting Eq. (LABEL:3pj3s1msbar) into the right side of Eq. (98) we find
that

$\displaystyle\frac{d}{d\log\mu_{\Lambda}}\langle 0|{\cal
O}_{0}^{Q\bar{Q}({}^{3}S_{1}^{[1]})}({}^{3}S_{1}^{[8]})|0\rangle_{\overline{\rm
MS}}^{(2)}$ $\displaystyle=$ $\displaystyle\frac{8\alpha_{s}}{3\pi
m^{2}}\,\frac{N_{c}^{2}-4}{4N_{c}}\langle 0|{\cal
O}_{0}^{Q\bar{Q}({}^{3}S_{1}^{[1]})}({}^{3}P^{[8]})|0\rangle_{\overline{\rm
MS}}^{(1)}.$ (99)

Equations (97) and (99) agree with the result in Eq. (B19b) of Ref. BBL  at
the leading nontrivial order in $v$ and with the corresponding result in Ref
z-g-he , but disagree with the corresponding result in Ref. Gremm:1997dq .

##  7 NRQCD matching and short-distance coefficients

In the discussion to follow, and in remainder of this paper, the short-
distance coefficients $d_{n}$ that appear are always in the
$\overline{\textrm{MS}}$ scheme. For brevity, we do not indicate the scheme
explicitly.

In relative order $v^{4}$, the NRQCD factorization equation that relates full
QCD and NRQCD for the fragmentation of a gluon into a ${}^{3}S_{1}^{[1]}$
$Q\bar{Q}$ pair is

$\displaystyle D_{4}[g\to Q\bar{Q}({}^{3}S_{1}^{[1]})]^{(3)}$ $\displaystyle=$
$\displaystyle\\{d_{4,1}[g\to Q\bar{Q}({}^{3}S_{1}^{[1]})]^{(3)}+d_{4,2}[g\to
Q\bar{Q}({}^{3}S_{1}^{[1]})]^{(3)}\\}$ (100) $\displaystyle\times\langle
0|\mathcal{O}_{4}^{Q\bar{Q}({}^{3}S_{1}^{[1]})}({}^{3}S_{1}^{[1]})|0\rangle^{(0)}$
$\displaystyle+$ $\displaystyle d_{0}[g\to
Q\bar{Q}({}^{3}P^{[8]})]^{(2)}\langle
0|\mathcal{O}_{0}^{Q\bar{Q}({}^{3}S_{1}^{[1]})}({}^{3}P^{[8]})|0\rangle^{(1)}$
$\displaystyle+$ $\displaystyle d_{0}[g\to
Q\bar{Q}({}^{3}S_{1}^{[8]})]^{(1)}\langle
0|\mathcal{O}_{0}^{Q\bar{Q}({}^{3}S_{1}^{[1]})}({}^{3}S_{1}^{[8]})|0\rangle^{(2)},$

where we have shown the contributions at the leading nontrivial order in
$\alpha_{s}$, namely, $\alpha_{s}^{3}$, and the LDME $\langle
0|\mathcal{O}_{4}^{Q\bar{Q}({}^{3}S_{1}^{[1]})}({}^{3}S_{1}^{[1]})|0\rangle^{(0)}$
means either $\langle
0|\mathcal{O}_{4,1}^{Q\bar{Q}({}^{3}S_{1}^{[1]})}({}^{3}S_{1}^{[1]})|0\rangle^{(0)}$
or $\langle 0|\mathcal{O}_{4,2}^{Q\bar{Q}({}^{3}S_{1}^{[1]})}$
$({}^{3}S_{1}^{[1]})|0\rangle^{(0)}$, since they are equal at the present
order of interest in $v$, as we have explained in Sec. 3. We wish to use this
matching equation to determine the sum of short-distance coefficients
$d_{4,1}[g\to Q\bar{Q}({}^{3}S_{1}^{[1]})]^{(3)}+d_{4,2}[g\to
Q\bar{Q}({}^{3}S_{1}^{[1]})]^{(3)}$. We have already computed the $Q\bar{Q}$
LDMEs on the right side of Eq. (100). In order to fix the short-distance
coefficients in the second and third terms on the right side of Eq. (100), we
make use of a matching equation at relative order $v^{2}$, which is

$\displaystyle D_{2}[g\to Q\bar{Q}({}^{3}P^{[8]})]^{(2)}$ $\displaystyle=$
$\displaystyle d_{0}[g\to Q\bar{Q}({}^{3}P^{[8]})]^{(2)}\langle
0|\mathcal{O}_{0}^{Q\bar{Q}({}^{3}P^{[8]})}({}^{3}P^{[8]})|0\rangle^{(0)}$
(101) $\displaystyle+$ $\displaystyle d_{0}[g\to
Q\bar{Q}({}^{3}S_{1}^{[8]})]^{(1)}\langle
0|\mathcal{O}_{0}^{Q\bar{Q}({}^{3}P^{[8]})}({}^{3}S_{1}^{[8]})|0\rangle^{(1)},$

and a matching equation at relative order $v^{0}$, which is

$\displaystyle D_{0}[g\to Q\bar{Q}({}^{3}S_{1}^{[8]})]^{(1)}$ $\displaystyle=$
$\displaystyle d_{0}[g\to Q\bar{Q}({}^{3}S_{1}^{[8]})]^{(1)}\langle
0|\mathcal{O}_{0}^{Q\bar{Q}({}^{3}S_{1}^{[8]})}({}^{3}S_{1}^{[8]})|0\rangle^{(0)}.$
(102)

First, we solve Eq. (102) for $d_{0}[g\to Q\bar{Q}({}^{3}S_{1}^{[8]})]^{(1)}$,
making use of Eq. (48) for $D_{0}[g\to Q\bar{Q}({}^{3}S_{1}^{[8]})]^{(1)}$ and
Eq. (72c) for $\langle
0|\mathcal{O}_{0}^{Q\bar{Q}({}^{3}S_{1}^{[8]})}({}^{3}S_{1}^{[8]})|0\rangle^{(0)}$.
The result is

$d_{0}[g\to
Q\bar{Q}({}^{3}S_{1}^{[8]})]^{(1)}=\frac{\pi\alpha_{s}}{(d-1)(N_{c}^{2}-1)m^{3}}\left(\frac{\mu^{2}_{\Lambda}}{4\pi}e^{{{\gamma}}_{\rm
E}}\right)^{\epsilon}\delta(1-z).$ (103)

The short-distance coefficient in Eq. (103) agrees with that in Refs.
BL:gfrag-NLO ; Lee:2005jw .

Next, we determine the short-distance coefficient $d_{0}[g\to
Q\bar{Q}({}^{3}P^{[8]})]^{(2)}$ by making use of the order-$v^{2}$ matching
equation (101). We substitute $D_{2}[g\to Q\bar{Q}({}^{3}P^{[8]})]^{(2)}$ in
Eq. (51), $d_{0}[g\to Q\bar{Q}({}^{3}S_{1}^{[8]})]^{(1)}$ in Eq. (103),
$\langle
0|\mathcal{O}_{0}^{Q\bar{Q}({}^{3}P^{[8]})}({}^{3}P^{[8]})|0\rangle^{(0)}$ in
Eq. (72f), and $\langle
0|\mathcal{O}_{0}^{Q\bar{Q}({}^{3}P^{[8]})}({}^{3}S_{1}^{[8]})|0\rangle^{(1)}$
in Eq. (86) into Eq. (101). Then, solving for $d_{0}[g\to
Q\bar{Q}({}^{3}P^{[8]})]^{(2)}$, we obtain

$\displaystyle d_{0}[g\to Q\bar{Q}({}^{3}P^{[8]})]^{(2)}$ $\displaystyle=$
$\displaystyle\frac{8\alpha_{s}^{2}}{3(d-1)(N_{c}^{2}-1)m^{5}}\left(\frac{N_{c}^{2}-4}{4N_{c}}\right)\left(\frac{\mu^{2}_{\Lambda}}{4\pi}e^{{{\gamma}}_{\rm
E}}\right)^{\epsilon}$ (104)
$\displaystyle\times\bigg{\\{}-\frac{\delta(1-z)}{2\epsilon_{\rm
IR}}\left[\frac{(1-\epsilon)\Gamma(1+\epsilon)}{1-\tfrac{2}{3}\epsilon}\left(\frac{\mu^{2}_{\Lambda}}{4m^{2}}e^{{{\gamma}}_{\rm
E}}\right)^{\epsilon}-1\right]$
$\displaystyle+\frac{(1-\epsilon)\Gamma(1+\epsilon)}{1-\tfrac{2}{3}\epsilon}\left(\frac{\mu^{2}_{\Lambda}}{4m^{2}}e^{{{\gamma}}_{\rm
E}}\right)^{\epsilon}f(z)\bigg{\\}},$

where $f(z)$ is given in Eq. (51b). In Ref. Braaten:1996rp , the color-singlet
short-distance coefficients $d_{0}[g\to Q\bar{Q}({}^{3}P_{J}^{[1]})]^{(2)}$
were computed. Summing the results for $d_{0}[g\to
Q\bar{Q}({}^{3}P_{J}^{[1]})]^{(2)}$ in Ref. Braaten:1996rp  over $J=0$, 1 and
2 and multiplying by
$[(N_{c}^{2}-1)/(4N_{c}^{2})]^{-1}[(N_{c}^{2}-4)/(4N_{c})]$ in order to obtain
the corresponding short-distance coefficient for the color-octet channel, we
find agreement with our result in Eq. (104).

The expression in Eq. (104) gives the exact $\epsilon$ dependence. Expanding
this expression to order $\epsilon^{1}$, using the expression for $f(z)$ in
Eq. (5.2), we find that

$\displaystyle d_{0}[g\to Q\bar{Q}({}^{3}P^{[8]})]^{(2)}$
$\displaystyle\quad=\frac{8\alpha_{s}^{2}}{3(d-1)(N_{c}^{2}-1)m^{5}}\,\frac{N_{c}^{2}-4}{4N_{c}}\left(\frac{\mu^{2}_{\Lambda}}{4\pi}e^{{{\gamma}}_{\rm
E}}\right)^{\epsilon}$
$\displaystyle\qquad\times\Bigg{\\{}\frac{1}{6}\delta(1-z)\bigg{[}1-6\log\frac{\mu_{\Lambda}}{2m}+\epsilon\left(\frac{2}{3}-\frac{\pi^{2}}{4}+2\log\frac{\mu_{\Lambda}}{2m}-6\log^{2}\frac{\mu_{\Lambda}}{2m}\right)\bigg{]}$
$\displaystyle\qquad\phantom{xx}+\left(\frac{1}{1-z}\right)_{+}\left[1-\epsilon\left(\frac{1}{3}-2\log\frac{\mu_{\Lambda}}{2m}\right)\right]-2\epsilon\left[\frac{\log(1-z)}{1-z}\right]_{+}$
$\displaystyle\qquad\phantom{xx}+\frac{13-7z}{4}\left[1-\epsilon\,\Big{(}\,\frac{1}{3}+\frac{3}{2}\log(1-z)-2\log\frac{\mu_{\Lambda}}{2m}\,\Big{)}\right]\log(1-z)$
$\displaystyle\qquad\phantom{xx}-\frac{1}{8}(1-2z)(8-5z)\left(1+2\epsilon\log\frac{\mu_{\Lambda}}{2m}\right)$
$\displaystyle\qquad\phantom{xx}+\frac{\epsilon}{24}\Big{[}8+48z-38z^{2}+15(8-11z+4z^{2})\log(1-z)\Big{]}\Bigg{\\}}+O(\epsilon^{2}).$
(105)

Finally, we determine the sum of short-distance coefficients $d_{4,1}[g\to
Q\bar{Q}({}^{3}S_{1}^{[1]})]^{(3)}+d_{4,2}[g\to
Q\bar{Q}({}^{3}S_{1}^{[1]})]^{(3)}$ by making use of the matching equation at
relative order $v^{4}$ [Eq. (100)]. We substitute $D_{4}[g\to
Q\bar{Q}({}^{3}S_{1}^{[1]})]^{(3)}$ in Eq. (68), $d_{0}[g\to
Q\bar{Q}({}^{3}P^{[8]})]^{(2)}$ in Eq. (104), $d_{0}[g\to
Q\bar{Q}({}^{3}S_{1}^{[8]})]^{(1)}$ in Eq. (103), $\langle
0|\mathcal{O}_{4}^{Q\bar{Q}({}^{3}S_{1}^{[1]})}({}^{3}S_{1}^{[1]})|0\rangle^{(0)}$
in Eq. (72e), $\langle
0|\mathcal{O}_{0}^{Q\bar{Q}({}^{3}S_{1}^{[1]})}({}^{3}P^{[8]})|0\rangle^{(1)}$
in Eq. (LABEL:3pj3s1msbar), and $\langle
0|\mathcal{O}_{0}^{Q\bar{Q}({}^{3}S_{1}^{[1]})}({}^{3}S_{1}^{[8]})|0\rangle^{(2)}$
in Eq. (94) into Eq. (100). Solving for $d_{4,1}[g\to
Q\bar{Q}({}^{3}S_{1}^{[1]})]^{(3)}+d_{4,2}[g\to
Q\bar{Q}({}^{3}S_{1}^{[1]})]^{(3)}$, we obtain

$\displaystyle d_{4,1}[g\to Q\bar{Q}({}^{3}S_{1}^{[1]})]^{(3)}+d_{4,2}[g\to
Q\bar{Q}({}^{3}S_{1}^{[1]})]^{(3)}$ $\displaystyle\qquad=\,\,d_{4}[g\to
Q\bar{Q}({}^{3}S_{1}^{[1]})]^{\rm
finite}\,\,+\,\,\frac{2\alpha_{s}^{3}(N_{c}^{2}-4)}{3\pi(d-1)^{3}N_{c}^{3}m^{7}}\,\,\bigg{\\{}\,\,\delta(1-z)\,\bigg{(}\frac{1}{24}-\frac{\pi^{2}}{6}$
$\displaystyle\phantom{xxxxx}-\frac{1}{3}\log\frac{\mu_{\Lambda}}{2m}+\log^{2}\frac{\mu_{\Lambda}}{2m}\bigg{)}+\left(\frac{1}{1-z}\right)_{\\!\\!+}\\!\\!\bigg{(}\frac{1}{3}-2\log\frac{\mu_{\Lambda}}{2m}\bigg{)}+2\left[\frac{\log(1-z)}{1-z}\right]_{\\!+}\\!$
$\displaystyle\phantom{xxxxx}-\frac{104-29z-10z^{2}}{24}+\frac{7[z+(1+z)\log(1-z)]}{2z^{2}}+\frac{(1-2z)(8-5z)}{4}\log\frac{\mu_{\Lambda}}{2m}$
$\displaystyle\phantom{xxxxx}+\frac{1+z}{4}\left(31-6z-\frac{36}{z}\right)\log(1-z)-\frac{z}{4}\left(\\!39-6z+\frac{8}{1-z}\\!\right)\log
z$
$\displaystyle\phantom{xxxxx}+\frac{13-7z}{2}\bigg{[}\,\bigg{(}\\!\log\frac{1-z}{z^{2}}-\log\frac{\mu_{\Lambda}}{2m}\bigg{)}\log(1-z)-\textrm{Li}_{2}(z)\,\bigg{]}\,\,\bigg{\\}}+O(\epsilon),$
(106)

where we have used the identity
$\textrm{Li}_{2}(1-z)=\pi^{2}/6-\textrm{Li}_{2}(z)-\log z\,\log(1-z)$.
$d_{4}[g\to Q\bar{Q}({}^{3}S_{1}^{[1]})]^{\rm finite}$ is defined by

$d_{4}[g\to Q\bar{Q}({}^{3}S_{1}^{[1]})]^{\rm finite}=\frac{D_{4}[g\to
Q\bar{Q}({}^{3}S_{1}^{[1]})]^{\rm finite}}{\langle
0|\mathcal{O}_{4}^{Q\bar{Q}({}^{3}S_{1}^{[1]})}({}^{3}S_{1}^{[1]})|0\rangle^{(0)}}=\frac{D_{4}[g\to
Q\bar{Q}({}^{3}S_{1}^{[1]})]^{\rm finite}}{2(d-1)N_{c}\bm{q}^{4}},$ (107)

where $D_{4}[g\to Q\bar{Q}({}^{3}S_{1}^{[1]})]^{\rm finite}$ and $\langle
0|\mathcal{O}_{4}^{Q\bar{Q}({}^{3}S_{1}^{[1]})}({}^{3}S_{1}^{[1]})|0\rangle^{(0)}$
are given in Eqs. (69) and (72e), respectively. The result in Eq. (7) is new.
As expected, both the double and single poles in $\epsilon_{\rm IR}$ have
cancelled in Eq. (7). These cancellations rely nontrivially on the correctness
of the infrared subtractions and on the NRQCD operator renormalizations in our
calculation.

Finally, we note that the color-singlet short-distance coefficient $d_{0}[g\to
Q\bar{Q}({}^{1}S_{0}^{[1]})]^{(2)}$ was computed in Ref. Braaten:1996rp . We
can obtain the corresponding color-octet short-distance coefficient by
multiplying by the ratio of the color-octet and color-singlet color factors,
namely, $[(N_{c}^{2}-1)/(4N_{c}^{2})]^{-1}_{N_{c}=3}[(N_{c}^{2}-4)/(4N_{c})]$
. Then, we find that

$d_{0}[g\to
Q\bar{Q}({}^{1}S_{0}^{[8]})]^{(2)}=\frac{\alpha_{s}^{2}}{8m^{3}}\frac{N_{c}^{2}-4}{4N_{c}}\left[3z-2z^{2}+2(1-z)\log(1-z)\right].$
(108)

##  8 Numerical results

In this section we describe the numerical results that derive from our
calculations.

We have evaluated $d_{4}[g\to Q\bar{Q}({}^{3}S_{1}^{[1]})]^{\rm finite}(z)$ in
Eq. (107) by carrying out the integrations over the phase space
$d\tilde{\Phi}_{2}$ in Eqs. (67), (70a), and (70b) numerically at $100$ points
in each of the ranges $z=0$ to $z=10^{-2}$, $z=10^{-2}$ to $z=1-10^{-2}$, and
$z=1-10^{-2}$ to $z=1$. We have then used the parametrization in Eq. (B) to
obtain a best fit to the numerical results, which leads to the parameters that
are given in Table 3. Details of this procedure are given in Appendix B. We
have also applied this procedure to $d_{0}[g\to
Q\bar{Q}({}^{3}S_{1}^{[1]})]^{(3)}(z)$ and $d_{2}[g\to
Q\bar{Q}({}^{3}S_{1}^{[1]})]^{(3)}(z)$. In Fig. 8, we plot the results of the
numerical calculations of the short-distance coefficients.

Figure 8: The color-singlet short-distance coefficients $d_{0}[g\to
Q\bar{Q}({}^{3}S_{1}^{[1]})]^{(3)}(z)$, $d_{2}[g\to
Q\bar{Q}({}^{3}S_{1}^{[1]})]^{(3)}(z)$ and $d_{4}[g\to
Q\bar{Q}({}^{3}S_{1}^{[1]})]^{\rm finite}(z)$, which are defined in Eqs. (21),
(22) and (7), respectively, as functions of $z$. The scaling factors are
$(\mathcal{N}_{0},\mathcal{N}_{2},\mathcal{N}_{4})=(10^{-3}\times\alpha_{s}^{3}/m^{3},10^{-2}\times\alpha_{s}^{3}/m^{5},10^{-2}\times\alpha_{s}^{3}/m^{7})$.

We would like to obtain estimates of the relative sizes of the contributions
of the various fragmentation functions to the cross section for $J/\psi$
production in hadron-hadron collisions. The contribution of a fragmentation
process to the $J/\psi$ production cross section, differential in the $J/\psi$
transverse momentum ${{p}}_{T}$, is

$\frac{d\sigma_{J/\psi}^{\rm
frag}}{d{{p}}_{T}}=\int_{0}^{1}\\!dz\,\frac{d\sigma_{g}}{d{{p}}_{T}}({{p}}_{T}/z)\,D(z),$
(109)

where $\frac{d\sigma_{g}}{d{{p}}_{T}}({{p}}_{T})$ is the cross section to
produce a gluon with transverse momentum ${{p}}_{T}$. Let us assume that
$\frac{d\sigma_{g}}{d{{p}}_{T}}({{p}}_{T})\propto 1/{{p}}_{T}^{\kappa}$, where
$\kappa$ is a fixed power. Then,

$\frac{d\sigma_{J/\psi}^{\rm frag}}{d{{p}}_{T}}\propto I_{\kappa}(D),$ (110)

where

$I_{\kappa}(D)=\int_{0}^{1}\\!dz\,z^{\kappa}\,D(z).$ (111)

Hence, we can obtain a rough estimate of the relative contribution of a
fragmentation process to the cross section by computing $I_{\kappa}(D)$.777The
precise calculation of a fragmentation contribution to a cross section would
require that one compute the convolution of the two-to-two partonic cross
sections that produce a final-state gluon with the fragmentation function and
with the appropriate parton distributions. Such a calculation is beyond the
scope of the present paper.

We can estimate $\kappa$ by taking advantage of the fact that $J/\psi$
production in the ${}^{3}S_{1}^{[8]}$ channel is dominated at LO in
$\alpha_{s}$ at large $p_{{}_{T}}$ by gluon fragmentation into a $J/\psi$ with
longitudinal-momentum fraction $z=1$ (Ref. Cho:1995vh ). Consequently, the
$p_{{}_{T}}$ dependence in this channel at high $p_{{}_{T}}$ is approximately
that of $d\sigma_{g}/dp_{{}_{T}}$. Furthermore, since the NLO $k$ factor for
this channel is essentially independent of $p_{{}_{T}}$ and amounts to a
correction of only about $14\%$ at the Tevatron Gong:2008ft , we do not expect
the NLO corrections to change the $p_{{}_{T}}$ dependence of this channel
significantly. Therefore, we estimate $\kappa$ by making use of the result for
the ${}^{3}S_{1}^{[8]}$ contribution to $d\sigma/dp_{{}_{T}}$ at NLO in
$\alpha_{s}$ that appears in Fig. 1(c) of Ref. Butenschoen:2010rq .
Specifically, we compute $(d/d\log p_{{}_{T}})\log[d\sigma/dp_{{}_{T}}]$ and
find that $\kappa\approx 5.2$ at $p_{{}_{T}}=20$ GeV.

We have computed $I_{\kappa}(d)$ for $\kappa=0$ and $5.2$, taking $d$ to be
the order-$v^{0}$, order-$v^{2}$ and order-$v^{4}$ fragmentation functions for
the fragmentation process $g\to Q\bar{Q}({}^{3}S_{1}^{[1]})$. The results are
shown in Table 1.

Table 1:  Results of computing $I_{\kappa}(d)$ for the order-$v^{0}$, order-$v^{2}$, and order-$v^{4}$ fragmentation functions for the fragmentation process $g\to Q\bar{Q}({}^{3}S_{1}^{[1]})$. The first two rows are the normalization integral $I_{0}(d)$ for the short-distance coefficients $d_{0}[g\to Q\bar{Q}({}^{3}S_{1}^{[1]})]^{(3)}$ [Eq. (21)], $d_{2}[g\to Q\bar{Q}({}^{3}S_{1}^{[1]})]^{(3)}$ [Eq. (22)] and $d_{4,1}[g\to Q\bar{Q}({}^{3}S_{1}^{[1]})]^{(3)}+d_{4,2}[g\to Q\bar{Q}({}^{3}S_{1}^{[1]})]^{(3)}$ [Eq. (7)]. The third and fourth rows are $I_{0}(d)$ times the ratio of LDMEs $R_{n}\equiv\langle 0|\mathcal{O}_{n}({}^{3}S_{1}^{[1]})|0\rangle/\langle 0|\mathcal{O}_{0}({}^{3}S_{1}^{[1]})|0\rangle$, as is described in the text. The fifth and sixth rows are the integral $I_{5.2}(d)$. The seventh and eighth rows are $I_{5.2}(d)$ times $R_{n}$. In the case of $d_{4,1}[g\to Q\bar{Q}({}^{3}S_{1}^{[1]})]^{(3)}+d_{4,2}[g\to Q\bar{Q}({}^{3}S_{1}^{[1]})]^{(3)}$, we have set $\mu_{\Lambda}=m$ and $2m$. The quantities in the table are divided by the factor $\mathcal{N}=10^{-4}\times\alpha_{s}^{3}/m^{3}$.  x | $I_{\kappa}(d_{n})\,\backslash\,d_{n}[g\to Q\bar{Q}({}^{3}S_{1}^{[1]})]^{(3)}$ | XXX | $d_{0}$ | XXX | $d_{2}$ | XX |  $d_{4,1}+d_{4,2}$ x
---|---|---|---|---|---|---|---
| $\phantom{{}_{.2}}I_{0}(d_{n})|_{\mu_{\Lambda}=m}$ |  | 8.29 |  | 20.3$/m^{2}$ |  | $-37.4\phantom{9}/m^{4}$
| $\phantom{{}_{.2}}I_{0}(d_{n})|_{\mu_{\Lambda}=2m}$ |  | 8.29 |  | 20.3$/m^{2}$ |  | $-6.54/m^{4}$
| $\phantom{{}_{.2}}I_{0}(d_{n})|_{\mu_{\Lambda}=m\phantom{2}}\times R_{n}$ |  | 8.29 |  | 3.95 |  | $-1.41$
| $\phantom{{}_{.2}}I_{0}(d_{n})|_{\mu_{\Lambda}=2m}\times R_{n}$ |  | 8.29 |  | 3.95 |  | $-0.247$
| $I_{5.2}(d_{n})|_{\mu_{\Lambda}=m}$ |  | 0.743 |  | 11.2$/m^{2}$ |  | $35.7\phantom{9}/m^{4}$
| $I_{5.2}(d_{n})|_{\mu_{\Lambda}=2m}$ |  | 0.743 |  | 11.2$/m^{2}$ |  | $84.5\phantom{9}/m^{4}$
| $I_{5.2}(d_{n})|_{\mu_{\Lambda}=m\phantom{2}}\times R_{n}$ |  | 0.743 |  | 2.18 |  |  -391.35
| $I_{5.2}(d_{n})|_{\mu_{\Lambda}=2m}\times R_{n}$ |  | 0.743 |  | 2.18 |  |  -393.19

In the first and second rows of Table 1, we give the normalization integrals
$I_{0}(d)$ for the short-distance coefficients. In the third and fourth rows,
we give $I_{0}(d)$ multiplied by the relative value of corresponding LDME. In
the fifth and sixth rows we give $I_{5.2}(d)$. In the seventh and eighth rows
we give $I_{5.2}(d)$ multiplied by the relative value of corresponding LDME.
We take the relative values of the LDMEs to be given by the generalized Gremm-
Kapustin relation Bodwin:2006dn :

$\langle 0|\mathcal{O}_{0}^{J/\psi}({}^{3}S_{1}^{[1]})|0\rangle:\langle
0|\mathcal{O}_{2}^{J/\psi}({}^{3}S_{1}^{[1]})|0\rangle:\langle
0|\mathcal{O}_{4}^{J/\psi}({}^{3}S_{1}^{[1]})|0\rangle=1:m^{2}\langle
v^{2}\rangle:m^{4}\langle v^{2}\rangle^{2}.$ (112)

We use the value of $m^{2}\langle v^{2}\rangle$ from Table I of Ref.
Bodwin:2007fz :

$m^{2}\langle
v^{2}\rangle=0.437\,\textrm{GeV}^{2}\,\,\,\textrm{at}\,\,\,m=1.5\,\textrm{GeV}.$
(113)

Examining the relative values in the fifth or sixth rows of Table 1, we see
that, if we exclude the factors from the NRQCD LDMEs, then the order-$v^{4}$
contribution is enhanced considerably relative to the order-$v^{0}$ and
order-$v^{2}$ contributions. We attribute this enhancement to the strong
peaking of $d_{4,1}[g\to Q\bar{Q}({}^{3}S_{1}^{[1]})]^{(3)}+d_{4,2}[g\to
Q\bar{Q}({}^{3}S_{1}^{[1]})]^{(3)}$ near $z=1$, which arises primarily from
the $\delta$ functions and $+$ distributions that it contains. We expect such
strong peaking to occur whenever the full-QCD process contains soft
divergences that correspond to divergences in the $Q\bar{Q}$ LDMEs.

From the seventh and eighth rows of Table 1, we see that the order-$v^{4}$
fragmentation contribution to the cross section is enhanced by approximately a
factor of 2 for $\mu_{\Lambda}=m$ and by approximately a factor of 4 for
$\mu_{\Lambda}=2m$, relative to the order-$v^{0}$ contribution. However, the
order-$v^{0}$ contribution to the $J/\psi$ cross section at ${{p}}_{T}=10$ GeV
lies about a factor of 30 below the measured cross section. Thus, while the
enhancement of the order-$v^{4}$ fragmentation contribution to the cross
section is substantial, it is not sufficient to make the order-$v^{4}$
fragmentation contribution to the cross section important for the
phenomenology of $J/\psi$ production at the current level of precision.

We can also compare the relative contributions to $d\sigma_{J/\psi}^{\rm
frag}$ of each of the three $Q\bar{Q}$ channels in that appear in $D[g\to
J/\psi]$ in order $v^{4}$ [Eq. (24)]. For each channel, we compute
$I_{\kappa}(5.2)$. Then, we multiply by the following LDMEs:

$\displaystyle\langle 0|\mathcal{O}_{4}^{J/\psi}({}^{3}S_{1}^{[1]})|0\rangle$
$\displaystyle=$ $\displaystyle m^{4}\langle v^{2}\rangle^{2}\times
1.32~{}{\rm GeV}^{3},$ (114a) $\displaystyle\langle
0|\mathcal{O}_{0}^{J/\psi}({}^{3}P^{[8]})|0\rangle$ $\displaystyle=$
$\displaystyle-0.109~{}{\rm GeV}^{5},$ (114b) $\displaystyle\langle
0|\mathcal{O}_{0}^{J/\psi}({}^{3}S_{1}^{[8]})|0\rangle$ $\displaystyle=$
$\displaystyle 3.12\times 10^{-3}~{}{\rm GeV}^{3},$ (114c)
$\displaystyle\langle 0|\mathcal{O}_{0}^{J/\psi}({}^{1}S_{0}^{[8]})|0\rangle$
$\displaystyle=$ $\displaystyle 4.50\times 10^{-2}~{}{\rm GeV}^{3},$ (114d)

where $m^{2}\langle v^{2}\rangle$ is given in Eq. (113), and the numerical
values on the right sides of Eq. (114) come from the fit at NLO in
$\alpha_{s}$ to the Tevatron and HERA data in Ref. Butenschoen:2010rq . The
results of this computation are shown in Table 2.

Table 2:  Relative contributions to $d\sigma_{J/\psi}^{\rm frag}$ in order $v^{4}$. The first and second rows give $I_{0}(d)$ for the short-distance coefficients in Eq. (24). The second and third rows give $I_{0}(d)$ times the LDMEs in Eq. (114). The fourth and fifth rows give $I_{5.2}(d)$. The fifth and sixth rows give $I_{5.2}(d)$ times the LDMEs in Eq. (114). We take $m=m_{c}=1.5\,\textrm{GeV}$. For compatibility with Ref. Butenschoen:2010rq , we take $\alpha_{s}=\alpha_{s}({{m}}_{T})$, where ${{m}}_{T}=\sqrt{{{p}}_{T}^{2}+4m_{c}^{2}}$. We choose the point ${{p}}_{T}=20$ GeV, which implies that $\alpha_{s}({{m}}_{T})=0.154$.  $I_{\kappa}(d)\,\backslash\,$channel | $\mathcal{O}_{0}^{J/\psi}({}^{1}S_{0}^{[8]})$ |  | $\mathcal{O}_{0}^{J/\psi}({}^{3}S_{1}^{[8]})$ | $\mathcal{O}_{0}^{J/\psi}({}^{3}P^{[8]})$ | $\mathcal{O}_{4}^{J/\psi}({}^{3}S_{1}^{[1]})$
---|---|---|---|---|---
$\phantom{{}_{.2}}I_{0}(d)|_{\mu_{\Lambda}=m\phantom{2}}\times\,10^{6}\,\textrm{GeV}^{3}$ | $122.$ |  | $5970$ | $-171.\,\textrm{GeV}^{-2}$ | $-0.799\,\textrm{GeV}^{-4}$
$\phantom{{}_{.2}}I_{0}(d)|_{\mu_{\Lambda}=2m}\times\,10^{6}\,\textrm{GeV}^{3}$ | $122.$ |  | $5970$ | $-271.\,\textrm{GeV}^{-2}$ | $-0.140\,\textrm{GeV}^{-4}$
$\phantom{{}_{.2}}I_{0}(d)|_{\mu_{\Lambda}=m\phantom{2}}\times\,10^{6}\,\textrm{LDME}$x | $\phantom{12}5.49$ |  | $\phantom{59}18.6$ | $\phantom{-1}18.6$ | $-0.201$
$\phantom{{}_{.2}}I_{0}(d)|_{\mu_{\Lambda}=2m}\times\,10^{6}\,\textrm{LDME}$ | $\phantom{12}5.49$ |  | $\phantom{59}18.6$ | $\phantom{-1}29.6$ | $-0.0353$
$I_{5.2}(d)|_{\mu_{\Lambda}=m\phantom{2}}\times\,10^{6}\,\textrm{GeV}^{3}$ | $\phantom{1}36.7$ |  | $5970$ | $-300.\,\textrm{GeV}^{-2}$ | $\phantom{-}0.763\,\textrm{GeV}^{-4}$
$I_{5.2}(d)|_{\mu_{\Lambda}=2m}\times\,10^{6}\,\textrm{GeV}^{3}$ | $\phantom{1}36.7$ |  | $5970$ | $-400.\,\textrm{GeV}^{-2}$ | $\phantom{-}1.81\phantom{2}\,\textrm{GeV}^{-4}$
$I_{5.2}(d)|_{\mu_{\Lambda}=m\phantom{2}}\times\,10^{6}\,\textrm{LDME}$ | $\phantom{12}1.65$ |  | $\phantom{59}18.6$ | $\phantom{-1}32.7$ | $\phantom{-}0.192$
$I_{5.2}(d)|_{\mu_{\Lambda}=2m}\times\,10^{6}\,\textrm{LDME}$ | $\phantom{12}1.65$ |  | $\phantom{59}18.6$ | $\phantom{-1}43.6$ | $\phantom{-}0.455$

We see from the second row of Table 2 that the $Q\bar{Q}({}^{3}S_{1}^{[1]})$
channel makes a small contribution to $D_{4}[g\to J/\psi]$ at ${p}_{T}=20$
GeV, confirming our previous conclusion that this channel is not important
phenomenologically at the current level of precision. We also see that the
$Q\bar{Q}({}^{3}S_{1}^{[8]})$ and $Q\bar{Q}({}^{3}P^{[8]})$ channels give
comparable contributions to $D_{4}[g\to J/\psi]$ at ${p}_{T}=20$ GeV. As we
have mentioned, at high $p_{{}_{T}}$ at LO in $\alpha_{s}$, the fragmentation
contribution gives the bulk of the contribution in the
$Q\bar{Q}({}^{3}S_{1}^{[8]})$ channel (Ref. Cho:1995vh ). As we have also
mentioned, the correction to the production rate in the
$Q\bar{Q}({}^{3}S_{1}^{[8]})$ channel at NLO in $\alpha_{s}$ is only about
$14\%$ at the Tevatron Gong:2008ft . Hence, we can estimate the fragmentation
contribution to $J/\psi$ production at the Tevatron in the
$Q\bar{Q}({}^{3}P^{[1]})$ channel by multiplying the ratio of the
$Q\bar{Q}({}^{3}P^{[8]})$ and $Q\bar{Q}({}^{3}S_{1}^{[8]})$ entries in the
seventh row of Table 2 by the value of $d\sigma/d{p}_{T}\times
B(J/\psi\to\mu\mu)$ at $20$ GeV from Fig. 1(c) of Ref. Butenschoen:2010rq  and
by dividing by $1.14$ to account for the NLO correction in the
$Q\bar{Q}({}^{3}S_{1}^{[8]})$ channel. Our estimate is that the fragmentation
contribution to $d\sigma/d{p}_{T}\times B(J/\psi\to\mu\mu)$ at ${p}_{T}=20$
GeV from the $Q\bar{Q}({}^{3}P^{[8]})$ channel is about $6\times 10^{-3}$
nb/GeV. We see that this is comparable to (about a factor of $2$ larger than)
the total NLO contribution in the $Q\bar{Q}({}^{3}P^{[8]})$ channel in Fig.
1(c) of Ref. Butenschoen:2010rq . A more precise calculation of the
fragmentation contribution in the $Q\bar{Q}({}^{3}P^{[8]})$ channel will be
necessary in order to determine whether it is the dominant contribution in
that channel at NLO in $\alpha_{s}$ at high $p_{{}_{T}}$.

##  9 Summary

We have calculated NRQCD short-distance coefficients for gluon fragmentation
into a ${}^{3}S_{1}$ heavy-quarkonium state through relative order $v^{4}$.
Our principal new result is the expression for the sum of the relative-
order-$v^{4}$ short-distance coefficients $d_{4,1}[g\to
Q\bar{Q}({}^{3}S_{1}^{[1]})]^{(3)}+d_{4,2}[g\to
Q\bar{Q}({}^{3}S_{1}^{[1]})]^{(3)}$ for gluon fragmentation through the
${}^{3}S_{1}$ color-singlet channel. This expression is given in Eq. (7) and
in the parametrization of $d_{4}[g\to Q\bar{Q}({}^{3}S_{1}^{[1]})]^{\rm
finite}$ in Eq. (B). As a byproduct of this calculation, we have computed the
short-distance coefficient for gluon fragmentation through the ${}^{3}S_{1}$
color-octet channel [Eq. (103)], finding agreement with the results in Refs.
BL:gfrag-NLO ; Lee:2005jw , and the sum of short-distance coefficients for
gluon fragmentation through the ${}^{3}P_{J}$ color-octet channel [Eq. (104)],
finding agreement with the result in Ref. Braaten:1996rp  for the
${}^{3}P_{J}$ color-singlet short-distance coefficients. We have also computed
the short-distance coefficients for gluon fragmentation through the
${}^{3}S_{1}$ color-singlet channel at leading order in $v$ and at relative-
order $v^{2}$ and find agreement with the results in Refs. Braaten:1993rw ;
Braaten:1995cj  and Ref. Bodwin:2003wh , respectively.

The analysis in this paper involves, for the first time in an NRQCD
factorization calculation, both single and double soft divergences in the
full-QCD process. These soft divergences manifest themselves as single and
double poles in $\epsilon=(4-d)/2$ in dimensional regularization. We have
dealt with the soft divergences by devising subtractions that remove both the
single and double poles in $\epsilon$. We have calculated the phase-space
integrals for the subtraction contributions analytically in $d=4-2\epsilon$
dimensions and have calculated the phase-space integrals for the finite
remainder contributions numerically in $d=4$ dimensions. We have also
calculated the perturbative NRQCD LDMEs that appear in the NRQCD matching
equations for the short-distance coefficients. These perturbative LDMEs
involve both one-loop and two-loop renormalizations of the NRQCD operators.
Our results for the renormalization-group evolution of the renormalized LDMEs
confirm some results in Refs. BBL ; Gremm:1997dq ; z-g-he , but disagree with
one of the results in Ref. Gremm:1997dq .

In Tables 1 and 2, we have given estimates of the relative sizes of the
contributions of the various channels to gluon fragmentation into $J/\psi$
through relative order $v^{4}$. As we have mentioned, the contribution of the
${}^{3}S_{1}$ color-octet channel is believed to be important
phenomenologically at large $p_{{}_{T}}$. Hence, one might expect the
order-$v^{4}$ contribution in the ${}^{3}S_{1}$ color-singlet channel to be
important as well, since the color-singlet channel mixes with the color-octet
channel through single and double logarithms of the NRQCD factorization scale
at order $v^{4}$. Indeed, we find that the contribution to the cross section
at $p_{{}_{T}}=20$ GeV of the ${}^{3}S_{1}$ color-singlet channel in order
$v^{4}$ is about a factor of $2$ larger than the contribution of the
${}^{3}S_{1}$ color-singlet channel at the leading order in $v$ when the
factorization scale is taken to be $\mu_{\Lambda}=m_{c}$ and is about a factor
of $4$ larger than the contribution of the ${}^{3}S_{1}$ color-singlet channel
at the leading order in $v$ when the factorization scale is taken to be
$\mu_{\Lambda}=2m_{c}$. This is a large enhancement, since one would nominally
expect the order-$v^{4}$ contribution to be about $\langle
v^{2}\rangle^{2}\approx 0.04$ times the order-$v^{0}$ contribution. In spite
of this large enhancement of the fragmentation contribution in order $v^{4}$,
the corresponding contribution to the $J/\psi$ production cross section at the
Tevatron or the LHC is not important at the current level of precision of the
phenomenology.

We attribute the large enhancement of the order-$v^{4}$ contribution to the
${}^{3}S_{1}$ color-singlet channel to the peaking of the sum of short-
distance coefficients near $z=1$. This peaking arises from terms that are
proportional to $\delta(1-z)$, $[1/(1-z)]_{+}$, and $[\log(1-z)/(1-z)]_{+}$.
These terms are remnants of the soft divergences that appear in the full-QCD
expression for the fragmentation process and that, ultimately, cancel in the
NRQCD matching equations for the short-distance coefficients. We expect such a
peaking, and the corresponding enhancement, to be present whenever soft
divergences appear in a full-QCD process and are cancelled in the NRQCD
matching equations.888There is also some peaking/enhancement in the
order-$v^{2}$ contribution to the ${}^{3}S_{1}$ color-singlet channel. The
order-$v^{2}$ contribution contains terms that are the product of an
order-$v^{0}$ amplitude with an order-$v^{2}$ amplitude. The order-$v^{2}$
amplitude is sufficiently singular that its square would produce a soft
divergence. However, in combination with the order-$v^{0}$ amplitude, it
produces only a peaking at $z=1$.

The signature of this peaking/enhancement in Table 2 is that the magnitude of
the $I_{5.2}$ entry is comparable to or greater than the magnitude of the
$I_{0}$ entry. We see this signature of peaking/enhancement in the
${}^{3}S_{1}$ and ${}^{3}P_{J}$ color-octet channels, as well as in the
${}^{3}S_{1}$ color-singlet channel in order $v^{4}$, but not in the
${}^{1}S_{0}$ color-octet channel. The ${}^{3}P_{J}$ color-octet channel
contains soft divergences in the full-QCD process, but the ${}^{1}S_{0}$
color-octet channel does not. (The fragmentation contribution in the
${}^{3}S_{1}$ color-octet channel is proportional to $\delta(1-z)$ at leading
order in $\alpha_{s}$.)

The fragmentation contribution to the ${}^{3}S_{1}$ color-octet channel
dominates the contribution of that channel to the $J/\psi$ production cross
section at hadron-hadron colliders at large $p_{{}_{T}}$, both at LO in
$\alpha_{s}$ (Ref. Cho:1995vh ) and at NLO in $\alpha_{s}$ (Ref. Gong:2008ft
). We estimate that the fragmentation contribution to the ${}^{3}P_{J}$ color-
octet channel gives a substantial part of the total contribution of that
channel at NLO in $\alpha_{s}$ to the $J/\psi$ production cross section at
hadron-hadron colliders at large $p_{{}_{T}}$. A more precise calculation of
the ${}^{3}P_{J}$ color-octet fragmentation contribution will be required in
order to determine if it is the dominant contribution at large $p_{{}_{T}}$ in
that channel.

###### Acknowledgements.

We thank Zhi-Guo He for providing us with his unpublished results for the
renormalization-group evolution of NRQCD LDMEs. We also thank Andrea Petrelli
for the use of some of his mathematica code. The work of G.T.B. in the High
Energy Physics Division at Argonne National Laboratory was supported by the U.
S. Department of Energy, Division of High Energy Physics, under Contract No.
DE-AC02-06CH11357. The work of U.R.K. and J.L. was supported by the MEST of
Korea under the NRF Grants No. 2011-0027559 and No. 2011-0003023,
respectively.

##  Appendix A Analytic calculation of phase-space integrals

In this appendix, we give some of the details of the analytic integrations
over the final-state phase space for $S_{12}$, $S_{1}$ and $S_{2}$.

We first carry out the average over the angles of the transverse components of
the final-state gluon momenta in the phase-space (32c). Under this angular
averaging we have

$\displaystyle\langle x\rangle_{\perp}$ $\displaystyle=$
$\displaystyle\frac{z_{1}z_{2}}{z}\left(\frac{e_{1}}{z_{1}}+\frac{e_{2}}{z_{2}}-\frac{1}{z}\right),$
(115a) $\displaystyle\langle x^{2}\rangle_{\perp}$ $\displaystyle=$
$\displaystyle\frac{z_{1}^{2}z_{2}^{2}}{z^{2}}\left(\frac{e_{1}}{z_{1}}+\frac{e_{2}}{z_{2}}-\frac{1}{z}\right)^{2}+\frac{z_{1}z_{2}}{(d-2)z^{2}}\left(2e_{1}-\frac{z_{1}}{z}\right)\left(2e_{2}-\frac{z_{2}}{z}\right).$
(115b)

After this angular averaging, the squared amplitude is independent of the
directions of $\bar{\bm{k}}_{1\perp}$ and $\bar{\bm{k}}_{2\perp}$.

In evaluating the integrals over $|\bar{\bm{k}}_{1\perp}|$ and
$|\bar{\bm{k}}_{2\perp}|$, it is convenient to define dimensionless variables
$u_{1}$ and $u_{2}$:

$\displaystyle u_{1}$ $\displaystyle=$
$\displaystyle\left(\frac{z}{z_{1}}\right)^{2}\bar{k}^{2}_{1\perp},$ (116a)
$\displaystyle u_{2}$ $\displaystyle=$
$\displaystyle\left(\frac{z}{z_{2}}\right)^{2}\bar{k}^{2}_{2\perp}.$ (116b)

Then $e_{1}$, $e_{2}$ and $x$ can be expressed as

$\displaystyle e_{1}$ $\displaystyle=$
$\displaystyle\frac{z_{1}}{2z}(1+u_{1}),$ (117a) $\displaystyle e_{2}$
$\displaystyle=$ $\displaystyle\frac{z_{2}}{2z}(1+u_{2}).$ (117b)

The integrals of $S_{12}$, $S_{1}$ and $S_{2}$ over the magnitudes of the
transverse components of the final-state gluon momenta can be expressed as
linear combinations of the following elementary integrals:

$\displaystyle J_{mn}^{i}$ $\displaystyle=$
$\displaystyle\int\frac{d^{d-2}\bm{\bar{k}}_{i\perp}}{(2\pi)^{d-2}}\frac{1}{e_{i}^{m}(1+2e_{i})^{n}}$
(118) $\displaystyle=$
$\displaystyle\frac{2^{m}}{(4\pi)^{1-\epsilon}\Gamma(1-\epsilon)}\left(\frac{z_{i}}{z}\right)^{2-m-n-2\epsilon}\int_{0}^{\infty}\frac{du_{i}}{u_{i}^{\epsilon}(1+u_{i})^{m}\left(1+\frac{z}{z_{i}}+u_{i}\right)^{n}},\phantom{xxxx}$

where $m$ takes on integer values equal to or greater than $-1$, and $n$ takes
on non-negative integer values. These integrals can be evaluated
straightforwardly to obtain

$\displaystyle J^{i}_{m0}$ $\displaystyle=$ $\displaystyle
2^{m}\left(\frac{z_{i}}{z}\right)^{2-m-2\epsilon}\frac{\Gamma(m-1+\epsilon)}{(4\pi)^{1-\epsilon}\Gamma(m)},$
(119a) $\displaystyle J^{i}_{0n}$ $\displaystyle=$
$\displaystyle\left(\frac{z_{i}}{z}\right)^{1-\epsilon}\left(1+\frac{z_{i}}{z}\right)^{1-n-\epsilon}\frac{\Gamma(n-1+\epsilon)}{(4\pi)^{1-\epsilon}\Gamma(n)},$
(119b) $\displaystyle J^{i}_{-1,2}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(J^{i}_{01}-J^{i}_{02}\right),$ (119c)
$\displaystyle J_{mn}^{i}$ $\displaystyle=$
$\displaystyle(-1)^{m+n}\left(\frac{z_{i}}{z}\right)^{1-\epsilon}\frac{\Gamma(\epsilon)}{(4\pi)^{1-\epsilon}\Gamma(m)\Gamma(n)}$
(119d) $\displaystyle\times\bigg{(}\frac{\partial}{\partial
a}\bigg{)}^{m-1}\bigg{(}\frac{\partial}{\partial
b}\bigg{)}^{n-1}\left\\{\frac{2}{b-2a}\bigg{[}\left(2a+\frac{z_{i}}{z}\right)^{-\epsilon}-\left(b+\frac{z_{i}}{z}\right)^{-\epsilon}\bigg{]}\right\\}\Bigg{|}_{a=0,\,b=1}.\phantom{xxxx}$

After we have carried out the integrations over the transverse components of
the final-state gluon momenta, the remaining integrals over $z_{1}$ and
$z_{2}$ for $S_{12}$ are simple. The integrations over $z_{1}$ and $z_{2}$ for
$S_{1}$ and $S_{2}$ can be expressed as linear combinations of the integrals

$\displaystyle A_{n}$ $\displaystyle=$
$\displaystyle\int_{z_{1}z_{2}}\frac{{z_{2}}^{n-\epsilon}}{{z_{1}}^{1+2\epsilon}(z+z_{2})^{1+\epsilon}},$
(120a) $\displaystyle B_{n}$ $\displaystyle=$
$\displaystyle\int_{z_{1}z_{2}}\frac{{z_{2}}^{n-\epsilon}}{{z_{1}}^{1+2\epsilon}(z+z_{2})^{\epsilon}},$
(120b) $\displaystyle C_{n}$ $\displaystyle=$
$\displaystyle\int_{z_{1}z_{2}}\frac{{z_{2}}^{n-2\epsilon}}{{z_{1}}^{1+2\epsilon}},$
(120c) $\displaystyle D_{n}$ $\displaystyle=$
$\displaystyle\int_{z_{1}z_{2}}\frac{{z_{2}}^{n-2\epsilon}}{{z_{1}}^{1+2\epsilon}(z+z_{2})},$
(120d)

where

$\int_{{z_{1}}{z_{2}}}\\!\\!\\!\\!F({z_{1}},{z_{2}})\,\,=\int_{0}^{1}\\!d{z_{1}}\int_{0}^{1}\\!d{z_{2}}\,\,F({z_{1}},{z_{2}})\,\delta(1-z-{z_{1}}-{z_{2}}).$
(121)

By making use of the identity $z_{2}/(z+z_{2})=1-z/(z+z_{2})$, we derive the
relations

$\displaystyle A_{n}$ $\displaystyle=$ $\displaystyle B_{n-1}-zA_{n-1},$
(122a) $\displaystyle D_{n}$ $\displaystyle=$ $\displaystyle
C_{n-1}-zD_{n-1}.$ (122b)

Applying the recursion relations (122) repeatedly, we can reduce all of the
integrals to the forms $A_{0}$, $B_{n}$ and $C_{n}$, with $1\leq n\leq 4$. It
turns out that the coefficient of $D_{0}$ vanishes. The coefficient of $A_{0}$
is of order $\epsilon^{0}$, while the coefficients of the $B_{n}$ are of order
$\epsilon^{-1}$ or $\epsilon^{0}$. Therefore, we evaluate $A_{0}$ only through
order $\epsilon^{0}$, and we evaluate the $B_{n}$ only through order
$\epsilon$. The expressions for $A_{0}$, $B_{n}$ and $C_{n}$ can be obtained
conveniently by making use of the identity

$\frac{1}{z_{1}^{1+2\epsilon}}=-\frac{(1-z)^{-2\epsilon}}{2\epsilon}\,\delta(z_{1})+\left[\frac{1}{z_{1}^{1+2\epsilon}}\right]_{1-z},$
(123)

which applies when the domain of integration is $0\leq z_{1}\leq 1-z$. The
distribution in the second term of Eq. (123) is defined by

$\int_{0}^{1-z}\\!\\!dz_{1}\,f(z_{1})\left[\frac{1}{z_{1}^{1+2\epsilon}}\right]_{1-z}=\int_{0}^{1-z}\\!\\!dz_{1}\,\frac{f(z_{1})-f(0)}{z_{1}^{1+2\epsilon}}.$
(124)

A straightforward evaluation of the integrals then gives

$\displaystyle A_{0}$ $\displaystyle=$
$\displaystyle-\frac{1}{2\epsilon(1-z)^{3\epsilon}}-\log z+O(\epsilon),$
(125a) $\displaystyle B_{n}$ $\displaystyle=$
$\displaystyle-\frac{(1-z)^{n-3\epsilon}}{2\epsilon}+(1-z)^{n-3\epsilon}\,\big{[}X_{n}+Y_{n}+Z_{n}+O(\epsilon^{2})\big{]},$
(125b) $\displaystyle C_{n}$ $\displaystyle=$
$\displaystyle(1-z)^{n-4\epsilon}\left\\{-\frac{1}{2\epsilon}+\frac{1}{2\epsilon}\left[1-\frac{\Gamma(n+1-2\epsilon)\Gamma(1-2\epsilon)}{\Gamma(n+1-4\epsilon)}\right]\right\\},$
(125c)

where

$\displaystyle X_{n}$ $\displaystyle=$
$\displaystyle\int_{0}^{1}\\!dt\,\,\frac{(1-t)^{n}-1}{t}=-\sum_{k=1}^{n}\frac{1}{k},$
(126a) $\displaystyle Y_{n}$ $\displaystyle=$
$\displaystyle-\epsilon\int_{0}^{1}\\!dt\,\,\frac{(1-t)^{n}}{t}\big{\\{}\log(1-t)+\log[1-(1-z)t]\,\big{\\}},$
(126b) $\displaystyle Z_{n}$ $\displaystyle=$
$\displaystyle-2\epsilon\int_{0}^{1}\\!dt\,\,\frac{(1-t)^{n}-1}{t}\log
t=-2\epsilon\sum_{k=1}^{n}\frac{1}{k}\sum_{\ell=1}^{k}\frac{1}{\ell}.$ (126c)

In Eq. (125c), the first term in the braces is the pole contribution, and the
remainder is finite. The results for $Y_{n}$ for $0\leq n\leq 4$ are

$\displaystyle Y_{0}$ $\displaystyle=$
$\displaystyle\epsilon\Bigg{[}\frac{\pi^{2}}{6}+{\rm Li}_{2}(1-z)\Bigg{]},$
(127a) $\displaystyle Y_{1}$ $\displaystyle=$
$\displaystyle\epsilon\Bigg{[}\frac{\pi^{2}}{6}+{\rm
Li}_{2}(1-z)-2-\frac{z}{1-z}\log z\Bigg{]},$ (127b) $\displaystyle Y_{2}$
$\displaystyle=$ $\displaystyle\epsilon\Bigg{[}\frac{\pi^{2}}{6}+{\rm
Li}_{2}(1-z)-\frac{5-6z}{2(1-z)}-\frac{z(2-3z)\log z}{2(1-z)^{2}}\Bigg{]},$
(127c) $\displaystyle Y_{3}$ $\displaystyle=$
$\displaystyle\epsilon\Bigg{[}\frac{\pi^{2}}{6}+{\rm
Li}_{2}(1-z)-\frac{49-110z+67z^{2}}{18(1-z)^{2}}-\frac{z(6-15z+11z^{2})\log
z}{6(1-z)^{3}}\Bigg{]},\phantom{xxx}$ (127d) $\displaystyle Y_{4}$
$\displaystyle=$ $\displaystyle\epsilon\Bigg{[}\frac{\pi^{2}}{6}+{\rm
Li}_{2}(1-z)-\frac{205-669z+756z^{2}-310z^{3}}{72(1-z)^{3}}$ (127e)
$\displaystyle\quad-\frac{z(12-42z+52z^{2}-25z^{3})\log
z}{12(1-z)^{4}}\Bigg{]}.\phantom{xxxxxx}$

##  Appendix B Parametrizations of the short-distance coefficients

In this appendix, we give parametrizations for the short-distance coefficients
$d_{0}[g\to Q\bar{Q}({}^{3}S_{1}^{[1]})]^{(3)}(z)$, $d_{2}[g\to
Q\bar{Q}({}^{3}S_{1}^{[1]})]^{(3)}(z)$, and $d_{4}[g\to
Q\bar{Q}({}^{3}S_{1}^{[1]})]^{\rm finite}(z)$, which are defined in Eqs. (21),
(22), and (107), respectively.

We observe that $d_{0}[g\to Q\bar{Q}({}^{3}S_{1}^{[1]})]^{(3)}(z)$,
$d_{2}[g\to Q\bar{Q}({}^{3}S_{1}^{[1]})]^{(3)}(z)$, and $d_{4}[g\to Q\bar{Q}$
$({}^{3}S_{1}^{[1]})]^{\rm finite}(z)$ are continuous functions of $z$ over
the whole range $0\leq z\leq 1$ and that they are not analytic at the
endpoints $z=0$ and $1$. One can prove that $d_{n}(0)=0$ for $n=0$, 2, and 4.
In addition, $d_{n}(1)=\alpha_{s}^{3}b_{n}/m^{3+n}$ is finite and calculable
analytically. The values of $b_{n}$ for $n=0$, 2, and 4 are

$\displaystyle b_{0}$ $\displaystyle=$ $\displaystyle 0,$ (128a)
$\displaystyle b_{2}$ $\displaystyle=$
$\displaystyle\frac{22\pi^{2}-15}{4374\pi},$ (128b) $\displaystyle b_{4}$
$\displaystyle=$ $\displaystyle-\frac{2922\pi^{2}-2485}{229635\pi}.$ (128c)

Using this information, we parametrize $d_{0}[g\to
Q\bar{Q}({}^{3}S_{1}^{[1]})]^{(3)}(z)$, $d_{2}[g\to
Q\bar{Q}({}^{3}S_{1}^{[1]})]^{(3)}(z)$, and $d_{4}[g\to
Q\bar{Q}({}^{3}S_{1}^{[1]})]^{\rm finite}(z)$ as follows:

$\displaystyle d_{n}^{\rm fit}(z)=\frac{\alpha_{s}^{3}}{m^{3+n}}$
$\displaystyle\Bigg{[}$ $\displaystyle
b_{n}z+\log(1-z)\sum_{k=1}^{n_{a}}\alpha_{nk}(1-z)^{k}+\log^{2}(1-z)\sum_{k=1}^{n_{b}}\beta_{nk}(1-z)^{k}$
$\displaystyle+\log
z\sum_{k=1}^{n_{c}}\mu_{nk}z^{k}+\log^{2}z\sum_{k=1}^{n_{d}}\nu_{nk}z^{k}+\sum_{k_{1}=1}^{n_{1}}\sum_{k_{2}=1}^{n_{2}}\omega_{nk_{1}k_{2}}z^{k_{1}}(1-z)^{k_{2}}\Bigg{]}.$

In order to fix the parameters in Eq. (B), we have computed $d_{0}[g\to
Q\bar{Q}({}^{3}S_{1}^{[1]})]^{(3)}(z)$, $d_{2}[g\to
Q\bar{Q}({}^{3}S_{1}^{[1]})]^{(3)}(z)$, and $d_{4}[g\to
Q\bar{Q}({}^{3}S_{1}^{[1]})]^{\rm finite}(z)$ by numerical integration for
$100$ points in each of the ranges $z=0$ to $z=10^{-2}$, $z=10^{-2}$ to
$z=1-10^{-2}$, and $z=1-10^{-2}$ to $z=1$. We have carried out the
integrations in two ways: (1) using the change of variables that was proposed
in Refs. Braaten:1993rw ; Braaten:1995cj  and is described in Appendix B of
Ref. Bodwin:2003wh  and (2) using the change of variables that is given in Eq.
(116). At each value of $z$, the numerical results from the two methods of
integration agree to better than $\Delta^{\rm int}_{n}=r_{n}d_{n}(z)$, where
$(r_{0},r_{2},r_{4})=(7.8\times 10^{-4},1.7\times 10^{-3},1.8\times 10^{-3})$.

Table 3: Results of fitting the parametrization in Eq. (B) to numerical values for the color-singlet short-distance coefficients $d_{0}[g\to Q\bar{Q}({}^{3}S_{1}^{[1]})]^{(3)}(z)$, $d_{2}[g\to Q\bar{Q}({}^{3}S_{1}^{[1]})]^{(3)}(z)$, and $d_{4}[g\to Q\bar{Q}({}^{3}S_{1}^{[1]})]^{\rm finite}(z)$, which are defined in Eqs. (21), (22), and (107), respectively. The quantity $\Delta_{\rm max}^{\rm fit}$ is described in the text. Parameters | $d_{0}(z)$ |  | $d_{2}(z)$ |  | $d_{4}(z)^{\rm finite}$
---|---|---|---|---|---
$b_{\phantom{11}}$ | $\phantom{-}0\phantom{11111111111}$ |  | $\phantom{-}1.4710\times 10^{-2}$ |  | $-3.6531\times 10^{-2}$
$\alpha_{1\phantom{1}}$ | $-4.9866\times 10^{-3}$ |  | $-1.8127\times 10^{-2}$ |  | $\phantom{-}5.2157\times 10^{-2}$
$\alpha_{2\phantom{1}}$ | $\phantom{-}9.8448\times 10^{-3}$ |  |  |  | $\phantom{-}2.4588\times 10^{-3}$
$\alpha_{3\phantom{1}}$ | $\phantom{-}1.9512\times 10^{-2}$ |  |  |  |
$\beta_{1\phantom{1}}$ | $-5.1697\times 10^{-6}$ |  | $-1.2282\times 10^{-2}$ |  | $\phantom{-}3.6020\times 10^{-2}$
$\beta_{2\phantom{1}}$ | $\phantom{-}1.0462\times 10^{-2}$ |  |  |  |
$\mu_{1\phantom{1}}$ | $-1.8921\times 10^{-3}$ |  | $\phantom{-}4.1069\times 10^{-3}$ |  | $-7.3565\times 10^{-3}$
$\mu_{2\phantom{1}}$ |  |  | $\phantom{-}1.8341\times 10^{-2}$ |  | $-2.5387\times 10^{-2}$
$\mu_{3\phantom{1}}$ |  |  |  |  | $-9.3477\times 10^{-3}$
$\mu_{4\phantom{1}}$ |  |  |  |  | $-3.6926\times 10^{-3}$
$\nu_{1\phantom{1}}$ | $\phantom{-}1.2154\times 10^{-3}$ |  | $-1.3653\times 10^{-3}$ |  | $\phantom{-}1.3839\times 10^{-3}$
$\nu_{2\phantom{1}}$ | $\phantom{-}1.3039\times 10^{-3}$ |  |  |  |
$\nu_{3\phantom{1}}$ | $-2.7246\times 10^{-3}$ |  |  |  |
$\nu_{4\phantom{1}}$ | $-1.4814\times 10^{-3}$ |  |  |  |
$\omega_{11}$ | $-1.6910\times 10^{-2}$ |  | $-2.1832\times 10^{-2}$ |  | $\phantom{-}7.7264\times 10^{-2}$
$\omega_{12}$ | $\phantom{-}3.8110\times 10^{-2}$ |  | $\phantom{-}4.1531\times 10^{-3}$ |  |
$\omega_{13}$ |  |  | $\phantom{-}8.4949\times 10^{-4}$ |  |
$\omega_{14}$ |  |  | $-3.8207\times 10^{-3}$ |  |
$\Delta_{\rm max}^{\rm fit}\times m^{3+n}/\alpha_{s}^{3}$ |  $3.26\times 10^{-8}\phantom{1}$ |  |  $1.12\times 10^{-6}\phantom{1}$ |  |  $2.05\times 10^{-6}\phantom{1}$

We fit the parametrizations to the numerical integration results, using
$\chi^{2}=\sum_{i}[d(z_{i})-d^{\rm fit}(z_{i})]^{2}/[\sigma(z_{i})]^{2}$ as
the criterion for goodness of fit, where $\sigma(z_{i})$ is the error in the
numerical integration that is given by the VEGAS integration program
Lepage:1977sw . In cases in which it is possible to reduce the number of
parameters in the fit without significantly affecting $\chi^{2}$, we have done
so.

In Table 3, we show the results of this fitting procedure, along with the
value for each fit of $\Delta_{\rm max}^{\rm fit}=\textrm{Max}|d(z_{i})-d^{\rm
fit}(z_{i})|$. In these fits, $\chi^{2}$ per degree of freedom is much larger
than one because the error in the fit is generally considerably larger than
the errors in the numerical integrations that are given by VEGAS. A
conservative estimate of the overall error in a parametrization at each value
of $z$ is $\sqrt{[\Delta_{\rm max}^{\rm fit}]^{2}+[\Delta_{n}^{\rm
int}]^{2}}$.

## References

  * (1) N. Brambilla, S. Eidelman, B.K. Heltsley, R. Vogt, G.T. Bodwin, E. Eichten, A.D. Frawley and A.B. Meyer et al., Heavy quarkonium: progress, puzzles, and opportunities, Eur. Phys. J. C 71 (2011) 1534 [arXiv:1010.5827] [InSPIRE].
  * (2) G.T. Bodwin, E. Braaten and G.P. Lepage, Rigorous QCD analysis of inclusive annihilation and production of heavy quarkonium, Phys. Rev. D 51 (1995) 1125 [ Erratum ibid. D 55 (1997) 5853] [hep-ph/9407339] [InSPIRE].
  * (3) G.C. Nayak, J.-W. Qiu and G.F. Sterman, Fragmentation, factorization and infrared poles in heavy quarkonium production, Phys. Lett. B 613 (2005) 45 [hep-ph/0501235] [InSPIRE].
  * (4) G.C. Nayak, J.-W. Qiu and G.F. Sterman, Fragmentation, nonrelativistic QCD, and NNLO factorization analysis in heavy quarkonium production, Phys. Rev. D 72 (2005) 114012 [hep-ph/0509021] [InSPIRE].
  * (5) M. Butenschoen and B.A. Kniehl, Complete Next-to-Leading-Order Corrections to $J/\psi$ Photoproduction in Nonrelativistic Quantum Chromodynamics, Phys. Rev. Lett. 104 (2010) 072001 [arXiv:0909.2798] [InSPIRE].
  * (6) Y.-Q. Ma, K. Wang and K.-T. Chao, $J/\psi$ ($\psi^{\prime}$) Production at the Tevatron and LHC at $\mathcal{O}(\alpha_{s}^{4}v^{4})$ in Nonrelativistic QCD, Phys. Rev. Lett. 106  (2011) 042002 [arXiv:1009.3655] [InSPIRE].
  * (7) M. Butenschoen and B.A. Kniehl,  Reconciling $J/\psi$ Production at HERA, RHIC, Tevatron, and LHC with nonrelativistic QCD Factorization at Next-to-Leading Order, Phys. Rev. Lett. 106  (2011) 022003 [arXiv:1009.5662] [InSPIRE].
  * (8) M. Butenschoen and B.A. Kniehl, $J/\psi$ production with NRQCD: HERA, Tevatron, RHIC and LHC, AIP. Conf. Proc.  1343 (2011) 409 [arXiv:1011.5619] [InSPIRE].
  * (9) Y.-Q. Ma, K. Wang and K.-T. Chao, Complete next-to-leading order calculation of the $J/\psi$ and $\psi^{\prime}$ production at hadron colliders, Phys. Rev.  D 84 (2011) 114001 [arXiv:1012.1030] [InSPIRE].
  * (10) M. Butenschoen and B.A. Kniehl, World data of $J/\psi$ production consolidate nonrelativistic QCD factorization at next-to-leading order, Phys. Rev.  D 84 (2011) 051501 [arXiv:1105.0820] [InSPIRE].
  * (11) M. Butenschoen and B.A. Kniehl, Probing Nonrelativistic QCD Factorization in Polarized $J/\psi$ Photoproduction at Next-to-Leading Order, Phys. Rev. Lett. 107 (2011) 232001 [arXiv:1109.1476] [InSPIRE].
  * (12) M. Butenschoen and B.A. Kniehl, $J/\psi$ Polarization at the Tevatron and the LHC: Nonrelativistic-QCD Factorization at the Crossroads, Phys. Rev. Lett. 108 (2012) 172002 [arXiv:1201.1872] [InSPIRE].
  * (13) K.-T. Chao, Y.-Q. Ma, H.-S. Shao, K. Wang and Y.-J. Zhang, $J/\psi$ Polarization at Hadron Colliders in Nonrelativistic QCD, Phys. Rev. Lett. 108 (2012) 242004 [arXiv:1201.2675] [InSPIRE].
  * (14) M. Butenschoen and B.A. Kniehl, $J/\psi$ production in NRQCD: A global analysis of yield and polarization, Nucl. Phys. Proc. Suppl. 222 (2012) 151 [arXiv:1201.3862] [InSPIRE].
  * (15) B. Gong, L. -P. Wan, J. -X. Wang and H. -F. Zhang, Polarization for Prompt $J/\psi$ $(\psi^{\prime})$ production at the Tevatron and LHC, [arXiv:1205.6682] [InSPIRE]
  * (16) P. Artoisenet, J.M. Campbell, J.P. Lansberg, F. Maltoni and F. Tramontano, $\Upsilon$ Production at Fermilab Tevatron and LHC Energies, Phys. Rev. Lett. 101 (2008) 152001 [arXiv:0806.3282] [InSPIRE].
  * (17) Z.-B. Kang, J.-W. Qiu and G. Sterman, Factorization and Quarkonium Production, Nucl. Phys. Proc. Suppl. 214 (2012) 39 [InSPIRE].
  * (18) Z.-B. Kang, J.-W. Qiu and G. Sterman, Heavy Quarkonium Production and Polarization, Phys. Rev. Lett. 108 (2012) 102002 [arXiv:1109.1520] [InSPIRE].
  * (19) E. Braaten and T.C. Yuan, Gluon fragmentation into heavy quarkonium, Phys. Rev. Lett. 71 (1993) 1673 [hep-ph/9303205] [InSPIRE].
  * (20) E. Braaten and T.C. Yuan, Gluon fragmentation into spin-triplet $S$-wave quarkonium, Phys. Rev. D 52 (1995) 6627 [hep-ph/9507398] [InSPIRE].
  * (21) G.T. Bodwin and J. Lee, Relativistic corrections to gluon fragmentation into spin-triplet $S$-wave quarkonium, Phys. Rev. D 69 (2004) 054003 [hep-ph/0308016] [InSPIRE].
  * (22) J.C. Collins and D.E. Soper, Parton Distribution and Decay Functions, Nucl. Phys. B 194 (1982) 445 [InSPIRE].
  * (23) M. Gremm and A. Kapustin, Annihilation of $S$-wave quarkonia and the measurement of $\alpha_{s}$, Phys. Lett. B 407 (1997) 323 [hep-ph/9701353] [InSPIRE].
  * (24) L.D. Landau, The moment of a 2-photon system, Dokl. Akad. Nawk., USSR 60 (1948) 207.
  * (25) C.-N. Yang, Selection Rules for the Dematerialization of a Particle into Two Photons, Phys. Rev. 77 (1950) 242 [InSPIRE].
  * (26) G.T. Bodwin and A. Petrelli, Order $v^{4}$ corrections to $S$-wave quarkonium decay, Phys. Rev. D 66 (2002) 094011 [hep-ph/0205210] [InSPIRE].
  * (27) N. Brambilla, E. Mereghetti and A. Vairo, Hadronic quarkonium decays at order $v^{7}$, Phys. Rev. D 79 (2009) 074002 [Erratum ibid. D 83 (2011) 079904] [arXiv:0810.2259] [InSPIRE].
  * (28) G.T. Bodwin, D. Kang and J. Lee, Potential-model calculation of an order-$v^{2}$ nonrelativistic QCD matrix element, Phys. Rev. D 74 (2006) 014014 [hep-ph/0603186] [InSPIRE].
  * (29) E. Braaten and Y.-Q. Chen, Dimensional regularization in quarkonium calculations, Phys. Rev. D 55 (1997) 2693 [hep-ph/9610401] [InSPIRE].
  * (30) E. Braaten and J. Lee, Next-to-leading order calculation of the color-octet ${}^{3}S_{1}$ gluon fragmentation function for heavy quarkonium, Nucl. Phys. B 586 (2000) 427 [hep-ph/0004228] [InSPIRE].
  * (31) J. Lee, Next-to-leading order calculation of a fragmentation function in a light-cone gauge, Phys. Rev. D 71 (2005) 094007 [hep-ph/0504285] [InSPIRE].
  * (32) A.C. Hearn, REDUCE User’s Manual v. 3.7 (The RAND Corporation, Santa Monica, 1999) [E-mail<EMAIL_ADDRESS>
  * (33) R. Mertig, M. Böhm and A. Denner, Feyn Calc – Computer-algebraic calculation of Feynman amplitudes, Comput. Phys. Commun. 64 (1991) 345 [InSPIRE].
  * (34) S. Wolfram, Mathematica, third edition (Wolfram Media, Champaign, Cambridge University Press, New York, 1996).
  * (35) R. Barbieri, R. Gatto and R. Kogerler, Calculation of the annihilation rate Of $P$ wave quark-antiquark bound states, Phys. Lett. B 60 (1976) 183 [InSPIRE].
  * (36) R. Barbieri, R. Gatto and E. Remiddi, Singular binding dependence in the hadronic widths of $1^{++}$ and $1^{+-}$ heavy quark antiquark bound states, Phys. Lett. B 61 (1976) 465 [InSPIRE].
  * (37) C.-H. Chang, Hadronic production of $J/\psi$ associated with a gluon, Nucl. Phys. B 172 (1980) 425 [InSPIRE].
  * (38) B. Guberina, J.H. Kühn, R.D. Peccei and R. Rückl, Rare decays of the $Z^{0}$, Nucl. Phys. B 174 (1980) 317 [InSPIRE].
  * (39) E.L. Berger and D. Jones, Inelastic photoproduction of $J/\psi$ And $\Upsilon$ by gluons, Phys. Rev. D 23 (1981) 1521 [InSPIRE].
  * (40) A. Petrelli, M. Cacciari, M. Greco, F. Maltoni and M.L. Mangano, NLO production and decay of quarkonium, Nucl. Phys. B 514 (1998) 245 [hep-ph/9707223] [InSPIRE].