#  A study of relativistic corrections to $J/\psi\rightarrow p\bar{p}$ decay

Nikolay Kivel
Physik-Department, Technische Universität München,
James-Franck-Str. 1, 85748 Garching, Germany

###### Abstract

We study relativistic corrections in exclusive $S$-wave charmonium decays into
proton-antiproton final state. We calculate the NRQCD corrections to the
dominant decay amplitude, which depend on the nucleon twist-3 light-cone
distribution amplitudes only. It is shown that in this case the collinear
factorisation is also valid beyond the leading-order approximation. Our
numerical estimates show that relativistic correction of relative order
$v^{2}$ provides large numerical impact.

##  1 Introduction

In recent years, the BESII and BESIII collaborations have obtained many
accurate data on baryon-antibaryon decays of $S$-wave charmonia, see e.g.
Refs.[1, 2, 3, 4, 5, 6, 7, 8]. These data allows one to get information about
various important hadronic parameters and provide an interesting possibility
to study the QCD dynamics.

There are different approaches for a description of the baryon-antybaryon
decays. In Refs.[9, 10] it was proposed to use a phenomenological model based
on an empirical effective Lagrangian with unknown parameters, which are
subject to restrictions from $SU(3)$ flavour symmetry and experimental data.

Another way is to use effective field theory framework, combining
nonrelativistic expansion and collinear factorisation [11, 12, 13, 14].
Appropriate effective Lagrangian is described by NRQCD [15, 16] and collinear
copies of QCD [11, 12, 14].

Such approach allows one to built a systematic expansion in powers of
$1/m_{c}$ describing decay amplitudes as a superposition of perturbative and
universal nonperturbative functions. The perturbative part can be computed in
pQCD at the scale of order $m_{c}$. The nonperturbative matrix elements are
defined at soft hadronic scale $\Lambda\ll m_{c}$ . They can be limited to
other processes or evaluated using non-perturbative methods such as QCD sum
rules or lattice QCD.

Since the real mass of the charm quark is not large enough, a systematic
description requires careful analysis of the various corrections. Existing
calculations of exclusive and inclusive reactions show that in many cases
radiative corrections in $\alpha_{s}(m_{c})$ and relativistic corrections can
have a significant numerical effect, see for instance, reviews [17, 18].

The first estimates for exclusive proton-antiproton decay based on the
leading-order (LO) approximation were obtained already long time ago in
Refs.[11, 14]. Recently, this analysis has been extended in Refs.[19, 20, 21],
where the power-law corrections $\sim\Lambda/m_{c}$ associated with the
collinear expansion were taken into account for different amplitudes. The
collinear power corrections are associated with the higher twist baryon light-
cone distribution amplitudes (LCDAs). It has been established that the effect
of such corrections is not very strong and one can describe various channels
of baryon-antibaryon decay with good accuracy. However, relativistic and QCD
radiative corrections for baryon–antibaryon decays have not yet been studied.

The available data for various baryons indicate the possible presence of
effects from relativistic corrections. For example, the so-called “12% rule”

$Q_{B}=\frac{Br[\psi(2S)\rightarrow B\bar{B}]}{Br[J/\psi\rightarrow
B\bar{B}]}\simeq\frac{Br[\psi(2S)\rightarrow e^{+}e^{-}]}{Br[J/\psi\rightarrow
e^{+}e^{-}]}\simeq 0.13,$ (1)

must hold if the width is dominated by the LO NRQCD approximation. The
available experimental data [22] show the following results:
$Q_{p}=0.1386(3)$, $Q_{n}=0.146(1)$, $Q_{\Lambda}=0.204(1)$,
$Q_{\Sigma^{0}}=0.21(3)$, $Q_{\Sigma^{+}}=0.072(5)$ and
$Q_{\Xi}^{+}=0.276(5)$. Consequently, the data for baryons other than nucleon
indicate a violation of this expectation.

For nucleon-antinucleon decays the rule (1) works quite well, but on the other
hand, the values of the polar angular distribution coefficient $\alpha_{p}$
for ground and excited $S$-wave charmonium ($\psi(2S)\equiv\psi^{\prime}$) are
very different [3, 7]

$\left.\alpha_{p}\right|_{J/\psi}=0.595\pm 0.012,\ \
\left.\alpha_{p}\right|_{\psi^{\prime}}=1.03\pm 0.06.$ (2)

This observation may also indicate a significant contribution of relativistic
corrections.

The purpose of this work is to study the effects of relativistic corrections
in proton-antiproton decays. The technique for calculating relativistic
corrections is well known. Such corrections were already studied in for
different exclusive decays as $J/\psi\rightarrow e^{+}e^{-}$ [23, 24],
$H\rightarrow\gamma J/\psi$ [25, 26], $\chi_{cJ}\rightarrow\gamma\gamma$ [27,
28] and for exclusive production $e^{+}e^{-}\rightarrow J/\psi\eta_{c}$ [29,
30]. However, relativistic corrections to exclusive hadronic decays have not
yet been considered.

The feature of hadronic decays is that the hard kernels are not simple numbers
but depend on light quark momentum fractions. Corresponding decay amplitude is
described by the convolution integral over the momentum fractions of the hard
kernel with LCDAs. In this case it is more convenient to perform the hard
matching at the amplitude level. In some cases the convolution integrals have
infrared divergencies that indicates about the violation of the collinear
factorisation. This is often the case for next-to-leading power contributions
in exclusive processes. However, for the dominant amplitude describing
$n^{3}S_{1}\rightarrow B\bar{B}$ decay, the collinear factorisation is still
applicable and this makes it possible to systematically study the relativistic
effects.

Our paper is organised as follows. In Sec. 2 we give important definitions and
kinematical notation. In Sec. 3 we calculate the hard kernels and perform the
resummation of some class of relativistic effects. After that we provide a
qualitative numerical analysis for different choices of the nucleon LCDAs. In
Sec. 4 we discuss the obtained results and make conclusions. In Appendix we
give important details about the nucleon twist-3 LCDAs and provide analytical
expressions for the hard kernels.

##  2 Kinematics and decay amplitudes

In this work we use notation from Ref.[21]. We define the kinematics of
$J/\psi(P)\rightarrow p(k)\bar{p}(k^{\prime})$ decay in the charmonium rest
frame

$P=M_{\psi}\omega,~{}\omega=(1,\vec{0}).\ $ (3)

The outgoing momenta $k$ and $k^{\prime}$ are directed along the $z$-axis and
read

$k=(M_{\psi}/2,0,0,M_{\psi}\beta/2),~{}k^{\prime}=(M_{\psi}/2,0,0,-M_{\psi}\beta/2),$
(4)

where $m_{N}$ is the nucleon mass and
$~{}\beta=\sqrt{1-4m_{N}^{2}/M_{\psi}^{2}}$. We also define auxiliary light-
cone vectors

$n=(1,0,0,-1),~{}\ \bar{n}=(1,0,0,1),$ (5)

so that any four-vector $V$ can be represented as

$V=V_{+}\frac{\bar{n}}{2}+V_{-}\frac{n}{2}+V_{\bot},~{}$ (6)

where $\ V_{+}\equiv(V\cdot n)=V_{0}+V_{3}$,
$V_{-}\equiv(V\cdot\bar{n})=V_{0}-V_{3}$.

The decay amplitude $J/\psi\rightarrow p\bar{p}$ is defined as

$\left\langle
k,k^{\prime}\right|~{}i\hat{T}~{}\left|P\right\rangle=(2\pi)^{4}\delta(P-k-k^{\prime})~{}iM,$
(7)

with

$M=\bar{N}(k)\left\\{A_{1}\hbox to0.0pt{\hbox
to5.00002pt{\hfil$\epsilon$\hfil}\hss}/_{\psi}+A_{2}\left(\epsilon_{\psi}\right)_{\mu}(k^{\prime}+k)_{\nu}\frac{i\sigma^{\mu\nu}}{2m_{N}}\right\\}V(k^{\prime})~{},$
(8)

where $\hbox to0.0pt{\hbox
to5.03125pt{\hfil/\hfil}\hss}p=p_{\mu}\gamma^{\mu}$. The nucleon $\bar{N}(k)$
and antinucleon $V(k^{\prime})$ spinors have standard normalisation
$\bar{N}N=2m_{N}$ and $\bar{V}V=-2m_{N}$. The charmonium polarisation vector
$\epsilon_{\psi}^{\mu}\equiv\epsilon_{\psi}^{\mu}(P,\lambda)$ satisfies

$\sum_{\lambda}\epsilon_{\psi}^{\mu}(P,\lambda)\epsilon_{\psi}^{\nu}(P,\lambda)=-g^{\mu\nu}+\frac{P^{\mu}P^{\nu}}{M_{\psi}^{2}}.$
(9)

The scalar amplitudes $A_{1}$ and $A_{2}$ describe the decay process. These
amplitudes are computed in the effective field framework by expansion with
respect to small relative heavy quark velocity $v$ and small ratio
$\lambda^{2}\sim\Lambda/m_{Q}$. The amplitude $A_{2}$ is suppressed as

$A_{2}/A_{1}\sim\lambda^{2}.$ (10)

The decay width is conveniently described in terms of two combinations

$\mathcal{G}_{M}=A_{1}+A_{2},~{}\
\mathcal{G}_{E}=A_{1}+\frac{M_{\psi}^{2}}{4m_{N}^{2}}A_{2},$ (11)

and reads

$\Gamma[J/\psi\rightarrow
p\bar{p}]=\frac{M_{\psi}\beta}{12\pi}\left(\left|\mathcal{G}_{M}\right|^{2}+\frac{2m_{N}^{2}}{M_{\psi}^{2}}\left|\mathcal{G}_{E}\right|^{2}\right).$
(12)

The ratio $\left|\mathcal{G}_{E}\right|/\left|\mathcal{G}_{M}\right|$ can be
measured through the angular behaviour of the cross section
$e^{+}e^{-}\rightarrow J/\psi\rightarrow p\bar{p}$, see e.g. Ref.[31].
Available data indicate that the contribution with the amplitude
$\mathcal{G}_{E}$ is quite small (about $10\%$) and the value of width is
dominated by the amplitude $\left|\mathcal{G}_{M}\right|^{2}$ .

To the leading-order approximation

$\mathcal{G}_{M}\simeq
A_{1}^{\text{lo}}+\mathcal{O}(\alpha_{s})+\mathcal{O}(v^{2})+\mathcal{O}(\lambda^{2}),$
(13)

where in the rhs we indicate various possible corrections. These are next-to-
leading (NLO) order contribution with respect to $\alpha_{s}$, relativistic
corrections of relative order $v^{2}$ and power corrections associated with
the higher twist light-cone distribution amplitudes, respectively. The NLO
radiative QCD corrections and relativistic corrections are not known yet. The
power suppressed contribution have been estimated in Ref.[21] and
corresponding numerical effect is moderate, about $25\%$. In this paper we
want to study the relativistic corrections to the amplitude
$A_{1}^{\text{lo}}$.

The amplitude $A_{1}^{\text{lo}}$ is described by the sum of two different
contributions: the hard and electromagnetic one, which correspond to the
annihilation into three gluons or one hard photon as shown in Fig.1. The
electromagnetic amplitude is relatively small, however it can provide
numerical impact through the interference with QCD amplitude Ref.[21]. In this
paper, for simplicity, the electromagnetic contribution is not considered.

Figure 1:  Hadronic $a)$ and electromagnetic $b)$ contributions to the decay
amplitudes.

The leading-order factorisation formula for the amplitude $A_{1}^{\text{lo}}$
reads [11, 12]

$A_{1}^{(0)}=\frac{\sqrt{2M_{\psi}}\left\langle
0\right|O\left|J/\psi\right\rangle}{~{}m_{c}^{2}}\frac{~{}f_{N}^{~{}2}}{m_{c}^{4}}~{}(\pi\alpha_{s})^{3}~{}\frac{10}{81}~{}J_{0},$
(14)

where the NRQCD matrix element is defined in the standard way

$\left\langle 0\right|O\left|J/\psi\right\rangle=\left\langle
0\right|\chi^{{\dagger}}\
\boldsymbol{\sigma}\cdot\boldsymbol{\epsilon}(\lambda)\psi\left|J/\psi(\lambda)\right\rangle=\sqrt{\frac{N_{c}}{2\pi}}R_{10}(0),$
(15)

