# Near Threshold Heavy Quarkonium Photoproduction at Large Momentum Transfer

Peng Sun  Department of Physics and Institute of Theoretical Physics, Nanjing
Normal University, Nanjing, Jiangsu 210023, China  Xuan-Bo Tong  School of
Science and Engineering, The Chinese University of Hong Kong, Shenzhen,
Shenzhen, Guangdong, 518172, P.R. China  University of Science and Technology
of China, Hefei, Anhui, 230026, P.R.China  Feng Yuan  Nuclear Science
Division, Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA

###### Abstract

Perturbative QCD is applied to investigate the near threshold heavy quarkonium
photoproduction at large momentum transfer. We take into account the
contributions from the leading three-quark Fock states of the nucleon. The
dominant contribution comes from the three-quark Fock state with one unit
quark orbital angular momentum (OAM) whereas that from zero quark OAM is
suppressed at the threshold. From our analysis, we also show that there is no
direct connection between the near threshold heavy quarkonium photoproduction
and the gluonic gravitational form factors of the nucleon. Based on the
comparison between our result and recent GlueX data of $J/\psi$
photoproduction, we make predictions for $\psi^{\prime}$ and $\Upsilon$
(1S,2S) states which can be tested in future experiments.

##  I Introduction

Exclusive heavy quarkonium production in high energy photon-proton scattering,

$\gamma^{(*)}+N\to V+N^{\prime}\ ,$ (1)

where the incoming photon can be real or virtual, has attracted great
attention in hadron physics community. This process is dominated by the two
gluon exchange Ryskin:1992ui ; Brodsky:1994kf  and can be formulated in the
generalized parton distribution (GPD) Ji:1996ek ; Ji:1996nm  framework
Collins:1996fb ; Hoodbhoy:1996zg ; Koempel:2011rc ; Cui:2018jha ;
Ivanov:2004vd . The theory advance has also pushed the perturbative QCD
computation of these processes to the next-to-leading order Ivanov:2004vd ;
Chen:2019uit ; Flett:2021ghh .

Recently, there has been a strong interest of this process at the lower end of
the energy range near the threshold Kharzeev:1995ij ; Kharzeev:1998bz ;
Gryniuk:2016mpk ; Hatta:2018ina ; Ali:2019lzf ; Hatta:2019lxo ;
Boussarie:2020vmu ; Mamo:2019mka ; Gryniuk:2020mlh ; Wang:2019mza ;
Zeng:2020coc ; Du:2020bqj ; Kharzeev:2021qkd ; Wang:2021dis ; Hatta:2021can ;
Mamo:2021krl ; Kou:2021bez ; Guo:2021ibg ; Mamo:2021tzd . In particular, it
was argued in Refs. Kharzeev:1995ij ; Kharzeev:1998bz  that this process can
provide a direct access to the so-called trace anomaly contribution to the
proton mass, while the origin of the proton mass is of fundamental in QCD
strong interaction theory Shifman:1978zn ; Ji:1994av ; Ji:1995sv ;
Hatta:2018sqd ; Metz:2020vxd ; Hatta:2020iin ; Ji:2021pys ; Ji:2021mtz ;
Lorce:2021xku .

In experiments, $J/\psi$ photo-production from the nuclear targets near the
threshold have been investigated before Gittelman:1975ix ; Camerini:1975cy .
More recently, high precision measurements have been carried out by the GlueX
collaboration at Jefferson Lab Ali:2019lzf . Future experiments will explore
both $J/\psi$ and $\Upsilon$ near threshold photo-production in great details
Joosten:2018gyo , including JLab-12GeV Dudek:2012vr ; Chen:2014psa  and
electron-ion colliders (EIC) Accardi:2012qut ; AbdulKhalek:2021gbh ;
Anderle:2021wcy .

In this paper, we will focus on one of the key aspects of the threshold
kinematics that the momentum transfer is relatively large: $-t\sim 2{\rm
GeV}^{2}$ and $10{\rm GeV}^{2}$ for $J/\psi$ and $\Upsilon$, respectively.
Here, $t$ is the momentum transfer squared from the nucleon target. Because of
this large momentum transfer, we can apply the QCD factorization argument to
compute the scattering amplitude. This factorization follows that of the
hadron form factor calculations in perturbative QCD Lepage:1979za ;
Brodsky:1981kj ; Efremov:1979qk ; Chernyak:1977as ; Chernyak:1980dj ;
Chernyak:1983ej ; Belitsky:2002kj ; Tong:2021ctu . For the heavy quarkonium
production in the final state, the non-relativistic QCD (NRQCD) Bodwin:1994jh
will be adopted and the associated color-singlet matrix element of the
quarkonium state is responsible for its production in the exclusive process.

In the perturbative calculations, the quark/gluon propagators in the
scattering amplitudes lead to the power behavior for the differential cross
section at large momentum transfer Brodsky:1973kr ; Matveev:1973ra ; Ji:2003yj
, which has been commonly assumed in the phenomenological studies of near
threshold heavy quarkonium production, see, e.g., Refs. Frankfurt:2002ka ;
Ali:2019lzf ; Kharzeev:2021qkd ; Wang:2021dis . In the following, we will
provide an explicit calculation to demonstrate this power behavior.

The hard exclusive processes at large momentum transfer depend on the non-
perturbative distribution amplitudes Lepage:1979za . In our derivations, we
take into account the contributions from the leading-twist and higher-twist
terms of the nucleon distribution amplitudes Braun:1999te ; Braun:2000kw .
They correspond to the three-quark Fock state light-cone wave functions of the
nucleon with zero orbital angular momentum (OAM) and one unit OAM components
Ji:2002xn , respectively. Their contributions lead to different power
behaviors at large $(-t)$, similar to the nucleon’s form factors
Belitsky:2002kj ; Tong:2021ctu .

We will also take the heavy quark mass limit and apply the following hierarchy
in scales:

$W_{\gamma p}^{2}\sim M_{V}^{2}\gg(-t)\gg\Lambda_{QCD}^{2}\ ,$ (2)

where $\Lambda_{QCD}$ for the non-perturbative scale. In addition, throughout
the following analysis, we take the threshold limit, i.e., $W_{\gamma p}\sim
M_{V}+M_{p}$, where $W_{\gamma p}$ represents the center of mass energy and
$M_{V}$ and $M_{p}$ for the heavy quarkonium and proton masses, respectively.
To determine the leading contribution, we introduce a parameter Brodsky:2000zc
,

$\chi=\frac{M_{V}^{2}+2M_{p}M_{V}}{W_{\gamma p}^{2}-M_{p}^{2}}\ ,$ (3)

which goes to 1 at the threshold. We will expand the amplitude in terms of
$(1-\chi)$. By applying this expansion, we find that the commonly used
$1/(-t)^{4}$ power term for the differential cross section is suppressed by
$(1-\chi)$. The dominant contribution at the threshold actually comes from the
higher-twist term with $1/(-t)^{5}$ power behavior.

As mentioned above, the exclusive heavy quarkonium production has been
extensively studied in the GPD framework and the scattering amplitude can be
written in terms of the gluon GPDs. In Refs. Boussarie:2020vmu ; Hatta:2021can
; Guo:2021ibg , the GPD formalism was applied in the the threshold kinematics,
where the connection to the gluonic gravitational form factors was explored.
One of the major objectives of this paper is to check the connection between
the near threshold heavy quarkonium photo-production and the gluonic
gravitational form factors. To do that, we compare the differential cross
section derived in this paper and those of the gluonic gravitational form
factors of the nucleon at large momentum transfer in Ref. Tong:2021ctu . We
will show explicitly that there is no direct connection between them.
Therefore, approximations have to be made to link the GPD formalism of this
process to the gluonic gravitational form factors Boussarie:2020vmu ;
Hatta:2021can ; Guo:2021ibg . A brief summary of our results has been
published in Ref. Sun:2021gmi . In the following, we provide more detailed
derivations.

The rest of the paper is organized as follows. In Sec. II, we will examine the
threshold kinematics and apply the expansion method mentioned above to
simplify the derivation. In Sec. III, we take the example of photon scattering
off a pion target. The leading Fock state of the pion contains quark and
antiquark and the derivation is much simpler compared to the nucleon case.
Sec. IV and V will be dedicated to the nucleon case. In Sec. IV, we study the
contribution from the leading component of the nucleon distribution amplitude
and show that its contribution is actually suppressed in the threshold limit.
In Sec. V, we perform the analysis of higher-twist component of the nucleon
distribution amplitude and show that its contribution to the differential
cross section does not vanish at the threshold. In Sec. VI, we discuss the
interpretation and consequence of our derivations. We conclude that there is
no direct connection between the near threshold photo-production of heavy
quarkonium and the gluonic gravitational form factors of the nucleon. In Sec.
VII, we provide phenomenological applications of our derivations. Predictions
on $\psi^{\prime}$ and $\Upsilon$ will be presented for future experiments
based on the comparison between our results and the GlueX data on near
threshold $J/\psi$ production at JLab. Finally, we summarize our paper in Sec.
VIII.

##  II Near Threshold Kinematics

The typical Feynman diagram of the two-gluon exchange contributions to the
near threshold heavy quarkonium photoproduction is shown in Fig. 1,

$\gamma(k_{\gamma})+N(p_{1})\to J/\psi(k_{\psi})+N^{\prime}(p_{2})\ ,$ (4)

where we have used $J/\psi$ as an example. In order to make the near threshold
expansion more evident, it is useful to examine the relevant kinematics for
the scattering amplitude. The center of mass energy and momentum transfer
squared can be written as,

$\displaystyle W_{\gamma
p}^{2}=(k_{\gamma}+p_{1})^{2}=(k_{\psi}+p_{2})^{2}\sim M_{V}^{2}\ ,$ (5)
$\displaystyle|t|=|(p_{2}-p_{1})^{2}|\ll M_{V}^{2}\ .$ (6)

Therefore, we will have the following approximations around the threshold
kinematics,

$\displaystyle p_{1}\cdot k_{\gamma}\sim p_{1}\cdot k_{\psi}\sim M_{V}^{2}\ ,$
(7) $\displaystyle p_{2}\cdot k_{\gamma}\sim p_{2}\cdot k_{\psi}\ll M_{V}^{2}\
.$ (8)

In addition, applying the heavy quark mass limit of $M_{V}^{2}\gg(-t)$, we
find that the invariant mass of the $t$-channel two gluons is much smaller
than heavy quarkonium mass. We will also take the approximation of
$M_{c}\approx M_{V}/2$ in the non-relativistic limit of the heavy quarkonium
system.

Figure 1: Schematics of two-gluon exchange contribution to the threshold heavy
quarkonium production.

The quark propagators in the upper part of the Feynman diagram of Fig. 1 are
all in order of $1/M_{V}$. For example, one of the quark propagators can be
simplified as

$\displaystyle\frac{1}{\left(k_{1}-k_{\psi}/2\right)^{2}-M_{c}^{2}}$
$\displaystyle=$ $\displaystyle\frac{1}{-k_{1}\cdot k_{\gamma}-k_{1}\cdot
k_{2}}$ (9) $\displaystyle\approx$ $\displaystyle\frac{1}{-k_{1}\cdot
k_{\gamma}}\ ,$

where we have applied $|k_{1}^{2}|\sim|k_{2}^{2}|\sim|k_{1}\cdot
k_{2}|\sim|t|\ll M_{V}^{2}$. Because $k_{1}$ carries certain momentum fraction
of the incoming nucleon, $k_{1}\cdot k_{\gamma}$ will be order of $M_{V}^{2}$.
Similarly, we have

$\displaystyle\frac{1}{\left(k_{2}-k_{\psi}/2\right)^{2}-M_{c}^{2}}$
$\displaystyle\approx$ $\displaystyle\frac{1}{-k_{2}\cdot k_{\gamma}}\ .$ (10)

The following propagator will also show up in some of the Feynman diagrams,

$\displaystyle\frac{1}{\left(k-k_{\psi}/2\right)^{2}-M_{c}^{2}}=\frac{1}{-k\cdot
k_{\gamma}}\approx\frac{2}{-M_{V}^{2}}\ ,$ (11)

where $k=k_{1}+k_{2}=p_{1}-p_{2}$. In the center of mass frame, $k$ is
dominated by $p_{1}$ because $p_{2}$ is soft.

