aainstitutetext: School of Physics, Beihang University, Beijing 100191,
Chinabbinstitutetext: Department of Physics, Liaoning University, Shenyang
110036 , China

# Relativistic Correction to $J/\psi$ and $\Upsilon$ Pair Production

Yi-Jie Li a ,Guang-Zhi Xu b ,Kui-Yong Liu a ,Yu-Jie Zhang<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

The relativistic corrections to the productions of double $J/\psi$ , double
$\Upsilon$, and $J/\psi+\Upsilon$ at the Tevatron and the LHC were
investigated within the frame of nonrelativistic QCD. The ratios of short
distance coefficients between relativistic correction and leading order result
for color singlet and color octet states in large $p_{T}$ limit are
approximately $-1,-11/3$, respectively, for pair production. And for
$J/\psi+\Upsilon$ process the ratio is $-11/6$ for the CO channel. The K
factors of relativistic corrections for color singlet and color octet double
$J/\psi$ production are $-0.23$ and $-0.84$ with adopting $v^{2}=0.23$. The
relativistic corrections significantly dilute the difference between the shape
of color-singlet differential cross sections and the color-octet’s at leading
order. Also our results show that the single parton scattering model may be
enough to explain the LHCb data for the $J/\psi$ pair production.

##  1 Introduction

As the most achievable effective field theory to study production cross
sections and decay rates for heavy quarkonium, the nonrelativistic QCD (NRQCD)
Bodwin:1994jh  has been widely accepted. In the NRQCD approach, a heavy
quarkonium can be described as a superposition of a short distance
coefficients $Q\bar{Q}_{n}(^{2s+1}L_{J})$ mutiplied by a long-distance NRQCD
matrix element $<O^{H}_{n}(^{2s+1}L_{J})>$, here, s, L, J denote the quantum
numbers for the spin, the orbital angular momentum, the total angular momentum
of the $Q\bar{Q}$ intermediate state, respectively, and n stands for the color
for the heavy-quark-antiquark pair state. The short distance coefficients
indicate the creation or annihilation of a heavy quark pair can be calculated
perturbatively with the expansions by the strong coupling constant
$\alpha_{s}$. And the long-distance NRQCD matrix element is the probability of
the intermediate heavy-quark-antiquark pair to evolve into an observed
quarkonium nonperturbatively. In general, the intermediate state can be in a
conventional color-singlet or a color-octet state. By introducing the color-
octet mechanism, one can get the infrared-safe calculations for the decay
rates of P-wave quarkonium statesHuang:1996fa ; Huang:1996sw ; Huang:1996cs ;
Bodwin:1992ye ; Brambilla:2008zg ; Brambilla:2006ph ; Brambilla:2002nu ;
Brambilla:2001xy  and D-waveHe:2008xb ; He:2009bf ; Fan:2009cj  decay widths.
Besides, one may also explain the Tevatron data on the surplus production of
$J/\psi$ and $\psi^{\prime}$ at large $p_{T}$ by fitting the data and
extracting the color-octet long-distance NRQCD matrix elements.

Leaving aside the successes NRQCD have achieved, there are still puzzles in
the comparisons between the NRQCD predictions and the experimental data. The
initial one is the transversely polarized quarkonium prediction at the
Tevatron failing in confronting with experimental measurements. This puzzle is
called a ”smoking gun” for the color-octet production mechanism. And recent
ones are the large discrepancies between the leading order calculations of
inclusiveCho:1996cg ; Yuan:1996ep ; Liu:2003jj  and exclusiveBraaten:2002fi ;
Liu:2002wq  charmonia productions and the experimental data Abe:2002rb ;
Abe:2004ww ; Aubert:2005tj  at B factories. Now the most achievable efforts to
resolve these puzzles seem to be higher order corrections. The next-to-leading
order (NLO) QCD correctionsZhang:2006ay ; Wang:2011qg ; Zhang:2008gp ;
Zhang:2009ym ; Gong:2009ng ; Gong:2009kp ; Ma:2008gq ; Dong:2011fb ;
Bodwin:2013ys , relativistic correctionsHe:2007te ; Jia:2009np ; He:2009uf ;
Fan:2012dy ; Fan:2012vw , and ${\mathcal{O}}(\alpha_{s}v^{2})$ correction
Dong:2012xx ; Li:2013qp  largely compensate the discrepancies between
theoretical values and experimental measurements at B factory. The B factories
analysis also show that the the color-octet matrix elements may be
overestimated previouslyChao:2012iv ; Zhang:2009ym ; Kniehl:1997fv . And the
contributions of higher order corrections for charmonium prediction and
polarization at hadronic colliders are still significant Campbell:2007ws ;
Gong:2008hk ; Gong:2010bk ; Gong:2008sn ; Gong:2012ah ; Gang:2012js ;
Chao:2012iv ; Ma:2010vd ; Butenschoen:2012px ; Gong:2012ug . There are more
information in the Ref.Brambilla:2010cs ; Brambilla:2004wf .

Most recently, the LHCb measured the $J/\psi$ pair production at a center-of-
mass energy of $\sqrt{s}=7TeV$, with an integrated luminosity of 35.2
$pb^{-1}$ and give the differential cross sections as a function of the
$J/\psi$ pair invariant massAaij:2011yc . These predictions are thought to be
useful to study the color-octet mechanism in quarkonium production by
distinguishing the color-siglet contribution and color-octet contribution at
different $p_{T}$ regionBarger:1995vx ; Qiao:2002rh ; Qiao:2009kg ; Li:2009ug
; Ko:2010xy . However, theoretical predictions within the single parton
scattering mechanism are small. So double parton scattering mechanism was
introduced in Ref.Berezhnoy:2011xy ; Kom:2011bd ; Baranov:2011ch . The
predictions combining the double parton scattering and the single parton
scattering lead a better agreement with data than that with the single parton
scattering alone.

However, as mentioned above, high order calculations including QCD and
relativistic corrections to both CS and CO states are essential for heavy
quarkonium production in the single parton scattering process. A complete next
leading order calculations to double quarkonium production would be helpful to
understand the multiple quarkonium production mechanism, COM mechanism, and
the double parton scattering effects, at hadron colliders.

Based on this motivation, in this paper, we studied the effects of
relativistic corrections to both the CS channel and CO channel of production
of double $J/\psi$ , production of double $\Upsilon$ at the LHC and the
Tevatron. And we also discussed the relativistic corrections to the
$J/\psi+\Upsilon$ process for a comparison.

This paper is organized as follows. In section II the NRQCD factorization and
the standard methods to calculate the short distance coefficients of
relativistic corrections are presented. In section III, we would give the
numerical results and make a discussion. Finally there is a summary.

##  2 The method of the calculations on relativistic corrections in NRQCD

The differential cross section of two heavy quarkonium production can be
discribed by convoluting the cross parton level sections in the parton model
as the following expression:

$\displaystyle
d\sigma\big{(}p+p(\bar{p}){\rightarrow}H_{1}+H_{2}+X\big{)}=\sum_{a,b}{\int}dx_{1}dx_{2}{f_{a/p}(x_{1})}{f_{b/p(\bar{p})}(x_{2})}d\hat{\sigma}(a+b{\rightarrow}H_{1}+H_{2}),$
(1)

where $f_{a(b)/p(\bar{p})}(x_{i})$ is the parton distribution function(PDF),
and $x_{i}$ is the parton momentum fraction denoted the fraction parton
carried from proton or antiproton.

Within the NRQCD framework, the parton level cross sections can be divided
into two parts: short distance coefficients and long distance matrix
elements:Bodwin:1994jh

$\hat{\sigma}(a+b{\rightarrow}H_{1}+H_{2})=\sum_{mn}\frac{F_{mn}(ab)}{m_{Q1}^{d_{m}-4}m_{Q2}^{d_{n}-4}}\langle
0|\mathcal{O}_{m}^{H_{1}}|0\rangle\langle 0|\mathcal{O}_{n}^{H_{2}}|0\rangle.$
(2)

On the right-hand side of the equation, the cross section is expanded to
sensible fock states. $F_{mn}$ noted by the subscript $m,n$, i.e. short
distance coefficients, which can be obtained matching with the perturbative
QCD, describe the process that produces intermediate $Q\bar{Q}$ pairs in short
range before harmonization to the physical meson state. $\langle
0|\mathcal{O}_{m,n}^{H}|0\rangle$, the long distance matrix elements,
represent the harmonization process that $Q\bar{Q}$ evolutes to the CS final
state by emitting soft gluons. The long matrix elements can be extracted from
the experiment dataFan:2009zq ; Guo:2011tz , determined by potential
modelBodwin:2007fz  or lattice calculationsBodwin:1996tg ; Bodwin:2001mk ;
Bodwin:2005gg . $\mathcal{O}_{m,n}^{H}$ are local four fermion operators. The
factor of $m_{Q}^{d_{m,n-4}}$ is introduced to make $F_{mn}$ dimensionless.

In this thesis, our calculations only refer to the S-wave spin-triplet
quarkonium ($J/\psi$,$\Upsilon$) and its wave function looks like

$|H\rangle=\mathcal{O}(1)|Q\bar{Q}(^{3}S_{1}^{[1]})\rangle+\mathcal{O}(v)|Q\bar{Q}(^{3}P_{J}^{[8]})g\rangle+\mathcal{O}(v^{2})|Q\bar{Q}(^{1}S_{0}^{[8]})gg\rangle+\mathcal{O}(v^{2})|Q\bar{Q}(^{3}S_{1}^{[8]})gg\rangle+\dots.$
(3)

So the fock state expansion of Eq.(2) may involve different combinations of
these states in above expressions.

For the computations of short-distance coefficients, there are different
processes contributing to the two meson production at parton level. In this
paper, we just consider the gluon fusion process which is the dominate one of
the parton level processes, i.e.
$g(k_{1})+g(k_{2}){\rightarrow}Q_{1}\bar{Q}_{1}+Q_{2}\bar{Q}_{2}{\rightarrow}H_{1}(p_{1})+H_{2}(p_{2})$.
Then the Lorentz-invariant Mandelstam variables at parton level are writen by
(Here we omit the hats on the top of the $s$, $t$, $u$ for simplification)

$\displaystyle s$ $\displaystyle=$
$\displaystyle(p_{1}+p_{2})^{2}=(k_{1}+k_{2})^{2},$ $\displaystyle t$
$\displaystyle=$ $\displaystyle(k_{1}-p_{1})^{2}=(k_{2}-p_{2})^{2},$
$\displaystyle u$ $\displaystyle=$
$\displaystyle(k_{2}-p_{1})^{2}=(k_{1}-p_{2})^{2},$ (4)

and they satisfy

$s+t+u=M_{H_{1}}^{2}+M_{H_{2}}^{2}.$ (5)

In the LO computation, the momenta of the heavy quark/antiquark in the
intermediate $Q\bar{Q}$ pair are set to half of the meson momentum. In an
orthodox method to deal with the relativistic corrections, the amplitudes
would be expanded by the rest relative momentum of the $Q\bar{Q}$ around the
quark/antiquark mass. In this paper, we also adhere to this strategy. We first
define the momentum of heavy quark/antiquark as $p_{Q}$/$p_{\overline{Q}}$. In
a arbitrary frame, they can be written as the following expressionsHe:2007te ;
Braaten:1996jt ; Ma:2012ex :

$\displaystyle{p_{Q}}$ $\displaystyle=$ $\displaystyle\frac{1}{2}p+q,$
$\displaystyle{p_{\overline{Q}}}$ $\displaystyle=$
$\displaystyle\frac{1}{2}p-q.$ (6)

where $p$ boosts from the momentum of the quarkonium in rest frame of the
meson, etc. $(2E_{q},\vec{0})$ and q boosts from $(0,\vec{q})$. Here,
$E_{q}=\sqrt{m_{q}^{2}+|\vec{q}|^{2}}$ is the rest energy of quark(antiquark).
$|\vec{q}|=m_{Q}v$ is a lorentz scalar and $v$ is the three relative velocity
of quark or the antiquark in the rest frame. We can estimate $v$ from the
Gremm-Kapustin relation in weak-coupling regimeGremm:1997dq ; Bodwin:2003wh

$v^{2}=v_{1}^{2}=v_{8}^{2}=\frac{M_{Q\overline{Q}}-2m_{Q}^{pole}}{2m_{Q}^{QCD}}.$
(7)

where subscript ”1” and ”8” represents the pair of quark and antiquark is in a
CS or CO state. And $m_{Q}^{QCD}$ is the mass of quark in NRQCD actions and
$m_{Q}^{pole}$ is the pole mass of quark. Taking $M_{Q\overline{Q}}$,
$m_{Q}^{pole}$, and $m_{Q}^{QCD}$ as the input parameters, we can estimate the
value of $v$. For example, in this paper, we select $M_{J/\psi}=3.1GeV$ and
$m_{c}^{QCD}=m_{c}^{pole}=1.39$ , and then get $v^{2}\sim 0.23$. And
$M_{\Upsilon}=9.4GeV$ and $m_{b}^{QCD}=m_{b}^{pole}=4.57$ ,and then get
$v^{2}\sim 0.1$. With the GK relation

$v^{2}=\frac{\langle
0|\mathcal{P}^{H(^{2s+1}L_{J})}|0\rangle}{m_{Q}^{2}\langle
0|\mathcal{O}^{H(^{2s+1}L_{J})}|0\rangle}[1+\mathcal{O}(v^{4})].$ (8)

one may use the above values of $v$ in estimation of the sizes of next leading
order LDMEs.