where $R_{10}(0)$ is the radial charmonium wave function at the origin. The
constant $f_{N}$ has the dimension of the square of the mass and describes the
normalization of the matrix element of the nucleon. The decay amplitude is
sensitive to nucleon structure that is described by the dimensionless
collinear convolution integral $J_{0}$, which reads

$\displaystyle J_{0}$ $\displaystyle=\frac{1}{4}\int Dx_{i}\int Dy_{i}\
\frac{1}{x_{1}x_{2}x_{3}}\frac{1}{y_{1}y_{2}y_{3}}\frac{y_{1}x_{3}}{D_{1}D_{3}}$
(16) $\displaystyle\ \ \ \ \ \ \ \ \ \ \ \
\times\left\\{\varphi_{3}(y_{123})~{}\varphi_{3}(x_{123})+\frac{1}{2}\left(\varphi_{3}(y_{123})+\varphi_{3}(y_{321})\right)\left(\varphi_{3}(x_{123})+\varphi_{3}(x_{321})\right)\right\\},$

where $D_{i}=x_{i}+y_{i}-2x_{i}y_{i}$ and we used short notation for the
twist-3 nucleon LCDA
$\varphi_{3}(x_{123})\equiv\varphi_{3}(x_{1},x_{2},x_{3})$. This
nonperturbative function depends on the quark momentum fractions $0<x_{i}<1$
and the factorisation scale, which is not shown explicitly. The light-cone
fractions $x_{i}$ satisfy the momentum conservation condition
$x_{1}+x_{2}+x_{3}=1$. Therefore the measure of the convolution integrals in
Eq.(16) includes the $\delta$-function

$Dx_{i}=dx_{1}dx_{2}dx_{3}\delta(1-x_{1}-x_{2}-x_{3}).$ (17)

The nucleon LCDA $\varphi_{3}$ is defined as the matrix element of a 3-quark
operator and is well known in the literature, see e.g. Refs.[12, 32]. For a
convenience, the definition and some important details are given in Appendix.

Our task is to calculate the relativistic corrections associated with the
higher order matrix elements in NRQCD

$\left\langle 0\right|O_{n}\left|J/\psi\right\rangle=\left\langle
0\right|\chi^{{\dagger}}\boldsymbol{\sigma}\cdot\boldsymbol{\epsilon}(\lambda)\left(-\frac{i}{2}\overleftrightarrow{\boldsymbol{D}}\right)^{2n}\psi\left|J/\psi(\lambda)\right\rangle,\
\ $ (18)

where $\left\langle 0\right|O_{0}\left|J/\psi\right\rangle\equiv\ \left\langle
0\right|O\left|J/\psi\right\rangle$. Usually, it is convenient to describe the
contribution of these matrix elements introducing the ratio

$\left\langle\boldsymbol{q}^{2n}\right\rangle=\frac{\left\langle
0\right|O_{n}\left|J/\psi\right\rangle}{\left\langle
0\right|O_{0}\left|J/\psi\right\rangle}.$ (19)

The set of the operators in Eq.(18) does not represent the complete set of the
all possible NRQCD operators. The operator $O_{2}$ provides the correct
description of the relativistic correction of relative order $v^{2}$. But to
describe relativistic corrections of a higher order, many more different
operators are used, for example, colour-octet ones. However, resumming the
contributions of all orders associated with the subset (19) is useful because
it allows one to study the convergence of the nonrelativistic expansion.

The value of the NLO matrix element (19) can be exactly found with the help of
Kapustin-Gremm relation Ref.[33]

$\left\langle 0\right|O_{1}\left|J/\psi\right\rangle=m_{c}E_{J/\psi}\
\left\langle 0\right|O_{0}\left|J/\psi\right\rangle,$ (20)

where $E_{J/\psi}$ denotes the binding energy. This gives

$\left\langle\boldsymbol{q}^{2}\right\rangle=m_{c}E_{J/\psi}.$ (21)

In Ref.[34] it was found that the higher order operators satisfy approximate
relation

$\left\langle\boldsymbol{q}^{2n}\right\rangle=\left\langle\boldsymbol{q}^{2}\right\rangle^{n}+\mathcal{O}(v^{2}).$
(22)

Neglecting the higher order corrections in the rhs of Eq.(22) allows one to
resumm the power series of the higher order terms with the matrix elements
(18) to all orders. Such resummation includes those relativistic corrections
that are contained in the quark-antiquark quarkonium wave function in the
leading potential model for the wave function [30].

The resummation of the relativistic corrections associated with the matrix
elements (19) is already considered for various processes in Refs.[23, 24, 30]
and this technique can be easily adapted for calculations of exclusive
hadronic decays.

##  3 Matching and all order resummation

In order to calculate the decay amplitude we use the covariant spin-projector
technique [35]. The generalisation of this technique to higher orders in
$v^{2}$ is considered in Refs.[23, 24] . Here we briefly describe the NRQCD
matching and clarify some features of the present calculation.

Our task is to compute the hard coefficient functions in front of higher order
operators (18). The matching can be conveniently done using $Q\bar{Q}$-state
as initial state instead of charmonium state. The corresponding heavy quark
and antiquark have the momenta $p$ and $p^{\prime}$, respectively

$p=(E,\boldsymbol{q}),\ p^{\prime}=(E,-\boldsymbol{q}),\ \ \
E=\sqrt{m_{Q}^{2}+\boldsymbol{q}^{2}}.$ (23)

Then total and relative momenta reads

$P=p+p^{\prime},\ \ q=\frac{1}{2}(p-p^{\prime}).$ (24)

In the rest frame

$P=(2E,\boldsymbol{0}),\ \ q=(0,\boldsymbol{q}).$ (25)

Let the perturbative amplitude is described as

$A_{QQ}=\bar{v}_{Q}(p^{\prime},s^{\prime})\ \hat{A}_{QQ}\ u_{Q}(p,s),$ (26)

where $\bar{v}_{Q}$ and $u_{Q}$ denote heavy quark spinors. The amplitude
projected onto $S$-wave state of heavy quark-antiquark $Q\bar{Q}(^{3}S_{1})$
reads

$A_{QQ}(^{3}S_{1})=\ \text{Tr}\left[\Pi_{1}\hat{A}_{QQ}(x_{i},y_{i})\right],$
(27)

where the spin-triplet projector $\Pi_{1}$ is given by [23]:

$\Pi_{1}=\frac{-1}{2\sqrt{2}E(E+m_{Q})}\left(\frac{1}{2}\hbox to0.0pt{\hbox
to7.80904pt{\hfil/\hfil}\hss}P+\hbox to0.0pt{\hbox
to5.00002pt{\hfil$q$\hfil}\hss}/+m_{Q}\right)\frac{\hbox to0.0pt{\hbox
to7.80904pt{\hfil/\hfil}\hss}P+2E}{4E}\hbox to0.0pt{\hbox
to5.00002pt{\hfil$\epsilon$\hfil}\hss}/\left(\frac{1}{2}\hbox to0.0pt{\hbox
to7.80904pt{\hfil/\hfil}\hss}P-\hbox to0.0pt{\hbox
to5.00002pt{\hfil$q$\hfil}\hss}/-m_{Q}\right)\otimes\frac{\mathbf{1}}{\sqrt{N_{c}}}.$
(28)

This progector is relativistically normalised

$\text{Tr[}\Pi_{1}\Pi_{1}^{{\dagger}}\text{]}=4E^{2}.$ (29)

In order to project the amplitude onto state with $L=0$ the expression in
Eq.(27) must be averaged over the angles of the momentum $\boldsymbol{q}$

$\left.A_{QQ}(\boldsymbol{q}^{2})\right|_{{}^{3}S_{1},\
L=0}=\bar{A}_{QQ}(\boldsymbol{q}^{2})=\frac{1}{4\pi}\int d\Omega\
\text{Tr}\left[\Pi_{1}\hat{A}_{QQ}(x_{i},y_{i})\right].$ (30)

In NRQCD the amplitude reads

$A_{\text{NRQCD}}=\sqrt{2M_{\psi}}\sum_{n\geq 0}c_{n}\ \left\langle
0\right|O_{n}\left|J/\psi\right\rangle,$ (31)

where $\sqrt{2M_{\psi}}$ arises from relativistic normalisation. The
coefficients $c_{n}$ can be obtained from the amplitude (30) using that

$\left\langle 0\right|O_{n}\left|Q\bar{Q}\right\rangle=\sqrt{2N_{c}}2E\
\boldsymbol{q}^{2n},$ (32)

and

$\bar{A}_{QQ}(\boldsymbol{q}^{2})=\sum_{n}c_{n}\left\langle
0\right|O_{n}\left|Q\bar{Q}\right\rangle=\sum_{n}c_{n}\sqrt{2N_{c}}2E\
\boldsymbol{q}^{2n}.$ (33)

Using this one finds

$c_{n}=\frac{1}{\sqrt{2N_{c}}}\frac{1}{n!}\frac{\partial}{\partial\boldsymbol{q}^{2n}}\frac{\bar{A}_{QQ}(\boldsymbol{q}^{2})}{2E}.$
(34)

Substituting this into Eq.(31) and using Eq.(22) one finds

$\displaystyle A_{\text{NRQCD}}\simeq\sqrt{2M_{\psi}}\left\langle
0\right|O\left|J/\psi\right\rangle\sum_{n}\frac{1}{n!}\left\langle\boldsymbol{q}^{2}\right\rangle^{n}\frac{\partial^{n}}{\partial\boldsymbol{q}^{2n}}\frac{\bar{A}_{QQ}(^{3}S_{1})}{\sqrt{2N_{c}}2E}$
(35) $\displaystyle=\sqrt{2M_{\psi}}\left\langle
0\right|O\left|J/\psi\right\rangle\frac{\bar{A}_{QQ}(\left\langle\boldsymbol{q}^{2}\right\rangle)}{\sqrt{2N_{c}}2E}.$
(36)

The specific of our calculation is that the amplitude
$\bar{A}_{QQ}(^{3}S_{1})$ also includes the nonperturbative matrix elements
describing the couplings with the nucleons. The amplitude $Q\bar{Q}\rightarrow
p\bar{p}$ is desribed by the diagrams in Fig. 2.

Figure 2:  The diagrams, which describe $Q\bar{Q}\rightarrow p\bar{p}$
subprocess.

Each diagram can be described as consisting of a heavy quark subdiagram
$\Gamma_{Q}$ describing the $Q\bar{Q}$ annihilation into three gluons, and a
residual part describing the transition of three virtual gluons to the
$p\bar{p}$ state. The expression for each diagram can be written in the
following form

$\displaystyle A_{QQ}^{icj}=\int Dx_{i}\int Dy_{i}\ \
\Delta_{1}^{\mu_{1}\mu_{1}^{\prime}}\Delta_{2}^{\mu_{2}\mu_{2}^{\prime}}\Delta_{3}^{\mu_{3}\mu_{3}^{\prime}}\
B_{\mu_{1}\mu_{2}\mu_{3}}(x_{i},y_{i})\
\bar{v}_{Q}(p^{\prime},s^{\prime})\Gamma_{Q}^{\mu_{i}^{\prime}\mu_{c}^{\prime}\mu_{j}^{\prime}}u_{Q}(p,s),$
(37)

where $\bar{v}_{Q}$ and $u_{Q}$ denote heavy quark spinors,
$\Delta_{i}^{\mu_{i}\mu_{i}^{\prime}}$ denotes the gluon propagator (we use
Feynman gauge and do not write the colour indices for simplicity)

$\Delta_{i}^{\mu_{i}\mu_{i}^{\prime}}=\frac{(-i)g^{\mu_{i}\mu_{i}^{\prime}}}{(k_{i}+k_{i}^{\prime})^{2}},\
k_{i}=x_{i}k,\ \ k_{i}^{\prime}=y_{i}k^{\prime},$ (38)

where the proton and antiproton momenta reads

$k\simeq E\ n,\ \ k^{\prime}\simeq E\ \bar{n},\ \ k^{2}=k^{\prime 2}=0.$ (39)

The heavy quark subdiagram reads