Applying the above approximations, we can simplify the photon-heavy quarkonium
transition amplitude. Let us define $\mu$ and $\nu$ for the polarization
indices for $k_{1}$ and $k_{2}$, respectively, and $\epsilon_{\gamma}$ and
$\epsilon_{\psi}$ for the photon polarization and $J/\psi$ polarization
vectors, respectively. To further simplify the derivation, we choose the
physical polarization for the incoming photon,

$\epsilon_{\gamma}\cdot k_{\gamma}=0,~{}~{}\epsilon_{\gamma}\cdot p_{1}=0\ .$
(12)

With this choice, we notice that the contributions from
$\epsilon_{\gamma}\cdot k_{1}$ and $\epsilon_{\gamma}\cdot k_{2}$ are also
suppressed in the heavy quark mass limit. Therefore, we will drop these terms
as well. We emphasize, all these approximations have been cross checked by a
full computation.

Finally, we have the following expression for the amplitude from the heavy
quarkonium side,

${\cal
M}_{\psi,ab}^{\mu\nu}=\frac{\delta^{ab}N_{\psi}\left[\epsilon_{\psi}^{*}\cdot\epsilon_{\gamma}{\cal
W}_{T}^{\mu\nu}+\epsilon_{\psi}^{*}\cdot k{\cal W}_{L}^{\mu\nu}+{\cal
W}_{S}^{\mu\nu}\right]}{k_{1}\cdot k_{\gamma}k_{2}\cdot k_{\gamma}}\ ,$ (13)

where $a$ and $b$ represent the color indices for the $t$-channel gluons. In
the above equation, $N_{\psi}$ is defined as

$N_{\psi}=-\frac{4e_{c}eg_{s}^{2}}{\sqrt{N_{c}M_{V}^{3}}}\psi_{J}(0)\ ,$ (14)

where $\psi_{J}(0)$ is the wave function of $J/\psi$ at the origin and is
related to the NRQCD matrix element Bodwin:1994jh . The tensor structures
${\cal W}_{T,L,S}^{\mu\nu}$ are defined as

$\displaystyle{\cal W}_{T}^{\mu\nu}$ $\displaystyle=$ $\displaystyle-
k_{1}\cdot k_{\gamma}k_{2}\cdot k_{\gamma}g^{\mu\nu}-k_{1}\cdot
k_{2}k_{\gamma}^{\mu}k_{\gamma}^{\nu}$ $\displaystyle+k_{1}\cdot
k_{\gamma}k_{2}^{\mu}k_{\gamma}^{\nu}+k_{2}\cdot
k_{\gamma}k_{1}^{\nu}k_{\gamma}^{\mu}$ $\displaystyle{\cal W}_{L}^{\mu\nu}$
$\displaystyle=$ $\displaystyle k_{1}\cdot
k_{\gamma}\epsilon_{\gamma}^{\nu}k_{2}^{\mu}+k_{2}\cdot
k_{\gamma}\epsilon_{\gamma}^{\mu}k_{1}^{\nu}$ $\displaystyle{\cal
W}_{S}^{\mu\nu}$ $\displaystyle=$ $\displaystyle-k_{1}\cdot
k_{2}\left(k_{1}\cdot
k_{\gamma}\epsilon_{\psi}^{*\mu}\epsilon_{\gamma}^{\nu}+k_{2}\cdot
k_{\gamma}\epsilon_{\psi}^{*\nu}\epsilon_{\gamma}^{\mu}\right.$ (15)
$\displaystyle~{}~{}\left.+k_{1}\cdot\epsilon_{\psi}^{*}k_{\gamma}^{\nu}\epsilon_{\gamma}^{\mu}+k_{2}\cdot\epsilon_{\psi}^{*}k_{\gamma}^{\mu}\epsilon_{\gamma}^{\nu}\right)\
.$

where ${\cal W}_{T}$ and ${\cal W}_{L}$ represent the amplitudes for a
transversely polarized and longitudinal polarized heavy quarkonium in the
final state, respectively, whereas ${\cal W}_{S}$ for a subleading term.

Clearly the above amplitude is symmetric under $k_{1},\mu\leftrightarrow
k_{2},\nu$. In the above equation, the first term is the leading contribution
in the heavy quark mass limit at the threshold. The second and third terms are
subleading contributions.

We have also carried out an important cross check for the above results. We
compute the full amplitude without any approximation. We then take the leading
contribution of the differential cross section (the amplitude squared) in the
heavy quark mass limit and threshold limit, and obtain the same result.

###  II.1 Vanishing of Three-gluon Exchange Contribution

Figure 2: Typical Feynman diagram from three-gluon exchange. These diagrams
vanish because of the $C$-parity conservation.

Before we start our derivations of the threshold scattering amplitudes, we
would like to comment on the three-gluon exchange contributions. The two-gluon
and three-gluon exchange diagrams were considered in Ref. Brodsky:2000zc  for
the threshold production of $J/\psi$ and it was argued that the three-gluon
exchange diagrams dominate the differential cross section contributions.

However, we find that the three-gluon exchange diagrams do not contribute in
our framework, due to the $C$-parity conservation. This is because the three
gluons from the nucleon side carry symmetric color structure (such as
$d_{abc}$) while those from the heavy quarkonium ($J/\psi$) side are
antisymmetric (such as $f_{abc}$), where $a$, $b$ and $c$ represent the color
indices for the three gluons in the $t$-channel, respectively. Explicitly,
from the nucleon side, we have, as shown in Fig. 2,

$\displaystyle\epsilon^{ijk}\epsilon^{lmn}T^{a}_{il}T^{b}_{jm}T^{c}_{kn}\propto
d^{abc}\ ,$ (16)

where $ijk$ and $lmn$ represent the color indices for the initial and final
three quarks, respectively. Here, we have applied the anti-symmetric color
structure for the three-quark Fock state wave function of the nucleon
Ji:2002xn . On the other hand, for the heavy quarkonium side, we have, instead

$\displaystyle{\rm
Tr}\left[T^{a}T^{b}T^{c}\right]=\frac{1}{4}\left(d^{abc}+if^{abc}\right)\ .$
(17)

However, because of $J/\psi$ is in the $1^{--}$ state, the photon-$J/\psi$
transition amplitude vanishes for the symmetric color configuration with three
gluons, i.e., $d^{abc}$ term from the above vanishes. Combining this with the
color structure from the nucleon side, we conclude the three-gluon exchange
diagrams do not contribute.

##  III Pion Case

In this section, we take the example of pion case to show the detailed of our
derivations. In this case, we have photon scatters on the pion target and
produces a $J/\psi$ in the final state close to the threshold,

$\gamma+\pi\to J/\psi+\pi\ ,$ (18)

where the dominant contribution is again a two-gluon exchange diagram. The two
gluons attach to the two quark lines from the pion target, as shown in Fig. 3.

Figure 3: The Feynman diagram contribution to the exclusive $\gamma\pi\to\pi
J/\psi$ at large momentum transfer. The two gluons attach to the quark and
antiquark lines, respectively.

Considering the leading Fock component of the pion, we have

$\displaystyle|\pi^{+}\rangle_{u\overline{d}}$ $\displaystyle=$
$\displaystyle\int
d[1]d[2]\psi_{u\overline{d}}(1,2)\frac{\delta_{ij}}{\sqrt{3}}\left[u^{\dagger}_{\uparrow
i}(1)\overline{d}^{\dagger}_{\downarrow j}(2)\right.$ (19)
$\displaystyle\left.-u^{\dagger}_{\downarrow
i}(1)\overline{d}^{\dagger}_{\uparrow j}(2)\right]|0\rangle\,$

where $i$ and $j=1,2,3$ are the color indices, and $\uparrow$ and $\downarrow$
label quark light-cone helicity $+1/2$ and $-1/2$, respectively. The color
factor $\delta_{ij}/\sqrt{3}$ is normalized to 1. The light-cone wave function
amplitude $\psi_{u\bar{d}}(1,2)$ is a function of quark momenta with argument
1 representing $x_{1}$ and $q_{1\perp}$ and so on. Since the momentum
conservation implies $\vec{q}_{1\perp}+\vec{q}_{2\perp}=0$ and
$x_{1}+x_{2}=1$, $\psi_{u\bar{d}}(1,2)$ depends on variables $x_{1}$ and
$q_{1\perp}$ only. The integration in the above equation is defined as,

$\displaystyle\int
d[1]d[2]=\int\frac{d^{2}q_{1\perp}}{(2\pi)^{3}}\frac{dx_{1}}{2\sqrt{x_{1}(1-x_{1})}}\
.$ (20)

From the light-cone wave function, we obtain the distribution amplitude,

$\displaystyle\phi(x)=\int\frac{d^{2}q_{1\perp}}{(2\pi)^{3}}\psi_{u\bar{d}}(1,2)\
.$ (21)

The final scattering amplitude of $\gamma+\pi^{+}\to J/\psi+\pi^{+}$ can be
computed in terms of the above distribution amplitude of pion,

$\displaystyle{\cal A}^{\pi}=$ $\displaystyle\int
dx_{1}dy_{1}\phi^{*}(y_{1})\phi(x_{1}){\cal
M}_{\psi}^{\mu\nu}(\epsilon_{\gamma},\epsilon_{\psi},x_{1},y_{1})$
$\displaystyle\times\frac{-g_{s}^{2}C_{F}}{2k_{1}^{2}k_{2}^{2}}{\rm
Tr}\left[\not{p}_{2}\gamma^{\mu}\not{p}_{1}\gamma^{\nu}\right]\ ,$ (22)

where ${\cal M}_{\psi}^{\mu\nu}$ has been given in the previous section.

###  III.1 Threshold Expansion

At the threshold, the amplitude squared can be further simplified as

$\displaystyle|\overline{{\cal
A}^{\pi}}|^{2}=G_{\psi}G_{\pi}(t)G_{\pi}^{*}(t)\ ,$ (23)

where the spin sum and average have been applied. Here, $G_{\psi}$ is defined
as

$G_{\psi}=|N_{\psi}|^{2}=\frac{384\pi^{2}e_{c}^{2}\alpha(4\pi\alpha_{s})^{2}}{N_{c}^{2}M_{\psi}^{3}}\langle
0|{\cal O}^{\psi}({}^{3}S_{1}^{(1)})|0\rangle\ ,$ (24)

where $\langle 0|{\cal O}({}^{3}S_{1})|0\rangle$ is the color-singlet NRQCD
matrix element for $J/\psi$. $G_{\pi}(t)$ is defined as

$G_{\pi}(t)=\frac{8\pi\alpha_{s}C_{F}}{t}\int
dx_{1}dy_{1}\phi^{*}(y_{1})\phi(x_{1})\frac{1}{x_{1}\bar{x}_{1}y_{1}\bar{y}_{1}}\
,$ (25)

where $C_{F}=(N_{c}^{2}-1)/2N_{c}$, $\bar{x}_{1}=1-x_{1}$ and
$\bar{y}_{1}=1-y_{1}$. Here, we have neglected high order corrections of
$t_{\psi}=-t/M_{V}^{2}$.

###  III.2 Compared to the Gravitational Form Factors

We now compare the above result to the gluonic gravitational form factors at
large momentum transfer, which have been computed in Ref. Tong:2021ctu . For
convenience, we list the results below. The gluonic gravitation form factors
of the pion are defined as

$\displaystyle\langle
p_{2}|T_{g}^{\mu\nu}|p_{1}\rangle=2\bar{P}^{\mu}\bar{P}^{\nu}A_{g}^{\pi}(t)+\frac{1}{2}(\Delta^{\mu}\Delta^{\nu}-g^{\mu\nu}\Delta^{2})C_{g}^{\pi}(t)$
$\displaystyle+2m^{2}g^{\mu\nu}\overline{C}_{g}^{\pi}(t)\ ,$ (26)

where $T_{g}^{\mu\nu}$ is the gluonic energy-momentum tensor in QCD and $m$
represents the pion mass. Here, $\bar{P}=(p_{1}+p_{2})/2$ is the average
momentum, $\Delta=p_{2}-p_{1}$ is the momentum transfer and hence
$t=\Delta^{2}$. From the results of Ref. Tong:2021ctu , we find