Taking advantage of the approximate relation $M_{H}=2E_{q}$ and Eqs.(2,5), the
Mandelstam variables may be expanded by the relative momentum $|\vec{q}|$. s
is the beam energy and is $|\vec{q}|$ independence. To expand $t,u$, we write
they in the frame of gluon mass. Then the expressions expanding to next order
may be given as

$\displaystyle t(|\vec{q}_{1}|,|\vec{q}_{2}|)$ $\displaystyle=$ $\displaystyle
2E_{q1}^{2}+2E_{q2}^{2}-\frac{s}{2}-\frac{1}{2}\lambda^{\frac{1}{2}}\big{(}s,(2E_{q1})^{2},(2E_{q2})^{2}\big{)}$
$\displaystyle=$ $\displaystyle
t(0,0)+Ft(2m_{Q1},2m_{Q2})|\vec{q}_{1}|^{2}+Ft(2m_{Q2},2m_{Q1})|\vec{q}_{2}|^{2}+\mathcal{O}(v^{4}),$
$\displaystyle u(|\vec{q}_{1}|,|\vec{q}_{2}|)$ $\displaystyle=$ $\displaystyle
2E_{q1}^{2}+2E_{q2}^{2}-\frac{s}{2}-\frac{1}{2}\lambda^{\frac{1}{2}}\big{(}s,(2E_{q2})^{2},(2E_{q2})^{2}\big{)}$
(9) $\displaystyle=$ $\displaystyle
u(0,0)+Fu(2m_{Q1},2m_{Q2})|\vec{q}_{1}|^{2}+Fu(2m_{Q2},2m_{Q1})|\vec{q}_{2}|^{2}+\mathcal{O}(v^{4}).$

where $\lambda(x,y,z)=x^{2}+y^{2}+z^{2}-2(xy+xz+yz)$ and

$\displaystyle Ft(x,y)$ $\displaystyle=$
$\displaystyle\frac{4\big{\\{}[t(0,0)+u(0,0)](t(0,0)+y^{2})-2y^{2}[t(0,0)+x^{2}]\big{\\}}}{s[s-2(x^{2}+y^{2})]+(x^{2}-y^{2})^{2}]},$
$\displaystyle Fu(x,y)$ $\displaystyle=$
$\displaystyle\frac{4\big{\\{}[t(0,0)+u(0,0)](u(0,0)+y^{2})-2y^{2}[u(0,0)+x^{2}]\big{\\}}}{s[s-2(x^{2}+y^{2})]+(x^{2}-y^{2})^{2}]}.$
(10)

For the case that the final two mesons are identical particals, the expanding
expressions would be simplified to

$\displaystyle
t(\vec{q}_{1},\vec{q}_{2})=t(0,0)-\frac{4\big{[}t(0,0)+4m_{Q}^{2}\big{]}}{s-16m_{Q}^{2}}(|\vec{q}_{1}|^{2}+|\vec{q}_{2}|^{2})+\mathcal{O}(v^{4}),$
$\displaystyle
u(\vec{q}_{1},\vec{q}_{2})=u(0,0)-\frac{4\big{[}u(0,0)+4m_{Q}^{2}\big{]}}{s-16m_{Q}^{2}}(|\vec{q}_{1}|^{2}+|\vec{q}_{2}|^{2})+\mathcal{O}(v^{4}).$
(11)

The short distance coefficients can be evaluated by matching the computations
of perturbative QCD and NRQCD:

$\displaystyle d\hat{\sigma}\Big{|}_{pert~{}QCD}$
$\displaystyle=\sum_{mn}\frac{F_{mn}}{m_{Q1}^{d_{m}-4}m_{Q2}^{d_{n}-4}}\langle
0|\mathcal{O}_{m}^{Q1\bar{Q1}}|0\rangle\langle
0|\mathcal{O}_{n}^{Q2\bar{Q2}}|0\rangle\Big{|}_{pert~{}NRQCD}.$ (12)

The covariant projection operator method should be adopted to compute the
expression on the left-hand side of the equation. Using this method, spin-
singlet and spin-triplet combinations of spinor bilinear in the amplitudes can
be written in covariant form. For instance, we give the spin-triplet
projection,

$\displaystyle\sum_{s\bar{s}}v(s)\bar{u}(\bar{s})\langle\frac{1}{2},s;\frac{1}{2},\bar{s}|1,S_{z}\rangle=\frac{1}{\sqrt{2}(E_{q}+m_{Q})}(\frac{\not{P}}{2}-\not{q}-m_{Q})\frac{\not{P}-2E_{q}}{4E_{q}}\not{\epsilon}(\frac{\not{P}}{2}+\not{q}+m_{Q}),$
(13)

where $\epsilon$ denotes the polarization vector of the spin-triplet state.
Furthermore in our calculation the Dirac spinors are normalized as
$\bar{u}u=-\bar{v}v=2m_{Q}$. Then the differential cross section satisfies:

$\displaystyle
d\hat{\sigma}(g+g{\rightarrow}Q_{1}\bar{Q}_{1}({}^{2s_{1}+1}L_{J1}^{[1,8]})+Q_{2}\bar{Q}_{2}({}^{2s_{2}+1}L_{J2}^{[1,8]}))$
$\displaystyle{\sim}\overline{\sum}|\mathcal{M}(g+g{\rightarrow}(Q_{1}\bar{Q}_{1})(^{2s_{1}+1}L_{J1}^{[1,8]})+(Q_{2}\bar{Q}_{2})(^{2s_{2}+1}L_{J2}^{[1,8]}))|^{2}$
$\displaystyle\times{\langle
0|\mathcal{O}^{H_{1}}(^{2s_{1}+1}L_{J1})|0\rangle\langle
0|\mathcal{O}^{H_{2}}(^{2s_{2}+1}L_{J2})|0\rangle}.$ (14)

where the amplitude $\mathcal{M}$ is calculated by the following expression

$\displaystyle\mathcal{M}(g+g{\rightarrow}(Q_{1}\bar{Q}_{1})(^{2s+1}L_{J}^{[1,8]})+(Q_{2}\bar{Q}_{2})(^{2s+1}L_{J}^{[1,8]}))$
$\displaystyle=\sum_{s_{1}\bar{s}_{1}s_{2}\bar{s}_{2}}\sum_{ijkl}\langle\frac{1}{2},s_{1};\frac{1}{2},\bar{s}_{1}|J_{1},J_{1z}\rangle\langle\frac{1}{2},s_{2};\frac{1}{2},\bar{s}_{2}|J_{2},J_{2z}\rangle\langle
3i;\bar{3j}|1,8a\rangle\langle 3k;\bar{3l}|1,8b\rangle$ $\displaystyle\hskip
28.45274pt\mathcal{A}(a+b{\rightarrow}Q_{1}^{i}+\bar{Q}_{1}^{j}+Q_{2}^{i}+\bar{Q}_{2}^{j}).$

Next, we first expand the amplitude in powers of $|\vec{q}_{1}|$ and
$|\vec{q}_{2}|$ to next order in $v^{2}$ (Here we only give the S-wave case
and it is the only case using in our consequent calculations)

$\displaystyle\mathcal{M}(g+g\rightarrow(Q_{1}\bar{Q}_{1})+(Q_{2}\bar{Q}_{2}))=(\frac{m_{Q1}}{E_{q1}})^{1/2}(\frac{m_{Q2}}{E_{q2}})^{1/2}\mathcal{M}(q_{1},q_{2},E_{q1},E_{q2})$
$\displaystyle=\mathcal{M}(0,0,m_{Q1},m_{Q2})+\big{\\{}\frac{1}{2}q_{1}^{\alpha}q_{1}^{\beta}\frac{\partial^{2}\mathcal{M}}{\partial
q_{1}^{\alpha}\partial
q_{1}^{\beta}}+\frac{1}{2}q_{2}^{\alpha}q_{2}^{\beta}\frac{\partial^{2}\mathcal{M}}{\partial
q_{2}^{\alpha}\partial q_{2}^{\beta}}$
$\displaystyle+\frac{|\vec{q}_{1}|^{2}}{2m_{Q1}}\frac{\partial\big{[}(m_{Q1}/E_{q1})^{1/2}\mathcal{M}\big{]}}{\partial
E_{q1}}+\frac{|\vec{q}_{2}|^{2}}{2m_{Q2}}\frac{\partial\big{[}(m_{Q2}/E_{q2})^{1/2}\mathcal{M}\big{]}}{\partial
E_{q2}}\big{\\}}\big{(}0,0,m_{Q1},m_{Q2}\big{)}$
$\displaystyle+\mathcal{O}(|\vec{q}_{1}|^{4})+\mathcal{O}(|\vec{q}_{2}|^{4})\equiv\mathcal{M}^{(0)}+\mathcal{M}^{(2)}+\mathcal{O}(|\vec{q}_{1}|^{4})+\mathcal{O}(|\vec{q}_{2}|^{4}).$
(15)

The factor $\displaystyle(\frac{m_{Q}}{E_{q}})^{1/2}$ comes from the
relativistic normalization of the Fock state. For the S-wave case,
$q^{\alpha}q^{\beta}=\frac{1}{3}|\vec{q}|^{2}(-g^{\alpha\beta}+\frac{p_{1}^{\alpha}p_{1}^{\beta}}{p_{1}^{2}})=\frac{1}{3}|\vec{q}|^{2}\Pi^{\alpha\beta}$
and the odd-power terms of $q$ vanish. Subsequently, by multiplying the
complex conjugate of the amplitude, the amplitude squared up to the next order
can be obtained:

$\displaystyle|\mathcal{M}|^{2}$ $\displaystyle=$
$\displaystyle\mathcal{M}^{0}\mathcal{M}^{0*}+(\mathcal{M}^{0}\mathcal{M}^{2*}+h.c.)+\mathcal{O}(v^{4})$
(16) $\displaystyle=$ $\displaystyle\mathcal{M}^{0}\mathcal{M}^{0*}$
$\displaystyle+$
$\displaystyle\frac{1}{6}|\vec{q}_{1}|^{2}\Pi^{\alpha\beta}[\frac{\partial^{2}\mathcal{M}}{\partial
q_{1}^{\alpha}\partial
q_{1}^{\beta}}(0,0,m_{Q1},m_{Q2})\mathcal{M}^{0\ast}+h.c.]+\frac{|\vec{q}_{1}|^{2}}{2m_{Q1}}\frac{\partial(m_{Q1}/E_{q1}\mathcal{M})}{\partial
E_{q1}}(0,0,m_{Q1},m_{Q2})$ $\displaystyle+$
$\displaystyle\frac{1}{6}|\vec{q}_{2}|^{2}\Pi^{\alpha\beta}[\frac{\partial^{2}\mathcal{M}}{\partial
q_{2}^{\alpha}\partial
q_{2}^{\beta}}(0,0,m_{Q1},m_{Q2})\mathcal{M}^{0\ast}+h.c.]+\frac{|\vec{q}_{2}|^{2}}{2m_{Q2}}\frac{\partial(m_{Q2}/E_{q2}\mathcal{M})}{\partial
E_{q2}}(0,0,m_{Q1},m_{Q2})$ $\displaystyle+$
$\displaystyle\mathcal{O}(v^{4}).$

Form the expressions, one can see the relativistic corrections contain four
parts to the two quarkonum production. To evaluate them, Eqs. (2 2) would be
used.

###  2.1 $J/\psi$($\Upsilon$) pair production

In Ref.Qiao:2010kn ; Ko:2010xy ; Ko:2010vh , the production of $J/\psi$ pair
was discussed in NRQCD. In the fock state expansion of the cross section, only
two combinations give the dominate contributions and other ones are suppressed
in all the $p_{T}$ region, i.e.
$|^{3}S_{1}^{[1]}\rangle|^{3}S_{1}^{[1]}gg\rangle$ state and
$|^{3}S_{1}^{[8]}g\rangle|^{3}S_{1}^{[8]}g\rangle$. The former item gives
dominate contribution at relative lower $p_{T}$ region and the latter becomes
important along with $p_{T}$ increasing. In this paper, our computations
restrict relativistic corrections of these two combinations. Therefore the
differential cross section in Eq.(2) up to next order in $v$ takes the
following form

$\displaystyle d\hat{\sigma}(g+g\rightarrow J/\psi+J/\psi)$ $\displaystyle=$
$\displaystyle\frac{F({}^{3}S_{1}^{[1]}{}^{3}S_{1}^{[1]})}{m_{c}^{4}}{\langle
0|\mathcal{O}^{J/\psi}({}^{3}S_{1}^{[1]})|0\rangle}{\langle
0|\mathcal{O}^{J/\psi}({}^{3}S_{1}^{[1]})|0\rangle}$ (17) $\displaystyle+$
$\displaystyle
2\times\frac{G({}^{3}S_{1}^{[1]}{}^{3}S_{1}^{[1]})}{m_{c}^{6}}{\langle
0|\mathcal{P}^{J/\psi}({}^{3}S_{1}^{[1]})|0\rangle}{\langle
0|\mathcal{O}^{J/\psi}({}^{3}S_{1}^{[1]})|0\rangle}$ $\displaystyle+$
$\displaystyle\frac{F({}^{3}S_{1}^{[8]}{}^{3}S_{1}^{[8]})}{m_{c}^{4}}{\langle
0|\mathcal{O}^{J/\psi}({}^{3}S_{1}^{[8]})|0\rangle}{\langle
0|\mathcal{O}^{J/\psi}({}^{3}S_{1}^{[8]})|0\rangle}$ $\displaystyle+$
$\displaystyle
2\times\frac{G({}^{3}S_{1}^{[8]}{}^{3}S_{1}^{[8]})}{m_{c}^{6}}{\langle
0|\mathcal{P}^{J/\psi}({}^{3}S_{1}^{[8]})|0\rangle}{\langle
0|\mathcal{O}^{J/\psi}({}^{3}S_{1}^{[8]})|0\rangle}.$