$\Gamma_{Q}^{\mu_{i}^{\prime}\mu_{c}^{\prime}\mu_{j}^{\prime}}(x_{i},y_{i})=i^{2}(ig)^{3}\frac{\gamma^{\mu_{i}}(-\hbox
to0.0pt{\hbox to7.80904pt{\hfil/\hfil}\hss}P/2+\hbox to0.0pt{\hbox
to5.00002pt{\hfil$q$\hfil}\hss}/+\hbox to0.0pt{\hbox
to5.52084pt{\hfil/\hfil}\hss}k_{i}+\hbox to0.0pt{\hbox
to5.52084pt{\hfil/\hfil}\hss}k_{i}^{\prime}+m_{Q})\gamma^{\mu_{c}}(\hbox
to0.0pt{\hbox to7.80904pt{\hfil/\hfil}\hss}P/2+\hbox to0.0pt{\hbox
to5.00002pt{\hfil$q$\hfil}\hss}/-\hbox to0.0pt{\hbox
to5.52084pt{\hfil/\hfil}\hss}k_{j}-\hbox to0.0pt{\hbox
to5.52084pt{\hfil/\hfil}\hss}k_{j}^{\prime}+m_{Q})\gamma^{\mu_{j}}}{\left[(k_{i}+k_{i}^{\prime}-P/2+q)^{2}-m_{Q}^{2}\right]\left[(P/2+q-k_{j}-k_{j}^{\prime})^{2}-m_{Q}^{2}\right]}.$
(40)

The set of indices $\\{\mu_{i},\mu_{c},\mu_{j}\\}$ describes six possible
transpositions of the indices $\\{\mu_{1},\mu_{2},\mu_{3}\\}$, which
correspond to the different diagrams: $\\{\mu_{1},\mu_{2},\mu_{3}\\}$,
$\\{\mu_{3},\mu_{2},\mu_{1}\\}$, $\\{\mu_{2},\mu_{1},\mu_{3}\\}$,
$\\{\mu_{3},\mu_{1},\mu_{2}\\}$, $\\{\mu_{1},\mu_{3},\mu_{2}\\}$,
$\\{\mu_{2},\mu_{3},\mu_{1}\\}$. Therefore each diagram can be coded by the
set $icj$ and the sum of these diagrams gives the total result.

The term $B_{\mu_{1}\mu_{2}\mu_{3}}(x_{i},y_{i})$ in Eq.(37) describes the
trace, which occur after contractions of the Dirac indices $\alpha_{i}$ and
$\alpha_{i}^{\prime}$ of the proton and antiproton light-cone matrix elements
with the $\gamma$-matrices of the quark-gluon vertices

$\displaystyle\left\langle
0\right|\left[O_{\text{tw3}}\right]_{\alpha_{1}^{\prime}\alpha_{2}^{\prime}\alpha_{3}^{\prime}}~{}\left|\bar{p}(k^{\prime})\right\rangle~{}\left[\left(ig\right)\gamma_{\mu_{1}}\right]_{\alpha_{1}^{\prime}\alpha_{1}}\left[\left(ig\right)\gamma_{\mu_{2}}\right]_{\alpha_{2}^{\prime}\alpha_{2}}\left[\left(ig\right)\gamma_{\mu_{3}}\right]_{\alpha_{3}^{\prime}\alpha_{3}}~{}\left\langle
0\right|\left[O_{\text{tw3}}\right]_{\alpha_{1}\alpha_{2}\alpha_{3}}~{}\left|p(k)\right\rangle$
(41) $\displaystyle=\
\text{FT}\left[B_{\mu_{1}\mu_{2}\mu_{3}\rho}(x_{i},y_{i})\right].$ (42)

The symbol “FT” denotes Fourier transformations, which occur from the light-
cone matrix elements (the definition of the light-cone operator
$O_{\text{tw3}}$ is given in Appendix A). The nucleon matrix elements are
described by the convenient combinations of LCDAs $V_{1}$,$\ A_{1}$ and
$T_{1}$, which are defined as

$\displaystyle V_{1}(x_{123})$
$\displaystyle=\frac{1}{2}\left(\varphi_{3}(x_{123})+\varphi_{3}(x_{213})\right),$
(43) $\displaystyle A_{1}(x_{123})$
$\displaystyle=\frac{1}{2}\left(\varphi_{3}(x_{123})-\varphi_{3}(x_{213})\right),$
(44)
$T_{1}(x_{123})=\frac{1}{2}\left(\varphi_{3}(x_{132})+\varphi_{3}(x_{231})\right).$
(45)

The $B$-term (42) is universal for all diagrams. The calculation of the Dirac
contractions yields

$B_{\mu_{1}\mu_{2}\mu_{3}}(x_{i},y_{i})=\bar{N}_{\bar{n}}\gamma_{\bot}^{\rho}{V}_{n}(ig)^{3}\frac{1}{4}(k\cdot
k^{\prime})\ ~{}T_{\mu_{1}\mu_{2}\mu_{3}\rho}(x_{i},y_{i}),$ (46)

with

$\displaystyle T_{\mu_{1}\mu_{2}\mu_{3}\rho}(x_{i},y_{i})$ $\displaystyle=\
g_{\mu_{1}\mu_{2}}^{\bot}g_{\rho\mu_{3}}^{\bot}\left\\{V_{1}(x_{123})V_{1}(y_{123})+A_{1}(x_{123})A_{1}(y_{123})\right\\}$
(47)
$\displaystyle-i\varepsilon_{\mu_{1}\mu_{2}}^{\bot}i\varepsilon_{\mu_{3}\rho}^{\bot}\left\\{V_{1}(y_{123})A_{1}(x_{123})+A_{1}(y_{123})V_{1}(x_{123})\right\\}$
$\displaystyle+G_{\mu_{1}\mu_{2}\alpha\sigma}^{\bot}G_{\rho\mu_{3}}^{\bot\alpha\sigma}~{}T_{1}(y_{123})T_{1}(x_{123}).$

where we used the following short notations

$g_{\alpha\beta}^{\bot}\
=g_{\alpha\beta}-\frac{1}{2}\left(n_{\alpha}\bar{n}_{\beta}+n_{\beta}\bar{n}_{\alpha}\right),\
\
i\varepsilon_{\alpha\beta}^{\bot}=\frac{1}{2}i\varepsilon_{\alpha\beta\sigma\lambda}n^{\sigma}\bar{n}^{\lambda},$
(48)
$G_{\alpha\beta\sigma\rho}^{\bot}=g_{\alpha\sigma}^{\bot}g_{\beta\rho}^{\bot}+g_{\alpha\rho}^{\bot}g_{\beta\sigma}^{\bot}-g_{\sigma\rho}^{\bot}g_{\alpha\beta}^{\bot}.\
$ (49)

For the projected nucleon spinors $\bar{N}(k)$ and are $V(k^{\prime})$ we used

$\bar{N}(k)\frac{\hbox to0.0pt{\hbox to6.00235pt{\hfil/\hfil}\hss}n\hbox
to0.0pt{\hbox to5.00002pt{\hfil$\bar{n}$\hfil}\hss}/}{4}=\bar{N}_{\bar{n}},\ \
\frac{\hbox to0.0pt{\hbox to6.00235pt{\hfil/\hfil}\hss}n\hbox to0.0pt{\hbox
to5.00002pt{\hfil$\bar{n}$\hfil}\hss}/}{4}V(k^{\prime})=V_{n}.$ (50)

Therefore the expression for the amplitude (30) reads

$i\bar{A}_{QQ}(\boldsymbol{q}^{2})=\frac{1}{2}\frac{1}{4\pi}\int
d\Omega\sum_{\\{icj\\}}\int Dx_{i}\int Dy_{i}\
T_{\mu_{1}\mu_{2}\mu_{3}}\Delta_{1}^{\mu_{1}\mu_{1}^{\prime}}\Delta_{2}^{\mu_{2}\mu_{2}^{\prime}}\Delta_{3}^{\mu_{3}\mu_{3}^{\prime}}\text{Tr}\left[\Pi_{1}\hat{A}_{QQ}^{\mu_{i}^{\prime}\mu_{c}^{\prime}\mu_{j}^{\prime}}\right],$
(51)

where we also show explicitly the symmetry factor $1/2$. Calculating the trace
and contracting the indices we obtain the following result111We used the
package Feyncalc [36].

$\bar{A}_{QQ}(\boldsymbol{q}^{2})=\sqrt{2N_{c}}\frac{2}{m_{Q}}\frac{f_{N}^{2}}{(m_{Q}^{2}+\boldsymbol{q}^{2})^{2}}\
(\pi\alpha_{s})^{3}\ \frac{10}{81}\frac{m_{Q}+E}{2m_{Q}}\
J(\boldsymbol{q}^{2}/m_{Q}^{2}),$ (52)

where the dimensionless collinear integral reads

$J(\boldsymbol{v}^{2})=\frac{1}{32}\frac{1}{f_{N}^{2}}\int\
\frac{Dx_{i}}{x_{1}x_{2}x_{3}}\int
Dy_{i}\frac{1}{y_{1}y_{2}y_{3}}\left\\{\frac{A(x_{i},y_{i})}{D_{1}D_{3}}+\frac{B(x_{i},y_{i})}{D_{1}D_{2}}+\frac{C(x_{i},y_{i})}{D_{2}D_{3}}\right\\}.\
$ (53)

The integrand in Eq.(53) has the following structure

$D_{i}=x_{i}+y_{i}-2x_{i}y_{i}=x_{i}(1-y_{i})+y_{i}(1-x_{i})\geq 0.$ (54)

The functions $A,B$ and $C$ also depend on
$\boldsymbol{v}^{2}=\boldsymbol{q}^{2}/m_{Q}^{2}$, which is not explicitly
shown for simplicity. They can be presented in the following way

$\displaystyle A(x_{123},y_{123})$
$\displaystyle=\sum_{i=0}^{4}A_{k}(x_{123},y_{123})\ I_{k}(1,3),\ $ (55)
$\displaystyle B(x_{123},y_{123})$
$\displaystyle=\sum_{i=0}^{4}B_{k}(x_{123},y_{123})\ I_{k}(1,2),$ (56)
$\displaystyle\ C(x_{123},y_{123})$
$\displaystyle=\sum_{i=0}^{4}C_{k}(x_{123},y_{123})\ I_{k}(2,3).$ (57)

These expressions include the angular integrals $I_{k}$ which appears due to
the integration over $d\Omega=d\varphi\ d\cos\theta$. The polar integration
can be easily computed and the remaining azimuth integrals read

$I_{k}[ij]=\frac{1}{2}\int_{-1}^{1}d\eta\
\frac{|\boldsymbol{v}|^{k}\eta^{k}}{\left[1+|\boldsymbol{v}|a_{i}\eta\right]\left[1-|\boldsymbol{v}|a_{j}\eta\right]},\
$ (58)

where we used $\ \eta\equiv\cos\theta,\ \
|\boldsymbol{v}|=|\boldsymbol{q}|/m_{Q}$ and convenient short notation

$a_{i}=\frac{m_{Q}}{E}\frac{x_{i}-y_{i}}{D_{i}}<1.$ (59)

The numerator in (58) arises from the trace and contractions in Eq.(51). The
denominators
$\left[1+|\boldsymbol{v}|a_{i}\eta\right]\left[1-|\boldsymbol{v}|a_{j}\eta\right]$
occur from the heavy quark propagators. The angular dependence occurs from the
scalar products

$(qk)=-Eq_{z}=-E|\boldsymbol{q}|\cos\theta,\ \
(qk^{\prime})=Eq_{z}=E|\boldsymbol{q}|\cos\theta.$ (60)

From the definition (58) it follows that

$I_{0}\sim\mathcal{O}(1),\ \ I_{1,2}\sim\mathcal{O}(\boldsymbol{v}^{2}),\ \ \
I_{3,4}\sim\mathcal{O}(\boldsymbol{v}^{4}).$ (61)

If one neglects by the contribution form the denominators setting
$a_{i}=a_{j}=0$ then

$I_{0}[ij]=1,\ I_{1}[ij]=0,\ I_{2}[ij]=\frac{\boldsymbol{v}^{2}}{3},\
I_{3}[ij]=0,\ I_{4}=\frac{\boldsymbol{v}^{2}}{5}.$ (62)