$\displaystyle
A^{\pi}_{g}(t)=C^{\pi}_{g}(t)=\frac{4m^{2}}{t}\overline{C}_{g}^{\pi}(t)$ (27)
$\displaystyle=\frac{4\pi\alpha_{s}C_{F}}{-t}\int
dx_{1}dy_{1}\phi^{*}(y_{1})\phi(x_{1})\left(\frac{1}{x_{1}\bar{x}_{1}}+\frac{1}{y_{1}\bar{y}_{1}}\right)\
.$

From the above results, we find that there is no direct connection between
$G_{\pi}(t)$ of Eq. (25) and any of the gravitational form factors of
$A_{g}^{\pi}(t)$, $C_{g}^{\pi}(t)$ or $\overline{C}_{g}^{\pi}(t)$ (Eq. (27)).
This indicates that we can not directly interpret the near threshold heavy
quarkonium photo-production in terms of the gluonic gravitational form
factors.

###  III.3 Compared to the GPD Formalism

As mentioned in the Introduction, the photo-production of heavy quarkonium has
been derived in the GPD framework. If we extend these derivations to the near
threshold kinematics, we obtain

$\displaystyle{\cal A}^{\pi}$ $\displaystyle=$ $\displaystyle
N_{\psi}\epsilon_{\psi}^{*}\cdot\epsilon_{\gamma}\int_{-1}^{1}dx\frac{H_{g}^{\pi}(x,\xi,t)}{(x+\xi-i\varepsilon)(x-\xi+i\varepsilon)}\
,$ (28)

for the pion target, where $N_{\psi}$ has been given in Eq. (14) and $\xi$ is
the skewness parameter. In the threshold limit we take $\xi=1$. In the above
equation, $H_{g}^{\pi}$ represents the GPD gluon distribution of the pion. The
GPD gluon distribution at large momentum transfer can be calculated in terms
of the distribution amplitudes as that of the quark GPD in Ref.
Hoodbhoy:2003uu , for which we list in Appendix A. If we substitute the result
of $H_{g}^{\pi}(x,\xi,t)$ from there, we will be able to reproduce the
scattering amplitude result from the direct computation in the above
subsection A. This provides a useful cross check for our derivations.

##  IV Nucleon Case: twist-three contributions

Now we turn to the proton cases. We show the typical Feynman diagram in Fig.
4. To compute these diagrams, we follow the factorization argument for the
hard exclusive processes Brodsky:1981kj , where the leading contributions come
from the three quark Fock state of the nucleon. The three-quark states can be
further classified into zero orbital angular momentum (OAM) and nonzero OAM
components Ji:2002xn . We will first examine the contribution from zero OAM
component. This corresponds to the twist-three contribution from the nucleon’s
distribution amplitude.

Figure 4: Typical Feynman diagram contributions to the threshold $J/\psi$
photoproduction at large momentum transfer from two-gluon exchange.

###  IV.1 Three-quark Fock State with Zero OAM

Because there is no quark OAM, the total quark spin equals to the nucleon
spin. The associated light-cone wave function amplitude is defined as

$\displaystyle|P\uparrow\rangle_{1/2}=\int
d[1]d[2]d[3]\left(\tilde{\psi}^{(1)}(1,2,3)\right)$
$\displaystyle\times\frac{\epsilon^{ijk}}{\sqrt{6}}u^{\dagger}_{i\uparrow}(1)\left(u^{\dagger}_{j\downarrow}(2)d^{\dagger}_{k\uparrow}(3)-d^{\dagger}_{j\downarrow}(2)u^{\dagger}_{k\uparrow}(3)\right)|0\rangle\
,$ (29)

where $ijk$ represent the color indices for the three quarks, respectively,
and the measure for the quark momentum is,

$\displaystyle
d[1]d[2]d[3]=\sqrt{2}\frac{dx_{1}dx_{2}dx_{3}}{\sqrt{2x_{1}2x_{2}2x_{3}}}\frac{d^{2}\vec{q}_{1\perp}d^{2}\vec{q}_{2\perp}d^{2}\vec{q}_{3\perp}}{(2\pi)^{9}}$
$\displaystyle\times(2\pi)^{3}\delta(1-x_{1}-x_{2}-x_{3})\delta^{(2)}(\vec{q}_{1\perp}+\vec{q}_{2\perp}+\vec{q}_{3\perp})\
.$ (30)

By integrating over the transverse momenta $q_{i\perp}$, we obtain the twist-
three distribution amplitude Braun:1999te

$\displaystyle\Phi_{3}(x_{1},x_{2},x_{3})$ $\displaystyle=$
$\displaystyle-2\sqrt{6}\int[dq_{\perp}]\tilde{\psi}^{(1)}(1,2,3)\ ,$ (31)

where
$[dq_{\perp}]=\frac{d^{2}\vec{q}_{1\perp}d^{2}\vec{q}_{2\perp}d^{2}\vec{q}_{3\perp}}{(2\pi)^{9}}\delta^{(2)}(\vec{q}_{1\perp}+\vec{q}_{2\perp}+\vec{q}_{3\perp})$.
In this configuration, the three quarks only carry longitudinal momenta to
form the nucleon state. The above parameterization applies to both initial and
final state nucleons. Of course, their momenta are different. In addition,
because the quark helicities are conserved, the nucleon helicity is also
conserved.

###  IV.2 Partonic Scattering Amplitude

Schematically, we can write the helicity-conserved amplitude as

$\displaystyle{\cal A}_{3}$ $\displaystyle=\langle
J/\psi(\epsilon_{\psi}),N^{\prime}_{\uparrow}|\gamma(\epsilon_{\gamma}),N_{\uparrow}\rangle$
$\displaystyle=\int[dx][dy]\Phi(x_{1},x_{2},x_{3})\Phi^{*}(y_{1},y_{2},y_{3})$
$\displaystyle~{}~{}\times{\cal
M}_{\psi}^{\mu\nu}(\epsilon_{\gamma},\epsilon_{\psi})\frac{1}{(-t)^{2}}{\cal
H}^{\mu\nu}(\\{x\\},\\{y\\})\ ,$ (32)

where $\\{x\\}=(x_{1},x_{2},x_{3})$ represent the momentum fractions carried
by the three quarks, $[dx]=dx_{1}dx_{2}dx_{3}\delta(1-x_{1}-x_{2}-x_{3})$, and
$\Phi_{3}(x_{1},x_{2},x_{3})$ is the twist-three distribution amplitude of the
proton Lepage:1980fj ; Braun:1999te . The partonic amplitude ${\cal
H}^{\mu\nu}$ is calculated from the lower part of Fig. 4, where the incoming
three quarks carry momenta of $x_{1}p_{1}$, $x_{2}p_{1}$ and $x_{3}p_{1}$ and
outgoing quarks with momenta of $y_{1}p_{2}$, $y_{2}p_{2}$ and $y_{3}p_{2}$,
respectively.

There are total of 12 diagrams (lower part) for the ${\cal H}^{\mu\nu}$.
However, all the diagrams can be generated by only two specific diagrams with
different helicity configurations and arrangement (permutation) of the momenta
for the quark lines. First, all these diagrams have the same color factor,

$\displaystyle C_{B}^{2}$
$\displaystyle\equiv\delta^{ac}\frac{1}{6}\epsilon_{ijk}\epsilon_{i^{\prime}j^{\prime}k^{\prime}}(T^{a})_{i^{\prime}i}(T^{c}T^{b})_{j^{\prime}j}(T^{b})_{k^{\prime}k}$
$\displaystyle=\left(\frac{2}{3}\right)^{2}.$ (33)

For these diagrams, a pair of quarks has zero total helicity. One can combine
these two fermion lines into a Dirac trace, by applying the following
identity,

$\displaystyle\bar{U}_{\uparrow/\downarrow}(p_{2})\Gamma
U_{\uparrow/\downarrow}(p_{1})=\bar{U}_{\downarrow/\uparrow}(p_{1})\Gamma_{R}U_{\downarrow/\uparrow}(p_{2})\
,$ (34)

where $\Gamma_{R}$ is a $\gamma$-matrix chain obtained by reversing the order
in $\Gamma$. This leads to the typical Dirac algebra for the partonic
amplitude ${\cal H}^{\mu\nu}$,

$\displaystyle\bar{U}_{\uparrow}(p_{2})\Gamma_{1}U_{\uparrow}(p_{1})\
\bar{U}_{\uparrow}(p_{1})\Gamma_{2R}U_{\uparrow}(p_{2})\
\bar{U}_{\uparrow}(p_{2})\Gamma_{3}U_{\uparrow}(p_{1}).$ (35)

It is easy to find out that the first two factors can be combined into a Dirac
trace, and we obtain the following expression,

$\displaystyle\text{Tr}\bigg{[}\frac{1+\gamma_{5}}{2}\not{p}_{2}\Gamma_{1}\frac{1+\gamma_{5}}{2}\not{p}_{1}\Gamma_{2R}\bigg{]}\
\bar{U}_{\uparrow}(p_{2})\Gamma_{3}U_{\uparrow}(p_{1})\ .$ (36)

We will apply the above simplification to all the diagrams.

Furthermore, by examining the two gluon kinematics, we realize that one of the
gluons’ kinematics is determined completely by the quark line that the gluon
attaches. Let us identify that gluon is “$k_{1}$”. Therefore,
$k_{1}=x_{i}p_{1}-y_{i}p_{2}$ where $i$ represents the quark line in the
diagram. With $k_{1}$ determined, we immediately deduce that
$k_{2}=\bar{x}_{i}p_{1}-\bar{y}_{i}p_{2}$.

Therefore, we can classify the partonic scattering amplitudes into two groups:
$k_{1}=x_{i}p_{1}-y_{i}p_{2}$ attaches to the helicity-up quark line (Type-I)
and $k_{1}$ attaches to the helcity-down quark line (Type-II). The derivations
for both types are similar but differ in some details.

Figure 5: Partonic scattering amplitude for the type I configuration: $k_{1}$
is determined by the quark line with helicity-up state.

The typical diagrams of Type-I are shown in Fig. 5, where we include all
possible attachments of $k_{2}$ and the additional gluon exchange between the
two quark lines. The contributions of all these four diagrams can be evaluated
at the same time and will be grouped together. For these diagrams, it is easy
to show that the amplitude can be written as,

$\bar{U}_{\uparrow}(p_{2})\gamma^{\mu}U_{\uparrow}(p_{1}){\rm
Tr}\left[\frac{1+\gamma_{5}}{2}\not{p}_{2}\cdots\gamma^{\nu}\cdots\frac{1+\gamma_{5}}{2}\not{p}_{1}\cdots\right]\
.$ (37)

Because there is no other vector than $p_{1}$, $p_{2}$ and $\nu$, we conclude
that the trace of the second factor is proportional to $p_{1}^{\nu}$ or
$p_{2}^{\nu}$. Explicitly, these four diagrams contribute,

$\displaystyle\frac{p_{2}^{\nu}}{x_{3}y_{3}\bar{x}_{1}}\
,~{}~{}\frac{p_{1}^{\nu}}{x_{3}y_{3}\bar{y}_{1}}\
,~{}~{}\frac{p_{2}^{\nu}}{x_{2}y_{2}\bar{x}_{1}}\
,~{}~{}\frac{p_{1}^{\nu}}{x_{2}y_{2}\bar{y}_{1}}\ .$ (38)

Adding them together, we have,

$\displaystyle\frac{1}{x_{1}y_{1}\bar{x}_{1}\bar{y}_{1}}\left(\frac{1}{x_{2}y_{2}}+\frac{1}{x_{3}y_{3}}\right)\frac{\bar{x}_{1}p_{1}^{\nu}+\bar{y}_{1}p_{2}^{\nu}}{\bar{x}_{1}\bar{y}_{1}}\
,$ (39)

where we have also included the $t$-channel gluon propagators.

Figure 6: Partonic scattering amplitude for the type II configuration: $k_{1}$
is determined by the quark line with helicity-down state.

Typical Type-II diagrams are shown in Fig. 6. The calculations are a little
bit involved. For example, the amplitude can be written in the following form,