where $F(_{\cdots})$, $G(_{\cdots})$ are the short distance coefficients and
for leading order and relativistic correction respectively. The factor 2 in
the relativistic terms is given because the contributions of relativistic
correction for identified particle are naturally equal. So we can calculate
only the contribution of one $J/\psi$, then multiply by 2 to get the total
correction. The four-Fermion operators $\mathcal{O}^{H}(\Lambda)$ with the
dimension-6 and the dimension-8 four-Fermion operators
$\mathcal{P}^{H}(\Lambda)$ for the relativistic correction operators are
defined asBodwin:1994jh

$\displaystyle\mathcal{O}^{H}({}^{3}S_{1}^{[1]})$ $\displaystyle=$
$\displaystyle\psi^{{\dagger}}\vec{\sigma}\chi\cdot\chi^{{\dagger}}\vec{\sigma}\psi$
$\displaystyle\mathcal{O}^{H}({}^{3}S_{1}^{[8]})$ $\displaystyle=$
$\displaystyle\psi^{{\dagger}}\vec{\sigma}T^{a}\chi\cdot\chi^{{\dagger}}\vec{\sigma}T^{a}\psi$
$\displaystyle\mathcal{P}^{H}({}^{3}S_{1}^{[1]})$ $\displaystyle=$
$\displaystyle\frac{1}{2}[\psi^{{\dagger}}\vec{\sigma}\chi\cdot\chi^{{\dagger}}\vec{\sigma}(-\frac{i}{2}\overleftrightarrow{\mathbf{D}})^{2}\psi+h.c.]$
$\displaystyle\mathcal{P}^{H}({}^{3}S_{1}^{[8]})$ $\displaystyle=$
$\displaystyle\frac{1}{2}[\psi^{{\dagger}}\vec{\sigma}T^{a}\chi\cdot\chi^{{\dagger}}\vec{\sigma}T^{a}(-\frac{i}{2}\overleftrightarrow{\mathbf{D}})^{2}\psi+h.c.]$
(18)

where $\chi$ and $\psi$ are the Pauli spinors describing anti quark creation
and quark annihilation respectively. $\sigma$ is the Pauli matrices and
$\mathbf{D}$ is the gauge-covariant derivative with
$\overleftrightarrow{\mathbf{D}}=\overrightarrow{\mathbf{D}}-\overleftarrow{\mathbf{D}}$.
$T$ is the $SU(3)$ color matrix.

The above analysis is also adequate to the production of $\Upsilon$ pair. We
use FeynArtsKublbeck:1990xc  to generate Feynman diagrams and amplitudes and
FeynCalcMertig:1990an  handle the amplitudes. For double $J/\psi^{\prime}s$
and double $\Upsilon^{\prime}s$ process there are 31 Feynman diagrams for
color-singlet channel and 72 Feynman diagrams for color-octet channel. We will
give the results of the short distance coefficients to leading order and
relativistic corrections in the Appendix. Our leading order results to color
singlet and color octet channel are consistent with the ref Ko:2010xy .

###  2.2 $J/\psi+\Upsilon$ channel

As the discussion in Refs.Ko:2010xy ; Ko:2010vh , the events to double
$J/\psi$ under the luminosity of LHC is 70 and it is difficult to study COM by
use of the process of $J/\psi$ pair production for the LHC experiments as well
as for the di-$\Upsilon$ pair production. They study the process
$pp\rightarrow{J/\psi+\Upsilon+X}$ and find it is nice probe to COM with
enough events at LHC. To study relativistic corrections to this process is a
meaningful work. At leading order of $O(\alpha_{s})$, the CS channel
$c\bar{c}(^{3}S_{1}^{[1]})+b\bar{b}(^{3}S_{1}^{[1]})$ have no contributions.
The CO channels include $c\bar{c}(^{3}S_{1}^{[1]})+b\bar{b}(^{3}S_{1}^{[8]})$,
$c\bar{c}(^{3}S_{1}^{[8]})+b\bar{b}(^{3}S_{1}^{[1]})$ and
$c\bar{c}(^{3}S_{1}^{[8]})+b\bar{b}(^{3}S_{1}^{[8]})$ combinations. In our
calculations, we neglect the relativistic corrections to $\Upsilon$ production
which is negligible compared with the corrections to $J/\psi$. Therefore we
have the obvious expression to next order in $v$

$\displaystyle d\hat{\sigma}(g+g\rightarrow J/\psi+\Upsilon)$ $\displaystyle=$
$\displaystyle\frac{F({}^{3}S_{1}^{[1]}{}^{3}S_{1}^{[8]})}{m_{c}^{2}m_{b}^{2}}{\langle
0|\mathcal{O}^{J/\psi}({}^{3}S_{1}^{[1]})|0\rangle}{\langle
0|\mathcal{O}^{\Upsilon}({}^{3}S_{1}^{[8]})|0\rangle}$ (19) $\displaystyle+$
$\displaystyle\frac{G({}^{3}S_{1}^{[1]}{}^{3}S_{1}^{[8]})}{m_{c}^{4}m_{b}^{2}}{\langle
0|\mathcal{P}^{J/\psi}({}^{3}S_{1}^{[1]})|0\rangle}{\langle
0|\mathcal{O}^{\Upsilon}({}^{3}S_{1}^{[8]})|0\rangle}$ $\displaystyle+$
$\displaystyle\frac{F({}^{3}S_{1}^{[8]}{}^{3}S_{1}^{[1]})}{m_{c}^{2}m_{b}^{2}}{\langle
0|\mathcal{O}^{J/\psi}({}^{3}S_{1}^{[8]})|0\rangle}{\langle
0|\mathcal{O}^{\Upsilon}({}^{3}S_{1}^{[1]})|0\rangle}$ $\displaystyle+$
$\displaystyle\frac{G({}^{3}S_{1}^{[8]}{}^{3}S_{1}^{[1]})}{m_{c}^{4}m_{b}^{2}}{\langle
0|\mathcal{P}^{J/\psi}({}^{3}S_{1}^{[8]})|0\rangle}{\langle
0|\mathcal{O}^{\Upsilon}({}^{3}S_{1}^{[1]})|0\rangle}$ $\displaystyle+$
$\displaystyle\frac{F({}^{3}S_{1}^{[8]}{}^{3}S_{1}^{[8]})}{m_{c}^{2}m_{b}^{2}}{\langle
0|\mathcal{O}^{J/\psi}({}^{3}S_{1}^{[8]})|0\rangle}{\langle
0|\mathcal{O}^{\Upsilon}({}^{3}S_{1}^{[8]})|0\rangle}$ $\displaystyle+$
$\displaystyle\frac{G({}^{3}S_{1}^{[8]}{}^{3}S_{1}^{[8]})}{m_{c}^{4}m_{b}^{2}}{\langle
0|\mathcal{P}^{J/\psi}({}^{3}S_{1}^{[8]})|0\rangle}{\langle
0|\mathcal{O}^{\Upsilon}({}^{3}S_{1}^{[8]})|0\rangle}.$

The operators are defined in Eq.(2.1) and the results of the short distance
coefficients can also be found in the Appendix.

##  3 The Numerical Result and Discussion

Using the method of phase space integration in Ref. Fan:2009zq , we can carry
out the numerical calculation the cross section at the LHC and Tevatron
respectively. Here we ignore the relativistic correction come from the phase
space. Fortran program was used to calculate the integration. We choose input
parameters: the mass of charm quark $m_{c}=1.5GeV$ and bottom quark
$m_{b}=4.7GeV$; the factorization scale $\mu$ taken as transverse mass
$\mu=m_{T}=(4{m_{Q}}^{2}+{p_{T}}^{2})$ ; specifically we use the parton
distribution function CTEQ6L1Pumplin:2002vw . The running coupling const
$\alpha_{s}$ is evaluated by the formula of CTEQ6L1. We employ the following
values for the NRQCD matrix elements for $J/\psi$ and $\Upsilon$ Bodwin:2007fz
; Braaten:1999qk ; Kang:2007uv ; Kramer:2001hh :

$\displaystyle\langle 0|\mathcal{O}^{J/\psi}({}^{3}S_{1}^{[1]})|0\rangle$
$\displaystyle=$ $\displaystyle 1.3GeV^{3},$ $\displaystyle\langle
0|\mathcal{O}^{J/\psi}({}^{3}S_{1}^{[8]})|0\rangle$ $\displaystyle=$
$\displaystyle 3.9{\times}10^{-3}GeV^{3}.$ (20) $\displaystyle\langle
0|\mathcal{O}^{\Upsilon}({}^{3}S_{1}^{[1]})|0\rangle$ $\displaystyle=$
$\displaystyle 9.2GeV^{3},$ $\displaystyle\langle
0|\mathcal{O}^{\Upsilon}({}^{3}S_{1}^{[8]})|0\rangle$ $\displaystyle=$
$\displaystyle 1.5{\times}10^{-1}GeV^{3}.$ (21)

As discussion in Sec.2, we use Eq.(8) to estimate the relativistic LDMES.

###  3.1 $J/\psi$($\Upsilon$) pair production

Under the large $p_{T}$ limit approximation,

$\displaystyle-\frac{m^{2}}{u}\sim-\frac{m^{2}}{t}<\frac{m^{2}}{p_{T}^{2}}\sim
0$ (22)

where $m$ is the mass of $J/\psi$ or $\Upsilon$, we can expand the short
distance coefficients with $m$. The ratios of the short distance coefficient
between LO (F) and its relativistic correction (G) in lowest order in the
expansion are:

$\displaystyle R({}^{3}S_{1}^{[1]}{}^{3}S_{1}^{[1]})\Big{|}_{p_{T}\gg
m}=\frac{G({}^{3}S_{1}^{[1]}{}^{3}S_{1}^{[1]})}{F({}^{3}S_{1}^{[1]}{}^{3}S_{1}^{[1]})}\Big{|}_{p_{T}\gg
m}$ $\displaystyle\sim$ $\displaystyle-1$ $\displaystyle
R({}^{3}S_{1}^{[8]}{}^{3}S_{1}^{[8]})\Big{|}_{p_{T}\gg
m}=\frac{G({}^{3}S_{1}^{[8]}{}^{3}S_{1}^{[8]})}{F({}^{3}S_{1}^{[8]}{}^{3}S_{1}^{[8]})}\Big{|}_{p_{T}\gg
m}$ $\displaystyle\sim$ $\displaystyle-\frac{11}{3}$ (23)

The ratio of ${}^{3}S_{1}^{[8]}$ is approximately $-11/3$ this ratio of color
octet state is consistent with Ref.Bodwin:2003wh ; Xu:2012am . In large
$p_{T}$ limit, the dominate contribution of this subprocess is $g^{*}\to
c\bar{c}({}^{3}S_{1}^{[8]})$. The propagator of virtual gluon $g^{*}$ is
proportional to $1/E_{q}^{2}$. This term offers a factor of $-2$ to the ratio
$R({}^{3}S_{1}^{[8]})$.

In Fig.1, we present the ratios at the Tevatron with $\sqrt{s}=1.96TeV$ and at
the LHC with $\sqrt{s}=7TeV$ or $\sqrt{s}=14TeV$ of $J/\psi$ and the curves
are very close at large $p_{T}$ and consistent with above analysis. In small
$p_{T}$ region, for CS, corrections are less than the values in large $p_{T}$
limit. On contrary, there is a sharp suppress in small $p_{T}$ region for CO.
That can obviously be shown in Fig.2, etc. differential cross sections of
$p_{T}$. The LO cross sections to CS and CO state suppressed by factor of
$23\%$ and $84\%$ when involving the effects of relativistic corrections. For
the LO case, the CS and CO curves cross at about $p_{T}=14GeV$. As $p_{T}$
increasing the contributions to CO state would surpass that to CS. And this
provide a probe to check the COM experimentally. But in our calculations as
shown in the diagrams, relativistic corrections reduce the CO contributions
largely and make the difference not apparently between the CSM results and the
COM results. Furthermore, we make a prediction of the cross sections at the
LHCb with $\sqrt{s}=8TeV$ as shown in Fig.3, COM with relativistic corrections
contribute little to total cross sections that is similar with the cases
mentioned above. For total cross sections, COM effects are also not obviously,
and this can be seen in Tab.1.

Figure 1:  The ratios of the short distance coefficient between LO $F$ and its
relativistic correction G of ${p+p}\rightarrow{J/\psi+J/\psi}$ process at the
Tevatron with $\sqrt{s}=1.96TeV$ and at the LHC with $\sqrt{s}=7TeV$ of
$\sqrt{s}=14TeV$.

Figure 2: The ${p+p}\rightarrow{J/\psi+J/\psi}$ process $p_{t}$ distribution
of deferential cross section at the LHC and Tevatron respectively.

Figure 3: The predicted ${p+p}\rightarrow{J/\psi+J/\psi}$ process $p_{t}$ distribution of deferential cross section at the LHC with $\sqrt{s}=8TeV$ of $|y|<2.4$ and $2.5<y<4.5$ , respectively. $m_{c}$ |  ${}^{3}{S_{1}}^{[1]}$ LO |  ${}^{3}{S_{1}}^{[1]}$ RC |  ${}^{3}{S_{1}}^{[8]}$ LO |  ${}^{3}{S_{1}}^{[8]}$ RC | total cross section
---|---|---|---|---|---
(GeV) | (nb) | (pb) | (pb) | (pb) | (nb)
1.60 | 5.18 | $-13.90$ | $1.63$ | $-1.94$ | 5.16
1.55 | 6.92 | $-22.12$ | $2.20$ | $-2.61$ | 6.89
1.50 | 9.31 | $-34.68$ | $3.00$ | $-3.55$ | 9.28
1.45 | 12.65 | $-53.09$ | $4.13$ | $-4.88$ | 12.60
1.40 | 17.34 | $-80.96$ | $5.73$ | $-6.77$ | 17.26
Table 1:  The total cross sections for color-singlet channel, color-octet
channel and their relativistic corrections at LHC with $\sqrt{s}=7TeV$,
integrating over $2<y<4.5$ and $p_{T}<10GeV$ and we choose charm quark mass
for input parameter from $1.4GeV$ to $1.6GeV$. ”RC” at the head of the table
denotes the relativistic contributions.