One can also easily find that

$I_{k}[ji]=(-1)^{k}I_{k}[ij].$ (63)

These integrals can be computed analytically that gives

$\displaystyle I_{0}[ij]$
$\displaystyle=\frac{1}{2}\frac{1}{a_{i}+a_{j}}\left(\ln\frac{1+|\boldsymbol{v}|a_{i}}{1-|\boldsymbol{v}|a_{i}}+\ln\frac{1+|\boldsymbol{v}|a_{j}}{1-|\boldsymbol{v}|a_{j}}\right)$
$\displaystyle=\frac{1}{a_{i}+a_{j}}\sum_{n=0}^{\infty}\frac{|\boldsymbol{v}|^{n}}{n+1}\left(a_{j}^{n+1}+(-1)^{n}a_{i}^{n+1}\right)\frac{1}{2}\left[1+(-1)^{n}\right],$
(64)

and for $k>0$

$I_{k}[ij]=\frac{1}{2}\frac{1}{a_{i}+a_{j}}\frac{(-1)^{k}}{a_{i}^{k-1}}\left(\frac{1}{|\boldsymbol{v}|a_{i}}\ln\frac{1+|\boldsymbol{v}|a_{i}}{1-|\boldsymbol{v}|a_{i}}-\sum_{l=0}^{k-1}|\boldsymbol{v}|^{l}a_{i}^{l}(-1)^{l}\frac{1+(-1)^{l}}{l+1}\right)$
(65)
$+\frac{1}{2}\frac{1}{a_{i}+a_{j}}\frac{1}{a_{j}^{k-1}}\left(\frac{1}{|\boldsymbol{v}|a_{i}}\ln\frac{1+|\boldsymbol{v}|a_{i}}{1-|\boldsymbol{v}|a_{i}}-\sum_{l=0}^{k-1}|\boldsymbol{v}|^{l}a_{j}^{l}\frac{1+(-1)^{l}}{l+1}\right)$
(66)
$=\frac{|\boldsymbol{v}|^{k}}{a_{i}+a_{j}}\sum_{n=0}^{\infty}\frac{|\boldsymbol{v}|^{n}}{n+1+k}\left(a_{j}^{n+1}+(-1)^{n}a_{i}^{n+1}\right)\frac{1}{2}\left[1+(-1)^{n+k}\right].$
(67)

The coefficients $A_{k},B_{k}$ and $C_{k}$ in Eqs.(55)-(57) depend on the
nucleon LCDAs and their general structure reads ($X=\\{A,B,C\\}$)

$\displaystyle X_{k}(x_{123},y_{123})$
$\displaystyle=\left[X_{k}\right]_{VV}\left\\{V_{1}(x_{123})V_{1}(y_{123})+A_{1}(x_{123})A_{1}(y_{123})\right\\}$
$\displaystyle+\left[X_{k}\right]_{AV}\left\\{A_{1}(x_{123})V_{1}(y_{123})+V_{1}(x_{123})A_{1}(y_{123})\right\\}$
$\displaystyle+\left[X_{k}\right]_{TT}\ T_{1}(x_{123})T_{1}(y_{123}),$ (68)

where coefficients $\left[X_{k}\right]_{VV,AV,TT}$ are polynomial functions of
the momentum fractions $x_{i}$ and $y_{i}$. Their analytical expressions are
presented in Appendix B.

The convolution integral $J$ in Eq.(53) is well defined that can be easily
understood taking into account the general behaviour of the LCDAs

$Z(x_{123})=x_{1}x_{2}x_{3}\times(\text{polynomial function in }x_{i}),\
Z=\left\\{V_{1},A_{1},T_{1}\right\\}.$ (69)

Using the expression in Eq.(52) one finds the NRQCD amplitude $A_{1}$ defined
in Eq.(8)

$\left(A_{1}\right)_{\text{NRQCD}}\ =\frac{\sqrt{2M_{\psi}}\left\langle
0\right|O\left|J/\psi\right\rangle}{m_{c}^{2}\sqrt{1+\left\langle\boldsymbol{v}^{2}\right\rangle}}\frac{f_{N}^{2}}{m_{c}^{4}(1+\left\langle\boldsymbol{v}^{2}\right\rangle)^{2}}\
(\pi\alpha_{s})^{3}\
\frac{10}{81}\left(1+\sqrt{1+\left\langle\boldsymbol{v}^{2}\right\rangle}\right)\frac{1}{2}\
J(\left\langle\boldsymbol{v}^{2}\right\rangle).$ (70)

In the limit $\left\langle\boldsymbol{v}^{2}\right\rangle\rightarrow 0$ this
result reproduces the well known leading-order approximation in Eq.(14). The
obtained expression depends on the power of heavy quark mass $m_{c}^{-6}$ as
required by the scale properties of the amplitude. The charm mass in the
factor $\sqrt{2M_{\psi}}$ is usually calculated as $M_{\psi}\simeq
2m_{c}\sqrt{1+\left\langle\boldsymbol{v}^{2}\right\rangle}$. The factor
$(1+\left\langle\boldsymbol{v}^{2}\right\rangle)^{-2}$ is naturally provided
by the hard propagators in the diagrams. The essential part of the
relativistic corrections is included in the convolution integral, which
depends on the twist-3 LCDAs, which describe the nonperturbtive overlaps with
outgoing nucleons. Therefore the total effect of the relativistic corrections
also depends on the nonpertubative structure of nucleon.

To see numerical impact of the obtained relativistic corrections we need to
compute the convolution integral $J$. This integral is not simple and, in the
general case, can only be calculated numerically. In order to understand
better the dependence of the integral on the model of the LCDA we consider two
different models. As an example of relatively simple case, consider the
asymptotic LCDA, which is defined as

$V_{1}^{as}=T_{1}^{as}=120\ f_{N}\ x_{1}x_{2}x_{3},\ \ \ A_{1}^{as}=0.\ \ $
(71)

Appropriate analysis allows better understanding numerical values of various
convolution integrals. For a more realistic description we use ABO-model [37],
which provides a reliable description of various data.

For our consideration, we use $m_{c}=1.4$ GeV for the pole heavy quark mass
and use the estimate for the binding energy $J/\psi$ from the Ref.[24]

$E_{J/\psi}=0.306\text{ GeV},$ (72)

that gives

$\left\langle\boldsymbol{v}^{2}\right\rangle_{J/\psi}\simeq 0.225\ .$ (73)

The value
$\left\langle\boldsymbol{v}^{2}\right\rangle_{\psi(2S)}\equiv\left\langle\boldsymbol{v}^{2}\right\rangle_{\psi^{\prime}}$
can be estimated as

$\left\langle\boldsymbol{v}^{2}\right\rangle_{\psi^{\prime}}=\frac{E_{\psi^{\prime}}}{m_{c}}\simeq\frac{M_{\psi^{\prime}}-M_{J/\psi}+E_{J/\psi}}{(M_{J/\psi}-E_{J/\psi})/2}=0.64,$
(74)

which is qute large due to the large difference
$M_{\psi^{\prime}}-M_{J/\psi}\simeq 589$ MeV.

###  3.1 Relativistic corrections with asymptotic nucleon LCDA

The asymptotic DA (71) is not a realistic model, but the use of this
approximation makes it possible to simplify analytical expressions and study
the properties of convolution integrals. At first step we also calculate the
correction of relative order $v^{2}$ only.

The results in Eqs. (64) and (67) yield

$I_{0}[ij]=1+\left\langle\boldsymbol{v}^{2}\right\rangle\frac{1}{3}\left(a_{i}^{2}-a_{i}a_{j}+a_{j}^{2}\right)+\mathcal{O}(v^{4})\simeq
1+\left\langle\boldsymbol{v}^{2}\right\rangle I_{0}^{(1)}[ij],$ (75)
$\displaystyle I_{1}[ij]$
$\displaystyle=\frac{\left\langle\boldsymbol{v}^{2}\right\rangle}{3}(a_{j}-a_{i})+\mathcal{O}(v^{4})\simeq
I_{1}^{(1)}[ij],\ \ $ (76) $\displaystyle\ \ I_{2}[ij]$
$\displaystyle=\frac{\left\langle\boldsymbol{v}^{2}\right\rangle}{3}+\mathcal{O}(v^{4}),\
\ I_{3,4}\sim\mathcal{O}(v^{4}),$ (77)

where $a_{i}$ is defined in Eq.(59). Therefore the contributions with
$I_{3,4}$ can be neglected. It is conveninet to define

$J_{as}(\left\langle\boldsymbol{v}^{2}\right\rangle)=\frac{120^{2}}{32}\sum_{k=0}^{4}J_{as}^{(k)}(\left\langle\boldsymbol{v}^{2}\right\rangle),$
(78)

where the coefficient in front of the sum is choosen in order to get a
convenient normalisation. In what follow we are going to discuss the values of
the integrals

$J_{as}^{(k)}(\left\langle\boldsymbol{v}^{2}\right\rangle)=\frac{1}{120^{2}f_{N}^{2}}\int\
\frac{Dx_{i}}{x_{1}x_{2}x_{3}}\int\frac{Dy_{i}}{y_{1}y_{2}y_{3}}\left\\{\frac{A_{k}^{as}I_{k}[13]}{D_{1}D_{3}}+\frac{B_{k}^{as}I_{k}[12]}{D_{1}D_{2}}+\frac{C_{k}^{as}I_{k}[23]}{D_{2}D_{3}}\right\\},$
(79)

and their dependence on the parameter
$\left\langle\boldsymbol{v}^{2}\right\rangle$.

For the coefficients $A_{k}^{as}$ we get

$\displaystyle A_{k}^{as}$
$\displaystyle=\left[A_{k}\right]_{VV}\left\\{V_{1}^{as}(x_{123})V_{1}^{as}(y_{123})\right\\}+\left[A_{k}\right]_{TT}\
T_{1}^{as}(x_{123})T_{1}^{as}(y_{123})$ (80) $\displaystyle=f_{N}^{2}\
120^{2}(x_{1}x_{2}x_{3})(y_{1}y_{2}y_{3})\
\left(\left[A_{k}\right]_{VV}+\left[A_{k}\right]_{TT}\right)=f_{N}^{2}\
120^{2}(x_{1}x_{2}x_{3})(y_{1}y_{2}y_{3})\ \bar{A}_{k}^{as},$ (81)

with

$\bar{A}_{k}^{as}=\left[A_{k}\right]_{VV}+\left[A_{k}\right]_{TT},$ (82)

and similarly for $B_{k}$ and $C_{k}$. Substituting this in Eq.(79) we obtain

$J_{as}^{(k)}\simeq\int Dx_{i}\int
Dy_{i}\left\\{\frac{\bar{A}_{k}^{as}I_{k}[13]}{D_{1}D_{3}}+\frac{\bar{B}_{k}^{as}I_{k}[12]}{D_{1}D_{2}}+\frac{\bar{C}_{k}^{as}I_{k}[23]}{D_{2}D_{3}}\right\\}$
(83) $=\int Dx_{i}\int
Dy_{i}\frac{I_{k}[12]}{D_{1}D_{2}}\left\\{\bar{B}_{k}^{as}+\hat{P}_{23}\bar{A}_{k}^{as}+(-1)^{k}\hat{P}_{13}\bar{C}_{k}^{as}\right\\},$
(84)

where the transposition operator $\hat{P}_{ij}$ intechanges the arguments
$\\{x_{i},y_{i}\\}\leftrightarrow\\{x_{j},y_{j}\\}$, for instance

$\hat{P}_{23}f(x_{123},y_{123})=f(x_{132},y_{132}).$ (85)

Substituting the explicit expressions for the coefficients $\bar{A}_{k}^{as}$,
$\bar{B}_{k}^{as}$ and $\bar{C}_{k}^{as}$ and neglecting the higher order
terms we get