$\bar{U}_{\uparrow}(p_{2})\gamma^{\rho}U_{\uparrow}(p_{1}){\rm
Tr}\left[\frac{1+\gamma_{5}}{2}\not{p}_{2}\cdots\gamma^{\nu}\cdots\frac{1+\gamma_{5}}{2}\not{p}_{1}\cdots\gamma_{\rho}\cdots\right]\
.$ (40)

Now the trace of the Gamma matrices can lead to a term like
$\epsilon^{p_{1}p_{2}\nu\rho}$, which can be simplified by applying the
following identity,

$\displaystyle\bar{U}_{\uparrow}(p_{2})\gamma^{\mu}U_{\uparrow}(p_{1})$
$\displaystyle=\frac{i\epsilon^{\mu\nu p_{1}p_{2}}}{p_{1}\cdot
p_{2}}\bar{U}_{\uparrow}(p_{2})\gamma_{\nu}U_{\uparrow}(p_{1})\ .$ (41)

In the end, we find that there is cancellation between different terms and the
Type-II diagrams vanish.

To derive the final result for the amplitude, we need to contract the ${\cal
M}^{\mu\nu}$ with ${\cal H}^{\mu\nu}$ in Eq. (32). The final results can be
summarized as

$\displaystyle{\cal A}_{3}$ $\displaystyle=\langle
J/\psi(\epsilon_{\psi}),N^{\prime}_{\uparrow}|\gamma(\epsilon_{\gamma}),N_{\uparrow}\rangle$
$\displaystyle=\int[dx][dy]\Phi(x_{1},x_{2},x_{3})\Phi^{*}(y_{1},y_{2},y_{3})\frac{1}{(-t)^{2}}$
$\displaystyle~{}~{}\times\bar{U}_{\uparrow}(p_{2})\not{k}_{\gamma}U_{\uparrow}(p_{1}){\cal
M}^{(3)}(\epsilon_{\gamma},\epsilon_{\psi},\\{x\\},\\{y\\})\ .$ (42)

The spinor structure in the above equation is a consequence of the leading-
twist amplitude which conserves the nucleon helicity. This is similar to the
$A$ form factor calculation in Ref. Tong:2021ctu .

###  IV.3 Threshold Expansion

In the threshold limit, we find that ${\cal M}^{(3)}$ can be further
simplified as

$\displaystyle{\cal
M}^{(3)}=\epsilon^{*}_{\psi}\cdot\epsilon_{\gamma}\frac{8e_{c}eg_{s}^{6}}{27\sqrt{3M_{\psi}^{7}}}\psi_{J}(0)\left(2{\cal
H}_{3}+{\cal H^{\prime}}_{3}\right)\ .\quad$ (43)

The coefficient ${\cal H}_{3}$ can be summarized as

$\displaystyle{\cal H}_{3}=I_{13}+I_{31}+I_{12}+I_{32},$ (44)

where

$I_{ij}=\frac{1}{x_{i}x_{j}y_{i}y_{j}\bar{x}_{i}^{2}\bar{y}_{i}}$ (45)

and ${\cal H}_{3}^{\prime}={\cal H}_{3}(y_{1}\leftrightarrow y_{3})$.

Similar to the pion case, we can reproduce the above result by applying the
GPD gluon distribution $H_{g}(x,\xi,t)$ at large momentum transfer in the GPD
formalism. For the reference, we list the GPD gluon distribution $H_{g}$ in
Appendix B.

The final result for the differential cross section will depend on the
threshold limit of the amplitude squared. In the limit of $\chi\to 1$ we find
the following result,

$|\overline{{\cal A}_{3}}|^{2}=(1-\chi)G_{\psi}G_{p3}(t)G_{p3}^{*}(t)\ ,$ (46)

which actually vanishes at the threshold. In the above, the spin sum and
average has been performed, and $G_{\psi}$ has been defined in Eq.(24).
$G_{p3}$ follows the form factor factorization and can be written as

$G_{p3}(t)=\frac{8\pi^{2}\alpha_{s}^{2}C_{B}^{2}}{3t^{2}}\int[dx][dy]\Phi_{3}(\\{x\\})\Phi_{3}^{*}(\\{y\\})\left[2{\cal
H}_{3}+{\cal H}_{3}^{\prime}\right]\ ,$ (47)

where ${\cal H}_{3}$ and ${\cal H}^{\prime}_{3}$ are given above, and
$C_{B}^{2}=(2/3)^{2}$ is the color factor related to partonic amplitudes.
Combining $G_{p3}$ and $G_{p3}^{*}$, this leads to $1/(-t)^{4}$ power behavior
for the amplitude squared, which is consistent with the conventional power
counting analysis. However, this contribution is suppressed at the threshold.

The suppression factor $(1-\chi)$ comes from the spinor structure in Eq. (42).
In order to obtain a nonvanishing contribution at the threshold, we have to go
beyond the leading-twist contributions. In the following section, we consider
the three-quark Fock states with one unit OAM, which are related to the twist-
four distribution amplitudes Ji:2002xn ; Braun:1999te .

##  V Nucleon case: Twist-four contributions

The twist-four contribution comes from the three-quark Fock state with one
unit quark OAM. Two important features emerge for nonzero OAM contributions.
First, the partonic scattering amplitudes conserve the quark helicities.
However, because of a nonzero OAM for one of the three-quark state, the
helicity of the nucleon states will be different. This contributes to the
hadron helicity-flip amplitude. Second, in order to get a nonzero
contribution, we have to perform the intrinsic transverse momentum expansion
for the hard partonic scattering amplitudes Belitsky:2002kj , which will
introduce an additional suppression factor of $1/(-t)$.

The twist-four distribution amplitudes are related to the three-quark Fock
states with one unit of OAM. This can comes from either the initial or final
state. For example, if we consider the contribution from the initial state of
spin-down nucleon, we can parameterize the Fock state as Ji:2002xn ,

$\displaystyle|p_{1}\downarrow\rangle_{1/2}=\int
d[1]d[2]d[3]\left((q_{1}^{x}-iq_{1}^{y})\tilde{\psi}^{(3)}(1,2,3)\right.$
$\displaystyle\left.+(q_{2}^{x}-iq_{2}^{y})\tilde{\psi}^{(4)}(1,2,3)\right)\frac{\epsilon^{ijk}}{\sqrt{6}}$
$\displaystyle\times\left(u^{\dagger}_{i\downarrow}(1)u^{\dagger}_{j\uparrow}(2)d^{\dagger}_{k\uparrow}(3)-d^{\dagger}_{i\downarrow}(1)u^{\dagger}_{j\uparrow}(2)u^{\dagger}_{k\uparrow}(3)\right)|0\rangle\
,$

where the total quark helicity equals to $+1/2$ with nucleon helicity $-1/2$.
With this choice, the final state nucleon’s Fock state can be taken as that in
the previous section.

An important step in the computation of twist-four contribution is to perform
the collinear expansion of the partonic scattering amplitude in terms of the
transverse momenta $q_{i\perp}$. In particular, the linear term of
$q_{i\perp}$ will lead to the twist-four distribution amplitudes when we
integrate over the $q_{i\perp}$ Belitsky:2002kj ,

$\displaystyle\Psi_{4}(x_{1},x_{2},x_{3})=-\frac{2\sqrt{6}}{x_{2}M}\int[dq_{\perp}]$
$\displaystyle~{}~{}\times\vec{q}_{2\perp}\cdot\left[\vec{q}_{1\perp}\tilde{\psi}^{(3)}(1,2,3)+\vec{q}_{2\perp}\tilde{\psi}^{(4)}(1,2,3)\right]\
,$ (48)
$\displaystyle\Phi_{4}(x_{2},x_{1},x_{3})=-\frac{2\sqrt{6}}{x_{3}M}\int[dq_{\perp}]$
$\displaystyle~{}~{}\times\vec{q}_{3\perp}\cdot\left[\vec{q}_{1\perp}\tilde{\psi}^{(3)}(1,2,3)+\vec{q}_{2\perp}\tilde{\psi}^{(4)}(1,2,3)\right]\
.$ (49)

To extract the linear dependence of the transverse momentum $q_{i\perp}$ from
the partonic amplitudes, one can first expand the spinor as

$\displaystyle U(x_{i}p_{1}+\vec{q}_{i\perp})$ $\displaystyle\approx
U(x_{i}p_{1})+\frac{\not{\vec{q}}_{i\perp}\not{p}_{2}}{2x_{i}p_{2}\cdot
p_{1}}U(x_{i}p_{1})\ .$ (50)

After the evaluation of the Dirac structures in the amplitudes following the
strategy in last section, all the linear dependence of $\vec{q}_{i}$ is
explicit and straightforward to find out. For the contributions associated
with the initial OAM, it will yield a structure like:

$\displaystyle{\Gamma}_{1}(\\{x\\},\\{y\\})(q_{1}^{x}+iq_{1}^{y})\bar{U}_{\uparrow}(p_{2})U_{\downarrow}(p_{1})$
$\displaystyle+{\Gamma}_{3}(\\{x\\},\\{y\\})(q_{3}^{x}+iq_{3}^{y})\bar{U}_{\uparrow}(p_{2})U_{\downarrow}(p_{1})\
.$ (51)

where the transverse momentum conservation
$\vec{q}_{2\perp}=-\vec{q}_{1\perp}-\vec{q}_{3\perp}$ is used, and the
identities $\gamma^{x}U_{\uparrow}(p)=U_{\downarrow}(p),\
\gamma^{y}U_{\uparrow}(p)=iU_{\downarrow}(p)$ have been applied.

Applying Eqs. (48,49) with the linear terms of $q_{i\perp}$ from the partonic
amplitudes, we obtain the twist-four contribution to the scattering process of
$\gamma p\to J/\psi p$ as

$\displaystyle{\cal A}_{4}$ $\displaystyle=$ $\displaystyle\langle
J/\psi(\epsilon_{\psi}),N^{\prime}_{\uparrow}|\gamma(\epsilon_{\gamma}),N_{\downarrow}\rangle$
(52) $\displaystyle=$
$\displaystyle\bar{U}_{\uparrow}(p_{2})U_{\downarrow}(p_{1})\frac{M_{p}}{(-t)^{3}}\int[dx][dy]\Phi_{3}^{*}(\\{y\\})$
$\displaystyle\times\left[\Psi_{4}(\\{x\\}){\cal
M}_{\Psi}^{(4)}+\Phi_{4}(\\{x\\}){\cal M}_{\Phi}^{(4)}\right]\ ,$

where $\Psi_{4}$ and $\Phi_{4}$ are the twist-four distributions introduced
above and ${\cal M}_{\Psi,\Phi}^{(4)}$ from the partonic amplitudes. From this
equation, we can clearly see that the nucleon helicity-flip is manifest in the
spinor structure. This amplitude is negligible at high energy, but will be
important at the threshold, because it is not suppressed in the limit of
$\chi\to 1$. The amplitude squared along with the associated spin sum and
average can be written as

$\displaystyle|\overline{{\cal
A}_{4}}|^{2}=\widetilde{m}_{t}^{2}G_{\psi}G_{p4}(t)G_{p4}^{*}(t)\ ,$ (53)

where $\widetilde{m}_{t}^{2}=M_{p}^{2}/(-t)$, $G_{\psi}$ is the same as above.
$G_{p4}$ depends on the twist-three and twist-four distribution amplitudes
Braun:1999te ; Braun:2000kw ,

$\displaystyle G_{p4}(t)$
$\displaystyle=\frac{C_{B}^{2}(4\pi\alpha_{s})^{2}}{12t^{2}}\int[dx][dy]\Phi_{3}(y_{1},y_{2},y_{3})$
$\displaystyle\times\left\\{x_{3}\Phi_{4}(x_{1},x_{2},x_{3}){T}_{4\Phi}(\\{x\\},\\{y\\})\right.$
$\displaystyle\left.+x_{1}\Psi_{4}(x_{2},x_{1},x_{3}){T}_{4\Psi}(\\{x\\},\\{y\\})\right\\}\
\ ,$ (54)

where the hard function has the following form

$\displaystyle{T}_{4\Psi}=2{\cal T}_{4\Psi}+{\cal T}^{\prime}_{4\Psi}\ ,$
$\displaystyle{T}_{4\Phi}=2{\cal T}_{4\Phi}+{\cal T}^{\prime}_{4\Phi}\ ,$ (55)