In comparison of process of single $J/\psi$ production in which QCD
corrections are essential than relativistic corrections to CO states and give
large contributions to total cross sections. It is expected for QCD
corrections to CS and CO states to give a large contributions.

To investigate the effects of relativistic corrections comparing with
experiment, we give the differential cross sections as functions of the
invariant mass of double $J/\psi$ and of the rapidity in Fig.4 to pair
production with the same constraints as the LHCb, i.e. $2<y<4.5$ and
$\sqrt{s}=7TeV$. From the figures, one can see the main part of the
theoretical calculations concentrates in the low region of the invariant mass
with a sharp peak at about $6.5GeV$, but LHCb data shows a peak at $7\sim
8GeV$. And theoretical value of the peak is nearby $5.5nb/GeV$ which is two
times large than the experimental peak value. In the low region less than
about $7.5GeV$, relativistic effects show a positive correction to LO curve
with a maximal enhanced factor 1.08. In the region large than $7.5GeV$,
relativistic effects suppress the LO curve but not obviously. Here we only
consider the contributions to direct pair production and the feed-down process
from excited state will contribute nearly $30\%$ to direct production
rateBerezhnoy:2011xy . Taking the feed-down corrections into account, the
theoretical prediction to differential cross section within the invariant mass
$6GeV\sim 8.5GeV$ would be large than the experimental values. This situation
is very different from the calculations in Refs.Berezhnoy:2011xy ; Kom:2011bd
. They showed that the single parton scattering contributions were too small
to fit the experimental data and introduced double parton scattering
mechanism. In our calculations, only LO contributions would be large enough
comparing with the experimental data. But as they discussed, single parton
scattering result peaks at too low a invariant mass value and with a too large
value, and the relativistic corrections do not relieve the puzzle but enhances
the peak value. Recently, an soft-collinear effective theory(SCET) has been
developed and seem to provide a possibility to suppress the end point peak
considerably and move it from the end pointBauer:2000yr . Besides, may NLO QCD
redia corrections give a large negative contribution especially for the peak?
It is expected the problem would be solved when involving these corrections.

Figure 4:  Left figure is the differential cross section to double $J/\psi$
production as a function of the invariant mass $m_{J/\psi J/\psi}$ at LHC with
$\sqrt{s}=7TeV$ and $2<y<4.5$, and dash and solid lines correspond to LO and
NLO($v^{2}$) results, respectively. Right figure is for rapidity distribution,
dash and solid lines also correspond to LO and NLO($v^{2}$), respectively. In
both two figures, it were set to $M_{J/\psi}=3.1$.

Figure 5:  The figure is about the differential cross section to double
$J/\psi$ production over rapidities of $J/\psi J/\psi$ at LHC with
$\sqrt{s}=7TeV$ and $2<y<4.5$,and $\sqrt{s}=8TeV$ and $2.5<y<4.5$
respectively.

Figure 6:  The deferential cross section of
${p+p}\rightarrow{\Upsilon+\Upsilon}$ process as a function of $p_{t}$ at the
LHC with $|\eta|<2.4$ and $\sqrt{s}=7TeV$ , $\sqrt{s}=14TeV$, respectively.

###  3.2 $J/\psi+\Upsilon$ channel

We give the differential cross sections as a function of $p_{t}$ at the LHC
with $|\eta|<2.4$ when choosing the beam energy of the LHC and the LHCb,
$\sqrt{s}=7TeV$ and $\sqrt{s}=14TeV$, respectively.

As discussed in Sec.2 we consider only the relativistic corrections of the
$J/\psi$ part. The radios of $F/G$ in the large $p_{T}$ limit are given below

$\displaystyle
R[{c\overline{c}}(^{3}S_{8}^{[8]}),{b\overline{b}}(^{3}S_{1}^{[8]})]\Big{|}_{p_{T}\gg
m}=\frac{G[{c\overline{c}}(^{3}S_{1}^{[8]}),{b\overline{b}}(^{3}S_{1}^{[8]})]}{F[{c\overline{c}}(^{3}S_{1}^{[8]}),{b\overline{b}}(^{3}S_{8}^{[8]})]}\Big{|}_{p_{T}\gg
m}$ $\displaystyle\sim$ $\displaystyle-\frac{11}{6}$ $\displaystyle
R[{c\overline{c}}(^{3}S_{8}^{[8]}),{b\overline{b}}(^{3}S_{1}^{[1]})]\Big{|}_{p_{T}\gg
m}=\frac{G[{c\overline{c}}(^{3}S_{1}^{[8]}),{b\overline{b}}(^{3}S_{1}^{[1]})]}{F[{c\overline{c}}(^{3}S_{1}^{[8]}),{b\overline{b}}(^{3}S_{1}^{[1]})]}\Big{|}_{p_{T}\gg
m}$ $\displaystyle\sim$ $\displaystyle-\frac{11}{6}$

For the channel ${c\bar{c}}(^{3}S_{8}^{[1]})+{b\bar{b}}(^{3}S_{1}^{[8]})$, the
ratio is not a const result and its value approximate 1.5 at the large $p_{T}$
region. In Fig.7, we give the $P_{T}$ distribution of the LO cross sections
and the next leading order ones for different channels. The relativistic
corrections for the production of $J/\psi+\Upsilon$ are limited to change the
LO CO contributions. This process may be still absorbing to research to be a
probe to COM.

Figure 7: The differential cross section of
${p+p}\rightarrow{J/\psi+\Upsilon}$ process as a function of $p_{t}$ at the
LHC with $\sqrt{s}=7TeV$ and $\sqrt{s}=14TeV$ with $|\eta|<2.4$ .

##  4 Summary

We calculate the relativistic corrections to the process
${p+p(\bar{p})}\rightarrow{J/\psi+J/\psi+X}$
,${p+p(\bar{p})}\rightarrow{\Upsilon+\Upsilon+X}$, and
${p+p(\bar{p})}\rightarrow{J/\psi+\Upsilon+X}$ at the Tevatron and at the LHC
based on NRQCD. The ratios of short distance coefficients between relativistic
corrections and LO results for CS and CO states are approximately
$-1,-\frac{11}{3}$ in large $p_{T}$ limit for the $J/\psi$ and $\Upsilon$ pair
productions. And for $J/\psi+\Upsilon$ the ratio is $-\frac{11}{6}$ for the CO
channel.

If fix the LO LDMEs and estimated the NLO($v^{2}$) LDMEs through the velocity
scaling rule with adopting $v^{2}=0.23$, the differential cross sections are
suppressed by $23\%$ and $84\%$ to LO cross sections for CS channel and CO
channel at large $p_{T}$ region, respectively. The relativistic corrections
significantly dilute the difference between the shape of color-singlet
differential cross sections and the color-octet’s at leading order. So, it is
difficult to study the color-octet mechanism in quarkonium production by
distinguishing the color-siglet contribution and color-octet contribution at
different $p_{T}$ regime, if the relativistic corrections are considered in
pair heavy quarkonia productions at hadron collider.

We also compute the differential cross sections as a function of the double
$J/\psi$ invariant mass and compare it with the recent LHCb experiments. Our
results show that the single parton scattering contributions are large enough
to fit experimental data and it would leave limit room to double parton
scattering mechanism.

Note: when this paper is preparing, another two thesises appeared which
discussed the relativistic corrections of the pair production of $J/\psi$ and
$\eta_{c}$ at the hadron colliderMartynenko:2012tf ; Martynenko:2013daa  but
they were computed within the relativistic quark model. Comparing with their
results, the corrections under the NRQCD frame are much smaller than that
given by relativistic quark model.

###### Acknowledgements.

The authors would like to thank Professor K.T. Chao for useful discussion.
This work was supported by the National Natural Science Foundation of China
(Grants No.10805002, No.10875055, and No.11075011), the Foundation for the
Author of National Excellent Doctoral Dissertation of China (Grants No.
2007B18 and No. 201020), the Project of Knowledge Innovation Program (PKIP) of
Chinese Academy of Sciences, Grant No. KJCX2.YW.W10, and the Education
Ministry of Liaoning Province.

##  5 Appendix

###  5.1  short distance coefficients of ${g+g}\rightarrow{J/\psi+J/\psi}$ and
${g+g}\rightarrow{\Upsilon+\Upsilon}$

The LO short distance coefficients and their relativistic corrections ones for
parton-level process
${g(k_{1})+g(k_{2})}\rightarrow{J/\psi(p_{1})+J/\psi(p_{2})}$ and
${g(k_{1})+g(k_{2})}\rightarrow{\Upsilon(p_{1})+\Upsilon(p_{2})}$ are given in
this appendix. Here $m=2m_{c}$ or $2m_{b}$ represents the $J/\psi$ mass or
$\Upsilon$ mass. (Here we adopt $t(u)$ to represent $t(0)(u(0))$ directly for
simplification)

The LO short distance coefficient for CS channel is:

$\displaystyle{F({}^{3}S_{1}^{[1]}{}^{3}S_{1}^{[1]})}=\frac{64\alpha_{s}^{4}\pi^{3}}{6561m^{2}(m^{2}-t)^{4}(m^{2}-u)^{4}(-2m^{2}+t+u)^{8}}\Big{[}2680m^{24}-14984m^{22}(t+u)$
$\displaystyle+m^{20}(31406t^{2}+89948tu+31406u^{2})-16m^{18}(1989t^{3}+12661t^{2}u$
$\displaystyle+12661tu^{2}+1989u^{3})+4m^{16}(4417t^{4}+57140t^{3}u+117714t^{2}u^{2}+57140tu^{3}+4417u^{4})$
$\displaystyle-4m^{14}(1793t^{5}+38340t^{4}u+134119t^{3}u^{2}+134119t^{2}u^{3}+38340tu^{4}+1793u^{5})+m^{12}(2956t^{6}$
$\displaystyle+76406t^{5}u+361624t^{4}u^{2}+571900t^{3}u^{3}+361624t^{2}u^{4}+76406tu^{5}+2956u^{6})-2m^{10}(397t^{7}$
$\displaystyle+15391t^{6}u+91227t^{5}u^{2}+167593t^{4}u^{3}+167593t^{3}u^{4}+91227t^{2}u^{5}+15391tu^{6}+397u^{7})$
$\displaystyle+m^{8}(47t^{8}+7642t^{7}u+73146t^{6}u^{2}+150334t^{5}u^{3}+132502t^{4}u^{4}+150334t^{3}u^{5}+73146t^{2}u^{6}$
$\displaystyle+7642tu^{7}+47u^{8})+2m^{6}(10t^{9}-411t^{8}u-8951t^{7}u^{2}-29063t^{6}u^{3}-17653t^{5}u^{4}-17653t^{4}u^{5}$
$\displaystyle-29063t^{3}u^{6}-8951t^{2}u^{7}-411tu^{8}+10u^{9})+m^{4}(t^{10}-66t^{9}u+2469t^{8}u^{2}+12874t^{7}u^{3}$
$\displaystyle+11928t^{6}u^{4}+1164t^{5}u^{5}+11928t^{4}u^{6}+12874t^{3}u^{7}+2469t^{2}u^{8}-66tu^{9}+u^{10})+4m^{2}t^{2}u^{2}(9t^{7}$
$\displaystyle-586t^{6}u-37t^{5}u^{2}-394t^{4}u^{3}-394t^{3}u^{4}-37t^{2}u^{5}-586tu^{6}+9u^{7})+2t^{4}u^{4}(349t^{4}-908t^{3}u$
$\displaystyle+1374t^{2}u^{2}-908tu^{3}+349u^{4})\Big{]}$ (25)

The LO short distance coefficient for CO channel is:

$\displaystyle{F({}^{3}S_{1}^{[8]}{}^{3}S_{1}^{[8]})}={\frac{\alpha_{s}^{4}\pi^{3}}{972m^{6}\left(m^{2}-t\right)^{4}\left(m^{2}-u\right)^{4}\left(-2m^{2}+t+u\right)^{8}}}$
$\displaystyle\Big{[}-104624m^{28}+522544m^{26}(t+u)-2m^{24}(137741t^{2}+1836350tu+137741u^{2})$
$\displaystyle-16m^{22}(80955t^{3}-345827t^{2}u-345827tu^{2}+80955u^{3})+4m^{20}t(667342t^{4}-701017t^{3}u$
$\displaystyle-1566126t^{2}u^{2}-701017tu^{3}+667342u^{4})-8m^{18}(316747t^{5}+1677t^{4}u+195158t^{3}u^{2}$
$\displaystyle+195158t^{2}u^{3}+1677tu^{4}+316747u^{5})+m^{16}(1455637t^{6}+281021t^{5}u+5910067t^{4}u^{2}$
$\displaystyle+14296294t^{3}u^{3}+5910067t^{2}u^{4}+281021tu^{5}+1455637u^{6})-2m^{14}(270334t^{7}-179453t^{6}u$
$\displaystyle+2260581t^{5}u^{2}+8407066t^{4}u^{3}+8407066t^{3}u^{4}+2260581t^{2}u^{5}-179453tu^{6}+270334u^{7})$
$\displaystyle+2m^{12}(66670t^{8}-239050t^{7}u+1059585t^{6}u^{2}+5792306t^{5}u^{3}+7997354t^{4}u^{4}+5792306t^{3}u^{5}$
$\displaystyle+1059585t^{2}u^{6}-239050tu^{7}+66670u^{8})-2m^{10}(10661t^{9}-110448t^{8}u+408662t^{7}u^{2}$
$\displaystyle+3075158t^{6}u^{3}+5105335t^{5}u^{4}+5105335t^{4}u^{5}+3075158t^{3}u^{6}+408662t^{2}u^{7}-110448tu^{8}$
$\displaystyle+10661u^{9})+2m^{8}(935t^{10}-25638t^{9}u+178326t^{8}u^{2}+1288892t^{7}u^{3}+2655522t^{6}u^{4}$
$\displaystyle+3136494t^{5}u^{5}+2655522t^{4}u^{6}+1288892t^{3}u^{7}+178326t^{2}u^{8}-25638tu^{9}+935u^{10})$
$\displaystyle+2m^{6}tu(2616t^{9}-65577t^{8}u-431588t^{7}u^{2}-1010108t^{6}u^{3}-1490351t^{5}u^{4}-1490351t^{4}u^{5}$
$\displaystyle-1010108t^{3}u^{6}-431588t^{2}u^{7}-65577tu^{8}+2616u^{9})+m^{4}tu(-243t^{10}+33354t^{9}u$
$\displaystyle+192771t^{8}u^{2}+555338t^{7}u^{3}+879860t^{6}u^{4}+1154016t^{5}u^{5}+879860t^{4}u^{6}+555338t^{3}u^{7}$
$\displaystyle+192771t^{2}u^{8}+33354tu^{9}-243u^{10})-162m^{2}t^{2}u^{2}(t+u)^{3}(27t^{6}+88t^{5}u+174t^{4}u^{2}$
$\displaystyle+208t^{3}u^{3}+174t^{2}u^{4}+88tu^{5}+27u^{6})+243t^{2}u^{2}(t+u)^{4}(t^{2}+tu+u^{2})^{3}\Big{]}$
(26)

The short distance coefficient to relativistic correction for CS channel is:

$\displaystyle{G({}^{3}S_{1}^{[1]}{}^{3}S_{1}^{[1]})}=\frac{256\alpha_{s}^{4}\pi^{3}}{19683m^{4}(m^{2}-t)^{5}(m^{2}-u)^{5}(-2m^{2}+t+u)^{10}(2m^{2}+t+u)}$
$\displaystyle\Big{[}370752m^{34}-2181536(t+u)m^{32}+64\big{[}73726(t^{2}+u^{2})+211649ut\big{]}m^{30}$
$\displaystyle-16\big{[}244805(t^{3}+u^{3})+1965309ut(t+u)\big{]}m^{28}-4\big{[}316463(t^{4}+u^{4})$
$\displaystyle-8160588ut(t^{2}+u^{2})-18782102u^{2}t^{2}\big{]}m^{26}+\big{[}5600590(t^{5}+u^{5})$
$\displaystyle-9749802ut(t^{3}+u^{3})-79767876u^{2}t^{2}(t+u)\big{]}m^{24}-2\big{[}2766605(t^{6}+u^{6})$
$\displaystyle+5407076ut(t^{4}+u^{4})-15844965u^{2}t^{2}(t^{2}+u^{2})-35085112u^{3}t^{3}\big{]}m^{22}$
$\displaystyle+2\big{[}1563020(t^{7}+u^{7})+6369573ut(t^{5}+u^{5})+2696929u^{2}t^{2}(t^{3}+u^{3})$
$\displaystyle-456546u^{3}t^{3}(t+u)\big{]}m^{20}-4\big{[}295064(t^{8}+u^{8})+1552733ut(t^{6}+u^{6})$
$\displaystyle+1230476u^{2}t^{2}(t^{4}+u^{4})+7465803u^{3}t^{3}(t^{2}+u^{2})+17306424u^{4}t^{4}\big{]}m^{18}$
$\displaystyle+2\big{[}157467(t^{9}+u^{9})+858538ut(t^{7}+u^{7})-1932037u^{2}t^{2}(t^{5}+u^{5})$
$\displaystyle+3904293u^{3}t^{3}(t^{3}+u^{3})+38310315u^{4}t^{4}(t+u)\big{]}m^{16}-2\big{[}30250(t^{10}+u^{10})$
$\displaystyle+134737ut(t^{8}+u^{8})-2493940u^{2}t^{2}(t^{6}+u^{6})-5184812u^{3}t^{3}(t^{4}+u^{4})$
$\displaystyle+15398042u^{4}t^{4}(t^{2}+u^{2})+35880758u^{5}t^{5}\big{]}m^{14}+\big{[}8417(t^{11}+u^{11})$
$\displaystyle+34279ut(t^{9}+u^{9})-2497569u^{2}t^{2}(t^{7}+u^{7})-10315939u^{3}t^{3}(t^{5}+u^{5})$
$\displaystyle+1057162u^{4}t^{4}(t^{3}+u^{3})+31198738u^{5}t^{5}(t+u)\big{]}m^{12}+\big{[}-817(t^{12}+u^{12})$
$\displaystyle-13616ut(t^{10}+u^{10})+672840u^{2}t^{2}(t^{8}+u^{8})+4910632u^{3}t^{3}(t^{6}+u^{6})$
$\displaystyle+5008157u^{4}t^{4}(t^{4}+u^{4})-7403400u^{5}t^{5}(t^{2}+u^{2})-13064072u^{6}t^{6}\big{]}m^{10}$
$\displaystyle-\big{[}21(t^{13}+u^{13})-4794ut(t^{11}+u^{11})+84325u^{2}t^{2}(t^{9}+u^{9})$
$\displaystyle+1346073u^{3}t^{3}(t^{7}+u^{7})+3171602u^{4}t^{4}(t^{5}+u^{5})-481617u^{5}t^{5}(t^{3}+u^{3})$
$\displaystyle-3414586u^{6}t^{6}(t+u)\big{]}m^{8}+\big{[}7(t^{14}+u^{14})-248ut(t^{12}+u^{12})$
$\displaystyle-2777u^{2}t^{2}(t^{10}+u^{10})+197002u^{3}t^{3}(t^{8}+u^{8})+973759u^{4}t^{4}(t^{6}+u^{6})$
$\displaystyle+633832u^{5}t^{5}(t^{4}+u^{4})-734333u^{6}t^{6}(t^{2}+u^{2})-747796u^{7}t^{7}\big{]}m^{6}$
$\displaystyle+tu\big{[}(t^{13}+u^{13})+137ut(t^{11}+u^{11})-1320u^{2}t^{2}(t^{9}+u^{9})$
$\displaystyle-171964u^{3}t^{3}(t^{7}+u^{7})-249123u^{4}t^{4}(t^{5}+u^{5})+36693u^{5}t^{5}(t^{3}+u^{3})$
$\displaystyle+125000u^{6}t^{6}(t+u)\big{]}m^{4}+2t^{3}u^{3}\big{[}2(t^{10}+u^{10})+761ut(t^{8}+u^{8})$
$\displaystyle+36064u^{2}t^{2}(t^{6}+u^{6})-14592u^{3}t^{3}(t^{4}+u^{4})+9630u^{4}t^{4}(t^{2}+u^{2})$
$\displaystyle-5362u^{5}t^{5}\big{]}m^{2}-2t^{5}u^{5}\big{[}227(t^{7}+u^{7})+6021ut(t^{5}+u^{5})$
$\displaystyle-9217u^{2}t^{2}(t^{3}+u^{3})+6041u^{3}t^{3}(t+u)\big{]}\Big{]}$
(27)

The short distance coefficient to relativistic correction for CO channel is:

$\displaystyle{G({}^{3}S_{1}^{[8]}{}^{3}S_{1}^{[8]})}=\frac{\alpha_{s}^{4}\pi^{3}}{729m^{8}s^{2}(m^{2}-t)^{5}(m^{2}-u)^{5}(-2m^{2}+t+u)^{8}(2m^{2}+t+u)}$
$\displaystyle\Big{[}-7840128m^{38}+58018240(t+u)m^{36}-16\big{[}4351865(t^{2}+u^{2})+36821254ut\big{]}m^{34}$
$\displaystyle-8\big{[}22135741(t^{3}+u^{3})-183792021ut(t+u)\big{]}m^{32}+4\big{[}173739553(t^{4}+u^{4})$
$\displaystyle-435034596ut(t^{2}+u^{2})-916724794u^{2}t^{2}\big{]}m^{30}-2\big{[}576450337(t^{5}+u^{5})$
$\displaystyle-575408499ut(t^{3}+u^{3})-2090528238u^{2}t^{2}(t+u)\big{]}m^{28}+2\big{[}612519223(t^{6}+u^{6})$
$\displaystyle-216110060ut(t^{4}+u^{4})-1597249863u^{2}t^{2}(t^{2}+u^{2})-1708173704u^{3}t^{3}\big{]}m^{26}$
$\displaystyle+2\big{[}-462600109(t^{7}+u^{7})+86546766ut(t^{5}+u^{5})+1187377312u^{2}t^{2}(t^{3}+u^{3})$
$\displaystyle+911415039u^{3}t^{3}(t+u)\big{]}m^{24}+\big{[}519274693(t^{8}+u^{8})-285877238ut(t^{6}+u^{6})$
$\displaystyle-2090011040u^{2}t^{2}(t^{4}+u^{4})-1710759498u^{3}t^{3}(t^{2}+u^{2})-1356254154u^{4}t^{4}\big{]}m^{22}$
$\displaystyle+\big{[}-221652624(t^{9}+u^{9})+375759263ut(t^{7}+u^{7})+1820975245u^{2}t^{2}(t^{5}+u^{5})$
$\displaystyle+1804776207u^{3}t^{3}(t^{3}+u^{3})+2287765845u^{4}t^{4}(t+u)\big{]}m^{20}+2\big{[}36158642(t^{10}+u^{10})$
$\displaystyle-146055034ut(t^{8}+u^{8})-613461077u^{2}t^{2}(t^{6}+u^{6})-576102238u^{3}t^{3}(t^{4}+u^{4})$
$\displaystyle-723890837u^{4}t^{4}(t^{2}+u^{2})-1150404320u^{5}t^{5}\big{]}m^{18}+\big{[}-17656681(t^{11}+u^{11})$
$\displaystyle+148210450ut(t^{9}+u^{9})+547766115u^{2}t^{2}(t^{7}+u^{7})+168319214u^{3}t^{3}(t^{5}+u^{5})$
$\displaystyle-601363505u^{4}t^{4}(t^{3}+u^{3})-515358665u^{5}t^{5}(t+u)\big{]}m^{16}+2\big{[}1515472(t^{12}+u^{12})$
$\displaystyle-25773835ut(t^{10}+u^{10})-63896082u^{2}t^{2}(t^{8}+u^{8})+200994011u^{3}t^{3}(t^{6}+u^{6})$
$\displaystyle+902427166u^{4}t^{4}(t^{4}+u^{4})+1542963240u^{5}t^{5}(t^{2}+u^{2})+1808885720u^{6}t^{6}\big{]}m^{14}$
$\displaystyle-\big{[}320244(t^{13}+u^{13})-12358347ut(t^{11}+u^{11})+10043059u^{2}t^{2}(t^{9}+u^{9})$
$\displaystyle+412260621u^{3}t^{3}(t^{7}+u^{7})+1621412381u^{4}t^{4}(t^{5}+u^{5})+3279793470u^{5}t^{5}(t^{3}+u^{3})$
$\displaystyle+4552155452u^{6}t^{6}(t+u)\big{]}m^{12}+\big{[}14662(t^{14}+u^{14})-1965047ut(t^{12}+u^{12})$
$\displaystyle+19157992u^{2}t^{2}(t^{10}+u^{10})+213718612u^{3}t^{3}(t^{8}+u^{8})+861766390u^{4}t^{4}(t^{6}+u^{6})$
$\displaystyle+1999127983u^{5}t^{5}(t^{4}+u^{4})+3170423452u^{6}t^{6}(t^{2}+u^{2})+3695142248u^{7}t^{7}\big{]}m^{10}$
$\displaystyle+tu\big{[}190300(t^{13}+u^{13})-7050997ut(t^{11}+u^{11})-69336231u^{2}t^{2}(t^{9}+u^{9})$
$\displaystyle-300672607u^{3}t^{3}(t^{7}+u^{7})-788642811u^{4}t^{4}(t^{5}+u^{5})-1422124344u^{5}t^{5}(t^{3}+u^{3})$
$\displaystyle-1887768238u^{6}t^{6}(t+u)\big{]}m^{8}+tu\big{[}-8019(t^{14}+u^{14})+1359759ut(t^{12}+u^{12})$
$\displaystyle+14495197u^{2}t^{2}(t^{10}+u^{10})+69109969u^{3}t^{3}(t^{8}+u^{8})+205737389u^{4}t^{4}(t^{6}+u^{6})$
$\displaystyle+417998865u^{5}t^{5}(t^{4}+u^{4})+630439593u^{6}t^{6}(t^{2}+u^{2})+718419166u^{7}t^{7}\big{]}m^{6}$
$\displaystyle-t^{2}u^{2}\big{[}143856(t^{13}+u^{13})+1824408ut(t^{11}+u^{11})+10108398u^{2}t^{2}(t^{9}+u^{9})$
$\displaystyle+33845791u^{3}t^{3}(t^{7}+u^{7})+78611535u^{4}t^{4}(t^{5}+u^{5})+132750949u^{5}t^{5}(t^{3}+u^{3})$
$\displaystyle+172324471u^{6}t^{6}(t+u)\big{]}m^{4}+81t^{2}u^{2}(t+u)^{6}\big{[}81(t^{8}+u^{8})+999ut(t^{6}+u^{6})$
$\displaystyle+2910u^{2}t^{2}(t^{4}+u^{4})+5459u^{3}t^{3}(t^{2}+u^{2})+6452u^{4}t^{4}\big{]}m^{2}$
$\displaystyle-2673t^{3}u^{3}(t+u)^{7}(t^{2}+ut+u^{2})^{3}\Big{]}$ (28)