$\displaystyle\sum_{k=0}^{2}J_{as}^{(k)}(\left\langle\boldsymbol{v}^{2}\right\rangle)$
$\displaystyle=\int Dx_{i}\int Dy_{i}\
\frac{1}{D_{1}D_{2}}\left\\{24x_{1}y_{2}-8\left\langle\boldsymbol{v}^{2}\right\rangle+24\left\langle\boldsymbol{v}^{2}\right\rangle
x_{1}y_{2}(I_{0}^{(1)}[12]+I_{1}^{(1)}[12])\right.$ $\displaystyle\ \ \ \ \ \
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\left.+\left\langle\boldsymbol{v}^{2}\right\rangle\left(2x_{1}y_{2}-\frac{1}{6}\left(x_{1}+x_{2}+y_{1}+y_{2}\right)\right)\right\\}.$
(86)

The integrals, which enter in these formulae can be easily calculated
numerically

$\text{LO:}\int Dx_{i}\int Dy_{i}\ \frac{\ x_{1}y_{2}}{D_{1}D_{2}}=0.140.$
(87) $\text{NLO}:\ 1^{st}=\int Dx_{i}\int Dy_{i}\ \frac{1}{D_{1}D_{2}}=2.029,$
(88) $\displaystyle\ \text{NLO: }2^{nd}$ $\displaystyle=\int Dx_{i}\int
Dy_{i}\ \frac{x_{1}y_{2}}{D_{1}D_{2}}(I_{0}^{(1)}[12]+I_{1}^{(1)}[12])=0.012,$
(89) $\displaystyle\text{NLO: }3^{d}$ $\displaystyle=\int Dx_{i}\int Dy_{i}\
\frac{1}{D_{1}D_{2}}\left(2x_{1}y_{2}-\frac{1}{6}\left(x_{1}+x_{2}+y_{1}+y_{2}\right)\right)=-0.121.$
(90)

This gives

$\displaystyle\sum J_{as}^{(k)}(\left\langle\boldsymbol{v}^{2}\right\rangle)$
$\displaystyle=24\ast 0.140\
-8\left\langle\boldsymbol{v}^{2}\right\rangle\ast\left(2.029\right)_{1st}+24\left\langle\boldsymbol{v}^{2}\right\rangle\left(0.012\right)_{2nd}+\left\langle\boldsymbol{v}^{2}\right\rangle\left(-0.121\right)_{3d}$
(91)
$\displaystyle=3.36_{\text{lo}}-16.232\left\langle\boldsymbol{v}^{2}\right\rangle_{1st}+0.167\left\langle\boldsymbol{v}^{2}\right\rangle_{2nd+3d}=3.36-16.40\left\langle\boldsymbol{v}^{2}\right\rangle.$
(92)

Substituting the numerical values for
$\left\langle\boldsymbol{v}^{2}\right\rangle_{J/\psi}$ (73) and
$\left\langle\boldsymbol{v}^{2}\right\rangle_{\psi^{\prime}}$ (74) we obtain

$\sum
J_{as}^{(k)}(\left\langle\boldsymbol{v}^{2}\right\rangle_{J/\psi})=-0.33,$
(93) $\sum
J_{as}^{(k)}(\left\langle\boldsymbol{v}^{2}\right\rangle_{\psi^{\prime}})=-7.14.$
(94)

These results show that the relativistic correction is negative and provides
about $100\%$ and $300\%$ numerical effect for $1S$ and $2S$-states,
respectively. To better understand the origin of such a large numerical
impact, consider the structure of the expression in the equation. (91).

From this result it follows that relativistic correction is dominated by the
terms which appear from the integral (88). This integral is about an order of
magnitude larger than the leading-order integral (87). This has the natural
mathematical explanation: the both integrands in Eqs.(87) and (88) are
positive, but the leading-order integrand in (87) has in the numerator small
factor $x_{1}y_{2}$ that reduces the value of the integral (87) comparing to
one in (88). The actual values of the NRQCD parameter
$\left\langle\boldsymbol{v}^{2}\right\rangle$ for charmonium states can not
compensate the numerical enhancement of the subleading integral. This
indicates about the relative smallness of the charm quark mass.

The mechanism of the enhancement of the subleading integral does not depend on
the model of LCDA of the nucleon. The structure of leading- and subleading
integrals is closely related to the structure of the numerators of the Feynman
diagrams and can be interpreted as a specific feature of the hard perturbative
subprocess. In order to see this consider the leading-order approximation with
$q=0$, then the trace in Eq.(51) reads

$\displaystyle\text{Tr}\left[\Pi_{1}^{(\text{lo})}\Gamma_{Q}^{\mu_{i}^{\prime}\mu_{c}^{\prime}\mu_{j}^{\prime}}(x_{i},y_{i})\right]\sim\text{Tr}\left[\left(\hbox
to0.0pt{\hbox to6.22453pt{\hfil/\hfil}\hss}\omega+1\right)\hbox to0.0pt{\hbox
to5.00002pt{\hfil$\epsilon$\hfil}\hss}/_{\psi}\gamma_{\bot}^{\mu_{i}}(-\hbox
to0.0pt{\hbox to7.80904pt{\hfil/\hfil}\hss}P/2+m_{Q}+\hbox to0.0pt{\hbox
to5.52084pt{\hfil/\hfil}\hss}k_{i}+\hbox to0.0pt{\hbox
to5.52084pt{\hfil/\hfil}\hss}k_{i}^{\prime})\gamma_{\bot}^{\mu_{c}}\right.$
$\displaystyle\left.\times(\hbox to0.0pt{\hbox
to7.80904pt{\hfil/\hfil}\hss}P/2+\hbox to0.0pt{\hbox
to5.52084pt{\hfil/\hfil}\hss}k_{j}+\hbox to0.0pt{\hbox
to5.52084pt{\hfil/\hfil}\hss}k_{j}^{\prime}+m_{Q})\gamma_{\bot}^{\mu_{j}}\right]$
(95)

where we used expression from Eq.(40) and

$\ \Pi_{1}^{(\text{lo})}=\frac{-1}{2\sqrt{2}}\left(\hbox to0.0pt{\hbox
to6.22453pt{\hfil/\hfil}\hss}\omega+1\right)\hbox to0.0pt{\hbox
to5.00002pt{\hfil$\epsilon$\hfil}\hss}/_{\psi}.$ (96)

Notice that all Dirac matrices with open Lorentz indices in Eq.(95) are
transvers, because of contraction with transverse tensor
$T_{\mu_{1}\mu_{2}\mu_{3}\rho}$ in Eq.(47). Remind, that the heavy qaurk
velocity

$\omega^{\mu}=\delta_{\mu 0}=\frac{1}{2}(n^{\mu}+\bar{n}^{\mu}).$ (97)

If one neglects collinear momenta $k_{i}$, $k_{i}^{\prime}$ or $k_{j}$,
$k_{j}^{\prime}$ in the trace Eq.(95) then the trace vanishes, for instance

$\displaystyle\text{Tr}\left[\left(\hbox to0.0pt{\hbox
to6.22453pt{\hfil/\hfil}\hss}\omega+1\right)\hbox to0.0pt{\hbox
to5.00002pt{\hfil$\epsilon$\hfil}\hss}/_{\psi}\gamma_{\bot}^{\mu_{i}}(-\hbox
to0.0pt{\hbox
to7.80904pt{\hfil/\hfil}\hss}P/2+m_{Q})\gamma_{\bot}^{\mu_{c}}(\hbox
to0.0pt{\hbox to7.80904pt{\hfil/\hfil}\hss}P/2+\hbox to0.0pt{\hbox
to5.52084pt{\hfil/\hfil}\hss}k_{j}+\hbox to0.0pt{\hbox
to5.52084pt{\hfil/\hfil}\hss}k_{j}^{\prime}+m_{Q})\gamma_{\bot}^{\mu_{j}}\right]$
(98) $\displaystyle=m_{Q}\text{Tr}\left[\left(\hbox to0.0pt{\hbox
to6.22453pt{\hfil/\hfil}\hss}\omega+1\right)\hbox to0.0pt{\hbox
to5.00002pt{\hfil$\epsilon$\hfil}\hss}/_{\psi}\ \gamma_{\bot}^{\mu_{i}}(-\hbox
to0.0pt{\hbox
to6.22453pt{\hfil/\hfil}\hss}\omega+1)\gamma_{\bot}^{\mu_{c}}(\ldots)\gamma_{\bot}^{\mu_{j}}\right]$
(99) $\displaystyle=m_{Q}\text{Tr}\left[\underline{\left(\hbox to0.0pt{\hbox
to6.22453pt{\hfil/\hfil}\hss}\omega+1\right)(-\hbox to0.0pt{\hbox
to6.22453pt{\hfil/\hfil}\hss}\omega+1)}\hbox to0.0pt{\hbox
to5.00002pt{\hfil$\epsilon$\hfil}\hss}/_{\psi}\
\gamma_{\bot}^{\mu_{i}}\gamma_{\bot}^{\mu_{c}}(\ldots)\gamma_{\bot}^{\mu_{j}}\right]=0,$
(100)

where we used $[\hbox to0.0pt{\hbox
to6.22453pt{\hfil/\hfil}\hss}\omega,\gamma_{\bot}]=[\hbox to0.0pt{\hbox
to6.22453pt{\hfil/\hfil}\hss}\omega,\hbox to0.0pt{\hbox
to5.00002pt{\hfil$\epsilon$\hfil}\hss}/_{\psi}]=0$. Therefore the nontrivial
result is obtained only from the term with collinear momenta in both quark
propagators

$\text{Tr}\left[\Pi_{1}^{(\text{lo})}\hat{A}_{QQ}^{\mu_{i}^{\prime}\mu_{c}^{\prime}\mu_{j}^{\prime}}(x_{i},y_{i})\right]\sim\text{Tr}\left[\left(\hbox
to0.0pt{\hbox to6.22453pt{\hfil/\hfil}\hss}\omega+1\right)\hbox to0.0pt{\hbox
to5.00002pt{\hfil$\epsilon$\hfil}\hss}/_{\psi}\gamma_{\bot}^{\mu_{i}}(\hbox
to0.0pt{\hbox to5.52084pt{\hfil/\hfil}\hss}k_{i}+\hbox to0.0pt{\hbox
to5.52084pt{\hfil/\hfil}\hss}k_{i}^{\prime})\gamma_{\bot}^{\mu_{c}}(\hbox
to0.0pt{\hbox to5.52084pt{\hfil/\hfil}\hss}k_{j}+\hbox to0.0pt{\hbox
to5.52084pt{\hfil/\hfil}\hss}k_{j}^{\prime})\gamma_{\bot}^{\mu_{j}}\right]\sim
x_{i}y_{j}+x_{j}y_{i}.$ (101)

This explains the structure of the leading-order integral in Eq.(87).

The calculation of the subleading correction involves the relative momentum
$q$ that introduces the traces with two insertions of $q$ instead of collinear
momenta, for instance

$\text{Tr}\left[\Pi_{1}\Gamma_{Q}^{\mu_{i}^{\prime}\mu_{c}^{\prime}\mu_{j}^{\prime}}(x_{i},y_{i})\right]\sim\text{Tr}\left[\left(\hbox
to0.0pt{\hbox to6.22453pt{\hfil/\hfil}\hss}\omega+1\right)\hbox to0.0pt{\hbox
to5.00002pt{\hfil$\epsilon$\hfil}\hss}/_{\psi}\gamma_{\bot}^{\mu_{i}}\hbox
to0.0pt{\hbox to5.00002pt{\hfil$q$\hfil}\hss}/\gamma_{\bot}^{\mu_{c}}\hbox
to0.0pt{\hbox to5.00002pt{\hfil$q$\hfil}\hss}/\gamma_{\bot}^{\mu_{j}}\right]+\
\ldots\ .$ (102)

Such subleading contributions give terms of relative order $v^{2}$ but the
corresponding integral has no momentum fractions in the numerator, see
Eq.(88). As a result, the large numerical contributions are generated by the
integrals $J_{as}^{(0)}$ and $J_{as}^{(2)}$ only, which can only get
appropriate contributions from the trace and contractions in Eq.(51).
Parametrically such integrals are suppressed by small velocity, however, in
reality such suppression is not sufficiently strong in order to provide a
small numerical correction with respect to the leading-order contribution.