and ${\cal T}^{\prime}_{4}$ is obtained from ${\cal T}_{4}$ by interchanging
$y_{1}$ and $y_{3}$. Then we have

$\displaystyle{\cal T}_{4\Psi}=$ $\displaystyle
x_{3}K_{1}(1+y_{2}/\bar{y}_{1})+2\bar{x}_{3}\tilde{K}_{1}$
$\displaystyle+2x_{3}(\tilde{K}_{2}-K_{2})-K_{3}/\bar{y}_{1}$
$\displaystyle+x_{3}(K_{4}+K_{5})/\bar{x}_{1}+2(\tilde{K}_{4}+\tilde{K}_{5})~{},$
$\displaystyle{\cal T}_{4\Phi}=$ $\displaystyle{\cal
T}_{4\Psi}(1\leftrightarrow 3)~{},$ (56)

where the functions $K_{i}$ and $\tilde{K}_{i}$ are defined as

$\displaystyle
K_{1}=\frac{1}{x_{1}x_{3}^{2}y_{1}y_{3}^{2}\bar{x}_{1}^{2}\bar{y}_{1}}~{},\quad
K_{2}=\frac{1}{x_{1}x_{2}x_{3}^{2}y_{2}y_{3}^{2}\bar{x}_{2}\bar{y}_{2}}~{},$
$\displaystyle
K_{3}=\frac{1}{x_{1}x_{2}y_{1}y_{2}\bar{x}_{1}^{2}\bar{y}_{1}}~{},\quad
K_{4}=\frac{1}{x_{1}x_{3}^{2}y_{1}y_{3}\bar{x}_{1}\bar{y}_{1}^{2}}~{},$
$\displaystyle
K_{5}=\frac{1}{x_{1}x_{2}x_{3}y_{1}y_{2}\bar{x}_{1}\bar{y}_{1}^{2}}~{},\quad\tilde{K}_{i}=K_{i}(1\leftrightarrow
3)~{}.$ (57)

As mentioned above, the twist-four distribution amplitudes can come from both
initial and final state nucleons. Because of the symmetric property of the
partonic scattering amplitudes, these two contributions are the same and have
been included in the above final result.

Eqs. (53) and (46) are the final results of our analysis. Comparing these two,
we find that the twist-four contribution is suppressed in $1/t$ but enhanced
at the threshold. These two features can be used to disentangle their
contributions in experiments. If we limit our discussions in the threshold
region, the only contribution comes from the twist-four term.

##  VI Interpretation in Terms of Gravitational Form Factor?

As mentioned in Introduction, the near threshold heavy quarkonium production
has been argued to provide a direct access to the gluonic gravitational form
factors of the nucleon. However, our explicit calculations for the pion case
have shown that there is no direct connection between them.

From the results in previous sections, we have calculated the near threshold
photo-production of heavy quarkonium on the nucleon target at large momentum
transfer. The gluonic form factors at large $(-t)$ have been recently
calculated in Ref. Tong:2021ctu . We conclude, again, we can not build a
direct connection between them.

###  VI.1 Construct the Gluonic Operator

The above conclusion can be understood from a detailed analysis of the photon-
quarkonium transition amplitude. As discussed in Sec. II, this amplitude can
be simplified in the heavy quark mass limit, $M_{V}^{2}\gg(-t)$,

$\displaystyle{\cal
M}_{\psi}^{\mu\nu}=N_{\psi}\epsilon_{\psi}^{*}\cdot\epsilon_{\gamma}\frac{k_{\gamma,\alpha}k_{\gamma,\beta}}{k_{1}\cdot
k_{\gamma}k_{2}\cdot k_{\gamma}}{\cal W}_{T}^{\alpha\beta\mu\nu}\ .$ (58)

Here, we only keep the leading term in this limit. For simplicity, we have
also dropped the associated color factors associated with the $t$-channel
gluons. In the above, ${\cal W}_{T}^{\alpha\beta\mu\nu}$ is defined as

$\displaystyle{\cal W}_{T}^{\alpha\beta\mu\nu}$ $\displaystyle=$
$\displaystyle-k_{1}^{\alpha}k_{2}^{\beta}g^{\mu\nu}-k_{1}\cdot
k_{2}g^{\alpha\mu}g^{\beta\nu}$ (59)
$\displaystyle+k_{1}^{\nu}k_{2}^{\beta}g^{\alpha\mu}+k_{2}^{\mu}k_{1}^{\alpha}g^{\beta\nu}\
,$

which can be identified as gluonic operator of
${F^{\alpha}}_{\rho}F^{\beta\rho}$ acting on the nucleon state. However, the
complete scattering amplitude involves the integral of the momenta $k_{1}$ and
$k_{2}$ with the associated propagators depending on them. In the end, the
$\gamma N\to J/\psi N^{\prime}$ amplitude can be schematically written as

$\displaystyle{\cal A}$ $\displaystyle=$ $\displaystyle
N_{\psi}\epsilon_{\psi}^{*}\cdot\epsilon_{\gamma}$ (60)
$\displaystyle\times\int
d^{4}k_{1}d^{4}k_{2}\frac{k_{\gamma,\alpha}k_{\gamma,\beta}}{(k_{1}\cdot
k_{\gamma}-i\varepsilon)(k_{2}\cdot k_{\gamma}-i\varepsilon)}$
$\displaystyle\times\int
d^{4}\eta_{1}d^{4}\eta_{2}e^{ik_{1}\cdot\eta_{1}+ik_{2}\cdot\eta_{2}}$
$\displaystyle~{}~{}~{}~{}\times\langle
N^{\prime}|{F^{a,\alpha}}_{\rho}(\eta_{1})F^{a,\beta\rho}(\eta_{2})|N\rangle\
.$

Clearly, if we neglect the $k_{1}$ and $k_{2}$ dependence in the pre-factor of
$\frac{1}{(k_{1}\cdot k_{\gamma}-i\varepsilon)(k_{2}\cdot
k_{\gamma}-i\varepsilon)}$, the above equation can be identified as a gluonic
gravitational form factor of the nucleon state. However, as discussed in Sec.
II, this pre-factor comes from the quark propagators in the photon-quarkonium
transition amplitude. The complete calculation will have a full dependence on
the momentum fractions of the incoming nucleon $p_{1}$ carried by the two
gluons $k_{1}$ and $k_{2}$.

We emphasize that the above discussions apply to all of the kinematics in
heavy quarkonium photo-production, including small and large $(-t)$.
Therefore, our conclusion is valid in the whole kinematics of this process
that there is no direct connection between the near threshold photo-prodution
of heavy quarkonium and the gluonic gravitational form factors of the nucleon.

###  VI.2 Compare to the GPD Formalism

It is interesting to find out that the above ${\cal W}_{T}^{\mu\nu}$ can be
directly compared to that for the gluon GPD calculations. Gluon GPD is defined
through the matrix element $\langle N^{\prime}|F^{+\alpha}F^{+}_{\
\alpha}|N\rangle$. The amplitude associated with this can be written as

$-n\cdot k_{1}n\cdot k_{2}g^{\mu\nu}-n^{\mu}n^{\nu}k_{1}\cdot
k_{2}+n^{\mu}k_{1}^{\nu}n\cdot k_{2}+n^{\nu}k_{2}^{\mu}n\cdot k_{1}\ ,$ (61)

where $n$ is the light-cone vector used in the GPD definition with $n\cdot
k=k^{+}$ for any momentum $k$. Here, $k_{1}$ and $k_{2}$ represent the gluon
momenta that couple to the nucleon state, and $\mu\nu$ for their polarization
indices. Clearly, this is the same structure as ${\cal W}_{T}^{\mu\nu}$ of
previous subsection if we identify $n\propto k_{\gamma}$.

Following this argument, the scattering amplitude of $\gamma N\to J/\psi
N^{\prime}$ can be formulated in terms of the gluon GPDs Hoodbhoy:1996zg ;
Ivanov:2004vd ; Koempel:2011rc ; Boussarie:2020vmu ; Hatta:2021can ;
Guo:2021ibg ,

$\displaystyle{\cal
A}=N_{\psi}\epsilon_{\psi}^{*}\cdot\epsilon_{\gamma}\int_{-1}^{1}dx\frac{1}{(x+\xi-i\varepsilon)(x-\xi+i\varepsilon)}$
$\displaystyle\times$
$\displaystyle\frac{1}{\bar{P}^{+}}\int\frac{d\eta^{-}}{2\pi}e^{ix\bar{P}^{+}\eta^{-}}\langle
N^{\prime}|{F^{a,+}}_{\alpha}(-\frac{\eta^{-}}{2})F^{a,\alpha+}(\frac{\eta^{-}}{2})|N\rangle\
,$

where the last factor defines the associated gluon GPDs. In previous sections,
we have given explicit examples that demonstrate the consistency between our
calculations with the GPD formalism.

Clearly, from the above GPD formalism, one can only link to the gluonic
gravitational form factors by making approximations of no $x$-dependence in
the pre-factor $\frac{1}{(x+\xi-i\varepsilon)(x-\xi+i\varepsilon)}$
Hatta:2021can ; Guo:2021ibg . This is the same as we discussed in the previous
subsection. Therefore, our conclusion of no direct connection between the near
threshold photo-production of heavy quarkonium state and the gluonic
gravitational form factors is consistent with the GPD formalism.

##  VII Phenomenology Applications

Taking into account the contributions derived in previous sections, we can
write down the differential cross section for the near threshold heavy
quarkonium photoproduction at large momentum transfer,

$\displaystyle\frac{d\sigma}{dt}|_{(-t)\gg\Lambda_{QCD}^{2}}$ $\displaystyle=$
$\displaystyle\frac{1}{16\pi(W_{\gamma
p}^{2}-M_{p}^{2})^{2}}\left(|\overline{{\cal A}_{3}}|^{2}+|\overline{{\cal
A}_{4}}|^{2}\right)$ (63) $\displaystyle\approx$
$\displaystyle\frac{1}{(-t)^{4}}\left[(1-\chi){\cal
N}_{3}+\widetilde{m}_{t}^{2}{\cal N}_{4}\right]\ ,$

where ${\cal N}_{3}$ and ${\cal N}_{4}$ represent the twist-three and twist-
four contributions, respectively. The most important consequence of our power
counting analysis is that the leading-twist contribution is suppressed at the
threshold. Away from the threshold point, it will start to contribute and may
dominate at large $(-t)$ because of the leading power feature. With high
precision future experiments Chen:2014psa ; AbdulKhalek:2021gbh ;
Joosten:2018gyo , we should be able to distinguish their contributions.

If we take the leading contribution of Eq. (63) at the threshold, i.e., the
${\cal N}_{4}$ term, the differential cross section only depends on the
momentum transfer $t$. This is an important signal from the perturbative QCD
analysis in this paper. Of course, away from the threshold region, we have to
take into account additional contribution from ${\cal N}_{3}$ and the
kinematic corrections in Eq. (63). In Ref. Sun:2021gmi , the twist-four
contribution has been applied to fit the GlueX data Ali:2019lzf  with the
following parameterization of the differential cross section,

$\displaystyle\frac{d\sigma}{dt}|^{twist-4}=\frac{N_{0}}{(-t+\Lambda^{2})^{5}}\
,$ (64)

where $\Lambda^{2}=1.41\pm 0.20~{}{\rm GeV}^{2}$ and $N_{0}=51\pm 22~{}{\rm
nb*GeV^{8}}$. The current data from the GlueX can be well described by the
above parameterization. In the following, we will apply this result to the
future experiments for the threshold photo-production of $\psi^{\prime}$ and
$\Upsilon$.