###  5.2  short distance coefficients of ${g+g}\rightarrow{J/\psi+\Upsilon}$

The LO short distance coefficient for
${g+g}\rightarrow{c\overline{c}(^{3}S_{1}^{[8]})+b\overline{b}(^{3}S_{1}^{[1]})}$
channel is:

$\displaystyle{F({c\overline{c}}(^{3}S_{1}^{[8]}){b\overline{b}}(^{3}S_{1}^{[1]}))}=\frac{10\alpha_{s}^{4}\pi^{3}}{243m^{3}{M}({M}^{2}-t)^{2}({M}^{2}-u)^{2}(-2m^{2}+t+u)^{2}s^{2}}$
$\displaystyle\Big{[}2m^{8}{M}^{2}+2m^{6}(7{M}^{4}-5{M}^{2}(t+u)+tu)+2m^{4}(7{M}^{6}-13{M}^{4}(t+u)$
$\displaystyle+{M}^{2}(7t^{2}+10tu+7u^{2})-tu(t+u))+m^{2}({M}^{8}-8{M}^{6}(t+u)$
$\displaystyle+{M}^{4}(15t^{2}+16tu+15u^{2})-4{M}^{2}(2t^{3}+3t^{2}u+3tu^{2}+2u^{3})+t^{2}u^{2})$
$\displaystyle+2{M}^{2}({M}^{4}(t^{2}+tu+u^{2})-{M}^{2}(2t^{3}+3t^{2}u+3tu^{2}+2u^{3})$
$\displaystyle+(t^{2}+tu+u^{2})^{2})\Big{]},$ (29)

where $m=m_{J/\psi}$, and $M=m_{\Upsilon}$.

The relativistic of short distance coefficient for
${g+g}\rightarrow{c\overline{c}(^{3}S_{1}^{[8]})+b\overline{b}(^{3}S_{1}^{[1]})}$
channel is:

$\displaystyle{G({c\overline{c}}(^{3}S_{1}^{[8]}){b\overline{b}}(^{3}S_{1}^{[1]}))}=-20\alpha_{s}^{4}\pi^{3}/(729m^{5}Ms^{2}({M}^{2}-t)^{3}({M}^{2}-u)^{3}(4m^{2}{M}^{2}-(t+u)^{2}))$
$\displaystyle/(-2m^{2}+t+u)^{3}\times\bigg{[}96{M}^{4}(2{M}^{2}-t-u)m^{14}+16{M}^{2}\big{[}61{M}^{6}-82(t+u){M}^{4}+(30(t^{2}+u^{2})$
$\displaystyle+37ut){M}^{2}+3tu(t+u)\big{]}m^{12}+4\big{[}52{M}^{10}-374(t+u){M}^{8}+(549(t^{2}+u^{2})+290ut){M}^{6}$
$\displaystyle+(-227(t^{3}+u^{3})+61ut(u+t)){M}^{4}-tu(85(t^{2}+u^{2})+214ut){M}^{2}+12t^{2}u^{2}(t+u)\big{]}m^{10}$
$\displaystyle-2\big{[}424{M}^{12}-916(t+u){M}^{10}+2(77(t^{2}+u^{2})+934ut){M}^{8}+(867(t^{3}+u^{3})$
$\displaystyle-1439ut(t+u)){M}^{6}+(-505(t^{4}+u^{4})+728ut(t^{2}+u^{2})+1778u^{2}t^{2}){M}^{4}-tu(241(t^{3}+u^{3})$
$\displaystyle+735ut(u+t)){M}^{2}+26t^{2}u^{2}(t+u)^{2}\big{]}m^{8}-2\big{[}44{M}^{14}-512(t+u){M}^{12}+2(545(t^{2}+u^{2})$
$\displaystyle+972ut){M}^{10}-(451(t^{3}+u^{3})+2589ut(t+u)){M}^{8}+(-572(t^{4}+u^{4})+1158ut(t^{2}+u^{2})$
$\displaystyle+2264u^{2}t^{2}){M}^{6}+(401(t^{5}+u^{5})-155ut(t^{3}+u^{3})-1018u^{2}t^{2}(t+u)){M}^{4}+2tu(41(t^{4}+u^{4})$
$\displaystyle+128ut(t^{2}+u^{2})+196u^{2}t^{2}){M}^{2}-3t^{2}u^{2}(t+u)^{3}\big{]}m^{6}+2\big{[}10(t+u){M}^{14}-(239(t^{2}+u^{2})$
$\displaystyle+390ut){M}^{12}+2(252(t^{3}+u^{3})+611ut(t+u)){M}^{10}-(117(t^{4}+u^{4})+1676u^{2}t^{2}$
$\displaystyle+1027ut(u^{2}+t^{2})){M}^{8}-(383(t^{5}+u^{5})+244ut(t^{3}+u^{3})-157u^{2}t^{2}(t+u)){M}^{6}$
$\displaystyle+(225(t^{6}+u^{6})+485ut(t^{4}+u^{4})+670u^{2}t^{2}(t^{2}+u^{2})+684u^{3}t^{3}){M}^{4}-tu(56(t^{5}+u^{5})$
$\displaystyle+205ut(t^{3}+u^{3})+325u^{2}t^{2}(t+u)){M}^{2}-t^{2}u^{2}(t+u)^{2}(t^{2}-9ut+u^{2})\big{]}m^{4}+\big{[}(71(t^{3}+u^{3})$
$\displaystyle+173ut(t+u)){M}^{12}-(t+u)^{2}(79(t^{2}+u^{2})+90ut){M}^{10}-(211(t^{5}+u^{5})+590ut(t^{3}+u^{3})$
$\displaystyle+859u^{2}t^{2}(t+u)){M}^{8}+(t+u)^{2}(375(t^{4}+u^{4})+634ut(t^{2}+u^{2})+1006u^{2}t^{2}){M}^{6}$
$\displaystyle-(156(t^{7}+u^{7})+827ut(t^{5}+u^{5})+2082u^{2}t^{2}(t^{3}+u^{3})+3115u^{3}t^{3}(t+u)){M}^{4}$
$\displaystyle+tu(t+u)^{2}(108(t^{4}+u^{4})+241ut(t^{2}+u^{2})+334u^{2}t^{2}){M}^{2}-5t^{3}u^{3}(t+u)^{3}\big{]}m^{2}$
$\displaystyle-22{M}^{2}({M}^{2}-t)({M}^{2}-u)(t+u)^{3}\big{[}(t^{2}+ut+u^{2}){M}^{4}-(2(t^{3}+u^{3})+3ut(t+u)){M}^{2}$
$\displaystyle+(t^{2}+ut+u^{2})^{2}\big{]}\bigg{]}$ (30)

The short distance coefficient of
${gg}\rightarrow{c\overline{c}(^{3}S_{1}^{[1]})+b\overline{b}(^{3}S_{1}^{[8]})}$
can be get through Exchanging $M$ and $m$. The relativistic correction of this
cross section is :

$\displaystyle{G({c\overline{c}}(^{3}S_{1}^{[1]}){b\overline{b}}(^{3}S_{1}^{[8]}))}=\frac{20\alpha_{s}^{4}\pi^{3}/(729m^{3}{M}^{3}s^{2})}{(m^{2}-t)^{3}(m^{2}-u)^{3}(-2{M}^{2}+t+u)^{4}(4m^{2}{M}^{2}-(t+u)^{2})}$
$\displaystyle\Big{[}(16{M}^{4}\big{[}14{M}^{2}-9(t+u)\big{]}m^{16}+4{M}^{2}\big{[}300{M}^{6}-452(t+u){M}^{4}+(161(t^{2}+u^{2})+274ut){M}^{2}$
$\displaystyle-9(t^{3}+u^{3})+13tu(t+u)\big{]}m^{14}-\big{[}288{M}^{10}+1648(t+u){M}^{8}-4(575(t^{2}+u^{2})+894ut){M}^{6}$
$\displaystyle+16(32(t^{3}+u^{3})+3ut(t+u)){M}^{4}+(t+u)^{2}(3(t^{2}+u^{2})+574ut){M}^{2}-4(t+u)^{3}(3(t^{2}+u^{2})$
$\displaystyle+2ut)\big{]}m^{12}-\big{[}800{M}^{12}-2080(t+u){M}^{10}+16(39(t^{2}+u^{2})+73ut){M}^{8}+4(121(t^{3}+u^{3})$
$\displaystyle-65ut(t+u)){M}^{6}+(34(t^{4}+u^{4})+2916ut(t^{2}+u^{2})+5444u^{2}t^{2}){M}^{4}+(17(t^{5}+u^{5})$
$\displaystyle-807ut(t^{3}+u^{3})-2346u^{2}t^{2}(t+u)){M}^{2}+2(t+u)^{4}(5(t^{2}+u^{2})-ut)\big{]}m^{10}$
$\displaystyle+\big{[}-768(t+u){M}^{12}+800{M}^{14}-48(t^{2}+32ut+u^{2}){M}^{10}-48(19(t^{3}+u^{3})-31tu(t+u)){M}^{8}$
$\displaystyle+2(787(t^{4}+u^{4})-122ut(t^{2}+u^{2})-1546u^{2}t^{2}){M}^{6}+(-932(t^{5}+u^{5})+688ut(t^{3}+u^{3})$
$\displaystyle+6004u^{2}t^{2}(t+u)){M}^{4}+(t+u)^{2}(347(t^{4}+u^{4})-133ut(t^{2}+u^{2})-592u^{2}t^{2}){M}^{2}$
$\displaystyle-2(t+u)^{3}(15(t^{4}+u^{4})+20ut(t^{2}+u^{2})+26u^{2}t^{2})\big{]}m^{8}+\big{[}256{M}^{16}-992(t+u){M}^{14}$
$\displaystyle+8(209(t^{2}+u^{2})+34ut){M}^{12}-48(55(t^{3}+u^{3})-29ut(t+u)){M}^{10}+2(1841(t^{4}+u^{4})$
$\displaystyle+704ut(t^{2}+u^{2})-1722u^{2}t^{2}){M}^{8}-8(395(t^{5}+u^{5})+478ut(t^{3}+u^{3})-183u^{2}t^{2}(t+u)){M}^{6}$
$\displaystyle+2(t+u)^{2}(782(t^{4}+u^{4})-808u^{2}t^{2}-181ut(t^{2}+u^{2})){M}^{4}-(435(t^{7}+u^{7})+1327ut(t^{5}+u^{5})$
$\displaystyle+1381u^{2}t^{2}(t^{3}+u^{3})+617u^{3}t^{3}(t+u)+28u^{2}t^{2}){M}^{2}+2(t+u)^{4}(21(t^{4}+u^{4})$
$\displaystyle+6ut(t^{2}+u^{2}))\big{]}m^{6}+\big{[}32(t^{2}-13ut+u^{2}){M}^{14}-8(41(t^{3}+u^{3})-229ut(t+u)){M}^{12}$
$\displaystyle+8(129(t^{4}+u^{4})-249ut(t^{2}+u^{2})-644u^{2}t^{2}){M}^{10}-2(765(t^{5}+u^{5})+85ut(t^{3}+u^{3})$
$\displaystyle-2282u^{2}t^{2}(t+u)){M}^{8}+2(t+u)^{2}(617(t^{4}+u^{4})-644u^{2}t^{2}-484ut(t^{2}+u^{2})){M}^{6}$
$\displaystyle-2(285(t^{7}+u^{7})+559ut(t^{5}+u^{5})-315u^{2}t^{2}(t^{3}+u^{3})-1753u^{3}t^{3}(t+u)){M}^{4}$
$\displaystyle+(t+u)^{2}(144(t^{6}+u^{6})+159ut(t^{4}+u^{4})-129u^{2}t^{2}(t^{2}+u^{2})-280u^{3}t^{3}){M}^{2}$
$\displaystyle-2(t+u)^{3}(7(t^{6}+u^{6})+6ut(t^{4}+u^{4})-u^{2}t^{2}(t^{2}+u^{2})-8u^{3}t^{3})\big{]}m^{4}$
$\displaystyle-
tu\big{[}-8(9t^{2}-10ut+9u^{2}){M}^{12}+8(31(t^{3}+u^{3})+25ut(t+u)){M}^{10}-2(141(t^{4}+u^{4})$
$\displaystyle+328ut(t^{2}+u^{2})+334u^{2}t^{2}){M}^{8}+4(25(t^{5}+u^{5})+84ut(t^{3}+u^{3})+107u^{2}t^{2}(t+u)){M}^{6}$
$\displaystyle+4(t+u)^{2}(3(t^{4}+u^{4})+20ut(t^{2}+u^{2})+35u^{2}t^{2}){M}^{4}-(t+u)^{3}(8(t^{4}+u^{4})$
$\displaystyle+27ut(t^{2}+u^{2})+34u^{2}t^{2}){M}^{2}+2(t+u)^{4}(t^{2}+ut+u^{2})^{2}\big{]}m^{2}$
$\displaystyle+{M}^{2}t^{2}u^{2}(t+u)^{2}(8{M}^{8}-16(t+u){M}^{6}+2(t^{2}+12ut+u^{2}){M}^{4}$
$\displaystyle+(6(t^{3}+u^{3})-2ut(t+u)){M}^{2}+tu(t+u)^{2}))\Big{]}$