The enhanced integral also takes place in the relative order $v^{4}$, but the
corresponding numerical effect is already suppressed by additional factor
$\left\langle\boldsymbol{v}^{2}\right\rangle$ and for $J/\psi$ such a
contribution can be estimated at about $20-30\%$ . But for the excited state
such a contribution still remains quite large because of large value of
$\left\langle\boldsymbol{v}^{2}\right\rangle_{\psi^{\prime}}$.

In order to study the more general situation and the possible numerical
effects of higher order terms in $v^{2}$, we present in Table 1 the results
for the integrals $J_{as}^{(k)}$, which are calculated with the resummed
relativistic corrections. Remind, that the leading-order approximation is
given by the integral
$J_{as}^{(0)}(\left\langle\boldsymbol{v}^{2}\right\rangle=0)=3.36.$

Table 1:  Numerical results for integrals $J_{as}^{(k)}$ for $J/\psi$ and
$\psi^{\prime}$, upper and bottom tables, respectively.

$\left\langle\boldsymbol{v}^{2}\right\rangle_{J/\psi}=0.225$ | $J_{as}^{(0)}$ | $J_{as}^{(1)}$ | $J_{as}^{(2)}$ | $J_{as}^{(3)}$ | $J_{as}^{(4)}$ | $\sum J_{as}^{(k)}$
---|---|---|---|---|---|---
NLO | $-1.84$ | $-0.18$ | $1.89$ | $-$ | $-$ | $-0.13$
all orders | $-1.00$ | $-0.19$ | $1.66$ | $-0.01$ | $0.13$ | $0.59$

$\left\langle\boldsymbol{v}^{2}\right\rangle_{\psi^{\prime}}=0.64$ | $J_{as}^{(0)}$ | $J_{as}^{(1)}$ | $J_{as}^{(2)}$ | $J_{as}^{(3)}$ | $J_{as}^{(4)}$ | $\sum J_{as}^{(k)}$
---|---|---|---|---|---|---
NLO | $-11.68$ | $-0.50$ | $5.37$ | $-$ | $-$ | $-6.81$
all orders | $-6.78$ | $-0.69$ | $4.53$ | $-0.14$ | $0.98$ | $-2.10$

The largest numerical effects are provided by the NLO corrections in integrals
$J_{as}^{(0,2)}$ for the same reason as discussed above. The numerical effect
from the the resummation is most visible for $J_{as}^{(0)}$ but corresponding
effects are not sufficiently large and do not change the qualitative picture
discussed for NLO approximation. Therefore we can conclude that relativistic
expansion for $J/\psi$ is well convergent. The large numerical impact is only
generated by the correction of relative order $v^{2}$ and associated with the
numerical enhancement of the NLO convolution integral.

In the Table 1 (bottom part) we also show the results for the $2S$-state,
which has much larger relativistic corrections. In this case the general
structure remains quite similar: the largest numerical effect is provided by
the NLO contribution in $J_{as}^{(0,2)}$ but the higher order power
corrections also make a significant numerical contribution, reducing the total
sum by a factor of three. This clearly illustrates the expected conclusion:
the relativistic expansion for excited charmonium has very large relativistic
corrections and converges rather slowly.

###  3.2 Relativistic corrections with realistic nucleon LCDA

The value of the collinear integral also depends on the model of the LCDA
$\varphi_{3}$, which describes a distribution of the quark longitudinal
momenta at zero transverse separation. To illustrate the possible effect of
the nucleon structure, we calculate in this section the convolution integrals
with a realistic LCDA model. For that purpose we use the ABOI-model from
Ref.[37]. This model gives a reliable description of the electromagnetic form
factor data within the light-cone sum rules [38] and also provides a good
description of $J/\psi\rightarrow p\bar{p}$ decay data in the leading-order
approximation [21]. The expression for the model reads

$\displaystyle\varphi_{3}^{\text{ABO}}(x_{i})=120x_{1}x_{2}x_{3}\left\\{1+\varphi_{10}\mathcal{P}_{10}(x_{i})+\varphi_{11}\mathcal{P}_{11}(x_{i})\right.$
$\displaystyle\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\left.\phantom{120x_{1}x_{2}x_{3}}+\varphi_{20}\mathcal{P}_{20}(x_{i})+\varphi_{21}\mathcal{P}_{21}(x_{i})+\varphi_{22}\mathcal{P}_{22}(x_{i})\right\\},$
(103)

where the orthogonal polynomials $\mathcal{P}_{ij}(x_{i})$ read

$\mathcal{P}_{10}(x_{i})=21(x_{1}-x_{3}),~{}~{}\mathcal{P}_{11}(x_{i})=7(x_{1}-2x_{2}+x_{3}),$
(104)
$\displaystyle\mathcal{P}_{20}(x_{i})=\frac{63}{10}\left[3(x_{1}-x_{3})^{2}-3x_{2}(x_{1}+x_{3})+2x_{2}^{2}\right],~{}$
(105)
$\displaystyle\mathcal{P}_{21}(x_{i})=\frac{63}{2}(x_{1}-3x_{2}+x_{3})(x_{1}-x_{3}),$
(106)
$\displaystyle\mathcal{P}_{22}(x_{i})=\frac{9}{5}\left[x_{1}^{2}+9x_{2}(x_{1}+x_{3})-12x_{1}x_{3}-6x_{2}^{2}+x_{3}^{2}\right].$
(107)

The moments $\varphi_{ij}\equiv$ $\varphi_{ij}(\mu)$ are multiplicatively
renormalisable, more details about properties of the polynomials
$\mathcal{P}_{ij}$ and about the evolution of the moments can be found in
Ref.[39]. In Appendix A we also provide some useful details. In our
calculation we fix the relatively low normalisation scale $\mu^{2}=1.5$ GeV
following to Ref.[21]. Then the values of the twist-3 moments read

$\varphi_{10}=0.051,\ \varphi_{11}=0.052,\ \ \varphi_{20}=0.078,\ \
\varphi_{21}=-0.028,\ \ \varphi_{22}=0.179.$ (108)

Again, we define the integral as the sum

$J_{\text{ABO}}(\left\langle\boldsymbol{v}^{2}\right\rangle)=\frac{120^{2}}{32}\sum_{k=0}^{4}J_{\text{ABO}}^{(k)}(\left\langle\boldsymbol{v}^{2}\right\rangle),$
(109)

where the integrals $J_{\text{ABO}}^{(k)}$ are defined as

$J_{\text{ABO}}^{(k)}(\left\langle\boldsymbol{v}^{2}\right\rangle)=\frac{1}{120^{2}f_{N}^{2}}\int\
\frac{Dx_{i}}{x_{1}x_{2}x_{3}}\int\frac{Dy_{i}}{y_{1}y_{2}y_{3}}\left\\{\frac{A_{k}^{\text{ABO}}I_{k}[13]}{D_{1}D_{3}}+\frac{B_{k}^{\text{ABO}}I_{k}[12]}{D_{1}D_{2}}+\frac{C_{k}^{\text{ABO}}I_{k}[23]}{D_{2}D_{3}}\right\\}.$
(110)

The value of the leading order integral reads

$J_{\text{ABO}}^{(0)}(\left\langle\boldsymbol{v}^{2}\right\rangle=0)=5.08.$
(111)

The numerical results for the different integrals $J_{\text{ABO}}^{(k)}$ for
$J/\psi$ and $\psi^{\prime}$ are presented in Table 2.

Table 2: Numerical results for the integrals $J_{\text{ABO}}^{(k)}$ for
$J/\psi$ and $\psi^{\prime}$, upper and bottom tables, respectively.

$\left\langle\boldsymbol{v}^{2}\right\rangle_{J/\psi}=0.225$ | $J_{\text{ABO}}^{(0)}$ | $J_{\text{ABO}}^{(1)}$ | $J_{\text{ABO}}^{(2)}$ | $J_{\text{ABO}}^{(3)}$ | $J_{\text{ABO}}^{(4)}$ | $\sum J_{\text{ABO}}^{(k)}$
---|---|---|---|---|---|---
NLO | $-5.18$ | $-0.27$ | $3.66$ | $-$ | $-$ | $-1.79$
all orders | $-3.65$ | $-0.35$ | $3.30$ | $-0.03$ | $0.22$ | $-0.51$

$\ \left\langle\boldsymbol{v}^{2}\right\rangle_{\psi^{\prime}}=0.64\ $ | $J_{\text{ABO}}^{(0)}$ | $J_{\text{ABO}}^{(1)}$ | $J_{\text{ABO}}^{(2)}$ | $J_{\text{ABO}}^{(3)}$ | $J_{\text{ABO}}^{(4)}$ | $\sum J_{\text{ABO}}^{(k)}$
---|---|---|---|---|---|---
NLO | $-24.09$ | $-0.75$ | $10.42$ | $-$ | $-$ | $-14.42$
all orders | $-15.69$ | $-1.39$ | $9.31$ | $-0.30$ | $1.78$ | $-6.29$

The qualitative picture remains the same as described above, the sum of the
integrals $J_{\text{ABO}}^{(0)}$+$J_{\text{ABO}}^{(2)}$ provides the largest
numerical impact but the values of the all integrals are larger. We can
conclude that the described mechanism of the large numerical effect also works
for the realistic model of LCDAs. In case of $J/\psi$ the total sum is
negative, which indicates that the negative relativistic correction is
somewhat larger than the LO contribution.

Additional numerical impact is also provided by
$\left\langle\boldsymbol{v}^{2}\right\rangle$-dependence of the coefficient in
front of of the convolution integral in Eq.(70). The dominant numerical effect
is provided by the factor
$1/(1+\left\langle\boldsymbol{v}^{2}\right\rangle)^{2}$, which occurs from the
hard propagators. It seems that this factor can be understood as an indication
that the charmonium mass $M_{\psi}^{2}\simeq
4m_{c}^{2}(1+\left\langle\boldsymbol{v}^{2}\right\rangle)$ is a more natural
scale for the hard gluon propagators in the Feynman diagrams. The
corresponding numerical effect is smaller in comparison with one in the
convolution integral discussed above. For $J/\psi$ the modification of the
coefficient in Eq.(70) gives a reduction about $37\%$. For $\psi^{\prime}$, a
similar effect is already $67\%$.

##  4 Discussion

We presented the first study of relativistic corrections in exclusive $S$-wave
charmonium decays into proton-antiproton final state. The relativistic
corrections are calculated within the NRQCD and collinear factorisation
frameworks. We only consider the helicity conserving amplitude $A_{1}$, which
provides the dominant numerical contribution to the decay width. In this case
the baryon non-perturbative light-cone matrix elements depend on the twist-3
LCDAs only and the formula for the amplitude with relativistic corrections
includes the collinear convolution integral, which is free from any infrared
singularities.

The result obtained provides a full correction of the relative order $v^{2}$,
and also includes the summation of all orders of relativistic corrections
associated with the quark-antiquark-quarkonia wave function in the potential
model [30].

The largest numerical impact is provided by the NLO relativistic correction
with the matrix element, which can be computed using equation of motion and is
proportional to the binding energy. Due to a specific structure of the LO and
NLO hard scattering kernels the value of the NLO collinear integral is about
an order of magnitude larger than the value of LO one. The charmonium
parameter $\left\langle\boldsymbol{v}^{2}\right\rangle$ is not small enough to
compensate for this effect, and the magnitude of the resulting relativistic
correction is numerically of the same order as the LO contribution. For the
case of $J/\psi$ the resulting NLO correction is negative and therefore, the
two contributions almost cancel out. The calculated higher order corrections
are relatively small and their resummation shows that relativistic expansion
converges sufficiently well. The numerical effect also depends on the model of
the nucleon LCDA however this dependence does not change the main qualitative
conclusions about the large relativistic corrections.

A description of excited state $\psi^{\prime}$ is more challenging because of
relatively large value of
$\left\langle\boldsymbol{v}^{2}\right\rangle_{\psi^{\prime}}$. In this case
the numerical effect from relativistic corrections is much larger than the LO
contribution and is more sensitive to the higher order contributions.
Therefore a description of $\psi^{\prime}$ baryonic decays suffers from large
uncertainties, which are associated with the higher-order contributions of the
relativistic expansions.

Taking into account these large relativistic corrections, one cannot expect
the $12\%$ rule in Eq.(1) to hold in the general case. The rather good
agreement for the nucleon channel is probably an accidental consequence of
different numerical cancellations.