We have also made an order of magnitude estimate of the differential cross
section by applying the twist-four contribution of Eq. (53) with model
assumptions for the twist-three and twist-four distribution amplitudes of the
nucleon Braun:1999te ; Braun:2000kw . There have been great efforts to compute
these distribution amplitudes from various methods Ioffe:1981kw ; Chung:1981cc
; Chernyak:1983ej ; Chernyak:1984bm ; King:1986wi ; Chernyak:1987nt ;
Chernyak:1987nu ; Bolz:1996sw ; Braun:2001tj ; Braun:2006hz ; Gockeler:2008xv
; Braun:2008ia ; QCDSF:2008qtn ; Passek-Kumericki:2008uqr ; Lenz:2009ar ;
Anikin:2013aka ; Braun:2014wpa ; Bali:2015ykx ; RQCD:2019hps , including the
lattice QCD, the light-cone sum rule and model calculations. The differential
cross sections calculated from the distribution amplitudes with realistic
model assumptions, e.g., those from Ref. Braun:2006hz , are consistent with
the experimental data around $-t=1.5\rm GeV^{2}$ and the fitted result of Eq.
(64) 111In the numeric calculation of the $G_{p4}(t)$ in Eq. (54), a lower
cutoff ($\sim(0.17{\rm GeV})^{2}/(-t)$) on the momentum fractions $x_{i}$ and
$y_{i}$ in the integral is imposed to avoid the end-point singularity. This is
similar to the Pauli form factor calculation at large momentum transfer in
Ref. Belitsky:2002kj ., whereas the results from the asymptotic distribution
amplitudes are an order of magnitude smaller.

###  VII.1 Predictions for $\psi^{\prime}$ and $\Upsilon$(nS) Production

Figure 7: Differential cross sections for $J/\psi$ and $\psi^{\prime}$ photo-
production as functions of the momentum transfer $t$ and the total cross
sections near the threshold as functions of $W_{\gamma p}$.

Extending our analysis to other heavy quarkonium states is straightforward and
similar formulas can be derived. As a first step, we take the differential
cross section from the twist-four contribution at the threshold,

$\displaystyle\frac{d\sigma(\gamma p\to
Vp)}{dt}|_{threshold}=\frac{N_{0}^{V}}{(-t+\Lambda^{2})^{5}}\ ,$ (65)

for a heavy quarkonium state $V$. In the heavy quark mass limit, the
$t$-dependence only comes from the nucleon side. Therefore, we will assume the
above $\Lambda$ parameter should be same for all heavy quarkonium states. On
the other hand, the normalization factor $N_{0}^{V}$ will depend on the
quarkonium state in the final state. From the derivations in previous section,
we know that the differential cross section is proportional to,

$\frac{d\sigma}{dt}\propto\frac{\alpha_{s}^{2}(M_{V})\langle 0|{\cal
O}^{V}({}^{3}S_{1}^{(1)})|0\rangle}{M_{V}^{7}}\ ,$ (66)

from which we derive the ratio between different heavy quarkonium states,

$\frac{N_{0}^{V}}{N_{0}}=\frac{\alpha_{s}^{2}(M_{V})\langle 0|{\cal
O}^{V}({}^{3}S_{1}^{(1)})|0\rangle/M_{V}^{7}}{\alpha_{s}^{2}(M_{\psi})\langle
0|{\cal O}^{\psi}({}^{3}S_{1}^{(1)})|0\rangle/M_{\psi}^{7}}\ .$ (67)

Substituting the associated NRQCD matrix elements for $J/\psi$,
$\psi^{\prime}$, and $\Upsilon$ (1S, 2S) from, e.g., Refs. Feng:2015sjx ;
Zhan:2021dlu , we find the following values for the normalization factors,

$\displaystyle{N_{0}^{\psi^{\prime}}}=0.20{N_{0}}\ ,$ (68)
$\displaystyle{N_{0}^{\Upsilon(1S)}}=5\times 10^{-3}{N_{0}}\ ,$ (69)
$\displaystyle{N_{0}^{\Upsilon(2S)}}=2.5\times 10^{-3}{N_{0}}\ .$ (70)

In Fig. 7, we show the threshold cross sections for $\gamma
p\to\psi^{\prime}p$. As comparison, we also show the results for $J/\psi$. For
Upsilon production, the results are plotted in Fig. 8.

Figure 8: The differential cross sections for $\Upsilon$(1S) and
$\Upsilon$(2S) photo-production as functions of the momentum transfer $t$ and
the total cross sections as functions $W_{\gamma p}$ near the threshold.

The comparison between different quarkonium states will provide an important
confirmation for the production mechanism. The comparison between Charmonium
and Bottomonium, in particular, will test the heavy quark limit we have
employed in this paper. Meanwhile, the momentum transfer range is much higher
for $\Upsilon$ as compared to $J/\psi$. This provides a unique opportunity to
explore the large momentum transfer region.

##  VIII Conclusion

In this paper, we have carried out a detailed derivation of near threshold
heavy quarkonium photoproduction at large momentum transfer. We have taken
into account the three-quark Fock state of the nucleon with zero and one unit
of quark OAM. We found that the contribution from the Fock state with zero
quark OAM is suppressed at threshold. The differential cross section is
dominated by the contribution from nonzero OAM Fock state and has a power
behavior of $1/(-t)^{5}$.

Our power counting predictions are consistent with recent experimental data of
near threshold photo-production of $J/\psi$ from the GlueX collaboration at
JLab. Based on the comparison between our derivation and the experimental
data, we have made predictions for $\psi^{\prime}$ and $\Upsilon$(1S,2S). All
these predictions can be tested at future facilities including the electron-
ion colliders Accardi:2012qut ; AbdulKhalek:2021gbh ; Anderle:2021wcy .

We have also shown that there is no direct connection between the near
threshold photo-production of heavy quarkonium state and the gluonic
gravitational form factors. The indirect connection can be built through GPD
gluon distributions. For example, we can parameterize the gluon GPDs and fit
to the experimental data, which, in return, can constrain the associated
gravitational form factors.

Acknowledgments: We thank Yoshitaka Hatta, Xiangdong Ji, and Nu Xu for
discussions and comments. This material is based upon work supported by the
LDRD program of Lawrence Berkeley National Laboratory, the U.S. Department of
Energy, Office of Science, Office of Nuclear Physics, under contract numbers
DE-AC02-05CH11231. P. Sun is supported by Natural Science Foundation of China
under grant No. 11975127 and No. 12061131006 as well as Jiangsu Specially
Appointed Professor Program. X. B. Tong is supported by the CUHK-Shenzhen
university development fund under grant No. UDF01001859.

##  Appendix A Gluon GPD for Pion

The gluon GPD for pion is defined as

$\displaystyle\int\frac{d\eta^{-}}{2\pi}e^{ix\bar{P}^{+}\eta^{-}}\Big{\langle}p_{2}\Big{|}F^{+\mu}_{a}(-\frac{\eta^{-}}{2}){\cal
L}_{ab}[-\frac{\eta^{-}}{2},\frac{\eta^{-}}{2}]$
$\displaystyle\quad\quad\times F^{\
+}_{\mu,b}(\frac{\eta^{-}}{2})\Big{|}p_{1}\Big{\rangle}$ $\displaystyle=$
$\displaystyle\bar{P}^{+}H^{(\pi)}_{g}(x,\xi,t)~{},$ (71)

where the gluon field strength tensor is
$F^{a}_{\mu\nu}=\partial_{\mu}A^{a}_{\nu}-\partial_{\nu}A^{a}_{\mu}-g_{s}f^{abc}A^{b}_{\mu}A^{c}_{\nu}$,
and the gauge link in the adjoint representation is

$\displaystyle{\cal L}_{ab}\left[z_{2},z_{1}\right]={\cal P}\
\text{exp}\left[g_{s}\int^{z_{2}}_{z_{1}}\ dz^{-}\
G^{+,c}(z^{-})f^{acb}\right]~{}.$ (72)

${\cal P}$ denote the path-ordering operation. In the definition of GPD,
$\bar{P}=(p_{1}+p_{2})/2$ is the average momentum, $\Delta=p_{2}-p_{1}$ is the
momentum transfer and $t=\Delta^{2}$. The skewness parameter $\xi$ is defined
as the projection of the momentum transfer $\Delta$ along $\bar{P}$ direction,
$\xi=-\frac{\Delta^{+}}{2\bar{P}^{+}}.$

In the large $(-t)$ limit, the gluon GPD of the pion has the following
factorization formula:

$\displaystyle H^{(\pi)}_{g}(x,\xi,t)=$ $\displaystyle\int
dx_{1}dy_{1}\phi(y_{1})\phi(x_{1}){\cal H}^{(\pi)}(x_{1},y_{1})~{}\ ,$ (73)

where $\phi$ represents the leading-twist distribution amplitude of pion. The
perturbative function at the leading order can be written as

$\displaystyle{\cal H}^{(\pi)}(x_{1},y_{1})=$
$\displaystyle\frac{g_{s}^{2}C_{F}}{-t}\left(\frac{(1-\xi)^{2}}{x_{1}\bar{x}_{1}}+\frac{(1+\xi)^{2}}{y_{1}\bar{y}_{1}}\right)$
$\displaystyle\times\bigg{(}\delta\big{[}x-(x_{1}-y_{1}+\xi(x_{1}+y_{1}-1))\big{]}$
$\displaystyle+\delta\big{[}x+(x_{1}-y_{1}+\xi(x_{1}+y_{1}-1))\big{]}\bigg{)}~{}.$
(74)

The above result is similar to that of the quark GPD calculated in Ref.
Hoodbhoy:2003uu  for the pion.

##  Appendix B Gluon GPD for Nucleon

The gluon GPD of nucleon is defined from

$\displaystyle\int\frac{d\eta^{-}}{2\pi}e^{ix\bar{P}^{+}\eta^{-}}\Big{\langle}p_{2},s^{\prime}\Big{|}F^{+\mu}_{a}(-\frac{\eta^{-}}{2}){\cal
L}_{ab}[-\frac{\eta^{-}}{2},\frac{\eta^{-}}{2}]$ $\displaystyle\quad\times
F^{\ +}_{\mu,b}(\frac{\eta^{-}}{2})\Big{|}p_{1},s\Big{\rangle}$
$\displaystyle=$
$\displaystyle\frac{1}{2}\bigg{(}H_{g}(x,\xi,t)\bar{U}(p_{2},s^{\prime})\gamma^{+}U(p_{1},s)$
$\displaystyle\quad+E_{g}(x,\xi,t)\bar{U}(p_{2},s^{\prime})\frac{i\sigma^{+\alpha}\Delta_{\alpha}}{2M_{p}}U(p_{1},s)\bigg{)}~{},$
(75)

where $s^{\mu}$ denote the covariant spin-vector of the proton.

Following the strategy in Hoodbhoy:2003uu ; Tong:2021ctu , the GPD $H_{g}$ can
be extracted from the helicity conserved amplitude, and one can show that at
the large momentum transfer, $H_{g}$ follows the following factorization
formula:

$\displaystyle H_{g}(x,\xi,t)=$
$\displaystyle\int[dx][dy]\Phi_{3}^{*}(y_{1},y_{2},y_{3})\Phi_{3}(x_{1},x_{2},x_{3})$
$\displaystyle\quad\quad\times{\cal H}(\\{x\\},\\{y\\})~{},$ (76)

where $\Phi_{3}$ is the twist-3 proton light-cone amplitude Braun:1999te , and
${\cal H}$ is the hard coefficient and perturbatively calculable. At the
leading order, we obtain

$\displaystyle{\cal H}(\\{x\\},\\{y\\})$ $\displaystyle=2\tilde{\cal
H}+\tilde{\cal H}(y_{1}\leftrightarrow y_{3})~{},$ (77)

where

$\displaystyle\tilde{\cal H}$
$\displaystyle=\frac{4\pi^{2}\alpha_{s}^{2}C_{B}^{2}}{3t^{2}}\times$ (78)
$\displaystyle\bigg{\\{}\bigg{(}\frac{x_{1}+y_{1}+\xi(x_{1}-y_{1})}{\bar{x}_{1}\bar{y}_{1}x_{1}x_{3}y_{1}y_{3}}+\frac{x_{1}+y_{1}+\xi(x_{1}-y_{1})}{\bar{x}_{1}\bar{y}_{1}x_{1}x_{2}y_{1}y_{2}}\bigg{)}$
$\displaystyle\times\bigg{(}\delta\big{[}x-(x_{1}-y_{1}+\xi(x_{1}+y_{1}-1))\big{]}$
$\displaystyle+\delta\big{[}x+(x_{1}-y_{1}+\xi(x_{1}+y_{1}-1))\big{]}\bigg{)}$
$\displaystyle+\bigg{(}\frac{x_{3}+y_{3}+\xi(x_{3}-y_{3})}{\bar{x}_{3}\bar{y}_{3}x_{3}x_{1}y_{3}y_{1}}+\frac{x_{3}+y_{3}+\xi(x_{3}-y_{3})}{\bar{x}_{3}\bar{y}_{3}x_{3}x_{2}y_{3}y_{2}}\bigg{)}$
$\displaystyle\times\bigg{(}\delta\big{[}x-(x_{3}-y_{3}+\xi(x_{3}+y_{3}-1))\big{]}$
$\displaystyle+\delta\big{[}x+(x_{3}-y_{3}+\xi(x_{3}+y_{3}-1))\big{]}\bigg{)}\bigg{\\}}~{}.$