For the process
${g+g}\rightarrow{c\overline{c}(^{3}S_{1}^{[8]})+b\overline{b}(^{3}S_{1}^{[8]})}$,
the express of relativistic correction term is too long to give here. We only
give the LO cross section.

$\displaystyle{F({c\overline{c}}(^{3}S_{1}^{[8]}){b\overline{b}}(^{3}S_{1}^{[8]}))}={\alpha_{s}^{4}\pi^{3}}/(108m^{3}{M}^{3}s^{2})/\Big{[}{(m^{2}-t)}(m^{2}-u)({M}^{2}-t)({M}^{2}-u)\Big{]}^{2}\times$
$\displaystyle
1/\Big{[}(t+u-2m^{2})(t+u-2{M}^{2})\Big{]}^{2}\times\Big{[}92{M}^{4}(t+u-2{M}^{2})^{2}m^{16}+4{M}^{2}(767{M}^{8}-1450(t+u){M}^{6}$
$\displaystyle+18(49(t^{2}+u^{2})+99ut){M}^{4}-2(86(t^{3}+u^{3})+279ut(t+u)){M}^{2}+tu(33ut+7(t^{2}+u^{2})))m^{14}$
$\displaystyle+2(2896{M}^{12}-9812(t+u){M}^{10}+4(5417ut+2899(t^{2}+u^{2}){M}^{8}-2(2847(t^{3}+u^{3})$
$\displaystyle+7877ut(t+u)){M}^{6}+(1003(t^{4}+u^{4})+6030u^{2}t^{2}+3710ut(u^{2}+t^{2})){M}^{4}-2tu(19(t^{3}+u^{3})$
$\displaystyle+211ut(t+u)){M}^{2}+76t^{2}u^{2}(t^{2}+ut+u^{2}))m^{12}+(3068{M}^{14}-19624(t+u){M}^{12}$
$\displaystyle+(41811(u^{2}+t^{2})+75706ut){M}^{10}-9(4421(t^{3}+u^{3})+11211ut(t+u)){M}^{8}+2(8729(t^{4}+u^{4})$
$\displaystyle+27487ut(t^{2}+u^{2})+40854u^{2}t^{2}){M}^{6}-(2924(t^{5}+u^{5})+10065ut(t^{3}+u^{3})+23123u^{2}t^{2}(t+u)){M}^{4}$
$\displaystyle+tu(2303ut(t^{2}+u^{2})-188(t^{4}+u^{4})+3378u^{2}t^{2}){M}^{2}-8t^{2}u^{2}(84(t^{3}+u^{3})+149ut(u+t)))m^{10}$
$\displaystyle+(368{M}^{16}-5800(t+u){M}^{14}+8(2899(t^{2}+u^{2})+5417ut){M}^{12}-9(4421(t^{3}+u^{3})$
$\displaystyle+11211ut(t+u)){M}^{10}+(33691(t^{4}+u^{4})+100610ut(t^{2}+u^{2})+147662u^{2}t^{2}){M}^{8}-2(6949(t^{5}+u^{5})$
$\displaystyle+21495ut(t^{3}+u^{3})+43712u^{2}t^{2}(t+u)){M}^{6}+(2274(t^{6}+u^{6})+4949ut(t^{4}+u^{4})+29902u^{3}t^{3}$
$\displaystyle+20878u^{2}t^{2}(t^{2}+u^{2})){M}^{4}+tu(718(t^{5}+u^{5})-4033ut(t^{3}+u^{3})-8973u^{2}t^{2}(t+u)){M}^{2}$
$\displaystyle+2t^{2}u^{2}(679(t^{4}+u^{4})+1669ut(t^{2}+u^{2})+2164u^{2}t^{2}))m^{8}-(368(t+u){M}^{16}-72(49(t^{2}+u^{2})$
$\displaystyle+99ut){M}^{14}+4(2847(t^{3}+u^{3})+7877ut(t+u)){M}^{12}-2(8729(t^{4}+u^{4})+27487ut(t^{2}+u^{2})$
$\displaystyle+40854u^{2}t^{2}){M}^{10}+2(6949(t^{5}+u^{5})+21495ut(t^{3}+u^{3})+43712u^{2}t^{2}(t+u)){M}^{8}-(5580(t^{6}+u^{6})$
$\displaystyle+12481ut(t^{4}+u^{4})+42144u^{2}t^{2}(t^{2}+u^{2})+58166u^{3}t^{3}){M}^{6}+(912(t^{7}+u^{7})-893ut(t^{5}+u^{5})$
$\displaystyle+11979u^{2}t^{2}(t^{3}+u^{3})+26506u^{3}t^{3}(t+u)){M}^{4}+tu(762(t^{6}+u^{6})-5413ut(t^{4}+u^{4})-21038u^{3}t^{3}$
$\displaystyle-15468u^{2}t^{2}(t^{2}+u^{2})){M}^{2}+2t^{2}u^{2}(789(t^{5}+u^{5})+2413ut(t^{3}+u^{3})+3910u^{2}t^{2}(t+u)))m^{6}$
$\displaystyle+(92(t+u)^{2}{M}^{16}-8(86(t^{3}+u^{3})+279ut(t+u)3){M}^{14}+2(1003(t^{4}+u^{4})+6030u^{2}t^{2}$
$\displaystyle+3710ut(t^{2}+u^{2})){M}^{12}-(2924(t^{5}+u^{5})+10065ut(t^{3}+u^{3})+23123u^{2}t^{2}(t+u)){M}^{10}$
$\displaystyle+(2274(t^{6}+u^{6})+4949ut(t^{4}+u^{4})+20878u^{2}t^{2}(t^{2}+u^{2})+29902u^{3}t^{3}){M}^{8}-(912(t^{7}+u^{7})$
$\displaystyle-893ut(t^{5}+u^{5})+11979u^{2}t^{2}(t^{3}+u^{3})+26506u^{3}t^{2}(t+u)){M}^{6}+(152(t^{8}+u^{8})$
$\displaystyle-1483ut(t^{6}+u^{6})+8243u^{2}t^{2}(t^{4}+u^{4})+24411u^{3}t^{3}(t^{2}+u^{2})+34058u^{4}t^{4}){M}^{4}$
$\displaystyle+tu(334(t^{7}+u^{7})-4775ut(t^{5}+u^{5})-15750u^{2}t^{2}(t^{3}+u^{3})-26417u^{3}t^{3}(t+u)){M}^{2}$
$\displaystyle+2t^{2}u^{2}(532(t^{6}+u^{6})+1939ut(t^{4}+u^{4})+3770u^{2}t^{2}(t^{2}+u^{2})+4618u^{3}t^{3}))m^{4}-tu(-4(7t^{2}$
$\displaystyle+33ut+7u^{2}){M}^{14}+(76(t^{3}+u^{3})+844ut(t+u)){M}^{12}+(188(t^{4}+u^{4})-2303ut(t^{2}+u^{2})$
$\displaystyle-3378u^{2}t^{2}){M}^{10}+(-718(t^{5}+u^{5})+4033ut(t^{3}+u^{3})+8973u^{2}t^{2}(t+u)){M}^{8}+(762(t^{6}+u^{6})$
$\displaystyle-5413ut(t^{4}+u^{4})-15468u^{2}t^{2}(t^{2}+u^{2})-21038u^{3}t^{3}){M}^{6}+(-334(t^{7}+u^{7})+4775ut(t^{5}+u^{5})$
$\displaystyle+15750u^{2}t^{2}(t^{3}+u^{3})+26417u^{3}t^{3}(t+u)){M}^{4}+(54(t^{8}+u^{8})-2182ut(t^{6}+u^{6})-8161u^{2}t^{2}(t^{4}+u^{4})$
$\displaystyle-16101u^{3}t^{3}(t^{2}+u^{2})-19812u^{4}t^{4}){M}^{2}+54tu(7(t^{7}+u^{7})+30ut(t^{5}+u^{5})+68u^{2}t^{2}(t^{3}+u^{3})$
$\displaystyle+99u^{3}t^{3}(t+u)))m^{2}+2t^{2}u^{2}(-2{M}^{2}+t+u)^{2}(19(t^{2}+ut+u^{2}){M}^{8}-(65(t^{3}+u^{3})$
$\displaystyle+111ut(t+u)){M}^{6}+(100(t^{4}+u^{4})+227ut(t^{2}+u^{2})+300u^{2}t^{2}){M}^{4}-27(3(t^{5}+u^{5})+8ut(t^{3}+u^{3})$
$\displaystyle+13u^{2}t^{2}(t+u)){M}^{2}+27(t^{2}+ut+u^{2})^{3})\Big{]}$ (32)

## References

  * (1) G. T. Bodwin, E. Braaten, and G. P. Lepage, Rigorous QCD analysis of inclusive annihilation and production of heavy quarkonium, Phys. Rev. D51 (1995) 1125–1171, [hep-ph/9407339].
  * (2) H.-w. Huang and K.-t. Chao, QCD radiative correction to the hadronic annihilation rate of $1^{+-}$ heavy quarkonium, Phys.Rev. D54 (1996) 3065–3072, [hep-ph/9601283].
  * (3) H.-W. Huang and K.-T. Chao, The cancellation of infrared divergences in the hadronic annihilation decay of heavy quarkonium, Phys. Rev. D55 (1997) 244–248, [hep-ph/9605362].
  * (4) H.-W. Huang and K.-T. Chao, QCD predictions for annihilation decays of P wave quarkonia to next-to-leading order in $\alpha_{s}$, Phys. Rev. D54 (1996) 6850–6854, [hep-ph/9606220].
  * (5) G. T. Bodwin, E. Braaten, and G. P. Lepage, Rigorous QCD predictions for decays of P wave quarkonia, Phys.Rev. D46 (1992) 1914–1918, [hep-lat/9205006].
  * (6) N. Brambilla, E. Mereghetti, and A. Vairo, Hadronic quarkonium decays at order $v^{7}$, Phys.Rev. D79 (2009) 074002, [arXiv:0810.2259].
  * (7) N. Brambilla, E. Mereghetti, and A. Vairo, Electromagnetic quarkonium decays at order $v^{7}$, JHEP 0608 (2006) 039, [hep-ph/0604190].
  * (8) N. Brambilla, D. Eiras, A. Pineda, J. Soto, and A. Vairo, Inclusive decays of heavy quarkonium to light particles, Phys.Rev. D67 (2003) 034018, [hep-ph/0208019].
  * (9) N. Brambilla, D. Eiras, A. Pineda, J. Soto, and A. Vairo, New predictions for inclusive heavy quarkonium P wave decays, Phys.Rev.Lett. 88 (2002) 012003, [hep-ph/0109130].
  * (10) Z.-G. He, Y. Fan, and K.-T. Chao, QCD prediction for the non-$D\bar{D}$ annihilation decay of $\psi(3770)$, Phys. Rev. Lett. 101 (2008) 112001, [arXiv:0802.1849].
  * (11) Z.-G. He, Y. Fan, and K.-T. Chao, NRQCD Predictions of D-Wave Quarkonia ${}^{3}D_{J}(J=1,2,3)$ Decay into Light Hadrons at Order $\alpha_{s}^{3}$, Phys.Rev. D81 (2010) 074032, [arXiv:0910.3939].
  * (12) Y. Fan, Z.-G. He, Y.-Q. Ma, and K.-T. Chao, Predictions of Light Hadronic Decays of Heavy Quarkonium ${}^{1}D_{2}$ States in NRQCD, Phys. Rev. D80 (2009) 014001, [arXiv:0903.4572].
  * (13) P. L. Cho and A. K. Leibovich, Color-singlet $\psi_{Q}$ production at $e^{+}e^{-}$ colliders, Phys. Rev. D54 (1996) 6690–6695, [hep-ph/9606229].
  * (14) F. Yuan, C.-F. Qiao, and K.-T. Chao, Prompt $J/\psi$ production at $e^{+}e^{-}$ colliders, Phys. Rev. D56 (1997) 321–328, [hep-ph/9703438].
  * (15) K.-Y. Liu, Z.-G. He, and K.-T. Chao, Inclusive charmonium production via double $c\bar{c}$ in $e^{+}e^{-}$ annihilation, Phys. Rev. D69 (2004) 094027, [hep-ph/0301218].
  * (16) E. Braaten and J. Lee, Exclusive double-charmonium production in $e^{+}e^{-}$ annihilation, Phys. Rev. D67 (2003) 054007, [hep-ph/0211085].
  * (17) K.-Y. Liu, Z.-G. He, and K.-T. Chao, Problems of double charm production in $e^{+}e^{-}$ annihilation at $\sqrt{s}=10.6~{}{GeV}$., Phys. Lett. B557 (2003) 45–54, [hep-ph/0211181].
  * (18) Belle Collaboration, K. Abe et al., Observation of double $c\bar{c}$ production in $e^{+}e^{-}$ annihilation at $\sqrt{s}$ approx. 10.6-GeV, Phys. Rev. Lett. 89 (2002) 142001, [hep-ex/0205104].
  * (19) Belle Collaboration, K. Abe et al., Study of double charmonium production in $e^{+}e^{-}$ annihilation at $\sqrt{s}$ approx. 10.6-GeV, Phys. Rev. D70 (2004) 071102, [hep-ex/0407009].
  * (20) BABAR Collaboration, B. Aubert et al., Measurement of double charmonium production in $e^{+}e^{-}$ annihilations at $\sqrt{s}=10.6$ GeV, Phys. Rev. D72 (2005) 031101, [hep-ex/0506062].
  * (21) Y.-J. Zhang and K.-T. Chao, Double charm production $e^{+}e^{-}\to J/\psi+c\bar{c}$ at B factories with next-to-leading order QCD correction, Phys.Rev.Lett. 98 (2007) 092003, [hep-ph/0611086].
  * (22) K. Wang, Y.-Q. Ma, and K.-T. Chao, QCD corrections to $e^{+}e^{-}toJ/\psi(\psi(2S))+\chi_{cJ}(J=0,1,2)$ at B Factories, Phys.Rev. D84 (2011) 034022, [arXiv:1107.2646].
  * (23) Y.-J. Zhang, Y.-Q. Ma, and K.-T. Chao, Factorization and NLO QCD correction in $e^{+}e^{-}\to J/\psi(\psi(2S))+\chi_{c0}$ at B Factories, Phys. Rev. D78 (2008) 054006, [arXiv:0802.3655].
  * (24) Y.-J. Zhang, Y.-Q. Ma, K. Wang, and K.-T. Chao, QCD radiative correction to color-octet J/psi inclusive production at B Factories, Phys.Rev. D81 (2010) 034015, [arXiv:0911.2166].
  * (25) B. Gong and J.-X. Wang, Next-to-leading-order QCD corrections to $e^{+}e^{-}\to J/\psi c\bar{c}$ at the B factories, Phys.Rev. D80 (2009) 054015, [arXiv:0904.1103].
  * (26) B. Gong and J.-X. Wang, Next-to-Leading-Order QCD Corrections to $e^{+}e^{-}\to J/\psi gg$ at the B Factories, Phys.Rev.Lett. 102 (2009) 162003, [arXiv:0901.0117].
  * (27) Y.-Q. Ma, Y.-J. Zhang, and K.-T. Chao, QCD correction to $e^{+}e^{-}\to J/\psi gg$ at B Factories, Phys.Rev.Lett. 102 (2009) 162002, [arXiv:0812.5106].
  * (28) H.-R. Dong, F. Feng, and Y. Jia, $O(\alpha_{s})$ corrections to $J/\psi+\chi_{cJ}$ production at $B$ factories, JHEP 1110 (2011) 141, [arXiv:1107.4351].
  * (29) G. T. Bodwin, H. S. Chung, and J. Lee, Endpoint Logarithms in $e^{+}e^{-}\to J/\psi+\eta_{c}$, arXiv:1301.3937.
  * (30) Z.-G. He, Y. Fan, and K.-T. Chao, Relativistic corrections to $J/\psi$ exclusive and inclusive double charm production at B factories, Phys.Rev. D75 (2007) 074011, [hep-ph/0702239].
  * (31) Y. Jia, Color-singlet relativistic correction to inclusive J/psi production associated with light hadrons at B factories, Phys.Rev. D82 (2010) 034017, [arXiv:0912.5498].
  * (32) Z.-G. He, Y. Fan, and K.-T. Chao, Relativistic correction to $e^{+}e^{-}\to J/\psi+gg$ at B factories and constraint on color-octet matrix elements, Phys.Rev. D81 (2010) 054036, [arXiv:0910.3636].
  * (33) Y. Fan, J. Lee, and C. Yu, Resummation of relativistic corrections to exclusive productions of charmonia in e+ e- collisions, arXiv:1211.4111.
  * (34) Y. Fan, J. Lee, and C. Yu, Higher-order corrections to exclusive production of charmonia at B factories, arXiv:1209.1875.
  * (35) H.-R. Dong, F. Feng, and Y. Jia, $O(\alpha_{s}v^{2})$ correction to $e^{+}e^{-}\to J/\psi+\eta_{c}$ at $B$ factories, Phys.Rev. D85 (2012) 114018, [arXiv:1204.4128].
  * (36) X.-H. Li and J.-X. Wang, ${\mathcal{O}}(\alpha_{s}v^{2})$ correction to $J/\psi$ plus $\eta_{c}$ production in $e^{+}e^{-}$ annihilation at $\sqrt{s}$=10.6GeV, arXiv:1301.0376.
  * (37) 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].
  * (38) B. A. Kniehl and G. Kramer, Color octet contributions to $J/\psi$ photoproduction via fragmentation at HERA, Phys.Lett. B413 (1997) 416–421, [hep-ph/9703280].
  * (39) J. M. Campbell, F. Maltoni, and F. Tramontano, QCD corrections to $J/\psi$ and $\Upsilon$ production at hadron colliders, Phys.Rev.Lett. 98 (2007) 252002, [hep-ph/0703113].
  * (40) B. Gong and J.-X. Wang, QCD corrections to polarization of $J/\psi$ and $\Upsilon$ at Tevatron and LHC, Phys. Rev. D78 (2008) 074011, [arXiv:0805.2469].
  * (41) B. Gong, J.-X. Wang, and H.-F. Zhang, QCD corrections to $\Upsilon$ production via color-octet states at the Tevatron and LHC, Phys.Rev. D83 (2011) 114021, [arXiv:1009.3839].
  * (42) B. Gong and J.-X. Wang, Next-to-leading-order QCD corrections to $J/\psi$ polarization at Tevatron and Large-Hadron-Collider energies, Phys. Rev. Lett. 100 (2008) 232001, [arXiv:0802.3727].
  * (43) B. Gong, J.-P. Lansberg, C. Lorce, and J. Wang, Next-to-leading-order QCD corrections to the yields and polarisations of $J/\psi$ and Upsilon directly produced in association with a Z boson at the LHC, arXiv:1210.2430.
  * (44) L. Gang, W. ShuangTe, S. Mao, and L. JiPing, Prompt heavy quarkonium production in association with a massive (anti)bottom quark at the LHC, Phys.Rev. D85 (2012) 074026, [arXiv:1203.0799].
  * (45) Y.-Q. Ma, K. Wang, and K.-T. Chao, QCD radiative corrections to $\chi_{cJ}$ production at hadron colliders, Phys.Rev. D83 (2011) 111503, [arXiv:1002.3987].
  * (46) M. Butenschoen and B. A. Kniehl, $J/\psi$ polarization at Tevatron and LHC: Nonrelativistic-QCD factorization at the crossroads, Phys.Rev.Lett. 108 (2012) 172002, [arXiv:1201.1872].
  * (47) B. Gong, L.-P. Wan, J.-X. Wang, and H.-F. Zhang, Polarization for Prompt $J/\psi$, $psi(2s)$ production at the Tevatron and LHC, arXiv:1205.6682.
  * (48) N. Brambilla, S. Eidelman, B. Heltsley, R. Vogt, G. Bodwin, et al., Heavy quarkonium: progress, puzzles, and opportunities, Eur.Phys.J. C71 (2011) 1534, [arXiv:1010.5827].
  * (49) Quarkonium Working Group Collaboration, N. Brambilla et al., Heavy quarkonium physics, hep-ph/0412158.
  * (50) LHCb Collaboration Collaboration, R. Aaij et al., Observation of $J/\psi$ pair production in $p\overline{p}$ collisions at $\sqrt{s}=7TeV$, Phys.Lett. B707 (2012) 52–59, [arXiv:1109.0963].
  * (51) V. D. Barger, S. Fleming, and R. J. N. Phillips, Double gluon fragmentation to $J/\psi$pairs at the Tevatron, Phys. Lett. B371 (1996) 111–116, [hep-ph/9510457].
  * (52) C.-F. Qiao, $J/\psi$ pair production at the Tevatron, Phys. Rev. D66 (2002) 057504, [hep-ph/0206093].
  * (53) C.-F. Qiao, L.-P. Sun, and P. Sun, Testing Charmonium Production Mechanism via Polarized $J/\psi$ Pair Production at the LHC, J. Phys. G37 (2010) 075019, [arXiv:0903.0954].
  * (54) R. Li, Y.-J. Zhang, and K.-T. Chao, Pair Production of Heavy Quarkonium and $B_{c}^{(*)}$ Mesons at Hadron Colliders, Phys.Rev. D80 (2009) 014020, [arXiv:0903.2250].
  * (55) P. Ko, C. Yu, and J. Lee, Inclusive double-quarkonium production at the Large Hadron Collider, JHEP 01 (2011) 070, [arXiv:1007.3095].
  * (56) A. V. Berezhnoy, A. K. Likhoded, A. V. Luchinsky, and A. A. Novoselov, Double $J/\psi$-meson Production at LHC and 4c-tetraquark state, Phys. Rev. D84 (2011) 094023, [arXiv:1101.5881].
  * (57) C. H. Kom, A. Kulesza, and W. J. Stirling, Pair production of $J/\psi$ as a probe of double parton scattering at LHCb, Phys. Rev. Lett. 107 (2011) 082002, [arXiv:1105.4186].
  * (58) S. Baranov, A. Snigirev, and N. Zotov, Double heavy meson production through double parton scattering in hadronic collisions, Phys.Lett. B705 (2011) 116–119, [arXiv:1105.6276].
  * (59) Y. Fan, Y.-Q. Ma, and K.-T. Chao, Relativistic Correction to $J/\psi$ Production at Hadron Colliders, arXiv:0904.4025.
  * (60) H.-K. Guo, Y.-Q. Ma, and K.-T. Chao, $O(\alpha_{s}v^{2})$ Corrections to Hadronic and Electromagnetic Decays of ${}^{1}S_{0}$ Heavy Quarkonium, Phys. Rev. D83 (2011) 114038, [arXiv:1104.3138].
  * (61) G. T. Bodwin, H. S. Chung, D. Kang, J. Lee, and C. Yu, Improved determination of color-singlet nonrelativistic QCD matrix elements for S-wave charmonium, Phys. Rev. D77 (2008) 094017, [arXiv:0710.0994].
  * (62) G. T. Bodwin, D. K. Sinclair, and S. Kim, Quarkonium decay matrix elements from quenched lattice QCD, Phys. Rev. Lett. 77 (1996) 2376–2379, [hep-lat/9605023].
  * (63) G. T. Bodwin, D. K. Sinclair, and S. Kim, Bottomonium decay matrix elements from lattice QCD with two light quarks, Phys. Rev. D65 (2002) 054504, [hep-lat/0107011].
  * (64) G. T. Bodwin, J. Lee, and D. K. Sinclair, Spin correlations and velocity-scaling in color-octet NRQCD matrix elements, Phys. Rev. D72 (2005) 014009, [hep-lat/0503032].
  * (65) E. Braaten and Y.-Q. Chen, Helicity Decomposition for Inclusive $J/\psi$ Production, Phys. Rev. D54 (1996) 3216–3227, [hep-ph/9604237].
  * (66) Y.-Q. Ma, General Form of $s$, $t$, $u$ Symmetric Polynomial and Heavy Quarkonium physics, arXiv:1207.3073.
  * (67) M. Gremm and A. Kapustin, Annihilation of S wave quarkonia and the measurement of alpha-s, Phys.Lett. B407 (1997) 323–330, [hep-ph/9701353].
  * (68) G. T. Bodwin and J. Lee, Relativistic corrections to gluon fragmentation into spin triplet S wave quarkonium, Phys.Rev. D69 (2004) 054003, [hep-ph/0308016].
  * (69) C.-F. Qiao, J. Wang, and Y.-h. Zheng, DJpsiFDC: an event generator for the process $gg\to J/\psi J/\psi$ at LHC, Chin.Phys. C35 (2011) 209–213, [arXiv:1007.0617].
  * (70) P. Ko, C. Yu, and J. Lee, $pp\to J/\psi+\Upsilon+X$ as a clean probe to the quarkonium production mechanism, Phys.Rev. D83 (2011) 054015, [arXiv:1006.3846].
  * (71) J. Kublbeck, M. Bohm, and A. Denner, Feyn arts: Computer algebraic generation of feynman graphs and amplitudes, Comput. Phys. Commun. 60 (1990) 165–180.
  * (72) R. Mertig, M. Bohm, and A. Denner, FEYN CALC: Computer algebraic calculation of Feynman amplitudes, Comput. Phys. Commun. 64 (1991) 345–359.
  * (73) J. Pumplin et al., New generation of parton distributions with uncertainties from global QCD analysis, JHEP 07 (2002) 012, [hep-ph/0201195].
  * (74) E. Braaten, B. A. Kniehl, and J. Lee, Polarization of prompt $J/\psi$ at the Tevatron, Phys. Rev. D62 (2000) 094005, [hep-ph/9911436].
  * (75) D. Kang, T. Kim, J. Lee, and C. Yu, Inclusive Charm Production in $\Upsilon(nS)$ Decay, Phys.Rev. D76 (2007) 114018, [arXiv:0707.4056].
  * (76) M. Kramer, Quarkonium production at high-energy colliders, Prog. Part. Nucl. Phys. 47 (2001) 141–201, [hep-ph/0106120].
  * (77) G.-Z. Xu, Y.-J. Li, K.-Y. Liu, and Y.-J. Zhang, Relativistic Correction to Color Octet $J/\psi$ Production at Hadron Colliders, arXiv:1203.0207.
  * (78) C. W. Bauer, S. Fleming, D. Pirjol, and I. W. Stewart, An Effective field theory for collinear and soft gluons: Heavy to light decays, Phys.Rev. D63 (2001) 114020, [hep-ph/0011336].
  * (79) A. Martynenko and A. Trunin, Relativistic corrections to double charmonium production in high energy proton-proton interaction, Phys.Rev. D86 (2012) 094003, [arXiv:1207.3245].
  * (80) A. Martynenko and A. Trunin, Relativistic corrections to $\eta_{c}$-pair production in high energy proton-proton collisions, arXiv:1302.6726.