The large effect from relativistic corrections raises the question about a
phenomenological description of baryon decays in the effective field theory
framework. Various existing phenomenological estimates for $J/\psi$ decay
width and angular behaviour are based on the LO approximation and provide a
qualitative reliable estimates [21]. However such qualitative picture is
violated by the large value of relativistic corrections. The large negative
correction almost cancel the LO contribution, which greatly reduces the
amplitude value and, therefore, makes the description problematic. A possible
solution of this problem might be associated with the large and positive NLO
radiative correction, which is not yet known. Cancellation of radiative and
relativistic corrections could resolve the situation. The large NLO radiative
corrections are already observed in the exclusive production
$e^{+}e^{-}\rightarrow J/\psi\eta_{c}$, see e.g. Ref.[40]. Probably a similar
situation also arises in the baryonic decays. Therefore, the calculation of
the NLO radiative correction is necessary in order to better understand the
hadronic decay dynamics.

## Acknowledgements

I am grateful to A. Vairo and V. Shtabovenko for very useful discussions. This
work was supported by the Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation) Project-ID 445769443.

## Appendix

##  Appendix A Definition and properties of nucleon LCDA

In this Appendix we provide definition of the light-cone matrix element for
proton state. The formulae for the antiproton state can be obtained by charge
conjugation. In order to simplify formulas, we also use the light-cone gauge

$n\cdot A(x)=0.$ (A.1)

In this case the light-cone three-quark operator which we need for the proton
state can be defined as ($i,j,k$ are the colour indices)

$O_{\text{tw3}}=\varepsilon^{ijk}\xi_{\alpha}^{i}\left(z_{1-}\right)\
\xi_{\beta}^{j}(z_{2-})\ \xi_{\sigma}^{k}(z_{3-}),$ (A.2)

where the projected quark field $\xi(z)\equiv\frac{\hbox to0.0pt{\hbox
to5.00002pt{\hfil/\hfil}\hss}\bar{n}\hbox to0.0pt{\hbox
to4.20164pt{\hfil/\hfil}\hss}n}{4}\psi(z)$, the arguments of the fields read

$z_{i-}=\left(\bar{n}\cdot z_{i}\right)\frac{n}{2}.$ (A.3)

The proton flavour structure implies that $O_{\text{tw3}}=u\ u\ d.$

The twist-3 light-cone matrix element is defined as

$\displaystyle\left\langle
0\left|O_{\text{tw3}}(z_{1},z_{2},z_{3})\right|p(k)\right\rangle=\frac{1}{4}~{}\left[\hbox
to0.0pt{\hbox
to5.52084pt{\hfil/\hfil}\hss}k~{}C\right]_{\alpha\beta}\left[\gamma_{5}N_{\bar{n}}\right]_{\sigma}~{}\text{FT}\left[V_{1}(y_{i})\right]$
$\displaystyle\mskip 100.0mu+\frac{1}{4}\left[\hbox to0.0pt{\hbox
to5.52084pt{\hfil/\hfil}\hss}k\gamma_{5}C\right]_{\alpha\beta}\left[N_{\bar{n}}\right]_{\sigma}~{}\text{FT}\left[A_{1}(y_{i})\right]+\frac{1}{4}~{}k^{\nu}\left[i\sigma_{\mu\nu}C\right]_{\alpha\beta}\left[\gamma_{\bot}^{\mu}\gamma_{5}N_{\bar{n}}\right]_{\sigma}\text{FT}\left[T_{1}(y_{i})\right],$
(A.4)

where the projected nucleon spinor reads

$N_{\bar{n}}=\frac{\hbox to0.0pt{\hbox
to5.00002pt{\hfil$\bar{n}$\hfil}\hss}/\hbox to0.0pt{\hbox
to6.00235pt{\hfil/\hfil}\hss}n}{4}N(k),$ (A.5)

and $C$ is charge conjugation matrix. The Fourier transformation “FT” is
defined as

$\text{FT}\left[F(y_{i})\right]=\int
Dy_{i}~{}e^{-iy_{1}k_{-}z_{1+}/2-iy_{2}k_{-}z_{2+}/2-iy_{3}k_{-}z_{3+}/2}F(y_{1,}y_{2},y_{3}),$
(A.6)

with the integration measure

$Dy_{i}=dy_{1}dy_{2}dy_{3}\delta(1-y_{1}-y_{2}-y_{3}).$ (A.7)

Three LCDAs $V_{1}$, $A_{1}$ and $T_{1}$ satisfy the following properties [32]

$V_{1}(y_{2},y_{1},y_{3})=V_{1}(y_{1},y_{2},y_{3}),$ (A.8)
$A_{1}(y_{2},y_{1},y_{3})=-A_{1}(y_{1},y_{2},y_{3}),$ (A.9)
$T_{1}(y_{2},y_{1},y_{3})=T_{1}(y_{1},y_{2},y_{3}).$ (A.10)

The isospin symmetry allows one to get the following relation, see details in
Refs. [32, 41]

$T_{1}(y_{1},y_{2},y_{3})=\frac{1}{2}\left(V_{1}-A_{1}\right)(y_{1},y_{3},y_{2})+\frac{1}{2}\left(V_{1}-A_{1}\right)(y_{2},y_{3},y_{1}).$
(A.11)

It is convenient to define the following combination

$f_{N}\
\varphi_{3}(y_{1},y_{2},y_{3})=V_{1}(y_{1},y_{2},y_{3})-A_{1}(y_{1},y_{2},y_{3}),$
(A.12)

which allows one to describe the set of three LCDAs $V_{1}$, $A_{1}$ and
$T_{1}$ in terms of the one function. The coupling $f_{N}$ describes the
normalisation so that $\varphi_{3}$ is dimensionless and normalised as

$\int Dy_{i}\text{ }\varphi_{3}(y_{1},y_{2},y_{3})=1.$ (A.13)

The LCDA $\varphi_{3}$ also depends on the factorisation scale, which is not
shown for simplicity

$\varphi_{3}(y_{1},y_{2},y_{3})\equiv\varphi_{3}(y_{1},y_{2},y_{3};\mu)\text{.}$
(A.14)

The evolution properties of this LCDA was studied in Ref.[39]. The moments in
Eq.(103) are multiplicatively renormalisable and can be calculated as

$\phi_{ij}(\mu)=\phi_{ij}(\mu_{0})\left(\frac{\alpha_{s}(\mu)}{\alpha_{s}(\mu_{0})}\right)^{\gamma_{ij}/\beta_{0}},$
(A.15)

where $\beta_{0}=11-2n_{f}/3$ and $\gamma_{ij}$ are the corresponding
anomalous dimensions

$\gamma_{10}=\frac{20}{9},\ \ \gamma_{11}=\frac{8}{3},\
\gamma_{20}=\frac{32}{9},\ \gamma_{21}=\frac{40}{9},\
\gamma_{22}=\frac{14}{3}.$ (A.16)

##  Appendix B Analytical expressions for the coefficients $A_{k},\ B_{k}$ and
$C_{k}$ defined in Eqs.(55)-(57)

Here we provide the explicit expressions for the coefficients
$\left[X_{k}\right]_{VV,AV,TT}$ defined in Eq.(68). These coefficients are
functions of the momentum fractions

$\displaystyle\left[X_{k}\right]_{VV,AV,TT}\equiv\left[X_{k}\right]_{VV,AV,TT}(x_{123},y_{123}).$
(B.1)

Below we also use following short notation

$\displaystyle\bar{m}=\frac{m_{c}}{E}=\frac{1}{\sqrt{1+\boldsymbol{v}^{2}}},\
\ C_{ij}=x_{i}y_{j}+x_{j}y_{i}.$ (B.2)

The coefficients $C_{k}$ can be obtained using the simple relations

$\displaystyle\left[C_{k}\right]_{VV,TT}(x_{123},y_{123})$ $\displaystyle=$
$\displaystyle\left[A_{k}\right]_{VV,TT}(x_{213},y_{213}),$ (B.3)
$\displaystyle\left[C_{k}\right]_{AV}(x_{123},y_{123})$ $\displaystyle=$
$\displaystyle-\left[A_{k}\right]_{AV}(x_{213},y_{213}).$ (B.4)

The other coefficients read

$\displaystyle\left[A_{0}\right]_{VV}$ $\displaystyle=$ $\displaystyle
2C_{13}-(1-\bar{m})(x_{3}+y_{3})-2\bar{m}(1-\bar{m}),$ (B.5)
$\displaystyle\left[A_{0}\right]_{AV}$ $\displaystyle=$
$\displaystyle\frac{2}{1+\bar{m}}(-2C_{13}+(1-\bar{m}^{2})(x_{1}-x_{2}+y_{1}-y_{2})),$
(B.6) $\displaystyle\left[A_{0}\right]_{TT}$ $\displaystyle=$
$\displaystyle-4(1-\bar{m})(\bar{m}+1-x_{3}-y_{3}).$ (B.7)
$\displaystyle\left[B_{0}\right]_{VV}$ $\displaystyle=$
$\displaystyle-(1-\bar{m})(2\bar{m}+x_{3}+y_{3}),$ (B.8)
$\displaystyle\left[B_{0}\right]_{AV}$ $\displaystyle=$ $\displaystyle
2(1-\bar{m})(x_{1}-x_{2}+y_{1}-y_{2}),$ (B.9)
$\displaystyle\left[B_{0}\right]_{TT}$ $\displaystyle=$ $\displaystyle
4(2C_{12}+(1-\bar{m})(x_{3}+y_{3})+\bar{m}^{2}-1).$ (B.10)
$\displaystyle\left[A_{1}\right]_{VV}$ $\displaystyle=$
$\displaystyle\frac{\bar{m}}{1+\bar{m}}(4(x_{1}y_{3}-x_{3}y_{1})+(1-\bar{m})(x_{3}-x_{1}+y_{1}-y_{3})),$
(B.11) $\displaystyle\left[A_{1}\right]_{AV}$ $\displaystyle=$
$\displaystyle-\left[A_{1}\right]_{VV},$ (B.12)
$\displaystyle\left[A_{1}\right]_{TT}$ $\displaystyle=$ $\displaystyle 0.$
(B.13) $\displaystyle\left[B_{1}\right]_{VV}$ $\displaystyle=$ $\displaystyle
0,$ (B.14) $\displaystyle\left[B_{1}\right]_{AV}$ $\displaystyle=$
$\displaystyle 0,$ (B.15) $\displaystyle\left[B_{1}\right]_{TT}$
$\displaystyle=$
$\displaystyle\frac{4\bar{m}}{1+\bar{m}}(4(x_{1}y_{2}-x_{2}y_{1})+(1-\bar{m})(x_{2}-x_{1}+y_{1}-y_{2})).$
(B.16) $\displaystyle\left[A_{2}\right]_{VV}$ $\displaystyle=$
$\displaystyle\frac{\bar{m}^{2}}{(1+\bar{m})^{2}}(2C_{13}+(1+\bar{m})(2+x_{3}+y_{3})-3(1-\bar{m}^{2})),$
(B.17) $\displaystyle\left[A_{2}\right]_{AV}$ $\displaystyle=$
$\displaystyle-\frac{\bar{m}^{2}}{1+\bar{m}}(4(x_{1}+y_{1}-1)+2(x_{3}+y_{3})+\bar{m}-1),$
(B.18) $\displaystyle\left[A_{2}\right]_{TT}$ $\displaystyle=$
$\displaystyle\frac{4\bar{m}^{2}}{1+\bar{m}}(1-x_{3}-y_{3}+\bar{m}).$ (B.19)
$\displaystyle\left[B_{2}\right]_{VV}$ $\displaystyle=$
$\displaystyle\frac{\bar{m}^{2}}{1+\bar{m}}(x_{3}+y_{3}+2\bar{m}),$ (B.20)
$\displaystyle\left[B_{2}\right]_{AV}$ $\displaystyle=$
$\displaystyle\frac{2\bar{m}^{2}}{1+\bar{m}}(x_{2}-x_{1}+y_{2}-y_{1}),$ (B.21)
$\displaystyle\left[B_{2}\right]_{TT}$ $\displaystyle=$
$\displaystyle\frac{4\bar{m}^{2}}{(1+\bar{m})^{2}}(2C_{12}+(1-\bar{m})(2-x_{3}-y_{3})-2(1-\bar{m}^{2})).$
(B.22) $\displaystyle\left[A_{3}\right]_{VV}$ $\displaystyle=$
$\displaystyle\frac{\bar{m}^{3}}{(1+\bar{m})^{2}}(2(x_{1}-y_{1})+y_{3}-x_{3}),$
(B.23) $\displaystyle\left[A_{3}\right]_{AV}$ $\displaystyle=$ $\displaystyle
0,$ (B.24) $\displaystyle\left[A_{3}\right]_{TT}$ $\displaystyle=$
$\displaystyle-\frac{4\bar{m}^{3}}{(1+\bar{m})^{2}}(x_{3}-y_{3}).$ (B.25)
$\displaystyle\left[B_{3}\right]_{VV}$ $\displaystyle=$
$\displaystyle\frac{\bar{m}^{3}}{(1+\bar{m})^{2}}(x_{1}-x_{2}-y_{1}+y_{2}),$
(B.26) $\displaystyle\left[B_{3}\right]_{AV}$ $\displaystyle=$ $\displaystyle
0,$ (B.27) $\displaystyle\left[B_{3}\right]_{TT}$ $\displaystyle=$
$\displaystyle\frac{4\bar{m}^{3}}{(1+\bar{m})^{2}}(x_{1}-x_{2}-y_{1}+y_{2}).$
(B.28) $\displaystyle\left[A_{4}\right]_{VV}$ $\displaystyle=$
$\displaystyle\left[B_{4}\right]_{VV}=\frac{2\bar{m}^{4}}{(1+\bar{m})^{2}},$
(B.29) $\displaystyle\left[A_{4}\right]_{AV}$ $\displaystyle=$
$\displaystyle\left[B_{4}\right]_{AV}=0,$ (B.30)
$\displaystyle\left[A_{4}\right]_{TT}$ $\displaystyle=$
$\displaystyle\left[B_{4}\right]_{TT}=\frac{4\bar{m}^{4}}{(1+\bar{m})^{2}}.$
(B.31)