Similar results for the quark GPDs $H_{q}$ of the nucleon have been calculated
in Ref. Hoodbhoy:2003uu . They share the same power behavior at large $(-t)$.

On the other hand, the GPD $E_{g}$ at large $(-t)$ is calculated from the
nucleon helicity-flip amplitude, and the related factorization formula can be
written as

$\displaystyle
E_{g}(x,\xi,t)=\int[dx][dy]\left\\{x_{3}\Phi_{4}(x_{1},x_{2},x_{3}){\cal
E}_{\Phi g}(\\{x\\},\\{y\\})\right.$
$\displaystyle\left.+x_{1}\Psi_{4}(x_{2},x_{1},x_{3}){\cal E}_{\Psi
g}(\\{x\\},\\{y\\})\right\\}\ \Phi_{3}(y_{1},y_{2},y_{3})\ ,$ (79)

where $\Psi_{4}$ and $\Phi_{4}$ are the twist-four distribution amplitude of
the proton Braun:2000kw . ${\cal E}_{g}$ can be written as,

$\displaystyle{\cal E}_{g}=2\tilde{\cal E}+\tilde{\cal E}^{\prime}~{},$ (80)

where $\tilde{\cal E}^{\prime}$ is obtained from $\tilde{\cal E}$ by
interchanging $y_{1}$ and $y_{3}$. The detailed calculation yields

$\displaystyle\tilde{\cal
E}_{\Psi}(\\{x\\},\\{y\\})=\frac{-C_{B}^{2}M^{2}_{p}}{12(-t)^{3}}(4\pi\alpha_{s})^{2}$
$\displaystyle\times$
$\displaystyle\Bigg{[}x_{3}K_{1}\tilde{\delta}[x_{1},y_{1}]\Bigl{(}(1+\xi)^{2}x_{1}\bar{x}_{1}+2\left(1-\xi^{2}\right)y_{3}\bar{x}_{1}$
$\displaystyle+(1-\xi)^{2}y_{1}y_{2}\Bigr{)}+\bar{x}_{3}\tilde{K}_{1}\tilde{\delta}[x_{3},y_{3}]\Bigl{(}(1+\xi)^{2}x_{3}\bar{x}_{3}$
$\displaystyle+(1-\xi)^{2}y_{3}\bar{y}_{3}\Bigr{)}+x_{3}(\tilde{K}_{2}-K_{2})\tilde{\delta}[x_{2},y_{2}]$
$\displaystyle\Bigl{(}(1+\xi)^{2}x_{2}\bar{x}_{2}+(1-\xi)^{2}y_{2}\bar{y}_{2}\Bigr{)}+K_{3}\tilde{\delta}[x_{1},y_{1}]$
$\displaystyle\Bigl{(}2\left(1-\xi^{2}\right)\bar{x}_{1}-(1-\xi)^{2}y_{1}\Bigr{)}+x_{3}(K_{4}+K_{5})$
$\displaystyle\tilde{\delta}[x_{1},y_{1}]\Bigl{(}2\left(1-\xi^{2}\right)\bar{y}_{1}+(1+\xi)^{2}x_{1}\Bigr{)}+(\tilde{K}_{4}+\tilde{K}_{5})$
$\displaystyle\tilde{\delta}[x_{3},y_{3}]\Bigl{(}(1+\xi)^{2}x_{3}\bar{x}_{3}+(1-\xi)^{2}y_{3}\bar{y}_{3}\Bigr{)}\Bigg{]}+(\xi\leftrightarrow-\xi)~{},$
(81)

and $\tilde{\cal E}_{\Phi}=\tilde{\cal E}_{\Psi}(1\leftrightarrow 3)$. Here we
have used the following notation to express the delta function of $x$:

$\displaystyle\tilde{\delta}[a,b]\equiv$
$\displaystyle\delta[x-(a-b+\xi(a+b-1))]$
$\displaystyle+\delta[x+(a-b+\xi(a+b-1))]~{}.$ (82)

The functions $K_{i}$ are the same as those defined in Eq.(57).