## References

  * [1] M. Ablikim et al. [BES], Phys. Lett. B 648 (2007), 149-155 doi:10.1016/j.physletb.2007.02.029 [arXiv:hep-ex/0610079 [hep-ex]].
  * [2] M. Ablikim et al. [BES], Phys. Rev. D 78 (2008), 092005 doi:10.1103/PhysRevD.78.092005 [arXiv:0810.1896 [hep-ex]].
  * [3] M. Ablikim et al. [BESIII], Phys. Rev. D 86 (2012), 032014 doi:10.1103/PhysRevD.86.032014 [arXiv:1205.1036 [hep-ex]].
  * [4] M. Ablikim et al. [BESIII], Phys. Rev. D 93 (2016) no.7, 072003 doi:10.1103/PhysRevD.93.072003 [arXiv:1602.06754 [hep-ex]].
  * [5] M. Ablikim et al. [BESIII], Phys. Lett. B 770 (2017), 217-225 doi:10.1016/j.physletb.2017.04.048 [arXiv:1612.08664 [hep-ex]].
  * [6] M. Ablikim et al. [BESIII], Phys. Rev. D 95 (2017) no.5, 052003 doi:10.1103/PhysRevD.95.052003 [arXiv:1701.07191 [hep-ex]].
  * [7] M. Ablikim et al. [BESIII], Phys. Rev. D 98 (2018) no.3, 032006 doi:10.1103/PhysRevD.98.032006 [arXiv:1803.02039 [hep-ex]].
  * [8] M. Ablikim et al. [BESIII], Phys. Rev. Lett. 125 (2020) no.5, 052004 doi:10.1103/PhysRevLett.125.052004 [arXiv:2004.07701 [hep-ex]].
  * [9] M. Alekseev, A. Amoroso, R. B. Ferroli, I. Balossino, M. Bertani, D. Bettoni, F. Bianchi, J. Chai, G. Cibinetto and F. Cossio, et al. Chin. Phys. C 43 (2019) no.2, 023103 doi:10.1088/1674-1137/43/2/023103 [arXiv:1809.04273 [hep-ph]].
  * [10] R. Baldini Ferroli, A. Mangoni, S. Pacetti and K. Zhu, Phys. Lett. B 799 (2019), 135041 doi:10.1016/j.physletb.2019.135041 [arXiv:1905.01069 [hep-ph]]; R. B. Ferroli, A. Mangoni, S. Pacetti and K. Zhu, Phys. Rev. D 103 (2021) no.1, 016005 doi:10.1103/PhysRevD.103.016005 [arXiv:2005.11265 [hep-ph]]; A. Mangoni, [arXiv:2202.08542 [hep-ph]].
  * [11] S. J. Brodsky and G. P. Lepage, Phys. Rev. D 24 (1981) 2848. doi:10.1103/PhysRevD.24.2848
  * [12] V. L. Chernyak and A. R. Zhitnitsky, Phys. Rept. 112 (1984) 173. doi:10.1016/0370-1573(84)90126-1
  * [13] N. Brambilla et al. [Quarkonium Working Group], hep-ph/0412158.
  * [14] V. L. Chernyak, A. A. Ogloblin and I. R. Zhitnitsky, Z. Phys. C 42 (1989) 583 [Yad. Fiz. 48 (1988) 1398] [Sov. J. Nucl. Phys. 48 (1988) 889]. doi:10.1007/BF01557664
  * [15] G. T. Bodwin, E. Braaten and G. P. Lepage, Phys. Rev. D 51 (1995) 1125 [Phys. Rev. D 55 (1997) 5853] [hep-ph/9407339].
  * [16] N. Brambilla, A. Pineda, J. Soto and A. Vairo, Rev. Mod. Phys. 77 (2005) 1423 [hep-ph/0410047].
  * [17] N. Brambilla et al. [Quarkonium Working Group], doi:10.5170/CERN-2005-005 [arXiv:hep-ph/0412158 [hep-ph]].
  * [18] M. B. Voloshin, Prog. Part. Nucl. Phys. 61 (2008), 455-511 doi:10.1016/j.ppnp.2008.02.001 [arXiv:0711.4556 [hep-ph]].
  * [19] N. Kivel, Eur. Phys. J. A 56 (2020) no.2, 64 [erratum: Eur. Phys. J. A 57 (2021) no.9, 271] doi:10.1140/epja/s10050-021-00575-9 [arXiv:1910.02850 [hep-ph]].
  * [20] N. Kivel, [arXiv:2109.05847 [hep-ph]].
  * [21] N. Kivel, Eur. Phys. J. A 58 (2022) no.7, 138 doi:10.1140/epja/s10050-022-00762-2 [arXiv:2203.00377 [hep-ph]].
  * [22] P.A. Zyla et al. [Particle Data Group], PTEP 2020 (2020) no.8, 083C01 doi:10.1093/ptep/ptaa104
  * [23] G. T. Bodwin and A. Petrelli, Phys. Rev. D 66 (2002), 094011 [erratum: Phys. Rev. D 87 (2013) no.3, 039902] doi:10.1103/PhysRevD.66.094011 [arXiv:hep-ph/0205210 [hep-ph]].
  * [24] G. T. Bodwin, H. S. Chung, D. Kang, J. Lee and C. Yu, Phys. Rev. D 77 (2008), 094017 doi:10.1103/PhysRevD.77.094017 [arXiv:0710.0994 [hep-ph]].
  * [25] G. T. Bodwin, H. S. Chung, J. H. Ee, J. Lee and F. Petriello, Phys. Rev. D 90 (2014) no.11, 113010 doi:10.1103/PhysRevD.90.113010 [arXiv:1407.6695 [hep-ph]].
  * [26] N. Brambilla, H. S. Chung, W. K. Lai, V. Shtabovenko and A. Vairo, Phys. Rev. D 100 (2019) no.5, 054038 doi:10.1103/PhysRevD.100.054038 [arXiv:1907.06473 [hep-ph]].
  * [27] J. P. Ma and Q. Wang, Phys. Lett. B 537 (2002), 233-240 doi:10.1016/S0370-2693(02)01937-8 [arXiv:hep-ph/0203082 [hep-ph]].
  * [28] N. Brambilla, E. Mereghetti and A. Vairo, JHEP 08 (2006), 039 [erratum: JHEP 04 (2011), 058] doi:10.1088/1126-6708/2006/08/039 [arXiv:hep-ph/0604190 [hep-ph]].
  * [29] E. Braaten and J. Lee, Phys. Rev. D 67 (2003), 054007 [erratum: Phys. Rev. D 72 (2005), 099901] doi:10.1103/PhysRevD.72.099901 [arXiv:hep-ph/0211085 [hep-ph]].
  * [30] G. T. Bodwin, J. Lee and C. Yu, Phys. Rev. D 77 (2008), 094018 doi:10.1103/PhysRevD.77.094018 [arXiv:0710.0995 [hep-ph]].
  * [31] T. Barnes, X. Li and W. Roberts, Phys. Rev. D 77 (2008), 056001 doi:10.1103/PhysRevD.77.056001 [arXiv:0709.4491 [hep-ph]].
  * [32] V. Braun, R. J. Fries, N. Mahnke and E. Stein, Nucl. Phys. B 589 (2000) 381 Erratum: [Nucl. Phys. B 607 (2001) 433] doi:10.1016/S0550-3213(00)00516-2, 10.1016/S0550-3213(01)00254-1 [hep-ph/0007279].
  * [33] M. Gremm and A. Kapustin, Phys. Lett. B 407 (1997), 323-330 doi:10.1016/S0370-2693(97)00744-2 [arXiv:hep-ph/9701353 [hep-ph]].
  * [34] G. T. Bodwin, D. Kang and J. Lee, Phys. Rev. D 74 (2006), 014014 doi:10.1103/PhysRevD.74.014014 [arXiv:hep-ph/0603186 [hep-ph]].
  * [35] J. H. Kuhn, J. Kaplan and E. G. O. Safiani, Nucl. Phys. B 157 (1979), 125-144 doi:10.1016/0550-3213(79)90055-5
  * [36] J. Kublbeck, H. Eck and R. Mertig, Nucl. Phys. B Proc. Suppl. 29 (1992), 204-208; V. Shtabovenko, R. Mertig and F. Orellana, Comput. Phys. Commun. 256 (2020), 107478 doi:10.1016/j.cpc.2020.107478 [arXiv:2001.04407 [hep-ph]].
  * [37] I. V. Anikin, V. M. Braun and N. Offen, Phys. Rev. D 88 (2013) 114021 doi:10.1103/PhysRevD.88.114021 [arXiv:1310.1375 [hep-ph]].
  * [38] I. I. Balitsky, V. M. Braun and A. V. Kolesnichenko, Sov. J. Nucl. Phys. 44 (1986), 1028; I. I. Balitsky, V. M. Braun and A. V. Kolesnichenko, Nucl. Phys. B 312 (1989), 509-550; V. L. Chernyak and I. R. Zhitnitsky, Nucl. Phys. B 345 (1990), 137-172 doi:10.1016/0550-3213(90)90612-H
  * [39] V. M. Braun, S. E. Derkachov, G. P. Korchemsky and A. N. Manashov, Nucl. Phys. B 553 (1999), 355-426 doi:10.1016/S0550-3213(99)00265-5 [arXiv:hep-ph/9902375 [hep-ph]]; V. M. Braun, A. N. Manashov and J. Rohrwild, Nucl. Phys. B 807 (2009) 89 doi:10.1016/j.nuclphysb.2008.08.012 [arXiv:0806.2531 [hep-ph]].
  * [40] Y. J. Zhang, Y. j. Gao and K. T. Chao, Phys. Rev. Lett. 96 (2006), 092001 doi:10.1103/PhysRevLett.96.092001 [arXiv:hep-ph/0506076 [hep-ph]]; B. Gong and J. X. Wang, Phys. Rev. D 77 (2008), 054028 doi:10.1103/PhysRevD.77.054028 [arXiv:0712.4220 [hep-ph]].
  * [41] V. L. Chernyak and I. R. Zhitnitsky, Nucl. Phys. B 246 (1984), 52-74 doi:10.1016/0550-3213(84)90114-7