## References

  * (1) M. G. Ryskin, Z. Phys. C 57, 89-92 (1993) doi:10.1007/BF01555742
  * (2) S. J. Brodsky, L. Frankfurt, J. F. Gunion, A. H. Mueller and M. Strikman, Phys. Rev. D 50, 3134-3144 (1994) doi:10.1103/PhysRevD.50.3134 [arXiv:hep-ph/9402283 [hep-ph]].
  * (3) X. D. Ji, Phys. Rev. Lett. 78, 610-613 (1997) doi:10.1103/PhysRevLett.78.610 [arXiv:hep-ph/9603249 [hep-ph]].
  * (4) X. D. Ji, Phys. Rev. D 55, 7114-7125 (1997) doi:10.1103/PhysRevD.55.7114 [arXiv:hep-ph/9609381 [hep-ph]].
  * (5) J. C. Collins, L. Frankfurt and M. Strikman, Phys. Rev. D 56, 2982-3006 (1997) doi:10.1103/PhysRevD.56.2982 [arXiv:hep-ph/9611433 [hep-ph]].
  * (6) P. Hoodbhoy, Phys. Rev. D 56, 388-393 (1997) doi:10.1103/PhysRevD.56.388 [arXiv:hep-ph/9611207 [hep-ph]].
  * (7) J. Koempel, P. Kroll, A. Metz and J. Zhou, Phys. Rev. D 85, 051502 (2012) doi:10.1103/PhysRevD.85.051502 [arXiv:1112.1334 [hep-ph]].
  * (8) Z. L. Cui, M. C. Hu and J. P. Ma, Eur. Phys. J. C 79 (2019) no.10, 812 doi:10.1140/epjc/s10052-019-7298-y [arXiv:1804.05293 [hep-ph]].
  * (9) D. Y. Ivanov, A. Schafer, L. Szymanowski and G. Krasnikov, Eur. Phys. J. C 34, no.3, 297-316 (2004) [erratum: Eur. Phys. J. C 75, no.2, 75 (2015)] doi:10.1140/epjc/s2004-01712-x [arXiv:hep-ph/0401131 [hep-ph]].
  * (10) Z. Q. Chen and C. F. Qiao, Phys. Lett. B 797, 134816 (2019) [erratum: Phys. Lett. B, 135759 (2020)] doi:10.1016/j.physletb.2019.134816 [arXiv:1903.00171 [hep-ph]].
  * (11) C. A. Flett, J. A. Gracey, S. P. Jones and T. Teubner, JHEP 08, 150 (2021) doi:10.1007/JHEP08(2021)150 [arXiv:2105.07657 [hep-ph]].
  * (12) D. Kharzeev, Proc. Int. Sch. Phys. Fermi 130, 105-131 (1996) doi:10.3254/978-1-61499-215-8-105 [arXiv:nucl-th/9601029 [nucl-th]].
  * (13) D. Kharzeev, H. Satz, A. Syamtomov and G. Zinovjev, Eur. Phys. J. C 9, 459-462 (1999) doi:10.1007/s100529900047 [arXiv:hep-ph/9901375 [hep-ph]].
  * (14) O. Gryniuk and M. Vanderhaeghen, Phys. Rev. D 94, no.7, 074001 (2016) doi:10.1103/PhysRevD.94.074001 [arXiv:1608.08205 [hep-ph]].
  * (15) Y. Hatta and D. L. Yang, Phys. Rev. D 98, no.7, 074003 (2018) doi:10.1103/PhysRevD.98.074003 [arXiv:1808.02163 [hep-ph]].
  * (16) A. Ali et al. [GlueX], Phys. Rev. Lett. 123, no.7, 072001 (2019) doi:10.1103/PhysRevLett.123.072001 [arXiv:1905.10811 [nucl-ex]].
  * (17) Y. Hatta, A. Rajan and D. L. Yang, Phys. Rev. D 100, no.1, 014032 (2019) doi:10.1103/PhysRevD.100.014032 [arXiv:1906.00894 [hep-ph]].
  * (18) R. Boussarie and Y. Hatta, Phys. Rev. D 101, no.11, 114004 (2020) doi:10.1103/PhysRevD.101.114004 [arXiv:2004.12715 [hep-ph]].
  * (19) K. A. Mamo and I. Zahed, Phys. Rev. D 101, no.8, 086003 (2020) doi:10.1103/PhysRevD.101.086003 [arXiv:1910.04707 [hep-ph]].
  * (20) O. Gryniuk, S. Joosten, Z. E. Meziani and M. Vanderhaeghen, Phys. Rev. D 102, no.1, 014016 (2020) doi:10.1103/PhysRevD.102.014016 [arXiv:2005.09293 [hep-ph]].
  * (21) R. Wang, J. Evslin and X. Chen, Eur. Phys. J. C 80, no.6, 507 (2020) doi:10.1140/epjc/s10052-020-8057-9 [arXiv:1912.12040 [hep-ph]].
  * (22) F. Zeng, X. Y. Wang, L. Zhang, Y. P. Xie, R. Wang and X. Chen, Eur. Phys. J. C 80, no.11, 1027 (2020) doi:10.1140/epjc/s10052-020-08584-6 [arXiv:2008.13439 [hep-ph]].
  * (23) M. L. Du, V. Baru, F. K. Guo, C. Hanhart, U. G. Meißner, A. Nefediev and I. Strakovsky, Eur. Phys. J. C 80, no.11, 1053 (2020) doi:10.1140/epjc/s10052-020-08620-5 [arXiv:2009.08345 [hep-ph]].
  * (24) D. E. Kharzeev, Phys. Rev. D 104, no.5, 054015 (2021) doi:10.1103/PhysRevD.104.054015 [arXiv:2102.00110 [hep-ph]].
  * (25) R. Wang, W. Kou, Y. P. Xie and X. Chen, Phys. Rev. D 103, no.9, L091501 (2021) doi:10.1103/PhysRevD.103.L091501 [arXiv:2102.01610 [hep-ph]].
  * (26) Y. Hatta and M. Strikman, Phys. Lett. B 817, 136295 (2021) doi:10.1016/j.physletb.2021.136295 [arXiv:2102.12631 [hep-ph]].
  * (27) K. A. Mamo and I. Zahed, Phys. Rev. D 103, no.9, 094010 (2021) doi:10.1103/PhysRevD.103.094010 [arXiv:2103.03186 [hep-ph]].
  * (28) W. Kou, R. Wang and X. Chen, [arXiv:2103.10017 [hep-ph]].
  * (29) Y. Guo, X. Ji and Y. Liu, Phys. Rev. D 103, no.9, 096010 (2021) doi:10.1103/PhysRevD.103.096010 [arXiv:2103.11506 [hep-ph]].
  * (30) K. A. Mamo and I. Zahed, Phys. Rev. D 104, no.6, 066023 (2021) doi:10.1103/PhysRevD.104.066023 [arXiv:2106.00722 [hep-ph]].
  * (31) M. A. Shifman, A. I. Vainshtein and V. I. Zakharov, Phys. Lett. B 78, 443-446 (1978) doi:10.1016/0370-2693(78)90481-1
  * (32) X. D. Ji, Phys. Rev. Lett. 74, 1071-1074 (1995) doi:10.1103/PhysRevLett.74.1071 [arXiv:hep-ph/9410274 [hep-ph]].
  * (33) X. D. Ji, Phys. Rev. D 52, 271-281 (1995) doi:10.1103/PhysRevD.52.271 [arXiv:hep-ph/9502213 [hep-ph]].
  * (34) Y. Hatta, A. Rajan and K. Tanaka, JHEP 12, 008 (2018) doi:10.1007/JHEP12(2018)008 [arXiv:1810.05116 [hep-ph]].
  * (35) A. Metz, B. Pasquini and S. Rodini, Phys. Rev. D 102, 114042 (2020) doi:10.1103/PhysRevD.102.114042 [arXiv:2006.11171 [hep-ph]].
  * (36) Y. Hatta and Y. Zhao, Phys. Rev. D 102, no.3, 034004 (2020) doi:10.1103/PhysRevD.102.034004 [arXiv:2006.02798 [hep-ph]].
  * (37) X. Ji and Y. Liu, Sci. China Phys. Mech. Astron. 64, no.8, 281012 (2021) doi:10.1007/s11433-021-1723-2 [arXiv:2101.04483 [hep-ph]].
  * (38) X. Ji, Front. Phys. (Beijing) 16, no.6, 64601 (2021) doi:10.1007/s11467-021-1065-x [arXiv:2102.07830 [hep-ph]].
  * (39) C. Lorcé, A. Metz, B. Pasquini and S. Rodini, [arXiv:2109.11785 [hep-ph]].
  * (40) B. Gittelman, K. M. Hanson, D. Larson, E. Loh, A. Silverman and G. Theodosiou, Phys. Rev. Lett. 35, 1616 (1975) doi:10.1103/PhysRevLett.35.1616
  * (41) U. Camerini, J. G. Learned, R. Prepost, C. M. Spencer, D. E. Wiser, W. Ash, R. L. Anderson, D. Ritson, D. Sherden and C. K. Sinclair, Phys. Rev. Lett. 35, 483 (1975) doi:10.1103/PhysRevLett.35.483
  * (42) S. Joosten and Z. E. Meziani, PoS QCDEV2017, 017 (2018) doi:10.22323/1.308.0017 [arXiv:1802.02616 [hep-ex]].
  * (43) J. Dudek, R. Ent, R. Essig, K. S. Kumar, C. Meyer, R. D. McKeown, Z. E. Meziani, G. A. Miller, M. Pennington and D. Richards, et al. Eur. Phys. J. A 48, 187 (2012) doi:10.1140/epja/i2012-12187-1 [arXiv:1208.1244 [hep-ex]].
  * (44) J. P. Chen et al. [SoLID], [arXiv:1409.7741 [nucl-ex]].
  * (45) A. Accardi, J. L. Albacete, M. Anselmino, N. Armesto, E. C. Aschenauer, A. Bacchetta, D. Boer, W. K. Brooks, T. Burton and N. B. Chang, et al. Eur. Phys. J. A 52, no.9, 268 (2016) doi:10.1140/epja/i2016-16268-9 [arXiv:1212.1701 [nucl-ex]].
  * (46) R. Abdul Khalek, A. Accardi, J. Adam, D. Adamiak, W. Akers, M. Albaladejo, A. Al-bataineh, M. G. Alexeev, F. Ameli and P. Antonioli, et al. [arXiv:2103.05419 [physics.ins-det]].
  * (47) D. P. Anderle, V. Bertone, X. Cao, L. Chang, N. Chang, G. Chen, X. Chen, Z. Chen, Z. Cui and L. Dai, et al. Front. Phys. (Beijing) 16, no.6, 64701 (2021) doi:10.1007/s11467-021-1062-0 [arXiv:2102.09222 [nucl-ex]].
  * (48) G. P. Lepage and S. J. Brodsky, Phys. Rev. Lett. 43, 545-549 (1979) [erratum: Phys. Rev. Lett. 43, 1625-1626 (1979)] doi:10.1103/PhysRevLett.43.545
  * (49) S. J. Brodsky and G. P. Lepage, Phys. Rev. D 24, 2848 (1981) doi:10.1103/PhysRevD.24.2848
  * (50) A. V. Efremov and A. V. Radyushkin, Phys. Lett. B 94, 245-250 (1980) doi:10.1016/0370-2693(80)90869-2
  * (51) V. L. Chernyak and A. R. Zhitnitsky, JETP Lett. 25, 510 (1977)
  * (52) V. L. Chernyak and A. R. Zhitnitsky, Sov. J. Nucl. Phys. 31, 544-552 (1980)
  * (53) V. L. Chernyak and A. R. Zhitnitsky, Phys. Rept. 112, 173 (1984) doi:10.1016/0370-1573(84)90126-1
  * (54) A. V. Belitsky, X. d. Ji and F. Yuan, Phys. Rev. Lett. 91, 092003 (2003) doi:10.1103/PhysRevLett.91.092003 [arXiv:hep-ph/0212351 [hep-ph]].
  * (55) X. B. Tong, J. P. Ma and F. Yuan, Phys. Lett. B 823, 136751 (2021) doi:10.1016/j.physletb.2021.136751 [arXiv:2101.02395 [hep-ph]].
  * (56) G. T. Bodwin, E. Braaten and G. P. Lepage, Phys. Rev. D 51, 1125-1171 (1995) [erratum: Phys. Rev. D 55, 5853 (1997)] doi:10.1103/PhysRevD.55.5853 [arXiv:hep-ph/9407339 [hep-ph]].
  * (57) S. J. Brodsky and G. R. Farrar, Phys. Rev. Lett. 31, 1153-1156 (1973) doi:10.1103/PhysRevLett.31.1153
  * (58) V. A. Matveev, R. M. Muradian and A. N. Tavkhelidze, Lett. Nuovo Cim. 7, 719-723 (1973) doi:10.1007/BF02728133
  * (59) X. d. Ji, J. P. Ma and F. Yuan, Eur. Phys. J. C 33, 75-90 (2004) doi:10.1140/epjc/s2003-01563-y [arXiv:hep-ph/0304107 [hep-ph]].
  * (60) L. Frankfurt and M. Strikman, Phys. Rev. D 66, 031502 (2002) doi:10.1103/PhysRevD.66.031502 [arXiv:hep-ph/0205223 [hep-ph]].
  * (61) V. M. Braun, S. E. Derkachov, G. P. Korchemsky and A. N. Manashov, Nucl. Phys. B 553, 355-426 (1999) doi:10.1016/S0550-3213(99)00265-5 [arXiv:hep-ph/9902375 [hep-ph]].
  * (62) V. Braun, R. J. Fries, N. Mahnke and E. Stein, Nucl. Phys. B 589, 381-409 (2000) [erratum: Nucl. Phys. B 607, 433-433 (2001)] doi:10.1016/S0550-3213(00)00516-2 [arXiv:hep-ph/0007279 [hep-ph]].
  * (63) X. d. Ji, J. P. Ma and F. Yuan, Nucl. Phys. B 652, 383-404 (2003) doi:10.1016/S0550-3213(03)00010-5 [arXiv:hep-ph/0210430 [hep-ph]].
  * (64) S. J. Brodsky, E. Chudakov, P. Hoyer and J. M. Laget, Phys. Lett. B 498, 23-28 (2001) doi:10.1016/S0370-2693(00)01373-3 [arXiv:hep-ph/0010343 [hep-ph]].
  * (65) P. Sun, X. B. Tong and F. Yuan, Phys. Lett. B 822, 136655 (2021) doi:10.1016/j.physletb.2021.136655 [arXiv:2103.12047 [hep-ph]].
  * (66) P. Hoodbhoy, X. d. Ji and F. Yuan, Phys. Rev. Lett. 92, 012003 (2004) doi:10.1103/PhysRevLett.92.012003 [arXiv:hep-ph/0309085 [hep-ph]].
  * (67) X. d. Ji, J. P. Ma and F. Yuan, Phys. Rev. Lett. 90, 241601 (2003) doi:10.1103/PhysRevLett.90.241601 [arXiv:hep-ph/0301141 [hep-ph]].
  * (68) G. P. Lepage and S. J. Brodsky, Phys. Rev. D 22, 2157 (1980) doi:10.1103/PhysRevD.22.2157
  * (69) B. L. Ioffe, Nucl. Phys. B 188, 317-341 (1981) [erratum: Nucl. Phys. B 191, 591-592 (1981)] doi:10.1016/0550-3213(81)90259-5
  * (70) Y. Chung, H. G. Dosch, M. Kremer and D. Schall, Nucl. Phys. B 197, 55-75 (1982) doi:10.1016/0550-3213(82)90154-7
  * (71) V. L. Chernyak and I. R. Zhitnitsky, Nucl. Phys. B 246, 52-74 (1984) doi:10.1016/0550-3213(84)90114-7
  * (72) I. D. King and C. T. Sachrajda, Nucl. Phys. B 279, 785-803 (1987) doi:10.1016/0550-3213(87)90019-8
  * (73) V. L. Chernyak, A. A. Ogloblin and I. R. Zhitnitsky, Sov. J. Nucl. Phys. 48, 536 (1988) IYF-87-135.
  * (74) V. L. Chernyak, A. A. Ogloblin and I. R. Zhitnitsky, Yad. Fiz. 48, 1410-1422 (1988) doi:10.1007/BF01557663
  * (75) J. Bolz and P. Kroll, Z. Phys. A 356, 327 (1996) doi:10.1007/s002180050186 [arXiv:hep-ph/9603289 [hep-ph]].
  * (76) V. M. Braun, A. Lenz, N. Mahnke and E. Stein, Phys. Rev. D 65, 074011 (2002) doi:10.1103/PhysRevD.65.074011 [arXiv:hep-ph/0112085 [hep-ph]].
  * (77) V. M. Braun, A. Lenz and M. Wittmann, Phys. Rev. D 73, 094019 (2006) doi:10.1103/PhysRevD.73.094019 [arXiv:hep-ph/0604050 [hep-ph]].
  * (78) M. Gockeler, R. Horsley, T. Kaltenbrunner, Y. Nakamura, D. Pleiter, P. E. L. Rakow, A. Schafer, G. Schierholz, H. Stuben and N. Warkentin, et al. Phys. Rev. Lett. 101, 112002 (2008) doi:10.1103/PhysRevLett.101.112002 [arXiv:0804.1877 [hep-lat]].
  * (79) V. M. Braun, A. N. Manashov and J. Rohrwild, Nucl. Phys. B 807, 89-137 (2009) doi:10.1016/j.nuclphysb.2008.08.012 [arXiv:0806.2531 [hep-ph]].
  * (80) V. M. Braun et al. [QCDSF], Phys. Rev. D 79, 034504 (2009) doi:10.1103/PhysRevD.79.034504 [arXiv:0811.2712 [hep-lat]].
  * (81) K. Passek-Kumericki and G. Peters, Phys. Rev. D 78, 033009 (2008) doi:10.1103/PhysRevD.78.033009 [arXiv:0805.1758 [hep-ph]].
  * (82) A. Lenz, M. Gockeler, T. Kaltenbrunner and N. Warkentin, Phys. Rev. D 79, 093007 (2009) doi:10.1103/PhysRevD.79.093007 [arXiv:0903.1723 [hep-ph]].
  * (83) I. V. Anikin, V. M. Braun and N. Offen, Phys. Rev. D 88, 114021 (2013) doi:10.1103/PhysRevD.88.114021 [arXiv:1310.1375 [hep-ph]].
  * (84) V. M. Braun, S. Collins, B. Gläßle, M. Göckeler, A. Schäfer, R. W. Schiel, W. Söldner, A. Sternbeck and P. Wein, Phys. Rev. D 89, 094511 (2014) doi:10.1103/PhysRevD.89.094511 [arXiv:1403.4189 [hep-lat]].
  * (85) G. S. Bali, V. M. Braun, M. Göckeler, M. Gruber, F. Hutzler, A. Schäfer, R. W. Schiel, J. Simeth, W. Söldner and A. Sternbeck, et al. JHEP 02, 070 (2016) doi:10.1007/JHEP02(2016)070 [arXiv:1512.02050 [hep-lat]].
  * (86) G. S. Bali et al. [RQCD], Eur. Phys. J. A 55, no.7, 116 (2019) doi:10.1140/epja/i2019-12803-6 [arXiv:1903.12590 [hep-lat]].
  * (87) Y. Feng and J. X. Wang, Adv. High Energy Phys. 2015, 726393 (2015) doi:10.1155/2015/726393 [arXiv:1510.05277 [hep-ph]].
  * (88) X. J. Zhan and J. X. Wang, Chin. Phys. C 45, no.2, 023112 (2021) doi:10.1088/1674-1137/abce11