# Photon-Odderon interference in exclusive $\chi_{c}$ charmonium production at
the Electron-Ion Collider

Sanjin Benić  Department of Physics, Faculty of Science, University of Zagreb,
Bijenička c. 32, 10000 Zagreb, Croatia  Adrian Dumitru  Department of Natural
Sciences, Baruch College, CUNY, 17 Lexington Avenue, New York, NY 10010, USA
The Graduate School and University Center, The City University of New York,
365 Fifth Avenue, New York, NY 10016, USA  Abhiram Kaushik  Centre for
Informatics and Computing, Rudjer Bošković Institute, HR-10002 Zagreb, Croatia
Leszek Motyka  Jagiellonian University, Institute of Theoretical Physics,
Łojasiewicza 11, 30-348 Kraków, Poland  Tomasz Stebel  Jagiellonian
University, Institute of Theoretical Physics, Łojasiewicza 11, 30-348 Kraków,
Poland

###### Abstract

Exclusive $C=+1$ scalar, axial-vector, and tensor quarkonium production in
high-energy electron-proton scattering requires a $C$-odd $t$-channel exchange
of a photon or a three gluon ladder. We derive the expressions for the
corresponding amplitudes. The relative phase of the photon vs. three gluon
exchange amplitudes is determined by the sign of the light-front matrix
element of the eikonal color current operator $d^{abc}J^{+a}J^{+b}J^{+c}$ at
moderate $x$, and is not affected by small-$x$ QCD evolution. Model
calculations predict constructive interference, which is particularly strong
for momentum transfer $|t|\sim 1$ GeV2 where the cross section for $\chi_{cJ}$
production exceeds that for pure photon exchange by up to a factor of 4. We
find that exclusive $\chi_{cJ}$ electroproduction at the Electron-Ion Collider
should occur with well measurable rates and measurements of these processes
should allow to find an evidence of the perturbative Odderon exchange. We also
compute the total electroproduction cross section as a function of energy and
provide first estimates of the number of $\chi_{cJ}$ events per month at the
Electron-Ion Collider design luminosity.

††preprint: ZTF-EP-24-02

###### Contents

  1. I Introduction
  2. II The production amplitude of $C$-even quarkonia
    1. II.1 Light-cone wave function of $C$-even quarkonia
    2. II.2 Final expressions for the amplitudes
    3. II.3 The $\gamma^{*}p\to\mathcal{H}p$ cross section
    4. II.4 Odderon amplitude and its evolution with $x$
  3. III The Primakoff contribution
    1. III.1 $|t|\to 0$ limit
    2. III.2 Adding the Pauli form factor
  4. IV Boosted Gaussian model for the $\chi_{cJ}$ wave functions
  5. V Numerical results
  6. VI Summary and Conclusions
  7. A Computation of the light-cone wave functions
  8. B Derivation of the amplitudes for axial and tensor quarkonia
  9. C The Primakoff contribution in specific kinematic limits
    1. C.1 The proof of formula (40)
    2. C.2 The $|t|\to 0$ limit for axial vector quarkonia and its connection to the Landau-Yang theorem
    3. C.3 The NRQCD limit of the Primakoff cross sections
  10. D The relative sign of the eikonal photon and Odderon exchange amplitudes
  11. E BJKP–BLV Odderon

##  I Introduction

A basic prediction of QCD, related to the non-zero cubic Casimir invariant, is
that there exists a contribution to high-energy scattering amplitudes from a
$C$-conjugation odd exchange driven mostly by gluons, with a weakly energy
dependent cross section. The possible existence of such an exchange, called
the Odderon, was predicted 50 years ago from general principles of quantum
field theory [1]. If the scattering involves hadrons of size much smaller than
the QCD confinement scale, or high momentum transfer, then the composition of
the exchanged “state” can be understood from perturbative QCD. At leading
order in perturbation theory, it corresponds to an exchange of three gluons
which is symmetric under permutations of their colors. At high energies the
three gluon exchange is dressed by higher order QCD corrections enhanced by
logarithms of energy that have to be resummed into a gluon ladder. It accounts
for interactions between the exchanged $t$-channel gluons, and for virtual
corrections. The resulting “reggeized” three gluon exchange is the hard
Odderon. In the weak field limit the corresponding Bartels-Jaroszewicz-
Kwieciński-Praszałowicz (BJKP) linear QCD evolution equation has been obtained
long ago [2, 3, 4], and it amounts to the reggeization of the exchanged gluons
by iteration of the Balitsky-Fadin-Kuraev-Lipatov (BFKL) interaction kernel
[5, 6, 7] for each pair of gluons in the ladder. The NLO corrections to the
BJKP equation were found in Ref. [8]. At asymptotically high energies the
leading solution to the BJKP equation corresponds to a configuration of gluons
in which two reggeized gluons are combined together [9] – the
Bartels–Lipatov–Vacca (BLV) solution. The BLV solution couples to color
dipoles so it is particularly important in $C$-even quarkonia (or meson)
photo- and electroproduction [9]. The intercept of this solution equals one,
in the leading logarithmic approximation [9], and it was argued that it
remains at this value to all orders in the perturbative expansion [10].

The non-linear Odderon evolution equation for dipole–proton scattering was
established in Refs. [11, 12]. This equation typically produces solutions
where the Odderon amplitude decreases with energy [13, 14, 15, 16], as will be
confirmed also in our present work. The seed for this evolution to small $x$
is provided by the cubic Casimir $d^{abc}J^{+a}J^{+b}J^{+c}$ in the effective
action describing color current ($J^{+a}$) fluctuations in the proton at
moderately small $x$ [17, 18], which allows the $C$-odd gluon ladder to couple
to the proton. The matrix element111We are here concerned with the off-forward
matrix element for non-zero momentum transfer $t$. In the $t\to 0$ limit,
instead, the “spin dependent Odderon” [19, 20, 21] is associated with a spin
flip of the proton in DIS. of this operator could also be evaluated directly
from the light-cone wave function of the proton [22, 23]. A key point for the
present analysis is that the Odderon evolution equation to small $x$ does not
alter the sign of the amplitude. Hence, the interference pattern with the
electromagnetic amplitude due to single photon exchange is determined by the
sign of the matrix element of the $d^{abc}J^{+a}J^{+b}J^{+c}$ operator at the
initial, moderately small $x$.

The TOTEM collaboration at the CERN-LHC has measured the differential cross
section for $pp$ elastic scattering at $\sqrt{s}=2.76$ GeV [24]. They observe
a significant difference to the cross section for $p\bar{p}$ scattering at
$\sqrt{s}=1.96$ GeV measured by D0 [25]. Assuming that the difference in
energy is negligible they conclude that these results provide evidence for the
exchange of a color singlet gluonic state with odd C-parity, the Odderon.
However, these measurements involve the scattering of hadrons with a size of
order the QCD confinement scale as well as low momentum transfers, $|t|\ll 1$
GeV2. Hence, an analysis of the nature of the $t$-channel exchange from
perturbative QCD would not appear reliable.

Exclusive production of pseudo-scalar $\eta_{c}$ quarkonia in DIS has
previously been proposed as a process suitable for the discovery of the hard
Odderon exchange in QCD [26, 27, 28, 29, 30, 31, 32, 16]. In practice, the
detection of this process is difficult as the $\eta_{c}$ has small branching
ratios to all relevant decay channels. For example, the radiative
$\eta_{c}\to\gamma\gamma$ mode has a tiny ($\sim 10^{-4}$) branching
ratio222We take all branching ratios and particle masses from the “Review of
Particle Physics” [33]., while the hadronic $\eta_{c}\to\rho\rho$ mode has a
branching ratio of $1.5\%$. Also, the cross section for exclusive production
of $J/\psi$ is expected to be far greater than that for $\eta_{c}$, and so
most $\eta_{c}$’s would originate from a $J/\psi\to\eta_{c}\gamma$
decay333This issue could be mitigated by considering $\eta_{c}(2S)$ instead,
since it lies between the $J/\psi$ and the $\psi(2S)$. The branching ratio
$\psi(2S)\to\eta_{c}+\gamma$ is favourably small ($(7\pm 5)\times 10^{-4}$)
but the detection of the $\eta_{c}(2S)$ is again problematic due to the small
branching ratios for many channels.. The experimental identification of the
soft photon ($M_{J/\psi}=3.097$ GeV, $M_{\eta_{c}}=2.984$ GeV) is difficult
[34, 35].

Alternative approaches for the discovery of the hard Odderon include
Pomeron–Odderon interference (with a background from Pomeron–photon
interference) which would manifest in asymmetries in exclusive production of
two charged particles in DIS [36, 37, 38]; for the $\gamma^{*}p\to\chi_{c}p$
process considered here, it is instead the interference of the Primakoff-like
amplitude with the QCD Odderon exchange amplitude which will be important,
i.e. photon–Odderon interference. Another process which involves a $C$-odd
$t$-channel gluonic exchange is exclusive production of a $J/\psi$ in double-
diffractive proton-proton scattering [39, 40]. At high-energy proton colliders
this requires instrumentation over a large range of rapidity, close to the
beams.

The GlueX Collaboration at Jefferson Lab has recently reported the observation
of $56.5\pm 8.2$ exclusive $\chi_{c1}(1P)$ and $12.7\pm 4.5$ $\chi_{c2}(1P)$
production events in DIS near threshold energy [41]. They noticed a “dramatic
difference” in the momentum transfer $t$-distribution of events as compared to
$J/\psi$ production in that the cross section for $\chi_{c}$ production
appears to drop off much less rapidly with increasing $|t|$. In other words,
the probability that the struck proton remains intact at high momentum
transfer is much greater in $\gamma^{*}p\to\chi_{c}p$ than in $\gamma^{*}p\to
J/\psi p$ events. This remarkable result illustrates the importance of the
underlying QCD dynamics as opposed to naive expectations based solely on the
mass of the produced quarkonia. A much harder $t$-dependence for exclusive
production of $\eta_{c}$ vs. $J/\psi$ in eikonal dipole–proton scattering is
predicted by simple light-front constituent quark models [32]. Inspired by the
GlueX measurement of exclusive $\chi_{c}$ production in DIS at Jefferson Lab,
here we consider the same process but at high energies appropriate for the
future Electron-Ion Collider (EIC) [42, 43] where dipole model factorization
[44, 45] applies, and where the $t$-channel exchange of a $C$-odd color
singlet state dominates.

Within the $\chi_{cJ}$ family ($J=0$, $1$, $2$), $\chi_{c1}$ has the largest
branching ratio for the radiative $\chi_{c1}\to J/\psi\gamma$ decay ($34.3\%$
as opposed to $1.4\%$ for $\chi_{c0}$ and $19.0\%$ for $\chi_{c2}$). Indeed,
this was the detection channel used by the GlueX measurement. Therefore, the
production of $\chi_{c1}$ axial-vector and $\chi_{c2}$ tensor quarkonia may
prove the most promising channels for discovery of the hard Odderon. This is
further corroborated by the fact that for $\chi_{c1}$ and $\chi_{c2}$ the
Odderon and photon exchange contributions become comparable at lower $|t|$
than for $\chi_{c0}$, within the acceptance of the proton spectrometer of the
EIC design detector [43]. However, for the $\chi_{cJ}$ quarkonia there is
again a feed-down channel from exclusive production of $\psi(2S)$ with
subsequent decay $\psi(2S)\to\chi_{cJ}\gamma$. Hence, the identification and
rejection of such feed down is required.

The paper is organized as follows. In Sec. II we start with the computation of
the amplitude and the corresponding $\gamma^{*}p$ cross section for exclusive
production of $C$-even scalar, axial vector and tensor quarkonia. We derive
the light-cone wave functions of quarkonia, whereas the amplitude is obtained
as an overlap with the photon light-cone wave function. We introduce the
Odderon exchange amplitude and discuss its evolution to small-$x$. Sec. III is
devoted to the computation of the Primakoff contribution, where we pay special
attention to the $t\to 0$ limit. In Sec. IV we perform a numerical fit of the
$\chi_{cJ}$ wave functions. The main results are shown in Sec. V, where we
numerically compute exclusive $\gamma^{*}p\to\chi_{cJ}p$ and
$ep\to\chi_{cJ}ep$ cross sections and the expected number of events at the
EIC. Our findings are summarized in the final Sec. VI. Several Appendices
follow where we explain the computational steps leading to some of the results
from the main text.

##  II The production amplitude of $C$-even quarkonia

We consider the process $\gamma^{*}(q)p(P)\to
p(P^{\prime})\mathcal{H}(\Delta)$ where $\mathcal{H}$ is the $C$-even
quarkonium state. In this paper, we specifically have in mind the case when
$\mathcal{H}$ are the $P$-wave charmonia $\chi_{cJ}$ states with $J=0$, $1$,
$2$. However our derivations are presented in a general form so that they
apply also to other $C=+1$ $q\bar{q}$ bound states such as bottomonia, or even
mesons with a light flavor content.

Our computation is performed in the dipole frame where
$q=(-Q^{2}/2q^{-},q^{-},0,0)$, $P=(P^{+},0,0,0)$ and the proton mass
$M_{p}\ll\sqrt{(q+P)^{2}}$ is neglected. We use the convention for the
components $(v^{+},v^{-},v^{1},v^{2})$ of a four vector $v$ where
$v^{\pm}=(v^{0}\pm v^{3})/\sqrt{2}$ are the $\pm$ light-cone components. The
amplitude is computed using light-cone gauge $A^{-}=A\cdot n=0$, that is,
$n^{\mu}=\delta^{\mu+}$. Following a similar procedure as in [32, 16] our
starting point is the amplitude

$\left\langle\mathcal{M}_{\lambda\bar{\lambda}}(\gamma^{*}p\to\mathcal{H}p)\right\rangle=2q^{-}N_{c}\int_{\boldsymbol{r}_{\perp}\boldsymbol{b}_{\perp}}\mathrm{e}^{-\mathrm{i}\boldsymbol{\Delta}_{\perp}\cdot\boldsymbol{b}_{\perp}}\mathrm{i}\mathcal{O}(\boldsymbol{r}_{\perp},\boldsymbol{b}_{\perp})\mathcal{A}_{\lambda\bar{\lambda}}(\boldsymbol{r}_{\perp},\boldsymbol{\Delta}_{\perp})\,,$
(1)

with $\lambda$ and $\bar{\lambda}$ the helicities of the incoming photon and
the outgoing quarkonium, respectively. For brevity, we will often times write
the amplitude simply as $\langle\mathcal{M}_{\lambda\bar{\lambda}}\rangle$
(omitting the round brackets). Here
$\mathcal{A}_{\lambda\bar{\lambda}}(\boldsymbol{r}_{\perp},\boldsymbol{\Delta}_{\perp})$
is the reduced amplitude, described physically in terms of the photon
($\Psi^{\gamma}_{\lambda,h\bar{h}}$) and the quarkonia
($\Psi^{\mathcal{H}}_{\lambda,h\bar{h}}$) wave function overlap. We have444We
have the following abbreviation for the transverse coordinate space integrals:
$\int_{\boldsymbol{r}_{\perp}}\equiv\int\mathrm{d}^{2}\boldsymbol{r}_{\perp}$.
For transverse and longitudinal momentum space integrals we use:
$\int_{\boldsymbol{l}_{\perp}}\equiv\int\mathrm{d}^{2}\boldsymbol{l}_{\perp}/(2\pi)^{2}$
and $\int_{z}\equiv\int_{0}^{1}\mathrm{d}z/(4\pi)$.

$\mathcal{A}_{\lambda\bar{\lambda}}(\boldsymbol{r}_{\perp},\boldsymbol{\Delta}_{\perp})=\int_{z}\int_{\boldsymbol{l}_{\perp}\boldsymbol{l}_{\perp}^{\prime}}\sum_{h\bar{h}}\Psi^{\gamma}_{\lambda,h\bar{h}}(\boldsymbol{l}_{\perp},z)\Psi^{\mathcal{H}*}_{\bar{\lambda},h\bar{h}}(\boldsymbol{l}_{\perp}^{\prime}-z\boldsymbol{\Delta}_{\perp},z)\mathrm{e}^{\mathrm{i}(\boldsymbol{l}_{\perp}-\boldsymbol{l}_{\perp}^{\prime}+\frac{1}{2}\boldsymbol{\Delta}_{\perp})\cdot\boldsymbol{r}_{\perp}}\,.$
(2)

A particular contribution to the amplitude (1) is shown in Fig. 1, where the
upper heavy quark loop represents the reduced amplitude
$\mathcal{A}_{\lambda\bar{\lambda}}(\boldsymbol{r}_{\perp},\boldsymbol{\Delta}_{\perp})$,
and the three-gluon contribution is the lowest order depiction of the Odderon
exchange amplitude
$\mathcal{O}(\boldsymbol{r}_{\perp},\boldsymbol{b}_{\perp})$ for the
particular case where the gluons connect to three different quarks in the
proton wave function, represented by the green blob.

The photon wave-function in momentum space is given by

$\Psi^{\gamma}_{\lambda,h\bar{h}}(\boldsymbol{k}_{\perp},z)\equiv\sqrt{z\bar{z}}\,\frac{\bar{u}_{h}(k)eq_{c}\not{\epsilon}(\lambda,q)v_{\bar{h}}(q-k)}{\boldsymbol{k}_{\perp}^{2}+\varepsilon^{2}}\,,$
(3)

where $q_{c}=+2/3$ is the fractional electric charge of the $c$-quark and
$u_{h}(k)$ ($v_{\bar{h}}(q-k)$) are particle (antiparticle) spinors with
$h=\pm 1$ indicating the sign of the helicity $h/2$. For explicit computations
we use Lepage-Brodsky (LB) spinors [46] \- see App. A for their expressions.
$\epsilon(\lambda,q)$ is the photon polarization vector: we have
$\epsilon(0,q)=(Q/q^{-},0,0,0)$ for longitudinal ($\lambda=0$) and
$\epsilon(\lambda,q)=(0,0,\boldsymbol{\epsilon}_{\perp}^{\lambda})$ for
transverse ($\lambda=\pm 1$) polarization. We are using the LB convention for
the 2D polarization vector:
$\boldsymbol{\epsilon}_{\perp}^{\lambda}\equiv(-\lambda,-\mathrm{i})/\sqrt{2}$,
and the common shorthands $\varepsilon\equiv\sqrt{m_{c}^{2}+z\bar{z}Q^{2}}$
and $\bar{z}\equiv 1-z$, with $z=k^{-}/q^{-}$. We compute the photon wave
function as described in App. A, with the result

$\begin{split}&\Psi^{\gamma}_{+,h\bar{h}}(\boldsymbol{k}_{\perp},z)=-\sqrt{2}eq_{c}\left[k_{\perp}\mathrm{e}^{\mathrm{i}\phi_{k}}\left(z\delta_{h+}\delta_{\bar{h}-}-\bar{z}\delta_{h-}\delta_{\bar{h}+}\right)+m_{c}\delta_{h+}\delta_{\bar{h}+}\right]\frac{1}{\boldsymbol{k}_{\perp}^{2}+\varepsilon^{2}}\,,\\\
&\Psi^{\gamma}_{-,h\bar{h}}(\boldsymbol{k}_{\perp},z)=-\sqrt{2}eq_{c}\left[k_{\perp}\mathrm{e}^{-\mathrm{i}\phi_{k}}\left(\bar{z}\delta_{h+}\delta_{\bar{h}-}-z\delta_{h-}\delta_{\bar{h}+}\right)+m_{c}\delta_{h-}\delta_{\bar{h}-}\right]\frac{1}{\boldsymbol{k}_{\perp}^{2}+\varepsilon^{2}}\,,\\\
&\Psi^{\gamma}_{0,h\bar{h}}(\boldsymbol{k}_{\perp},z)=\frac{eq_{c}2Qz\bar{z}}{\boldsymbol{k}_{\perp}^{2}+\varepsilon^{2}}\delta_{h,-\bar{h}}\,,\end{split}$
(4)

where $k_{\perp}\mathrm{e}^{\pm\mathrm{i}\phi_{k}}=k^{1}\pm\mathrm{i}k^{2}$.
Eq. (4) coincides with the result of Ref. [47] for both $\lambda=0$ and
$\lambda=\pm 1$ after adjusting for an overall factor of $-\sqrt{N_{c}}$ due
to the difference in conventions. Comparing to Ref. [48] the result coincides
for $\lambda=0$ (again up to a factor $-\sqrt{N_{c}}$) but there is an
opposite overall sign for $\lambda=\pm 1$ (the relative factor being
$+\sqrt{N_{c}}$), as also explicitly noted in [48]. While such overall signs
play no role for the processes involving overlaps of photons and vector
quarkonia, the sign is important for tensor quarkonia, since the tensor wave
function can be understood as a linear superposition of the vector wave
functions, see Eq. (12) below. The computation described in App. A takes the
LB spinors as a starting point from which we proceed to explicitly compute the
photon wave function as well as all the subsequent quarkonia wave functions in
a systematic and self-contained fashion.

Figure 1: An example Feynman diagram illustrating the amplitude for exclusive
production of $\chi_{cJ}$’s via three gluon (Odderon) exchange from the proton
in $\gamma^{*}p$ collision. The green blobs represent the proton wave function
that is shown here in terms of three valence quarks. The magenta blob is the
$\chi_{cJ}$ wave function. In the dipole approach the two gluons coupling to
the same quark line in the top part of the diagram occupy the same point in
transverse position space.

The explicit expressions for the quarkonia wave functions will be given in the
following subsection. Here we want to emphasize that thanks to the $C=+1$
parity of the quarkonium wave function (see Eq. (10)) below) and to the $C=-1$
parity of the photon, the amplitude is strictly proportional to the $C=-1$
Odderon amplitude
$\mathcal{O}(\boldsymbol{r}_{\perp},\boldsymbol{b}_{\perp})$. In the high
energy limit of eikonal dipole-proton scattering we have [11, 12]

$\mathcal{O}(\boldsymbol{r}_{\perp},\boldsymbol{b}_{\perp})=-\frac{1}{2\mathrm{i}N_{c}}{\rm
tr}\left\langle
V(\boldsymbol{x}_{\perp})V^{\dagger}(\boldsymbol{y}_{\perp})-V(\boldsymbol{y}_{\perp})V^{\dagger}(\boldsymbol{x}_{\perp})\right\rangle\,,$
(5)

where $\boldsymbol{r}_{\perp}$ is the dipole size and $\boldsymbol{b}_{\perp}$
the impact parameter

$\boldsymbol{r}_{\perp}\equiv\boldsymbol{x}_{\perp}-\boldsymbol{y}_{\perp}\,,\qquad\boldsymbol{b}_{\perp}\equiv\frac{\boldsymbol{x}_{\perp}+\boldsymbol{y}_{\perp}}{2}\,,$
(6)

respectively 555The impact parameter of the dipole-proton collision is
actually
$\tilde{\boldsymbol{b}}_{\perp}=z\boldsymbol{x}_{\perp}+\bar{z}\boldsymbol{y}_{\perp}=\boldsymbol{b}_{\perp}-(\frac{1}{2}-z)\boldsymbol{r}_{\perp}$,
and $\tilde{\boldsymbol{b}}_{\perp}$ is the Fourier conjugate of the
transverse momentum transfer $\boldsymbol{\Delta}_{\perp}$. This is the origin
of the “off-forward phase”
$-i\boldsymbol{\delta}_{\perp}\cdot\boldsymbol{r}_{\perp}$ in Eq. (18)..
$V(\boldsymbol{x}_{\perp})$ and $V^{\dagger}(\boldsymbol{x}_{\perp})$
represent Wilson lines at transverse coordinate $\boldsymbol{x}_{\perp}$ which
describe the propagation of the quark and anti-quark, respectively, through
the color field of the target proton in covariant gauge, and
$\langle\cdots\rangle$ denotes an average over the configurations of that
field, see App. D for the conventions used in this work. In Fourier space, the
Odderon amplitude is given by

$\mathcal{O}(\boldsymbol{r}_{\perp},\boldsymbol{\Delta}_{\perp})\equiv\int_{\boldsymbol{b}_{\perp}}\mathrm{e}^{-\mathrm{i}\boldsymbol{\Delta}_{\perp}\cdot\boldsymbol{b}_{\perp}}\mathcal{O}(\boldsymbol{r}_{\perp},\boldsymbol{b}_{\perp})\,.$
(7)

The $C$-invariance of the full amplitude (1) is now easily verified. We can
conventionally start by first exchanging the quark-antiquark coordinates
$\boldsymbol{x}_{\perp}\leftrightarrow\boldsymbol{y}_{\perp}$ so that
$\boldsymbol{r}_{\perp}\to-\boldsymbol{r}_{\perp}$ in the Odderon amplitude.
Because of the phases in (2) this needs to be followed by an exchange of quark
and antiquark transverse momenta
$\boldsymbol{l}_{\perp}\leftrightarrow\boldsymbol{l}_{\perp}^{\prime}-z\boldsymbol{\Delta}_{\perp}$
and $z\leftrightarrow\bar{z}$. By this transformation only the photon wave
function picks up a sign which cancels with the sign in the transformation of
the Odderon amplitude as
$\mathcal{O}(-\boldsymbol{r}_{\perp},\boldsymbol{b}_{\perp})=-\mathcal{O}(\boldsymbol{r}_{\perp},\boldsymbol{b}_{\perp})$.

###  II.1 Light-cone wave function of $C$-even quarkonia

The light-cone wave function of the $C$-even quarkonia,
$\Psi^{\mathcal{H}}_{\bar{\lambda},h\bar{h}}$, is modeled by the following
covariant ansatz

$\Psi^{\mathcal{H}}_{\bar{\lambda},h\bar{h}}(\boldsymbol{k}_{\perp},z)\equiv\frac{1}{\sqrt{z\bar{z}}}\bar{u}_{h}(k)\Gamma_{\bar{\lambda}}^{\mathcal{H}}(k,k^{\prime})v_{\bar{h}}(k^{\prime})\phi_{\mathcal{H}}(\boldsymbol{k}_{\perp},z)\,.$
(8)

Here $\bar{\lambda}$ is the quarkonia helicity and $k$
($k^{\prime}=\Delta_{0}-k$) is the quark (antiquark) momenta, with
$\Delta_{0}$ being the invariant four-momentum of the $c\bar{c}$ pair.
$\Gamma^{\mathcal{H}}_{\bar{\lambda}}$ is the appropriate Dirac matrix vertex
function for either a scalar ($\mathcal{H}=\mathcal{S}$), axial vector
($\mathcal{A}$) or a tensor ($\mathcal{T}$) quarkonia which we take as666Our
convention is
$\gamma^{5}\equiv+\mathrm{i}\gamma^{0}\gamma^{1}\gamma^{2}\gamma^{3}$ and
$\epsilon^{0123}\equiv+1=-\epsilon^{+-12}$.

$\Gamma^{\mathcal{H}}_{\bar{\lambda}}(k,k^{\prime})=\begin{cases}1~{},&\mathcal{H}=\mathcal{S}\,,\\\
\mathrm{i}\gamma_{5}\not{E}(\bar{\lambda},\Delta_{0})~{},&\mathcal{H}=\mathcal{A}\,,\\\
\frac{1}{2}\left(\gamma_{\mu}(k_{\nu}-k^{\prime}_{\nu})+\gamma_{\nu}(k_{\mu}-k^{\prime}_{\mu})\right)E^{\mu\nu}(\bar{\lambda},\Delta_{0})~{},&\mathcal{H}=\mathcal{T}\,.\end{cases}$
(9)

$\phi_{\mathcal{H}}(\boldsymbol{k}_{\perp},z)$ is a non-perturbative scalar
function that we model later (see Sec. IV).

The structure of (9) is motivated in part to ensure the correct $C$-even
property of the quarkonia wave function in (8). Namely, by exchanging the
quark and the antiquark momenta ($k\leftrightarrow k^{\prime}$) and helicity
($h\leftrightarrow\bar{h}$) we have

$\bar{u}_{\bar{h}}(k^{\prime})\Gamma_{\bar{\lambda}}^{\mathcal{H}}(k^{\prime},k)v_{h}(k)=+\bar{u}_{h}(k)\Gamma^{\mathcal{H}}_{\bar{\lambda}}(k,k^{\prime})v_{\bar{h}}(k^{\prime})\,,$
(10)

which holds thanks to the relation
$C(\Gamma^{\mathcal{H}}_{\bar{\lambda}}(k^{\prime},k))^{T}C^{-1}=+\Gamma_{\bar{\lambda}}^{\mathcal{H}}(k,k^{\prime})$
with $C=\mathrm{i}\gamma^{0}\gamma^{2}$. Furthermore, the
$\Gamma^{\mathcal{T}}_{\bar{\lambda}}(k,k^{\prime})$ vertex is modelled as a
coupling to the fermionic energy-momentum tensor, see e. g. [49, 50, 51, 52].

In (9) $E(\bar{\lambda},\Delta_{0})$ describe the polarization state of the
axial quarkonia. To ensure that the axial vector carries the correct quantum
numbers, we require the transversality condition

$\Delta_{0}\cdot E(\bar{\lambda},\Delta_{0})=0\,.$ (11)

This way, in the $c\bar{c}$ rest frame, $E^{\mu}(\bar{\lambda},\Delta_{0})$
reduces to its spatial components corresponding to the total spin vector of
the $c\bar{c}$ state. Recall that in the light-cone formulation the “$-$” and
$\perp$ components of $\Delta_{0}$ are exactly equal to the corresponding
components of the quarkonia four-momentum $\Delta$, but that the “$+$”
components differ:
$\Delta^{+}=(M_{\mathcal{H}}^{2}+\boldsymbol{\Delta}_{\perp}^{2})/2\Delta^{-}$
and $\Delta_{0}^{+}=({M_{0}}^{2}+\boldsymbol{\Delta}_{\perp}^{2})/2\Delta^{-}$
with $M_{0}^{2}=\Delta_{0}^{2}$ being the invariant mass of the $c\bar{c}$
pair and $M^{2}_{\mathcal{H}}=\Delta^{2}$ the quarkonia mass. This is because
the $c\bar{c}$ pair is allowed to be virtual in the “$+$” component, conjugate
to the light-cone time. Explicitly, we have

$E^{\mu}(\bar{\lambda}=\pm
1,\Delta_{0})=\left(\frac{\boldsymbol{\epsilon}_{\perp}^{\bar{\lambda}}\cdot\boldsymbol{\Delta}_{\perp}}{\Delta^{-}},0,\boldsymbol{\epsilon}_{\perp}^{\bar{\lambda}}\right)\,,\quad
E^{\mu}(\bar{\lambda}=0,\Delta_{0})=\left(\frac{\boldsymbol{\Delta}_{\perp}^{2}-M_{0}^{2}}{2{M_{0}}\Delta^{-}},\frac{\Delta^{-}}{{M_{0}}},\frac{\boldsymbol{\Delta}_{\perp}}{{M_{0}}}\right)=\frac{1}{{M_{0}}}\Delta_{0}^{\mu}-\frac{{M_{0}}}{\Delta^{-}}n^{\mu}\,.$
(12)

Using $\Delta$ instead of $\Delta_{0}$ has no effect on the $\bar{\lambda}=\pm
1$ components, namely $E^{\mu}(\bar{\lambda}=\pm
1,\Delta_{0})=E^{\mu}(\bar{\lambda}=\pm 1,\Delta)$. However, the
$\bar{\lambda}=0$ case differs as

$E(\bar{\lambda}=0,\Delta)=\frac{1}{M_{\mathcal{A}}}\Delta^{\mu}-\frac{M_{\mathcal{A}}}{\Delta^{-}}n^{\mu}\,.$
(13)

We have checked that the projector
$\mathrm{i}\gamma_{5}\not{E}(\bar{\lambda}=0,\Delta)$ (instead of
$\mathrm{i}\gamma_{5}\not{E}(\bar{\lambda}=0,\Delta_{0})$) leads to incorrect
results for the $c\bar{c}$ state as it leaves an admixture with spin different
from 1. This happens because in the light-cone formulation the rest frames of
the quarkonia and of the $c\bar{c}$ partonic state in general move with
different velocities, which is a consequence of the difference in their “+”
components of the four-momenta. The necessary condition for the correct
projector, $\Delta_{0}\cdot E(\bar{\lambda},\Delta)=0$, is fulfilled if and
only if ${M_{0}}=M_{\mathcal{A}}$. In other words, we will assume that the
amplitude for the transition from $E^{\mu}(\bar{\lambda}=0,\Delta_{0})$ to
$E^{\mu}(\bar{\lambda}=0,\Delta)$ is equal to one. A similar approach was used
in Ref. [53]. The underlying principle of parton–hadron duality is also the
basis for the successful approach to exclusive $\rho^{0}$ electroproduction
[54, 55].

The tensor quarkonia polarizations $E^{\mu\nu}(\bar{\lambda},\Delta_{0})$ can
be obtained in terms of $E^{\mu}(\bar{\lambda},\Delta_{0})$, using Clebsch-
Gordan coefficients. We have [49]

$\begin{split}&E^{\mu\nu}(\pm 2,\Delta_{0})=E^{\mu}(\pm
1,\Delta_{0})E^{\nu}(\pm 1,\Delta_{0})\,,\\\ &E^{\mu\nu}(\pm
1,\Delta_{0})=\frac{1}{\sqrt{2}}\big{(}E^{\mu}(\pm
1,\Delta_{0})E^{\nu}(0,\Delta_{0})+E^{\mu}(0,\Delta_{0})E^{\nu}(\pm
1,\Delta_{0})\big{)}\,,\\\
&E^{\mu\nu}(0,\Delta_{0})=\frac{1}{\sqrt{6}}\big{(}E^{\mu}(+1,\Delta_{0})E^{\nu}(-1,\Delta_{0})+E^{\mu}(-1,\Delta_{0})E^{\nu}(+1,\Delta_{0})+2E^{\mu}(0,\Delta_{0})E^{\nu}(0,\Delta_{0})\big{)}\,.\\\
\end{split}$ (14)

These polarization tensors are traceless,
$g_{\mu\nu}E^{\mu\nu}(\bar{\lambda},\Delta_{0})=0$, and symmetric,
$E^{\mu\nu}(\bar{\lambda},\Delta_{0})=E^{\nu\mu}(\bar{\lambda},\Delta_{0})$.
We also have transversality conditions,
$\Delta_{0}^{\mu}E_{\mu\nu}(\bar{\lambda},\Delta_{0})=0=\Delta_{0}^{\nu}E_{\mu\nu}(\bar{\lambda},\Delta_{0})$,
as a consequence of (11). We have checked that Eqs. (14) are consistent with
the polarization tensors written in [51], up to overall signs.

It is instructive to explicitly compute the quarkonia wave functions in this
approach. Using LB spinors [46] we find for the scalar quarkonia

$\Psi^{\mathcal{S}}_{h\bar{h}}(\boldsymbol{k}_{\perp},z)=\frac{1}{z\bar{z}}\left[hk_{\perp}\mathrm{e}^{-\mathrm{i}h\phi_{k}}\delta_{h\bar{h}}-m_{c}(z-\bar{z})\delta_{h,-\bar{h}}\right]\phi_{\mathcal{S}}(\boldsymbol{k}_{\perp},z)\,.$
(15)

The helicity structure of the wave function agrees completely with the result
in [56], see (A.18) there. For the axial quarkonia we have

$\begin{split}&\Psi^{\mathcal{A}}_{+1,h\bar{h}}(\boldsymbol{k}_{\perp},z)=-\frac{\sqrt{2}\mathrm{i}}{z\bar{z}}\left[k_{\perp}\mathrm{e}^{\mathrm{i}\phi_{k}}(\bar{z}\delta_{h+}\delta_{\bar{h}-}+z\delta_{h-}\delta_{\bar{h}+})+m_{c}(z-\bar{z})\delta_{h+}\delta_{\bar{h}+}\right]\phi_{\mathcal{A},T}(\boldsymbol{k}_{\perp},z)\,,\\\
&\Psi^{\mathcal{A}}_{-1,h\bar{h}}(\boldsymbol{k}_{\perp},z)=-\frac{\sqrt{2}\mathrm{i}}{z\bar{z}}\left[k_{\perp}\mathrm{e}^{-\mathrm{i}\phi_{k}}(\bar{z}\delta_{h+}\delta_{\bar{h}-}+z\delta_{h-}\delta_{\bar{h}+})-m_{c}(z-\bar{z})\delta_{h-}\delta_{\bar{h}-}\right]\phi_{\mathcal{A},T}(\boldsymbol{k}_{\perp},z)\,,\\\
&\Psi^{\mathcal{A}}_{0,h\bar{h}}(\boldsymbol{k}_{\perp},z)=-\frac{\mathrm{i}}{z\bar{z}}\frac{1}{M_{\mathcal{A}}}\left(2\boldsymbol{k}_{\perp}^{2}h\delta_{h,-\bar{h}}+2m_{c}k_{\perp}\mathrm{e}^{-\mathrm{i}h\phi_{k}}\delta_{h\bar{h}}\right)\phi_{\mathcal{A},L}(\boldsymbol{k}_{\perp},z)\,,\end{split}$
(16)

where $k_{\perp}\mathrm{e}^{\pm\mathrm{i}\phi_{k}}\equiv
k^{1}\pm\mathrm{i}k^{2}$. At this point we have introduced different scalar
functions $\phi_{\mathcal{A},T}(r_{\perp},z)$ and
$\phi_{\mathcal{A},L}(r_{\perp},z)$ for transversely and longitudinally
polarized quarkonia777For the sake of simplicity, in the case of
$\bar{\lambda}=0$, we have redefined
$\phi_{\mathcal{A},L}(\boldsymbol{k}_{\perp},z)/M_{0}\to\phi_{\mathcal{A},L}(\boldsymbol{k}_{\perp},z)/M_{\mathcal{A}}$,
which introduces only subleading effects in the heavy quark limit.. The
helicity structure of (16) coincides with the axial part of the $Z$-boson wave
function from Refs. [57, 58]. In Ref. [59] the axial quarkonia wave function
was computed starting from the quarkonia rest frame followed by a Melosh
transformation. The resulting wave functions, their Eqs. (A.7) and (A.9),
agree up to an overall normalization constant with (16).

For the tensor quarkonia we obtain the following results

$\begin{split}&\Psi^{\mathcal{T}}_{+2,h\bar{h}}(\boldsymbol{k}_{\perp},z)=-\frac{2}{z\bar{z}}k_{\perp}\mathrm{e}^{\mathrm{i}\phi_{k}}\left[k_{\perp}\mathrm{e}^{\mathrm{i}\phi_{k}}(z\delta_{h+}\delta_{\bar{h}-}-\bar{z}\delta_{h-}\delta_{\bar{h}+})+m_{c}\delta_{h+}\delta_{\bar{h}+}\right]\phi_{\mathcal{T},T2}(\boldsymbol{k}_{\perp},z)\,,\\\
&\Psi^{\mathcal{T}}_{-2,h\bar{h}}(\boldsymbol{k}_{\perp},z)=\frac{2}{z\bar{z}}k_{\perp}\mathrm{e}^{-\mathrm{i}\phi_{k}}\left[k_{\perp}\mathrm{e}^{-\mathrm{i}\phi_{k}}(\bar{z}\delta_{h+}\delta_{\bar{h}-}-z\delta_{h-}\delta_{\bar{h}+})+m_{c}\delta_{h-}\delta_{\bar{h}-}\right]\phi_{\mathcal{T},T2}(\boldsymbol{k}_{\perp},z)\,,\\\
&\Psi^{\mathcal{T}}_{+1,h\bar{h}}(\boldsymbol{k}_{\perp},z)=\frac{M_{\mathcal{T}}}{z\bar{z}}\left[-k_{\perp}\mathrm{e}^{\mathrm{i}\phi_{k}}\left((3z-4z^{2})\delta_{h+}\delta_{\bar{h}-}+(3\bar{z}-4\bar{z}^{2})\delta_{h-}\delta_{\bar{h}+}\right)+m_{c}(z-\bar{z})\delta_{h+}\delta_{\bar{h}+}\right]\phi_{\mathcal{T},T}(\boldsymbol{k}_{\perp},z)\,,\\\
&\Psi^{\mathcal{T}}_{-1,h\bar{h}}(\boldsymbol{k}_{\perp},z)=\frac{M_{\mathcal{T}}}{z\bar{z}}\left[k_{\perp}\mathrm{e}^{-\mathrm{i}\phi_{k}}\left((3\bar{z}-4\bar{z}^{2})\delta_{h+}\delta_{\bar{h}-}+(3z-4z^{2})\delta_{h-}\delta_{\bar{h}+}\right)+m_{c}(z-\bar{z})\delta_{h-}\delta_{\bar{h}-}\right]\phi_{\mathcal{T},T}(\boldsymbol{k}_{\perp},z)\,,\\\
&\Psi^{\mathcal{T}}_{0,h\bar{h}}(\boldsymbol{k}_{\perp},z)=\frac{\sqrt{2}}{\sqrt{3}}\frac{1}{z\bar{z}}\left[(3\boldsymbol{k}_{\perp}^{2}+2m_{c}^{2})(z-\bar{z})\delta_{h,-\bar{h}}+m_{c}(k_{\perp}\mathrm{e}^{-\mathrm{i}\phi_{k}}\delta_{h+}\delta_{\bar{h}+}-k_{\perp}\mathrm{e}^{\mathrm{i}\phi_{k}}\delta_{h-}\delta_{\bar{h}-})\right]\phi_{\mathcal{T},L}(\boldsymbol{k}_{\perp},z)\,.\\\
\end{split}$ (17)

Here we have three different scalar functions
$\phi_{\mathcal{T},T2}(\boldsymbol{k}_{\perp},z)$,
$\phi_{\mathcal{T},T}(\boldsymbol{k}_{\perp},z)$ and
$\phi_{\mathcal{T},L}(\boldsymbol{k}_{\perp},z)$ accounting for the two
possible transverse and longitudinal polarizations of the tensor
quarkonia888In the case $\bar{\lambda}=\pm 1$ we have redefined
$M_{0}\phi_{\mathcal{T},T}(\boldsymbol{k}_{\perp},z)\to
M_{\mathcal{T}}\phi_{\mathcal{T},T}(\boldsymbol{k}_{\perp},z)$..

At this point we make a comment on the spin-momentum structure of the obtained
wave functions in the rest frame. Since the $C$ parity of a $q\bar{q}$ bound
state is $C=(-1)^{L+S}$, the quarkonia with $J=0$, $1$ and $2$, are all
$P$-waves ($L=1$), and have $S=1$. Based on our model wave function (8) we
therefore explicitly checked that, after replacing the light-cone LB spinors
with Dirac spinors ($u_{h}(k)\to u_{s}(k)$ and $v_{h}(k)\to v_{s}(k)$), the
wave function in the quarkonia rest frame have their expected non-relativistic
${}^{3}P_{J}$ spin-momentum structure. The scalar wave function is
proportional to
$\xi^{\dagger}_{\bar{s}}(\boldsymbol{\sigma}\cdot\boldsymbol{k})\tilde{\xi}_{s}$,
while for the axial case we find
$\xi^{\dagger}_{\bar{s}}\boldsymbol{\sigma}\cdot(\boldsymbol{k}\times\boldsymbol{E})\tilde{\xi}_{s}$
[60, 53, 56, 59]. Here $\xi_{s}$ are the standard two-component Pauli spinors.
In case of the tensor wave function we obtain
$\xi^{\dagger}_{\bar{s}}E^{ij}\sigma_{i}k_{j}\tilde{\xi}_{s}$ [61].

###  II.2 Final expressions for the amplitudes

For explicit computations it is convenient to write the reduced amplitude (2)
in the following equivalent form, see e.g. Refs. [47, 62, 16]

$\mathcal{A}_{\lambda\bar{\lambda}}(\boldsymbol{r}_{\perp},\boldsymbol{\Delta}_{\perp})=eq_{c}\int_{z}\mathrm{e}^{-\mathrm{i}\boldsymbol{\delta}_{\perp}\cdot\boldsymbol{r}_{\perp}}\int_{\boldsymbol{l}_{\perp}}\frac{\mathrm{e}^{\mathrm{i}\boldsymbol{l}_{\perp}\cdot\boldsymbol{r}_{\perp}}}{\boldsymbol{l}_{\perp}^{2}+\varepsilon^{2}}\int_{\boldsymbol{l}_{\perp}^{\prime}}\mathrm{e}^{-\mathrm{i}(\boldsymbol{l}_{\perp}^{\prime}-z\boldsymbol{\Delta}_{\perp})\cdot\boldsymbol{r}_{\perp}}\phi_{\mathcal{H}}(\boldsymbol{l}_{\perp}^{\prime}-z\boldsymbol{\Delta}_{\perp},z)\frac{1}{z\bar{z}}A_{\lambda\bar{\lambda}}(\boldsymbol{l}_{\perp},\boldsymbol{l}_{\perp}^{\prime}-z\boldsymbol{\Delta}_{\perp},z)\,,$
(18)

where
$\boldsymbol{\delta}_{\perp}=\frac{1}{2}(z-\bar{z})\boldsymbol{\Delta}_{\perp}=(z-\frac{1}{2})\boldsymbol{\Delta}_{\perp}$.
The helicity sum in (2) is turned into a covariant Dirac trace

$A_{\lambda\bar{\lambda}}(\boldsymbol{l}_{\perp},\boldsymbol{l}_{\perp}^{\prime}-z\boldsymbol{\Delta}_{\perp},z)=\frac{1}{(2q^{-})^{2}}{\rm
tr}\left[(\not{l}+m_{c})\not{\epsilon}(\lambda,q)(\not{l}-\not{q}+m_{c})\gamma^{-}(\not{l}^{\prime}-\not{\Delta}+m_{c})\Gamma_{\bar{\lambda}}^{\mathcal{H}*}(l^{\prime},\Delta-l^{\prime})(\not{l}^{\prime}+m_{c})\gamma^{-}\right]\,,$
(19)

containing the physical information on the polarization dependent part of the
photon-quarkonia wave function overlap. Together with Eq. (1), Eqs. (18) and
(19) comprise the main formulas that will be used below to write the
$\gamma^{*}p\to\mathcal{H}p$ amplitudes in a form suitable for numerical
computations. For explicit evaluation of the traces we have used FeynCalc
[63].

We start with the computation of the amplitude for the scalar quarkonia.
Inserting the first line of (9) into (19) yields

$\begin{split}&A_{\lambda=0}(\boldsymbol{l}_{\perp},\boldsymbol{l}_{1\perp},z)=-4m_{c}Qz\bar{z}(z-\bar{z})\,,\\\
&A_{\lambda=\pm
1}(\boldsymbol{l}_{\perp},\boldsymbol{l}_{1\perp},z)=-2m_{c}\left[(z-\bar{z})^{2}(\boldsymbol{\epsilon}_{\perp}^{\lambda}\cdot\boldsymbol{l}_{\perp})-(\boldsymbol{\epsilon}_{\perp}^{\lambda}\cdot\boldsymbol{l}_{1\perp})\right]\,,\\\
\end{split}$ (20)

where
$\boldsymbol{l}_{1\perp}\equiv\boldsymbol{l}_{\perp}^{\prime}-z\boldsymbol{\Delta}_{\perp}$.
As a cross check, note that the traces (20) are odd under joint
$\boldsymbol{l}_{\perp}\to-\boldsymbol{l}_{\perp}$,
$\boldsymbol{l}_{1\perp}\to-\boldsymbol{l}_{1\perp}$ and $z\to\bar{z}$
transformation. In other words, the reduced amplitude
$\mathcal{A}_{\lambda}(\boldsymbol{r}_{\perp},\boldsymbol{\Delta}_{\perp})$ in
(18) is odd under $\boldsymbol{r}_{\perp}\to-\boldsymbol{r}_{\perp}$. This is
consistent with
$\mathcal{O}(\boldsymbol{r}_{\perp},\boldsymbol{\Delta}_{\perp})$ being odd
under $\boldsymbol{r}_{\perp}\to-\boldsymbol{r}_{\perp}$ in the full amplitude
(1).

We now plug (20) into (18) and Fourier transform to coordinate space. After
separating the explicit polarization dependence, the reduced amplitudes are
found to be

$\begin{split}&\mathcal{A}_{0}(\boldsymbol{r}_{\perp},\boldsymbol{\Delta}_{\perp})=eq_{c}\int_{z}\mathrm{e}^{-\mathrm{i}\boldsymbol{\delta}_{\perp}\cdot\boldsymbol{r}_{\perp}}\mathcal{A}_{L}(r_{\perp})\,,\\\
&\mathcal{A}_{\lambda=\pm
1}(\boldsymbol{r}_{\perp},\boldsymbol{\Delta}_{\perp})=eq_{c}\lambda\mathrm{e}^{\mathrm{i}\lambda\phi_{r}}\int_{z}\mathrm{e}^{-\mathrm{i}\boldsymbol{\delta}_{\perp}\cdot\boldsymbol{r}_{\perp}}\mathcal{A}_{T}(r_{\perp})\,,\\\
\end{split}$ (21)

where

$\begin{split}&\mathcal{A}_{L}(r_{\perp})\equiv-\frac{2}{\pi}m_{c}Q(z-\bar{z})K_{0}(\varepsilon
r_{\perp})\phi_{\mathcal{S}}(r_{\perp},z)\,,\\\
&\mathcal{A}_{T}(r_{\perp})\equiv\frac{\mathrm{i}\sqrt{2}}{2\pi}\frac{m_{c}}{z\bar{z}}\left[(z-\bar{z})^{2}\varepsilon
K_{1}(\varepsilon r_{\perp})\phi_{\mathcal{S}}(r_{\perp},z)-K_{0}(\varepsilon
r_{\perp})\frac{\partial\phi_{\mathcal{S}}}{\partial
r_{\perp}}\right]\,.\end{split}$ (22)

In the amplitude (1) we are left with the $\boldsymbol{r}_{\perp}$ integral
and the $z$ integral. We use the Jacobi-Anger expansion to compute the
$\phi_{r}$ integral and find

$\begin{split}\left\langle\mathcal{M}_{\lambda=0}\right\rangle&\equiv
q^{-}\widetilde{\mathcal{M}}_{\lambda=0}\equiv
q^{-}\widetilde{\mathcal{M}}_{L}\,,\\\ \left\langle\mathcal{M}_{\lambda=\pm
1}\right\rangle&\equiv q^{-}\widetilde{\mathcal{M}}_{\lambda=\pm 1}\equiv
q^{-}\lambda\mathrm{e}^{\mathrm{i}\lambda\phi_{\Delta}}\widetilde{\mathcal{M}}_{T}\,,\\\
\end{split}$ (23)

where in the first equality we have separated the flux factor and in the
second equality the explicit polarization dependence (together with the
overall phase). The remainder consists of two different scalar functions

$\begin{split}&\widetilde{\mathcal{M}}_{L}=8\pi
N_{c}eq_{c}\sum_{k=0}^{\infty}(-1)^{k}\int_{z}\int_{0}^{\infty}r_{\perp}\mathrm{d}r_{\perp}\mathcal{O}_{2k+1}(r_{\perp},\Delta_{\perp})\mathcal{A}_{L}(r_{\perp})\,{\rm
sgn}(z-\bar{z})J_{2k+1}(r_{\perp}\delta_{\perp})\,,\\\
&\widetilde{\mathcal{M}}_{T}=4\pi\mathrm{i}N_{c}eq_{c}\sum_{k=0}^{\infty}(-1)^{k}\int_{z}\int_{0}^{\infty}r_{\perp}\mathrm{d}r_{\perp}\mathcal{O}_{2k+1}(r_{\perp},\Delta_{\perp})\mathcal{A}_{T}(r_{\perp})\left[J_{2k}(r_{\perp}\delta_{\perp})-J_{2k+2}(r_{\perp}\delta_{\perp})\right]\,,\end{split}$
(24)

which are our main results to be used in the numerical computations in Sec. V.
The sign function appears due to the definition of $\delta_{\perp}$ as the
modulus of the vector $\boldsymbol{\delta}_{\perp}$:
$\delta_{\perp}=\frac{1}{2}|z-\bar{z}|\Delta_{\perp}$.

The quantities $\mathcal{O}_{2k+1}(r_{\perp},\Delta_{\perp})$ appearing in Eq.
(24) are the azimuthal harmonics of the Odderon amplitude,

$\mathcal{O}_{2k+1}(r_{\perp},\Delta_{\perp})=2\pi\mathrm{i}(-1)^{k+1}\int_{0}^{\infty}b_{\perp}\mathrm{d}b_{\perp}\mathcal{O}_{2k+1}(r_{\perp},b_{\perp})J_{2k+1}(\Delta_{\perp}b_{\perp})\,,$
(25)

where $\mathcal{O}_{2k+1}(r_{\perp},b_{\perp})$ is extracted from
$\mathcal{O}(\boldsymbol{r}_{\perp},\boldsymbol{b}_{\perp})$ as its Fourier
series coefficients

$\mathcal{O}(\boldsymbol{r}_{\perp},\boldsymbol{b}_{\perp})=2\sum_{k=0}^{\infty}\mathcal{O}_{2k+1}(r_{\perp},b_{\perp})\cos((2k+1)\phi_{rb})\,,\qquad\mathcal{O}_{2k+1}(r_{\perp},b_{\perp})=\frac{1}{2\pi}\int_{0}^{2\pi}\mathrm{d}\phi_{rb}\mathcal{O}(\boldsymbol{r}_{\perp},\boldsymbol{b}_{\perp})\cos((2k+1)\phi_{rb})\,,$
(26)

with $\phi_{rb}=\phi_{r}-\phi_{b}$.

The calculation for axial vector quarkonia follows the same steps as for
scalar case. Detailed derivation of the formulas can be found in App. B. We
end up with three different scalar functions:

$\begin{split}&\widetilde{\mathcal{M}}_{B}=4\pi\mathrm{i}N_{c}eq_{c}\sum_{k=0}^{\infty}(-1)^{k}\int_{z}\int_{0}^{\infty}r_{\perp}\mathrm{d}r_{\perp}\mathcal{O}_{2k+1}(r_{\perp},\Delta_{\perp})\mathcal{A}_{B}(r_{\perp})\left[J_{2k}(r_{\perp}\delta_{\perp})-J_{2k+2}(r_{\perp}\delta_{\perp})\right]\,,\qquad
B=TL,LT\\\ &\widetilde{\mathcal{M}}_{TT}=8\pi
N_{c}eq_{c}\sum_{k=0}^{\infty}(-1)^{k}\int_{z}\int_{0}^{\infty}r_{\perp}\mathrm{d}r_{\perp}\mathcal{O}_{2k+1}(r_{\perp},\Delta_{\perp})\mathcal{A}_{TT}(r_{\perp})\,{\rm
sgn}(z-\bar{z})J_{2k+1}(r_{\perp}\delta_{\perp})\,,\end{split}$ (27)

with respective reduced amplitudes as

$\begin{split}&\mathcal{A}_{LT}(r_{\perp})\equiv\frac{\sqrt{2}}{\pi}QK_{0}(\varepsilon
r_{\perp})\frac{\partial\phi_{\mathcal{A},T}}{\partial r_{\perp}}\,,\\\
&\mathcal{A}_{TL}(r_{\perp})\equiv\frac{\sqrt{2}}{\pi}\frac{1}{z\bar{z}}\frac{1}{M_{\mathcal{A}}}\left[-m_{c}^{2}K_{0}(\varepsilon
r_{\perp})\frac{\partial\phi_{\mathcal{A},L}}{\partial r_{\perp}}+\varepsilon
K_{1}(\varepsilon
r_{\perp})\boldsymbol{\nabla}^{2}_{\perp}\phi_{\mathcal{A},L}\right]\,,\\\
&\mathcal{A}_{TT}(r_{\perp})\equiv-\frac{\mathrm{i}}{\pi}\frac{z-\bar{z}}{z\bar{z}}\left[\frac{\partial\phi_{\mathcal{A},T}}{\partial
r_{\perp}}\varepsilon K_{1}(\varepsilon r_{\perp})-m_{c}^{2}K_{0}(\varepsilon
r_{\perp})\phi_{\mathcal{A},T}\right]\,.\end{split}$ (28)

Here the notation $\widetilde{\mathcal{M}}_{TT}$ stands for amplitude with
transition from transversely polarized photon to transversely polarized axial
quarkonia. In principle, both the polarization preserving and polarization
flipping transition would be possible. However we find that the polarization
flipping transition vanishes, and so $\widetilde{\mathcal{M}}_{TT}$ describes
only a polarization preserving transition. Similar notation is used for other
contributing amplitudes. We find the amplitude when both the photon and the
axial quarkonia would be longitudinal also vanishes, see App. B for more
details.

The procedure is again similar for tensor quarkonia, with detailed steps to be
found in App. B. The final results for the amplitudes are

$\begin{split}\widetilde{\mathcal{M}}_{B}=-4\pi
N_{c}eq_{c}\sum_{k=0}^{\infty}(-1)^{k}&\int_{z}\int_{0}^{\infty}r_{\perp}\mathrm{d}r_{\perp}\mathcal{O}_{2k+1}(r_{\perp},\Delta_{\perp})\mathcal{A}_{B}(r_{\perp})\\\
&\times{\rm
sgn}(z-\bar{z})\left[J_{2k+3}(r_{\perp}\delta_{\perp})+J_{2k-1}(r_{\perp}\delta_{\perp})\right]\,,\quad
B=LT2,TTf\,,\\\
\widetilde{\mathcal{M}}_{B}=4\pi\mathrm{i}N_{c}eq_{c}\sum_{k=0}^{\infty}(-1)^{k}&\int_{z}\int_{0}^{\infty}r_{\perp}\mathrm{d}r_{\perp}\mathcal{O}_{2k+1}(r_{\perp},\Delta_{\perp})\mathcal{A}_{B}(r_{\perp})\\\
&\times\left[J_{2k}(r_{\perp}\delta_{\perp})-J_{2k+2}(r_{\perp}\delta_{\perp})\right]\,,\quad
B=TT2p,LT,TL\,,\\\ \widetilde{\mathcal{M}}_{B}=8\pi
N_{c}eq_{c}\sum_{k=0}^{\infty}(-1)^{k}&\int_{z}\int_{0}^{\infty}r_{\perp}\mathrm{d}r_{\perp}\mathcal{O}_{2k+1}(r_{\perp},\Delta_{\perp})\mathcal{A}_{B}(r_{\perp}){\rm
sgn}(z-\bar{z})J_{2k+1}(r_{\perp}\delta_{\perp})\,,\quad B=TTp,LL\,,\\\
\widetilde{\mathcal{M}}_{TT2f}=4\pi\mathrm{i}N_{c}eq_{c}\sum_{k=0}^{\infty}(-1)^{k}&\int_{z}\int_{0}^{\infty}r_{\perp}\mathrm{d}r_{\perp}\mathcal{O}_{2k+1}(r_{\perp},\Delta_{\perp})\mathcal{A}_{TT2f}(r_{\perp})\left[J_{2k+4}(r_{\perp}\delta_{\perp})-J_{2k-2}(r_{\perp}\delta_{\perp})\right]\,,\\\
\end{split}$ (29)

where

$\begin{split}&\mathcal{A}_{LT2}(r_{\perp})\equiv\frac{2}{\pi}(z-\bar{z})QK_{0}(\varepsilon
r_{\perp})\left(\frac{\partial^{2}\phi_{\mathcal{T},T2}}{\partial
r_{\perp}^{2}}-\frac{1}{r_{\perp}}\frac{\partial\phi_{\mathcal{T},T2}}{\partial
r_{\perp}}\right)\,,\\\
&\mathcal{A}_{TT2p}(r_{\perp})\equiv-\frac{\mathrm{i}\sqrt{2}}{\pi}\frac{1}{z\bar{z}}\left((z^{2}+\bar{z}^{2})\varepsilon
K_{1}(\varepsilon
r_{\perp})\left(\frac{\partial^{2}\phi_{\mathcal{T},T2}}{\partial
r_{\perp}^{2}}-\frac{1}{r_{\perp}}\frac{\partial\phi_{\mathcal{T},T2}}{\partial
r_{\perp}}\right)+m_{c}^{2}K_{0}(\varepsilon
r_{\perp})\frac{\partial\phi_{\mathcal{T},T2}}{\partial
r_{\perp}}\right)\,,\\\
&\mathcal{A}_{TT2f}(r_{\perp})\equiv\frac{\mathrm{i}2\sqrt{2}}{\pi}\varepsilon
K_{1}(\varepsilon
r_{\perp})\left(\frac{\partial^{2}\phi_{\mathcal{T},T2}}{\partial
r_{\perp}^{2}}-\frac{1}{r_{\perp}}\frac{\partial\phi_{\mathcal{T},T2}}{\partial
r_{\perp}}\right)\,,\\\
&\mathcal{A}_{LT}(r_{\perp})\equiv-\frac{\mathrm{i}}{\pi}QM_{\mathcal{T}}\left(3-4(z^{2}+\bar{z}^{2})\right)K_{0}(\varepsilon
r_{\perp})\frac{\partial\phi_{\mathcal{T},T}}{\partial r_{\perp}}\,,\\\
&\mathcal{A}_{TTp}(r_{\perp})=-\frac{M_{\mathcal{T}}}{\sqrt{2}\pi}\frac{z-\bar{z}}{z\bar{z}}\left[m_{c}^{2}K_{0}(\varepsilon
r_{\perp})\phi_{\mathcal{T},T}(r_{\perp},z)-(z-\bar{z})^{2}\varepsilon
K_{1}(\varepsilon r_{\perp})\frac{\partial\phi_{\mathcal{T},T}}{\partial
r_{\perp}}\right]\,,\\\
&\mathcal{A}_{TTf}(r_{\perp})=-\frac{2\sqrt{2}M_{\mathcal{T}}}{\pi}(z-\bar{z})\varepsilon
K_{1}(\varepsilon r_{\perp})\frac{\partial\phi_{\mathcal{T},T}}{\partial
r_{\perp}}\,,\\\
&\mathcal{A}_{LL}(r_{\perp})=-\frac{2\sqrt{2}Q}{\sqrt{3}\pi}(z-\bar{z})K_{0}(\varepsilon
r_{\perp})\left(3\boldsymbol{\nabla}_{\perp}^{2}-2m_{c}^{2}\right)\phi_{\mathcal{T},L}(r_{\perp},z)\,,\\\
&\mathcal{A}_{TL}(r_{\perp})=\frac{\mathrm{i}}{\pi\sqrt{3}}\frac{1}{z\bar{z}}\left[\varepsilon
K_{1}(\varepsilon
r_{\perp})(z-\bar{z})^{2}\left(3\boldsymbol{\nabla}_{\perp}^{2}-2m_{c}^{2}\right)\phi_{\mathcal{T},L}-m_{c}^{2}K_{0}(\varepsilon
r_{\perp})\frac{\partial\phi_{\mathcal{T},L}}{\partial
r_{\perp}}\right]\,.\end{split}$ (30)

In this case the $LL$-type transition is allowed. Also, both the polarization
preserving and the polarization flipped contributions are allowed, thus
explaining the above used notation: $p$ = preserving, $f$ = flipped.

###  II.3 The $\gamma^{*}p\to\mathcal{H}p$ cross section

Using the $\widetilde{\mathcal{M}}_{\lambda\bar{\lambda}}$ amplitudes obtained
in the previous Sec. II.2 for the process
$\gamma^{*}(q)p(P)\to\mathcal{H}(\Delta)p(P^{\prime})$, the cross section for
longitudinal ($L$) and transverse ($T$) photons is given by

$\begin{split}&\frac{\mathrm{d}\sigma_{L}(\gamma^{*}p\to\mathcal{H}p)}{\mathrm{d}|t|}=\frac{1}{16\pi}\sum_{\bar{\lambda}=-J}^{J}\left|\widetilde{\mathcal{M}}_{\lambda=0,\bar{\lambda}}(\gamma^{*}p\to\mathcal{H}p)\right|^{2}\,,\\\
&\frac{\mathrm{d}\sigma_{T}(\gamma^{*}p\to\mathcal{H}p)}{\mathrm{d}|t|}=\frac{1}{16\pi}\frac{1}{2}\sum_{\lambda=\pm
1}\sum_{\bar{\lambda}=-J}^{J}\left|\widetilde{\mathcal{M}}_{\lambda\bar{\lambda}}(\gamma^{*}p\to\mathcal{H}p)\right|^{2}\,.\end{split}$
(31)

Inserting the scalar quarkonia amplitudes (23) into (31) gives

$\begin{split}&\frac{\mathrm{d}\sigma_{L}(\gamma^{*}p\to\mathcal{S}p)}{\mathrm{d}|t|}=\frac{1}{16\pi}\left|\widetilde{\mathcal{M}}_{L}\right|^{2}\,,\\\
&\frac{\mathrm{d}\sigma_{T}(\gamma^{*}p\to\mathcal{S}p)}{\mathrm{d}|t|}=\frac{1}{16\pi}\left|\widetilde{\mathcal{M}}_{T}\right|^{2}\,.\end{split}$
(32)

where $\widetilde{\mathcal{M}}_{B}$ ($B=L$, $T$) are given in Eqs. (24). For
the axial vector quarkonia we find

$\begin{split}&\frac{\mathrm{d}\sigma_{L}(\gamma^{*}p\to\mathcal{A}p)}{\mathrm{d}|t|}=\frac{1}{8\pi}\left|\widetilde{\mathcal{M}}_{LT}\right|^{2}\,,\\\
&\frac{\mathrm{d}\sigma_{T}(\gamma^{*}p\to\mathcal{A}p)}{\mathrm{d}|t|}=\frac{1}{16\pi}\left(\left|\widetilde{\mathcal{M}}_{TL}\right|^{2}+\left|\widetilde{\mathcal{M}}_{TT}\right|^{2}\right)\,,\end{split}$
(33)

with $\widetilde{\mathcal{M}}_{B}$ ($B=LT$, $TL$, $TT$) given in (27). We have
used (72) (see Appendix B) to relate
$\widetilde{\mathcal{M}}_{\lambda\bar{\lambda}}$ to
$\widetilde{\mathcal{M}}_{A}$. For the tensor quarkonia, the cross sections
are

$\begin{split}&\frac{\mathrm{d}\sigma_{L}(\gamma^{*}p\to\mathcal{T}p)}{\mathrm{d}|t|}=\frac{1}{16\pi}\left(\left|\widetilde{\mathcal{M}}_{LL}\right|^{2}+2\left|\widetilde{\mathcal{M}}_{LT}\right|^{2}+2\left|\widetilde{\mathcal{M}}_{LT2}\right|^{2}\right)\,,\\\
&\frac{\mathrm{d}\sigma_{T}(\gamma^{*}p\to\mathcal{T}p)}{\mathrm{d}|t|}=\frac{1}{16\pi}\left(\left|\widetilde{\mathcal{M}}_{TL}\right|^{2}+\left|\widetilde{\mathcal{M}}_{TTp}\right|^{2}+\left|\widetilde{\mathcal{M}}_{TTf}\right|^{2}+\left|\widetilde{\mathcal{M}}_{TT2p}\right|^{2}+\left|\widetilde{\mathcal{M}}_{TT2f}\right|^{2}\right)\,,\end{split}$
(34)

with $\widetilde{\mathcal{M}}_{B}$ ($B$ = $LL$, $LT$, $LT2$, $TL$, $TTp$,
$TTf$, $TT2p$, $TT2f$) collected in (29) and the relationship to
$\widetilde{\mathcal{M}}_{\lambda\bar{\lambda}}$ in (76).

###  II.4 Odderon amplitude and its evolution with $x$

Here we discuss the Odderon amplitude
$\mathcal{O}(\boldsymbol{r}_{\perp},\boldsymbol{b}_{\perp})$ used later in
Sec. V for numerical predictions of the cross section for exclusive
$\chi_{cJ}$ production. The initial
$\mathcal{O}(\boldsymbol{r}_{\perp},\boldsymbol{b}_{\perp})$ at
$x_{\mathcal{P}}=x_{0}=0.01$ is approximated by the matrix element of the
$C$-odd three gluon operator. This matrix element has been evaluated
numerically and tabulated in Ref. [64] using a phenomenological non-
perturbative three quark light-cone wave function [65, 66] supplemented by the
first correction of perturbative QCD [23]. We refer to these references for
further details. Specifically, we use the table for the first azimuthal
harmonic $a_{1}(r_{\perp},b_{\perp})$ of the
$\mathcal{O}(\boldsymbol{r}_{\perp},\boldsymbol{b}_{\perp})$ defined in Eq.
(5) of Ref. [64]. As explained in the App. D the definition of
$\mathcal{O}(\boldsymbol{r}_{\perp},\boldsymbol{b}_{\perp})$ in [64] is the
same as in this work999However, our own numerical investigation and private
communication with the authors of Ref. [64] revealed that their plots of
$\mathcal{O}(\boldsymbol{r}_{\perp},\boldsymbol{b}_{\perp})$ and their
tabulated values actually correspond to minus the function defined in their
Eq. (3) and in our Eq. (105). We therefore reversed the signs of their table
entries for $a_{1}(r_{\perp},b_{\perp})$.. The first azimuthal moment is
defined in [64] through
$\mathcal{O}(\boldsymbol{r}_{\perp},\boldsymbol{b}_{\perp})=\cos(\phi_{rb})a_{1}(r_{\perp},b_{\perp})+\dots$.
Thus, from (26) we have the relationship

$\mathcal{O}_{1}(r_{\perp},b_{\perp})=\frac{1}{2}a_{1}(r_{\perp},b_{\perp})\,.$
(35)

To determine $\mathcal{O}(\boldsymbol{r}_{\perp},\boldsymbol{b}_{\perp})$ for
$x_{\mathcal{P}}<x_{0}$ we solve the impact parameter dependent [67, 68]
extension of the Balitsky-Kovchegov (BK) equation (with running coupling
kernel – rcBK) to the coupled non-linear evolution equations for the Pomeron
and the Odderon derived in Refs. [11, 12] and analyzed in Refs. [13, 14, 15].
These equations describe the evolution of the real and imaginary parts of the
dipole $S$-matrix, and the initial condition is evolved in rapidity
$Y=\ln(x_{0}/x_{\mathcal{P}})$. For the real part (the Pomeron), we are using
the initial condition from [69]. The coupled BK equations for the Pomeron-
Odderon system are solved in the “local approximation” [70, 69], where the
impact parameter $\boldsymbol{b}_{\perp}$ becomes an external parameter. This
approximation may be less justified for the Odderon than for the Pomeron since
the former amplitude peaks at smaller $b_{\perp}$ [64]. The $t$-dependence of
three gluon exchange is then obtained via Fourier transform, Eq. (7).
Technical details of the implementation and additional numerical results can
be found in Sec. III of Ref. [16].

Figure 2: The first azimuthal harmonic of the Odderon amplitude
$\mathrm{i}\mathcal{O}_{1}(r_{\perp},\Delta_{\perp})$ as a function of
$\Delta_{\perp}^{2}$ for $r_{\perp}=0.3$ fm and $\alpha_{S}=0.25$. Different
line styles correspond to different values of $x_{\mathcal{P}}$.

Fig. 2 shows the first azimuthal harmonic of the odderon amplitude
$\mathrm{i}\mathcal{O}_{1}(r_{\perp},\Delta_{\perp})$ as a function of
$\Delta_{\perp}^{2}$ for a dipole of size $r_{\perp}=0.3$ fm. We have divided
by the parameter $\sigma_{0}/2=16.36$ mb = 42.027 GeV-2 (area of the proton)
[69] to compensate for the dimensionality of
$\mathrm{i}\mathcal{O}_{1}(r_{\perp},\Delta_{\perp})$. At
$x_{\mathcal{P}}=10^{-2}$ this corresponds to the perturbative exchange of
three gluons (with negative $C$-parity) while for smaller $x_{\mathcal{P}}$
the BK resummation has been performed as described above. We note, first of
all, the small magnitude of
$\mathcal{O}_{1}(r_{\perp},\Delta_{\perp})<10^{-2}\ll 1$, i. e. that
scattering of perturbative (small) dipoles via $C$-odd exchanges is predicted
to be very weak. However, we also observe that the $C$-odd amplitude peaks at
fairly large $\Delta_{\perp}^{2}\simeq 1$ GeV2 and that it decreases by merely
a factor of $\sim 2$ as $\Delta_{\perp}^{2}$ increases further to $\simeq 3$
GeV2. This is a manifestation of the “Landshoff mechanism” [71], see also [72,
73, 74] whereby large-$t$ scattering via three gluon exchange is less likely
to break up the proton than scattering via two gluon exchange101010The proton
model light-cone wave function employed here, however, includes a perturbative
correction [23] due to the emission of a gluon from one of the valence quarks,
or the internal exchange of a gluon, not considered in the quoted papers by
Donnachie and Landshoff.. In simple terms this is due to the fact that a large
momentum transfer can be shared by the exchanged gluons, and so a comparable
transverse momentum can be transmitted to up to three partons in the proton,
resulting in a smaller increase of their invariant mass than in case of two
gluon exchange. Such a contribution has been illustrated on Fig. 1. The light-
cone Fock space amplitudes of the proton are strongly suppressed if the
invariant mass of the parton system is far from $\sim
N_{c}\Lambda_{\mathrm{QCD}}$. On the other hand, an amplitude which depends
weakly on momentum transfer, in impact parameter space would peak at small
$b_{\perp}$. Amplitudes with negative parity must vanish at $b_{\perp}=0$,
however, and hence are expected to have small magnitude.

Fig. 2 also shows the decrease of the Odderon exchange amplitude with
$x_{\mathcal{P}}$ by approximately a factor of 2 from
$x_{\mathcal{P}}=10^{-2}$ to $x_{\mathcal{P}}=10^{-4}$. The initial decrease
of the Odderon amplitude has been understood to originate from gluon Regge
trajectory suppressions in the BJKP equation [11, 13, 15]. In App. E we have
checked that by omitting the unitarity corrections, the Odderon still
decreases but with a smaller rate. Importantly, the sign of the Odderon
amplitude is preserved by the evolution equation, which is simply a
consequence of the fact that its evolution equation is linear in
$\mathcal{O}(\boldsymbol{r}_{\perp},\boldsymbol{b}_{\perp})$.

##  III The Primakoff contribution

In this section we focus on the Primakoff process, wherein the exclusive
$\gamma^{*}p\to\mathcal{H}p$ production proceeds via photon rather than
Odderon exchange. Typically, the cross section for a QCD process far exceeds
its QED counterparts, and so the latter can be safely neglected. This is,
however, not so in case of the Odderon because its QCD cross section starts at
order $\alpha_{S}^{6}$ and so the QED contribution becomes a competitive
background in Odderon searches, in particular for low momentum transfer.

We treat the Primakoff process in the high-energy approximation of eikonal
scattering, replacing the Odderon amplitude
$\mathcal{O}(\boldsymbol{r}_{\perp},\boldsymbol{\Delta}_{\perp})$ in Eq. (1)
with [32, 16]

$\mathcal{O}(\boldsymbol{r}_{\perp},\boldsymbol{\Delta}_{\perp})\to
8\pi\mathrm{i}q_{c}\alpha\sin\left(\frac{\boldsymbol{\Delta}_{\perp}\cdot\boldsymbol{r}_{\perp}}{2}\right)\frac{F_{1}(\Delta_{\perp})}{\boldsymbol{\Delta}_{\perp}^{2}}\,,$
(36)

at leading order in the electromagnetic coupling $\alpha=e^{2}/4\pi$. The
amplitude can be written as

$\left\langle\mathcal{M}_{\lambda\bar{\lambda}}(\gamma^{*}p\to\mathcal{H}p)\right\rangle\equiv-\frac{eF_{1}(\ell_{\perp})}{\boldsymbol{\ell}_{\perp}^{2}}n_{\mu}\mathcal{M}^{\mu}_{\lambda\bar{\lambda}}(\gamma^{*}\gamma^{*}\to\mathcal{H})\,,$
(37)

where
$\mathcal{M}^{\mu}_{\lambda\bar{\lambda}}(\gamma^{*}\gamma^{*}\to\mathcal{H})$
is the
$\gamma^{*}(\lambda,q)\gamma^{*}(\mu,\ell)\to\mathcal{H}(\bar{\lambda},\Delta)$
amplitude (in the high energy limit) that is obtained by inserting (36) into
(1). In (37), $\ell\equiv P-P^{\prime}$ is the momentum transfer from the
proton, with $\boldsymbol{\ell}_{\perp}=\boldsymbol{\Delta}_{\perp}$ in the
dipole frame and $F_{1}(\ell_{\perp})$ is the Dirac charge form factor of the
proton with $F_{1}(0)=1$, see App. D for more detail. Since (36) represents
single photon exchange it is more convenient to start from the expression in
Eq. (18) for the reduced amplitude and perform the $\boldsymbol{r}_{\perp}$
integral. This leads to

$\begin{split}\frac{1}{q^{-}}n_{\mu}\mathcal{M}^{\mu}_{\lambda\bar{\lambda}}(\gamma^{*}\gamma^{*}\to\mathcal{H})&=-8\pi\mathrm{i}q_{c}^{2}\alpha
N_{c}\int_{z}\int_{\boldsymbol{l}_{\perp}}\left[\frac{A_{\lambda\bar{\lambda}}(\boldsymbol{l}_{\perp}-\bar{z}\boldsymbol{\Delta}_{\perp},\boldsymbol{l}_{\perp},z)}{(\boldsymbol{l}_{\perp}-\bar{z}\boldsymbol{\Delta}_{\perp})^{2}+\varepsilon^{2}}-\frac{A_{\lambda\bar{\lambda}}(\boldsymbol{l}_{\perp}+z\boldsymbol{\Delta}_{\perp},\boldsymbol{l}_{\perp},z)}{(\boldsymbol{l}_{\perp}+z\boldsymbol{\Delta}_{\perp})^{2}+\varepsilon^{2}}\right]\frac{\phi_{\mathcal{H}}(\boldsymbol{l}_{\perp},z)}{z\bar{z}}\,,\end{split}$
(38)

The term $eF_{1}(\ell_{\perp})/\boldsymbol{\ell}_{\perp}^{2}$ in Eq. (37) is
understood as a part of the usual Weizsäcker-Williams photon field

$A^{\mu}(\ell)=(2\pi)\delta(\ell^{-})\frac{1}{\boldsymbol{\ell}_{\perp}^{2}}eF_{1}(\ell_{\perp})n^{\mu}\,,$
(39)

where $n^{\mu}$ is the light-cone gauge vector.

###  III.1 $|t|\to 0$ limit

In the high energy approximation it is possible to express the $|t|\to 0$
limit of the Primakoff cross section for a $C$-even quarkonia $\mathcal{H}$
with spin $J\neq 1$ in terms of the two-photon decay width
$\Gamma(\mathcal{H}\to\gamma\gamma)$

$\lim_{|t|\to 0}|t|\frac{\mathrm{d}\sigma(\gamma
p\to\mathcal{H}p)}{\mathrm{d}|t|}=\frac{8\pi\alpha(2J+1)\Gamma(\mathcal{H}\to\gamma\gamma)}{M_{\mathcal{H}}^{3}}\,.$
(40)

Eq. (40) is a model independent result and not affected by QCD corrections.
Hence it represents a stringent constraint on the model predictions and a
useful cross check. The key point for obtaining (40) is the relation of the
high energy impact factor for the quarkonia photoproduction, to the
quarkonia’s two-photon decay amplitude. We provide a derivation of (40) in
App. C. For axial vector quarkonia, like $\chi_{c1}$,
$\Gamma(\mathcal{A}\to\gamma\gamma)$ vanishes and Eq. (40) implies only that
$\lim_{|t|\to 0}|t|\mathrm{d}\sigma(\gamma p\to\mathcal{A}p)/\mathrm{d}|t|=0$,
as is consistent with the Landau-Yang (LY) theorem [75, 76], but does not
constrain the value of the cross section at $|t|\to 0$. In App. C we thus
provide a separate and general analysis of the $|t|\to 0$ limit for the axial
quarkonia.

In practice, we will use (40) to deduce $\Gamma(\mathcal{H}\to\gamma\gamma)$
from the computation of its left hand side. This will be required in Sec. IV
where we perform a fit of the quarkonia wave functions. Therefore, we now
focus on obtaining explicit expressions of the $\gamma^{*}p\to\mathcal{H}p$
amplitudes in the $|t|\to 0$ (and $Q^{2}\to 0$) limit.

Starting with the scalar quarkonia we use (38) and insert the appropriate
expressions for the reduced amplitudes
$A_{\lambda,\bar{\lambda}}(\boldsymbol{l}_{\perp},\boldsymbol{l}_{1\perp},z)$
found in Sec. II.2. The relevant reduced amplitude in the $Q^{2}\to 0$ limit
concerns only the second line of (20). The expression in the square bracket in
(38) is conveniently factored as

$\frac{A_{\lambda=\pm
1}(\boldsymbol{l}_{\perp}-\bar{z}\boldsymbol{\Delta}_{\perp},\boldsymbol{l}_{\perp},z)}{(\boldsymbol{l}_{\perp}-\bar{z}\boldsymbol{\Delta}_{\perp})^{2}+\varepsilon^{2}}-\frac{A_{\lambda=\pm
1}(\boldsymbol{l}_{\perp}+z\boldsymbol{\Delta}_{\perp},\boldsymbol{l}_{\perp},z)}{(\boldsymbol{l}_{\perp}+z\boldsymbol{\Delta}_{\perp})^{2}+\varepsilon^{2}}\to
2m_{c}\frac{(\boldsymbol{\epsilon}_{\perp}^{\lambda}\cdot\boldsymbol{\Delta}_{\perp})}{(\boldsymbol{l}_{\perp}^{2}+m_{c}^{2})^{2}}\left[\boldsymbol{l}_{\perp}^{2}+(z-\bar{z})^{2}m_{c}^{2}\right]\,,$
(41)

where we have set $Q^{2}=0$, performed an expansion around $\Delta_{\perp}\to
0$ to linear order, as well as the angular average according to
$l_{\perp}^{i}l_{\perp}^{j}\to\boldsymbol{l}_{\perp}^{2}\delta^{ij}/2$. We now
insert (41) into (38) and use second line of (23) to extract
$\widetilde{\mathcal{M}}_{T}$ as

$\widetilde{\mathcal{M}}_{T}\to 8\sqrt{2}\pi\mathrm{i}q_{c}^{2}\alpha
N_{c}m_{c}\frac{eF_{1}(\Delta_{\perp})}{\Delta_{\perp}}\int_{z}\int_{\boldsymbol{l}_{\perp}}\frac{\boldsymbol{l}_{\perp}^{2}+(z-\bar{z})^{2}m_{c}^{2}}{(\boldsymbol{l}_{\perp}^{2}+m_{c}^{2})^{2}}\frac{\phi_{\mathcal{S}}(\boldsymbol{l}_{\perp},z)}{z\bar{z}}\,.$
(42)

Performing similar steps for axial quarkonia we find that for the square
brackets in (38) the $O(\Delta_{\perp})$ contributions vanish for
$\bar{\lambda}=0$ and for $\bar{\lambda}=\pm 1$. A non-vanishing contribution
is found for $\bar{\lambda}=\pm 1$ at $O(\Delta_{\perp}^{2})$ leading to

$\widetilde{\mathcal{M}}_{TT}\to-16\pi q_{c}^{2}\alpha
N_{c}eF_{1}(\Delta_{\perp})\int_{z}\int_{\boldsymbol{l}_{\perp}}\frac{(z-\bar{z})^{2}m_{c}^{2}}{(\boldsymbol{l}_{\perp}^{2}+m_{c}^{2})^{2}}\frac{\phi_{\mathcal{A}}(\boldsymbol{l}_{\perp},z)}{z\bar{z}}\,,$
(43)

while the $\bar{\lambda}=0$ contribution is $O(\Delta_{\perp}^{3})$. Such a
special $\Delta_{\perp}$ dependence of the axial quarkonia amplitude is
explained in App. C.2 from general considerations.

For the calculation of the $|t|\to 0$ limit for the tensor quarkonia we also
need

$l_{\perp}^{i}l_{\perp}^{j}l_{\perp}^{m}l_{\perp}^{n}\to\frac{1}{8}\boldsymbol{l}_{\perp}^{4}(\delta^{ij}\delta^{mn}+\delta^{im}\delta^{jn}+\delta^{in}\delta^{jm})\,.$
(44)

Computing for various tensor polarizations $\bar{\lambda}$ we find that
$\bar{\lambda}=\pm 1$ does not contribute in the $|t|\to 0$ limit. As for the
$\bar{\lambda}=\pm 2$ and $\bar{\lambda}=0$ cases, we extract the amplitudes
(29) as

$\begin{split}&\widetilde{\mathcal{M}}_{TT2,p}\to
16\sqrt{2}\pi\mathrm{i}q_{c}^{2}\alpha
N_{c}\frac{eF_{1}(\Delta_{\perp})}{\Delta_{\perp}}\int_{z}\int_{\boldsymbol{l}_{\perp}}\frac{\boldsymbol{l}_{\perp}^{2}}{(\boldsymbol{l}_{\perp}^{2}+m_{c}^{2})^{2}}\left[(z^{2}+\bar{z}^{2})\boldsymbol{l}_{\perp}^{2}+m_{c}^{2}\right]\frac{\phi_{\mathcal{T},T2}(\boldsymbol{l}_{\perp},z)}{z\bar{z}}\,,\\\
&\widetilde{\mathcal{M}}_{TT2,f}\to 0\,,\\\
&\widetilde{\mathcal{M}}_{TL}\to-\frac{16}{\sqrt{3}}\pi\mathrm{i}q_{c}^{2}\alpha
N_{c}\frac{eF_{1}(\Delta_{\perp})}{\Delta_{\perp}}\int_{z}\int_{\boldsymbol{l}_{\perp}}\frac{m_{c}^{2}}{(\boldsymbol{l}_{\perp}^{2}+m_{c}^{2})^{2}}\left[\boldsymbol{l}_{\perp}^{2}-2(z-\bar{z}^{2})\left(\frac{3}{2}\boldsymbol{l}_{\perp}^{2}+m_{c}^{2}\right)\right]\frac{\phi_{\mathcal{T},L}(\boldsymbol{l}_{\perp},z)}{z\bar{z}}\,.\end{split}$
(45)

###  III.2 Adding the Pauli form factor

In general, the $\gamma^{*}p\to\mathcal{H}p$ amplitude can be related to the
$\gamma^{*}\gamma^{*}\to\mathcal{H}$ amplitude by separating out the Dirac
current with the Dirac ($F_{1}$) and Pauli ($F_{2}$) form factors

$\mathcal{M}_{\lambda\bar{\lambda}}(\gamma^{*}p\to\mathcal{H}p)=\mathcal{M}_{\lambda\bar{\lambda}}^{\mu}(\gamma^{*}\gamma^{*}\to\mathcal{H})\frac{g_{\mu\nu}}{\ell^{2}}\bar{u}_{h^{\prime}}(P^{\prime})\left[eF_{1}(\ell_{\perp})\gamma^{\nu}+\mathrm{i}\sigma^{\nu\rho}\ell_{\rho}\frac{eF_{2}(\ell_{\perp})}{2m_{N}}\right]u_{h}(P)\,,$
(46)

where $u_{h}(P)$ is a proton spinor with helicity $h$ and
$\sigma^{\nu\rho}=\mathrm{i}[\gamma^{\nu},\gamma^{\rho}]/2$. The dominant
contribution of the first term in the high energy limit is with $\gamma^{+}$,
while in the second term the dominant contribution is with $\sigma^{+i}$.
Thus, we can write in a covariant way

$\mathcal{M}_{\lambda\bar{\lambda}}(\gamma^{*}p\to\mathcal{H}p)\approx
P_{\mu}\mathcal{M}_{\lambda\bar{\lambda}}^{\mu}(\gamma^{*}\gamma^{*}\to\mathcal{H})\frac{1}{q\cdot
P}\bar{u}_{h^{\prime}}(P^{\prime})\left[\frac{eF_{1}(\ell_{\perp})}{\ell^{2}}\not{q}+\mathrm{i}\sigma^{\nu\rho}q_{\nu}\ell_{\rho}\frac{eF_{2}(\ell_{\perp})}{2m_{N}\ell^{2}}\right]u_{h}(P)\,.$
(47)

The above spinor vertices in the high energy limit become111111The calculation
steps are very similar to the one performed in App. A and also available in
[46].

$\begin{split}&\bar{u}_{h^{\prime}}(P^{\prime})\not{q}u_{h}(P)\approx
W^{2}\delta_{hh^{\prime}}\,,\\\
&\bar{u}_{h^{\prime}}(P^{\prime})\mathrm{i}\sigma^{\nu\rho}q_{\nu}\ell_{\rho}u_{h}(P)\approx\ell_{\perp}h\mathrm{e}^{\mathrm{i}h\phi_{\ell}}W^{2}\delta_{h,-h^{\prime}}\,.\end{split}$
(48)

Taking into account the connection to the usual notion of the color averaged
amplitude
$\left\langle\mathcal{M}_{\lambda\bar{\lambda}}\right\rangle=q^{-}\widetilde{\mathcal{M}}_{\lambda\bar{\lambda}}=q^{-}\mathcal{M}_{\lambda\bar{\lambda}}/W^{2}$,
we have

$\left\langle\mathcal{M}_{\lambda\bar{\lambda}}(\gamma^{*}p\to\mathcal{H}p)\right\rangle=-n_{\mu}\mathcal{M}_{\lambda\bar{\lambda}}^{\mu}(\gamma^{*}\gamma^{*}\to\mathcal{H})\left[\frac{eF_{1}(\ell_{\perp})}{\boldsymbol{\ell}_{\perp}^{2}}\delta_{hh^{\prime}}+\frac{eF_{2}(\ell_{\perp})}{\boldsymbol{\ell}_{\perp}^{2}}\frac{\ell_{\perp}}{2m_{N}}h\mathrm{e}^{\mathrm{i}h\phi_{\ell}}\delta_{h,-h^{\prime}}\right]\,,$
(49)

where the first term reproduces (37). Thanks to the difference in the helicity
structure in the high energy limit there is no interference between the
$F_{1}$ and $F_{2}$ contributions in (49). For the same reason there is no
intereference of the $F_{2}$ term with the Odderon exchange amplitude. The
addition of $F_{2}$ to the Primakoff cross section therefore amounts to a
simple replacement

$F_{1}^{2}(\ell_{\perp})\to
F_{1}^{2}(\ell_{\perp})+\frac{\boldsymbol{\ell}_{\perp}^{2}}{4m_{N}^{2}}F_{2}^{2}(\ell_{\perp})\,.$
(50)

Thanks to the additional $\ell_{\perp}$ dependence, the $F_{2}$ contribution
is negligible at $|t|=\boldsymbol{\ell}_{\perp}^{2}\to 0$ but in practice
becomes up to about $50\%$ correction at finite $|t|$ [77].

##  IV Boosted Gaussian model for the $\chi_{cJ}$ wave functions

In this section we perform a fit of the scalar part of the $C$-even quarkonia
wave function $\phi_{\mathcal{H},i}(\boldsymbol{k}_{\perp},z)$. We adopt a
Boosted Gaussian ansatz [78] originally written in coordinate space as

$\phi_{\mathcal{H},B}(r_{\perp},z)=\mathcal{N}_{\mathcal{H},B}z\bar{z}\exp\left(-\frac{m_{c}^{2}\mathcal{R}_{\mathcal{H}}^{2}}{8z\bar{z}}-\frac{2z\bar{z}r_{\perp}^{2}}{\mathcal{R}_{\mathcal{H}}^{2}}+\frac{1}{2}m_{c}^{2}\mathcal{R}_{\mathcal{H}}^{2}\right)\,,$
(51)

where $\mathcal{R}_{\mathcal{H}}$ is the radii parameter and
$\mathcal{N}_{\mathcal{H},i}$ are normalization parameters for different
quarkonia species and polarizations that are to be determined below. For
$\mathcal{H}=\mathcal{S}$ we have only one normalization,
$\mathcal{N}_{\mathcal{S}}$, for $\mathcal{H}=\mathcal{A}$, we have $B=T$,
$L$, denoting transverse and longitudinal polarizations. For
$\mathcal{H}=\mathcal{T}$, the index $B$ spans over $B=T2$, $T$, $L$. The wave
function is normalized as

$1=N_{c}\sum_{h\bar{h}}\int_{z}\int_{\boldsymbol{r}_{\perp}}\left|\Psi^{\mathcal{H}}_{\bar{\lambda},h\bar{h}}(\boldsymbol{r}_{\perp},z)\right|^{2}\,.$
(52)

In practice, we convert (52) to momentum space and express the helicity sum as
a Dirac trace as

$1=N_{c}\int_{\boldsymbol{k}_{\perp}}\int_{z}\frac{1}{z\bar{z}}{\rm
tr}\left[(\not{k}-m_{c})\gamma^{0}\Gamma^{\mathcal{H},{\dagger}}_{\bar{\lambda}}(k,k^{\prime})\gamma^{0}(\not{k}+m_{c})\Gamma^{\mathcal{H}}_{\bar{\lambda}}(k,k^{\prime})\right]\left|\phi_{\mathcal{H},B}(\boldsymbol{k}_{\perp},z)\right|^{2}\,.$
(53)

To constrain the parameters $\mathcal{R}_{\mathcal{H}}$ and
$\mathcal{N}_{\mathcal{H},i}$, we are using (53) and the
$\mathcal{H}\to\gamma\gamma$ decay width for $\chi_{c0}$ and $\chi_{c2}$. The
$\chi_{c1}\to\gamma\gamma$ decay is forbidden due to the LY theorem. In this
case our assumption is that the $\mathcal{R}$ parameter of $\chi_{c1}$ is
equal to that for $\chi_{c2}$.

The decay rate for $\mathcal{H}=\mathcal{S},$ $\mathcal{T}$ can be obtained
from the correspondence in Eq. (40). For scalar quarkonia we first insert the
amplitude (42) into (32) to find the $|t|\to 0$ limit of the cross section. We
then extract the decay rate via (40) as

$\begin{split}&\Gamma(\mathcal{S}\to\gamma\gamma)=\frac{\pi\alpha^{2}}{4}M_{\mathcal{S}}^{3}F_{\mathcal{S}}^{2}\,,\\\
&F_{\mathcal{S}}\equiv
4q_{c}^{2}m_{c}N_{c}\int_{z}\int_{\boldsymbol{k}_{\perp}}\frac{\boldsymbol{k}_{\perp}^{2}+(z-\bar{z})^{2}m_{c}^{2}}{(\boldsymbol{k}_{\perp}^{2}+m_{c}^{2})^{2}}\frac{\phi_{\mathcal{S}}(\boldsymbol{k}_{\perp},z)}{z\bar{z}}\,,\end{split}$
(54)

which agrees perfectly with Eq. (13) in Ref. [79] after adjusting the
conventions. We have additionally confirmed that in the NRQCD limit (54)
agrees with the known result in [80] (see App. C).

For the actual fit, the LO result in (54) is supplemented by the NLO QCD
corrections [81, 80]. This amounts to a replacement

$\Gamma(\mathcal{S}\to\gamma\gamma)\to\left(1+\frac{\alpha_{S}}{\pi}\left(\frac{\pi^{2}}{3}-\frac{28}{9}\right)\right)\Gamma(\mathcal{S}\to\gamma\gamma)\approx(1+0.06\alpha_{S})\Gamma(\mathcal{S}\to\gamma\gamma)\,.$
(55)

For the experimental value of $\chi_{c0}\to\gamma\gamma$ decay width we use
the most recent PDG value: $\Gamma(\chi_{c0}\to\gamma\gamma)=2.203\times
10^{-6}$ GeV [33]. Using $m_{c}=1.4$ GeV and $M_{\chi_{c0}}=3.414$ GeV and
$\alpha_{S}(2m_{c})\approx 0.25$, we determine $\mathcal{N}_{\chi_{c0}}=1.148$
and $\mathcal{R}_{\chi_{c0}}=1.539$ GeV-1. As expected, the radius parameter
$\mathcal{R}$ of the $\chi_{c0}$ is similar to that of the $J/\psi$ [78].

For $\mathcal{T}$ we insert the amplitude (45) into (34) and extract the decay
rate via (40) as

$\begin{split}&\Gamma(\mathcal{T}\to\gamma\gamma)=\frac{\pi\alpha^{2}}{20}M_{\mathcal{T}}^{3}\left(F_{\mathcal{T},L}^{2}+F_{\mathcal{T},T2}^{2}\right)\,,\\\
&F_{\mathcal{T},L}=\frac{4\sqrt{2}}{\sqrt{3}}q_{c}^{2}N_{c}\int_{z}\int_{\boldsymbol{k}_{\perp}}\frac{m_{c}^{2}\left[\boldsymbol{k}_{\perp}^{2}-2(z-\bar{z})^{2}\left(\frac{3}{2}\boldsymbol{k}_{\perp}^{2}+m_{c}^{2}\right)\right]}{(\boldsymbol{k}_{\perp}^{2}+m_{c}^{2})^{2}}\frac{\phi_{\mathcal{T},L}(\boldsymbol{k}_{\perp},z)}{z\bar{z}}\,,\\\
&F_{\mathcal{T},T2}=8q_{c}^{2}N_{c}\int_{z}\int_{\boldsymbol{k}_{\perp}}\frac{\boldsymbol{k}_{\perp}^{2}\left[(z^{2}+\bar{z}^{2})\boldsymbol{k}_{\perp}^{2}+m_{c}^{2}\right]}{(\boldsymbol{k}_{\perp}^{2}+m_{c}^{2})^{2}}\frac{\phi_{\mathcal{T},T2}(\boldsymbol{k}_{\perp},z)}{z\bar{z}}\,,\end{split}$
(56)

where $F_{\mathcal{T},L}$ ($F_{\mathcal{T},T2}$) originates from the
$\bar{\lambda}=0$ ($\bar{\lambda}=\pm 2$) polarizations. We find only a
partial agreement of our result (56) and Eqs. (22)-(24) of [79]. The third
line of (56) can be brought in agreement with Eq. (24) in [79] after
appropriate adjustments in conventions and also after a judicious
identification of $M_{\mathcal{T}}$ with the pair invariant mass $M_{0}$.
However, the square brackets of the second line in (56) contain a term
$\boldsymbol{k}_{\perp}^{2}$ while in (23) of [79] they rather have
$-m_{c}^{2}$. Using (56) in the NRQCD limit we find the $\bar{\lambda}=0$
contribution vanishes with the $\bar{\lambda}=\pm 2$ contribution saturating
the NRQCD limit completely and in accordance with the known result [80].
Taking the NRQCD limit of (22)-(24) in [79] gives an incorrect result as it
leads to a finite contribution from the $\bar{\lambda}=0$ polarization.

The QCD corrections [81, 80] turn out to be sizeable in this particular
channel [81, 80]

$\Gamma(\mathcal{T}\to\gamma\gamma)\to\left(1-\frac{16\alpha_{S}}{3\pi}\right)\Gamma(\mathcal{T}\to\gamma\gamma)\approx(1-1.7\alpha_{S})\Gamma(\mathcal{T}\to\gamma\gamma)\,.$
(57)

Using $M_{\chi_{c2}}=3.556$ GeV and the PDG value
$\Gamma(\chi_{c2}\to\gamma\gamma)=5.614\times 10^{-7}$ GeV [33], we find
$\mathcal{N}_{\chi_{c2},T2}=0.609$ GeV-1 $=\mathcal{N}_{\chi_{c2},L}$,
$\mathcal{N}_{\chi_{c2},T}=0.592$ GeV-1 and $\mathcal{R}_{\chi_{c2}}=1.482$
GeV-1. In these numerical estimates we used $\alpha_{S}=0.25$ in the NLO
correction factors.

Finally, for $\chi_{c1}$ we assume
$\mathcal{R}_{\chi_{c1}}=\mathcal{R}_{\chi_{c2}}$, allowing us to fix the
parameters $\mathcal{N}_{\chi_{c1},T,L}$ via normalization. We find
$\mathcal{N}_{\chi_{c1},T}=1.388$ and $\mathcal{N}_{\chi_{c1},L}=1.401$.

From the above fits we may compute charge radii of the $\chi_{cJ}$ states via
the following formula [82]:

$\left\langle
r_{c}^{2}\right\rangle=N_{c}\sum_{h\bar{h}}\int_{z}\int_{\boldsymbol{r}_{\perp}}\left(A_{\bar{\lambda}}^{\mathcal{H}}\bar{z}^{2}\boldsymbol{r}_{\perp}^{2}\right)\left|\Psi^{\mathcal{H}}_{\bar{\lambda},h\bar{h}}(\boldsymbol{r}_{\perp},z)\right|^{2}\,.$
(58)

$(A_{\bar{\lambda}}^{\mathcal{H}})^{-1}$ denotes the expectation value of
$\boldsymbol{r}_{\perp}^{2}/\boldsymbol{r}^{2}=\sin^{2}\theta$ in the state

$|J,\bar{\lambda}\rangle=\sum_{m,m_{s}}C(J,\bar{\lambda}|L,m;S,m_{s})\,|L,m\rangle\,|S,m_{s}\rangle~{},$
(59)

where $C(J,\bar{\lambda}|L,m;S,m_{s})$ denotes the Clebsch-Gordan coefficient
for $L=S=1$. This gives $A^{\chi_{c0}}=3/2$, $A^{\chi_{c1}}_{0}=5/4$,
$A^{\chi_{c1}}_{1}=5/3$, $A^{\chi_{c2}}_{0}=15/8$, $A^{\chi_{c2}}_{1}=5/3$,
$A^{\chi_{c2}}_{2}=5/4$. For $\chi_{c0}$ we then find $\sqrt{\langle
r^{2}_{c}\rangle}=0.270$ fm, which is close to the value obtained in Ref. [82]
from potential models. For $\chi_{c1}$ we find $\sqrt{\langle
r^{2}_{c}\rangle}=0.277$ fm for $\bar{\lambda}=\pm 1$, while $\sqrt{\langle
r^{2}_{c}\rangle}=0.236$ fm was obtained for $\bar{\lambda}=0$. For
$\chi_{c2}$ we have $\sqrt{\langle r^{2}_{c}\rangle}=0.263$ fm for
$\bar{\lambda}=\pm 2$, $\sqrt{\langle r^{2}_{c}\rangle}=0.262$ fm for
$\bar{\lambda}=\pm 1$ and $\sqrt{\langle r^{2}_{c}\rangle}=0.237$ fm for
$\bar{\lambda}=0$. Thus, all charge radii are of comparable magnitude and less
than 0.3 fm, which appears reasonable.

##  V Numerical results

In this Section we will firstly show the numerical results for the
$\gamma^{*}p\to\chi_{cJ}p$ cross section based on the formulas in Secs. II.2
and III and relevant for the EIC kinematics. Using the
$\gamma^{*}p\to\chi_{cJ}p$ cross sections, we will also compute the electro-
production cross section
$e(k)p(P)\to\chi_{cJ}(\Delta)e(k^{\prime})p(P^{\prime})$ through the standard
formula [83, 84, 85, 86]

$\frac{\mathrm{d}\sigma_{ep}}{\mathrm{d}x_{\mathcal{P}}\mathrm{d}Q^{2}\mathrm{d}|t|}=\frac{\alpha}{2\pi
Q^{2}x_{\mathcal{P}}}\left\\{2(1-y)\frac{\mathrm{d}\sigma_{L}}{\mathrm{d}|t|}+\left(1+(1-y)^{2}-2(1-y)\frac{Q^{2}_{\rm
min}}{Q^{2}}\right)\frac{\mathrm{d}\sigma_{T}}{\mathrm{d}|t|}\right\\}\,,$
(60)

where $Q^{2}_{\rm min}=m_{e}^{2}y^{2}/(1-y)$ accounts for projectile mass
corrections, with $m_{e}$ being the electron mass. $y$ is the inelasticity,
given by

$y\equiv\frac{q\cdot P}{k\cdot
P}=\frac{W^{2}+Q^{2}-m_{N}^{2}}{S-m_{N}^{2}}\,,$ (61)

where $S\equiv(k+P)^{2}$ is the $ep$ collision energy, $W^{2}\equiv(q+P)^{2}$
is the $\gamma^{*}p$ collision energy and $m_{N}$ is the proton mass. Recall
that the Odderon amplitude evolves with Pomeron-$x_{\mathcal{P}}$, which is
related to the conventional Bjorken-$x$ variable $x_{B}=Q^{2}/(2q\cdot
P)=Q^{2}/(W^{2}+Q^{2}-m_{N}^{2})$ via

$x_{\mathcal{P}}\equiv\frac{x_{B}}{\beta}\,,\qquad\beta\equiv\frac{Q^{2}}{2q\cdot\ell}=\frac{Q^{2}}{Q^{2}+M_{\mathcal{H}}^{2}+|t|}\,.$
(62)

In the numerical computations, we augment the Primakoff cross section taking
into account the Pauli form factor $F_{2}$ as explained in Sec. III.2. For
$\chi_{c0}$ and $\chi_{c2}$ we take into account the available QCD corrections
to the Primakoff cross section via (40). The total cross section is based on
taking the coherent sum of the Primakoff and the Odderon exchange amplitudes
in which case their relative sign becomes crucial. Our careful analysis in
App. D reveals that the Primakoff and Odderon amplitudes are in-phase, that
is, they interfere constructively for each value of the Odderon evolution
parameter $x_{\mathcal{P}}$, as we have also explained in Sec. II.4. The
contribution of the $F_{2}$ does not interfere with the coherent sum of the
Odderon and Primakoff amplitudes. The interference between the QCD correction
to the Primakoff amplitude (available for $\chi_{c0,2}$) and the Odderon
exchange is determined through the relative phase of the Odderon and the
Primakoff amplitude at tree level. For the Odderon component we are using the
solutions of the rcBK evolution, keeping only the $k=0$ harmonic of the
Fourier series (26), as explained in Sec. II.4. For the Primakoff component we
use the recent fits of $F_{1}$ and $F_{2}$ from [77]. Whenever it is not
explicitly stated, our results are based on the standard choice for the
numerical value of $\alpha_{S}$: $\alpha_{S}(2m_{c})\approx 0.25$.

Figure 3: The $\gamma^{*}p\to\chi_{cJ}p$ cross sections as functions of $|t|$
for the transverse photon. The label “sum” stands for the coherent sum of the
odderon and the Primakoff contributions. We have set $Q^{2}=0.01$ GeV2.
Different line styles in the odderon and the summed cross section correspond
to different values of $x_{\mathcal{P}}$. The comparison of the Primakoff
contribution with the full and dashed orange lines emphasizes the impact of
$F_{2}$.

We start with the numerical results for the $\gamma^{*}p\to\chi_{cJ}p$ cross
sections shown in Fig. 3 as functions of $|t|$, for different values of
$x_{\mathcal{P}}$. We have set $Q^{2}=0.01$ GeV2. Since $Q^{2}$ is low, we
focus on the transverse cross section $\mathrm{d}\sigma_{T}/\mathrm{d}|t|$. At
small $|t|$ the cross section is, of course, dominated by the Primakoff
process (photon exchange). However, we note that our predictions for the
Primakoff cross sections are substantially lower than those shown in Fig. 4c
of Ref. [87]. This emphasizes the importance of the constraint (40) from the
two photon decay rate.

The QCD Odderon exchange amplitude reaches a comparable magnitude at
$|t|\approx 1$ GeV2, depending on $x_{\mathcal{P}}$ and $J$. The lowest cross-
over from the Primakoff dominated to the Odderon dominated regime at $|t|<1$
GeV2 is seen for $\chi_{c1}$ and $\chi_{c2}$ where the Primakoff-like
background is lower than for $\chi_{c0}$. In the regime of $|t|$ where
individual Primakoff and Odderon cross sections are of similar magnitude,
thanks to their constructive interference, the coherently summed cross section
is four times greater than the Primakoff cross section alone. At high momentum
transfers Odderon exchange dominates due to its slower fall-off with $|t|$.

Figure 4: The $\gamma^{*}p\to\chi_{cJ}p$ cross sections as functions of $W$
for the transverse photon. The label “sum” stands for the coherent sum of the
Odderon and Primakoff contributions. We have set $Q^{2}=0.01$ GeV2 and
integrated in $|t|$ over $0.5<|t|<1.5$ GeV2.

For Fig. 4 we have integrated $\mathrm{d}\sigma_{T}/\mathrm{d}|t|$ over the
range $0.5<|t|<1.5$ GeV2 as appropriate for the EIC design detector [43]. We
show the result as a function of of $W$. As a consequence of small-$x$
evolution, the Odderon cross section drops with increasing $W$. With the
Primakoff cross section alone being constant in our framework, at the lower
end of $W$ the coherent sum is about five times greater than the Primakoff
component alone. Interestingly, $\sigma_{T}$ displays a negative slope due to
the decreasing Odderon amplitude towards smaller $x$.

We find that the decrease of the Odderon cross section with $W$ is driven
mostly by the nonlinear corrections in the unitarized equation for the
Odderon. In App. E (see Fig. 9, right) we have computed $\sigma_{T}$ based on
linear evolution of the BLV Odderon and found a slower decrease with $W$ as
one would anticipate from the asymptotics of the BLV solution with the
intercept equal to one [9].

Figure 5: The $ep\to\chi_{cJ}ep$ ($J=0$, $1$, $2$) cross sections as functions
of $|t|$ at the top EIC energy of $\sqrt{S}=140$ GeV.

We also calculate the differential $ep\to\chi_{cJ}ep$ cross section for top
EIC energy, $\sqrt{S}=140$ GeV, and use the following kinematic cuts:
$0.01<y<0.95$ [43], $x_{\cal P}<0.01$, and $Q^{2}_{\rm min}<Q^{2}<(2\times
3.5\,{\rm GeV})^{2}$. The resulting $|t|$-dependence of the cross section is
shown in Fig. 5. With these cuts the crossover from photon to Odderon exchange
occurs at about $|t|\approx 1$ GeV2 where the total cross section is several
times greater than due to the Primakoff process alone.

Figure 6: Total electroproduction cross section of the $\chi_{cJ}$ quarkonia
as a function of the $ep$ center of mass energy $\sqrt{S}$. Kinematic cuts are
explained in the text.

The total $\chi_{cJ}$ cross sections are shown in Fig. 6. They have been
integrated over $|t|$ over the range $0.5$ GeV${}^{2}<|t|<1.5$ GeV2 where
according to Fig. 5 the Odderon contribution is appreciable. Both the photon
and Odderon exchange contributions level off towards top EIC energy since our
kinematic cuts are energy independent, and neither exchange involves a
positive intercept.

$\chi_{c1}$ has the highest branching ratio ${\rm BR}(\chi_{c1}\to
J/\psi\gamma)=34.3\%$ [33]. As an example, lets estimate the total number of
$\chi_{c1}$’s per month at the EIC. From Fig. 6 the total cross section is
about $\sigma_{ep}\approx 20$ fb at $\sqrt{S}=140$ GeV. Multiplying by the
expected luminosity at the EIC ($L=10^{34}$ cm-2s${}^{-1}=10^{-5}$ fb-1 s-1)
gives about $2\times 10^{-4}$ events/second or $\approx 518$ events/month.
After taking into account ${\rm BR}(\chi_{c1}\to J/\psi\gamma)$ results in
about 177 events/month. $J/\psi$’s are detected through the $J/\psi\to
e^{+}e^{-}$ or $\to\mu^{+}\mu^{-}$ decays, and the combined corresponding
branching ratio is about ${\rm BR}(J/\psi\to\ell^{+}\ell^{-})=12\%$, after
which we end up with about 21 events/month.

Figure 7: A summary of expected Primakoff (orange) vs. total (magenta) number
of exclusive $\chi_{cJ}\to J/\psi\gamma\to l^{+}l^{-}\gamma$ events per month
at the top EIC energy of $\sqrt{S}=140$ GeV, for three different $|t|$ cuts.
The remaining kinematic cuts are the same as on Fig. 6. The two different bins
for the total number of events correspond to $\alpha_{S}=0.25$ and
$\alpha_{S}=0.2$, respectively.

In the final Fig. 7, we summarize the expected number of exclusive $\chi_{cJ}$
events per month, $N_{\chi_{cJ}}$, at the EIC design luminosity for two ranges
of momentum transfer $|t|$. The result is obtained according to:

$N_{\chi_{cJ}}=L\times\sigma_{ep}(ep\to\chi_{cJ}ep)\times{\rm BR}(\chi_{cJ}\to
J/\psi\gamma)\times{\rm BR}(J/\psi\to\ell^{+}\ell^{-})\,.$ (63)

The most statistics is expected for $\chi_{c1}$ and $\chi_{c2}$ tensor
quarkonia with about a factor of 2-3 excess events over the Primakoff process
in the interval $0.5$ GeV${}^{2}<|t|<1.5$ GeV2. For $\chi_{c0}$, which has the
largest cross section from Fig. 6 is compensated by a small branching ratio
${\rm BR}(\chi_{c0}\to J/\psi\gamma)=1.4\%$ [33].

The presented estimates of $\chi_{cJ}$ production cross sections carry some
theoretical uncertainties. For the photon exchange contribution they are
smaller than for the Odderon exchange. The coupling of the photon to the
proton is well constrained experimentally [77], and so the main sources of
uncertainty are the $\gamma^{*}\gamma^{*}\to\chi_{cJ}$ amplitudes. They depend
on the details of the quarkonia wave functions and are sensitive to unknown
higher order QCD corrections. For $\chi_{c0}$ and $\chi_{c2}$ the differential
cross section $\mathrm{d}\sigma/\mathrm{d}|t|$ for the Primakoff process obeys
a stringent constraint at small $|t|\ll 1$ GeV2 imposed by Eq. (40) that holds
at all orders in QCD. Hence, for these charmonia the uncertainties mentioned
above affect mostly the details of the $t$-shape of
$\mathrm{d}\sigma/\mathrm{d}|t|$. They are expected to be small as the
$\chi_{cJ}$ wave functions are probed mostly at short distances $\sim
1/m_{c}$, where they are well constrained by the $\gamma\gamma$ decay width.
The value of the coupling $\alpha_{S}(2m_{c})\simeq 0.25$ is not large and the
$t\to 0$ limit of the QCD corrections is known, so the uncertainties from QCD
corrections at $t\neq 0$ should be small as well. For $\chi_{c1}$ the
theoretical uncertainty of the Primakoff cross section is larger than for
$\chi_{c0}$ and $\chi_{c2}$, as the constraint (40) cannot be imposed due to
the LY theorem and, in consequence, the vanishing $\gamma\gamma$ decay width
of $\chi_{c1}$. In this case measurements of the differential cross sections
$\mathrm{d}\sigma/\mathrm{d}|t|$ for $|t|\ll 1$ GeV2 where the photon exchange
dominates should greatly reduce the uncertainty associated with this
contribution at $|t|\sim 1$ GeV2, where we expect to isolate a significant
Odderon signal.

In the Odderon exchange, in addition to uncertainties from unknown details of
the quarkonia wave function and from higher order QCD corrections, there are
uncertainties associated with the model for the proton wave function, and from
the value of $\alpha_{S}$. The model of the proton employed here obeys general
constraints coming from the measurements of the proton size, exclusive
$J/\psi$ production off the proton [32] and open charm electroproduction at
HERA [64]. Those constraints, however, originate from measurements in the
$C$-even sector, and the emerging $C$-odd correlators have never been probed
experimentally. Furthermore, the perturbative Odderon is strongly sensitive to
the numerical value of $\alpha_{S}$ – the amplitude being $O(\alpha_{S}^{3})$
at the lowest order. To quantify the uncertainty in the value of $\alpha_{S}$
we take an interval $0.2<\alpha_{S}<0.25$. The lower value $\alpha_{S}=0.2$
comes from a fit of the proton and $C$-even exchange models to open charm
electroproduction at HERA [64]. Therein the average value of the hard scale
$\mu=\sqrt{Q^{2}+4m_{c}^{2}}$ is about $5$ GeV, so the obtained value
$\alpha_{S}=0.2$ is fully consistent with the running of $\alpha_{S}(\mu)$. In
our case the bulk of exclusive $\chi_{cJ}$ production in $ep$ collisions
occurs at small photon virtualities and moderate $|t|$, and so
$\alpha_{S}(2m_{c})\approx 0.25$ is more appropriate. The range
$0.2<\alpha_{S}<0.25$ results in an uncertainty of about a factor of 2.5
between the minimal and maximal values of the Odderon amplitude, when NLO
effects in the proton wave function are included.

For the combined Primakoff and Odderon contributions, the $\alpha_{S}$-induced
uncertainty is negligible for $|t|\ll 1$ GeV2, where the Primakoff process
dominates; whereas it about a factor of $2$ when Primakoff and Odderon (with
$\alpha_{S}=0.25$) amplitudes are close to each other, and up to a factor of
about $6$ when the Odderon strongly dominates over the Primakoff channel.
Since the uncertainty coming from $\alpha_{S}$ is sizeable, we assume that it
is the dominant uncertainty for the Odderon exchange and take it as our
estimate of the theoretical uncertainty. We consider the choice of
$\alpha_{S}=0.25$ to be realistic, as the scale $2m_{c}$ is already fairly
high. Another choice is $\mu=m_{c}$, leading to $\alpha_{S}\approx 0.35$, as
obtained for exclusive $J/\psi$ photoproduction using the same model for the
proton at leading order [32]. Thus, we consider the lower value of
$\alpha_{S}=0.2$ to be a conservative choice.

##  VI Summary and Conclusions

In this paper we have derived amplitudes for exclusive scalar, axial-vector
and tensor quarkonium production, Eqs. (22, 24), (28, 27), and (30, 29),
respectively, in electron-proton scattering. We have provided first estimates
of the cross sections for exclusive production of $\chi_{cJ}$ quarkonia with
positive $C$-parity at the EIC. This process requires a $C$-odd exchange in
the $t$-channel. In the limit of heavy quarks and/or high transverse momentum
transfer this could be the exchange of a photon, a Primakoff process, or the
exchange of a color-symmetric three gluon ladder, the Odderon. Our estimates
suggest that for $|t|\approx 1$ GeV2 the two amplitudes are of similar
magnitude and that there are strong interference effects. Importantly, the
relative phase is not affected by QCD evolution of the Odderon towards small
$x$, and so it is determined by the three gluon exchange amplitude at
moderately small $x_{0}$. In turn, the matrix element of the eikonal color
current operator $d^{abc}J^{+a}J^{+b}J^{+c}$ at $x_{0}$, for transverse
momenta of order $\mathrm{max}(R^{-1}_{\chi_{c}},|t|^{1/2})$ ($R_{\chi_{c}}$
being here roughly the size of $\chi_{c}$), can be computed in a truncated
Fock space for the proton that encompasses the states that are relevant in
that kinematic regime.

We find that photon and Odderon exchanges interfere constructively, leading to
an enhancement of the differential cross section for $\chi_{cJ}$ production
around $|t|\approx 1$ GeV2 by up to a factor of 4 over the pure photon
exchange contribution. Given that both the normalization and the
$t$-dependence of the Primakoff process are reasonably well determined, this
presents a very exciting opportunity to potentially discover at the EIC the
hard $C$-odd exchange predicted by QCD. Furthermore, we find that towards top
EIC energy, $\sqrt{S}=140$ GeV, the total electroproduction cross section of
$\chi_{cJ}$ quarkonia (with kinematic cuts specified in the previous section)
levels off as neither photon nor Odderon exchange involve a positive
intercept. Therefore, it would be important to measure the energy dependence
[88] of the cross section from $\sqrt{S}=50$ GeV up to top EIC energy
$\sqrt{S}=140$ GeV, where photon and Odderon amplitudes are similar, and where
constructive interference of amplitudes leads to a total cross section for
production of $\chi_{c1}$ and $\chi_{c2}$ quarkonia which exceeds the
Primakoff component by a factor of $\approx 3-4$.

The predictions for the Odderon exchange process of course involve a number of
uncertainties such as the matrix element in the proton of
$d^{abc}J^{+a}J^{+b}J^{+c}$ or the value of the strong coupling $\alpha_{S}$.
The associated cross section scales approximately like $\alpha_{S}^{6}$. Most
importantly, the discovery of the hard QCD $C$-odd exchange requires fixing
the normalization of the Primakoff background from measurements at low
$\sqrt{|t|}<0.5$ GeV; the relation (40) of the $t\to 0$ limit of the cross
section to the $\gamma\gamma$ decay width provides an important constraint on
theoretical predictions for this background. The $t$-dependence of the
Primakoff cross section is then determined by the Dirac and Pauli
electromagnetic form factors of the proton which are well known. Hence, one
may then look for excess events above this background, and for a change of
slope of the differential cross section $\mathrm{d}\sigma/\mathrm{d}t$, at
higher $|t|\sim 1$ GeV2 and beyond.

Our estimates indicate very weak dipole-proton (hard) scattering with $C$-odd
exchange; this is largely since these amplitudes are also parity odd and
vanish for impact parameter $\boldsymbol{b}_{\perp}\to 0$, unlike parity even
amplitudes, as well as due to the fact that they fall off more rapidly towards
large $\boldsymbol{b}_{\perp}$. The total electro-production cross section of
$\chi_{cJ}$ quarkonia at top EIC energy is estimated to have a magnitude
$10-20$ fb (for $J=1$) and $60$ fb (for $J=2$). Thus, it was not possible to
observe these processes at the HERA accelerator. However, in view of the
projected high luminosity of the EIC, data collection over a time span of
several months to a year may be sufficient for the discovery of the hard
Odderon. A promising alternative would be to allow low-mass excitations of the
proton, while requiring a large rapidity gap to the $C=+1$ quarkonia. Such
rapidity gap, diffractive processes have greater cross-sections than exclusive
ones. Furthermore, they would extend the reach to higher $|t|$ where Odderon
exchange would more clearly dominate over photon exchange.

###### Acknowledgements.

We thank Y. Hatta, H. Mäntysaari, L. Pentchev, and R. Venugopalan for useful
discussions. The work of S. B. and A. K. is supported by the Croatian Science
Foundation (HRZZ) no. 5332 (UIP-2019-04). A. D. acknowledges support by the
DOE Office of Nuclear Physics through Grant DE-SC0002307, and The City
University of New York for PSC-CUNY Research grant 65079-00 53. He also thanks
the EIC Theory Institute at Brookhaven National Laboratory for their
hospitality during a visit in July 2023 when initial work on this project was
done. L. M. acknowledges the support of the Polish National Science Centre
(NCN) grant no. 2017/27/B/ST2/02755. T. S. kindly acknowledges the support of
the Polish National Science Center (NCN) Grants No. 2019/32/C/ST2/00202 and
2021/43/D/ST2/03375. S. B. thanks for the hospitality at the Jagiellonian
University where part of this work was done.

##  Appendix A Computation of the light-cone wave functions

In order to compute the light-cone wave functions, the key element is a vertex
contraction between spinors, e. g. $\bar{u}_{h}(k)\Gamma
v_{\bar{h}}(k^{\prime})$ ($\Gamma$ is some general Dirac vertex). For the
spinors we use the LB basis, defined through

$\begin{split}&u_{h}(k)=\frac{1}{2^{1/4}\sqrt{k^{-}}}\left(\sqrt{2}k^{-}+\gamma^{0}m+\boldsymbol{\alpha}_{\perp}\cdot\boldsymbol{k}_{\perp}\right)\chi_{h}=\frac{1}{2^{1/4}\sqrt{k^{-}}}\left(\sqrt{2}k^{-}+\gamma^{0}m+\not{\boldsymbol{k}}_{\perp}\gamma^{0}\right)\chi_{h}\,,\\\
&v_{h}(k)=\frac{1}{2^{1/4}\sqrt{k^{-}}}\left(\sqrt{2}k^{-}-\gamma^{0}m+\boldsymbol{\alpha}_{\perp}\cdot\boldsymbol{k}_{\perp}\right)\chi_{-h}=\frac{1}{2^{1/4}\sqrt{k^{-}}}\left(\sqrt{2}k^{-}-\gamma^{0}m+\not{\boldsymbol{k}}_{\perp}\gamma^{0}\right)\chi_{-h}\,,\end{split}$
(64)

where $\alpha^{i}=\gamma^{0}\gamma^{i}$. In accordance with the conventions
used in this work, we have replaced plus and minus light-cone coordinates in
the above expression with respect to the original LB convention [46]. Thus the
spinors $\chi_{h}$ are eigenstates of $\gamma^{+}\gamma^{-}$, namely

$\frac{1}{2}\gamma^{+}\gamma^{-}\chi_{h}=\chi_{h}\,,\qquad\frac{1}{2}\gamma^{-}\gamma^{+}\chi_{h}=0\,.$
(65)

Explicitly, we have $\chi_{+}=(-1,0,1,0)/\sqrt{2}$,
$\chi_{-}=(0,1,0,1)/\sqrt{2}$. The spinors (64) can now be written as

$\begin{split}&u_{h}(k)=\frac{1}{2^{1/4}\sqrt{k^{-}}}(\not{k}+m)\gamma^{0}\chi_{h}\,,\\\
&v_{h}(k)=\frac{1}{2^{1/4}\sqrt{k^{-}}}(\not{k}-m)\gamma^{0}\chi_{-h}\,.\end{split}$
(66)

By defining a projection matrix

$\chi_{h\bar{h}}\equiv\gamma^{0}\chi_{-\bar{h}}\bar{\chi}_{h}\gamma^{0}\,,$
(67)

one obtains

$\bar{u}_{h}(k)\Gamma v_{\bar{h}}(k^{\prime})={\rm
tr}\left[(\not{k}+m)\Gamma(\not{k}^{\prime}-m)\chi_{h\bar{h}}\right]\,,$ (68)

and so the computation of the light-cone wave function comes down to the
computation of the above Dirac trace. We use this method to compute all the
wave functions considered in this work. Using (67), the explicit forms of the
projection matrices are found to be

$\begin{split}&\chi_{\pm\pm}=\pm\frac{1}{2\sqrt{2}}\gamma^{-}(\gamma^{1}\mp\mathrm{i}\gamma^{2})\,,\\\
&\chi_{\pm\mp}=\frac{1}{2\sqrt{2}}\gamma^{-}(1\mp\gamma_{5})\,.\end{split}$
(69)

##  Appendix B Derivation of the amplitudes for axial and tensor quarkonia

In this appendix we present the steps of the derivation of the amplitudes from
Sec. II.2 for the axial vector and tensor quarkonia.

In the case of axial quarkonia we plug the second line of (9) into the traces
in (19) to obtain

$\begin{split}&A_{\lambda=0,\bar{\lambda}=0}(\boldsymbol{l}_{\perp},\boldsymbol{l}_{1\perp},z)=0\,,\\\
&A_{\lambda=0,\bar{\lambda}=\pm}(\boldsymbol{l}_{\perp},\boldsymbol{l}_{1\perp},z)=4Qz\bar{z}(\boldsymbol{\epsilon}_{\perp}^{\bar{\lambda}*}\times\boldsymbol{l}_{1\perp})\,,\\\
&A_{\lambda=\pm,\bar{\lambda}=0}(\boldsymbol{l}_{\perp},\boldsymbol{l}_{1\perp},z)=\frac{4}{M_{\mathcal{A}}}\left(m_{c}^{2}(\boldsymbol{\epsilon}_{\perp}^{\lambda}\times\boldsymbol{l}_{1\perp})+\boldsymbol{l}_{1\perp}^{2}(\boldsymbol{\epsilon}_{\perp}^{\lambda}\times\boldsymbol{l}_{\perp})\right)\,,\\\
&A_{\lambda=\pm,\bar{\lambda}=\pm}(\boldsymbol{l}_{\perp},\boldsymbol{l}_{1\perp},z)=2(z-\bar{z})\left((\boldsymbol{\epsilon}_{\perp}^{\lambda}\cdot\boldsymbol{\epsilon}_{\perp}^{\bar{\lambda}*})(\boldsymbol{l}_{\perp}\times\boldsymbol{l}_{1\perp})-(\boldsymbol{l}_{\perp}\cdot\boldsymbol{l}_{1\perp}+m_{c}^{2})(\boldsymbol{\epsilon}_{\perp}^{\lambda}\times\boldsymbol{\epsilon}_{\perp}^{\bar{\lambda}*})\right)\,,\end{split}$
(70)

The last line is understood to contain four different combinations of photon
and quarkonia transverse polarizations. We used the following notations for
the 2D cross products:
$\boldsymbol{x}_{\perp}\times\boldsymbol{y}_{\perp}\equiv
x^{1}y^{2}-x^{2}y^{1}=-\epsilon^{+-\boldsymbol{x}_{\perp}\boldsymbol{y}_{\perp}}$.
The non-zero reduced amplitudes are

$\begin{split}&\mathcal{A}_{\lambda=0,\bar{\lambda}=\pm}(\boldsymbol{r}_{\perp},\boldsymbol{\Delta}_{\perp})=eq_{c}\mathrm{e}^{-\mathrm{i}\bar{\lambda}\phi_{r}}\int_{z}\mathrm{e}^{-\mathrm{i}\boldsymbol{\delta}_{\perp}\cdot\boldsymbol{r}_{\perp}}\mathcal{A}_{LT}(r_{\perp})\,,\\\
&\mathcal{A}_{\lambda=\pm,\bar{\lambda}=0}(\boldsymbol{r}_{\perp},\boldsymbol{\Delta}_{\perp})=eq_{c}\mathrm{e}^{\mathrm{i}\lambda\phi_{r}}\int_{z}\mathrm{e}^{-\mathrm{i}\boldsymbol{\delta}_{\perp}\cdot\boldsymbol{r}_{\perp}}\mathcal{A}_{TL}(r_{\perp})\,,\\\
&\mathcal{A}_{\lambda=\pm,\bar{\lambda}=\pm}(\boldsymbol{r}_{\perp},\boldsymbol{\Delta}_{\perp})=eq_{c}\lambda\delta_{\lambda\bar{\lambda}}\int_{z}\mathrm{e}^{-\mathrm{i}\boldsymbol{\delta}_{\perp}\cdot\boldsymbol{r}_{\perp}}\mathcal{A}_{TT}(r_{\perp})\,,\\\
\end{split}$ (71)

where $\mathcal{A}_{A}$ are given in Eqs. (28). Notice that the first term
from the fourth line of (70) is proportional to
$\boldsymbol{r}_{\perp}\times\boldsymbol{r}_{\perp}=0$ in coordinate space and
so it does not contribute to $\mathcal{A}_{TT}(r_{\perp})$. In the last line
of Eq. (71), out of four possible combinations of photon and axial quarkonia
transverse polarizations, only the two polarization preserving transitions
survive since this amplitude is proportional to the 2D cross product,
$\boldsymbol{\epsilon}_{\perp}^{\lambda}\times\boldsymbol{\epsilon}_{\perp}^{\bar{\lambda}*}=(-\mathrm{i}/2)(\lambda+\bar{\lambda})=-\mathrm{i}\lambda\delta_{\lambda\bar{\lambda}}$
of the polarization vectors of the incoming photon and outgoing axial vector
quarkonia. Accordingly, the two polarization flipping transitions vanish,
unlike in case of vector quarkonia.

Computing now the $\phi_{r}$ integral we separate the helicity dependence and
obtain the amplitudes

$\begin{split}&\left\langle\mathcal{M}_{\lambda=0,\bar{\lambda}=\pm
1}\right\rangle\equiv q^{-}\widetilde{\mathcal{M}}_{0,\bar{\lambda}=\pm
1}\equiv
q^{-}\mathrm{e}^{-\mathrm{i}\bar{\lambda}\phi_{\Delta}}\widetilde{\mathcal{M}}_{LT}\,,\\\
&\left\langle\mathcal{M}_{\lambda=\pm 1,\lambda=0}\right\rangle\equiv
q^{-}\widetilde{\mathcal{M}}_{\lambda=\pm 1,0}\equiv
q^{-}\mathrm{e}^{\mathrm{i}\lambda\phi_{\Delta}}\widetilde{\mathcal{M}}_{TL}\,,\\\
&\left\langle\mathcal{M}_{\lambda=\pm 1,\bar{\lambda}=\pm
1}\right\rangle\equiv q^{-}\widetilde{\mathcal{M}}_{\lambda=\pm
1,\bar{\lambda}=\pm 1}\equiv
q^{-}\lambda\delta_{\lambda\bar{\lambda}}\widetilde{\mathcal{M}}_{TT}\,,\end{split}$
(72)

where $\widetilde{\mathcal{M}}_{B}$ are given in Eqs. (27).

For the tensor quarkonia, we proceed along similar steps as above, only this
time there is a total of 15 overlaps that need to be computed. Starting from
the traces (19), we plug in the third line of (9) and find

$\begin{split}A_{\lambda=0,\bar{\lambda}=\pm
2}(\boldsymbol{l}_{\perp},\boldsymbol{l}_{1\perp},z)&=-8Qz\bar{z}(z-\bar{z})(\boldsymbol{\epsilon}_{\perp}^{\pm*}\cdot\boldsymbol{l}_{1\perp})^{2}\,,\\\
A_{\lambda=\pm 1,\bar{\lambda}=\pm
2}(\boldsymbol{l}_{\perp},\boldsymbol{l}_{1\perp},z)&=-4(\boldsymbol{\epsilon}_{\perp}^{\pm*}\cdot\boldsymbol{l}_{1\perp})\big{[}(\boldsymbol{\epsilon}_{\perp}^{\lambda}\cdot\boldsymbol{\epsilon}_{\perp}^{\pm*})(\boldsymbol{l}_{\perp}\cdot\boldsymbol{l}_{1\perp}+m_{c}^{2})\\\
&+(z-\bar{z})^{2}(\boldsymbol{\epsilon}_{\perp}^{\lambda}\cdot\boldsymbol{l}_{\perp})(\boldsymbol{\epsilon}_{\perp}^{\pm*}\cdot\boldsymbol{l}_{1\perp})-(\boldsymbol{\epsilon}_{\perp}^{\lambda}\cdot\boldsymbol{l}_{1\perp})(\boldsymbol{\epsilon}_{\perp}^{\pm*}\cdot\boldsymbol{l}_{\perp})\big{]}\,,\\\
A_{\lambda=0,\bar{\lambda}=\pm
1}(\boldsymbol{l}_{\perp},\boldsymbol{l}_{1\perp},z)&=2\sqrt{2}QM_{\mathcal{T}}z\bar{z}(3-4(z^{2}+\bar{z}^{2}))(\boldsymbol{\epsilon}_{\perp}^{\bar{\lambda}*}\cdot\boldsymbol{l}_{1\perp})\,,\\\
A_{\lambda=\pm 1,\bar{\lambda}=\pm
1}(\boldsymbol{l}_{\perp},\boldsymbol{l}_{1\perp},z)&=-\sqrt{2}M_{\mathcal{T}}(z-\bar{z})\big{[}(\boldsymbol{\epsilon}_{\perp}^{\lambda}\cdot\boldsymbol{\epsilon}_{\perp}^{\bar{\lambda}*})(\boldsymbol{l}_{\perp}\cdot\boldsymbol{l}_{1\perp}+m_{c}^{2})\\\
&-(3-4(z^{2}+\bar{z}^{2}))(\boldsymbol{\epsilon}_{\perp}^{\lambda}\cdot\boldsymbol{l}_{\perp})(\boldsymbol{\epsilon}_{\perp}^{\bar{\lambda}*}\cdot\boldsymbol{l}_{1\perp})-(\boldsymbol{\epsilon}_{\perp}^{\lambda}\cdot\boldsymbol{l}_{1\perp})(\boldsymbol{\epsilon}_{\perp}^{\bar{\lambda}*}\cdot\boldsymbol{l}_{\perp})\big{]}\,,\\\
A_{\lambda=0,\bar{\lambda}=0}(\boldsymbol{l}_{\perp},\boldsymbol{l}_{1\perp},z)&=\frac{8\sqrt{2}Q}{\sqrt{3}}z\bar{z}(z-\bar{z})\left(\frac{3}{2}\boldsymbol{l}_{1\perp}^{2}+m_{c}^{2}\right)\,,\\\
A_{\lambda=\pm
1,\bar{\lambda}=0}(\boldsymbol{l}_{\perp},\boldsymbol{l}_{1\perp},z)&=\frac{2\sqrt{2}}{\sqrt{3}}\Bigg{[}2(z-\bar{z})^{2}(\boldsymbol{\epsilon}_{\perp}^{\lambda}\cdot\boldsymbol{l}_{\perp})\left(\frac{3}{2}\boldsymbol{l}_{1\perp}^{2}+m_{c}^{2}\right)+(\boldsymbol{\epsilon}_{\perp}^{\lambda}\cdot\boldsymbol{l}_{1\perp})m_{c}^{2}\Bigg{]}\,.\end{split}$
(73)

The results for $\bar{\lambda}=0$ have been simplified using

$(\boldsymbol{\epsilon}_{\perp}^{+*}\cdot\boldsymbol{u}_{\perp})(\boldsymbol{\epsilon}_{\perp}^{-*}\cdot\boldsymbol{v}_{\perp})+(\boldsymbol{\epsilon}_{\perp}^{-*}\cdot\boldsymbol{u}_{\perp})(\boldsymbol{\epsilon}_{\perp}^{+*}\cdot\boldsymbol{v}_{\perp})=-\boldsymbol{u}_{\perp}\cdot\boldsymbol{v}_{\perp}\,,\\\
$ (74)

which holds for general 2D vectors $\boldsymbol{u}_{\perp}$ and
$\boldsymbol{v}_{\perp}$. Inserting (73) in (18) we obtain the reduced
amplitudes

$\begin{split}&\mathcal{A}_{\lambda=0,\bar{\lambda}=\pm
2}(\boldsymbol{r}_{\perp},\boldsymbol{\Delta}_{\perp})=eq_{c}\mathrm{e}^{-\mathrm{i}\bar{\lambda}\phi_{r}}\int_{z}\mathrm{e}^{-\mathrm{i}\boldsymbol{\delta}_{\perp}\cdot\boldsymbol{r}_{\perp}}\mathcal{A}_{LT2}(r_{\perp})\,,\\\
&\mathcal{A}_{\lambda=\pm 1,\bar{\lambda}=\pm
2}(\boldsymbol{r}_{\perp},\boldsymbol{\Delta}_{\perp})=eq_{c}\lambda\mathrm{e}^{\mathrm{i}(\lambda-\bar{\lambda})\phi_{r}}\int_{z}\mathrm{e}^{-\mathrm{i}\boldsymbol{\delta}_{\perp}\cdot\boldsymbol{r}_{\perp}}\left(\delta_{\lambda,\bar{\lambda}/2}\mathcal{A}_{TT2p}(r_{\perp})+\delta_{\lambda,-\bar{\lambda}/2}\mathcal{A}_{TT2f}(r_{\perp})\right)\,,\\\
&\mathcal{A}_{\lambda=0,\bar{\lambda}=\pm
1}(\boldsymbol{r}_{\perp},\boldsymbol{\Delta}_{\perp})=eq_{c}\bar{\lambda}\mathrm{e}^{-\mathrm{i}\bar{\lambda}\phi_{r}}\int_{z}\mathrm{e}^{-\mathrm{i}\boldsymbol{\delta}_{\perp}\cdot\boldsymbol{r}_{\perp}}\mathcal{A}_{LT}(r_{\perp})\,,\\\
&\mathcal{A}_{\lambda=\pm 1,\bar{\lambda}=\pm
1}(\boldsymbol{r}_{\perp},\boldsymbol{\Delta}_{\perp})=eq_{c}\mathrm{e}^{\mathrm{i}(\lambda-\bar{\lambda})\phi_{r}}\int_{z}\mathrm{e}^{-\mathrm{i}\boldsymbol{\delta}_{\perp}\cdot\boldsymbol{r}_{\perp}}\left(\delta_{\lambda\bar{\lambda}}\mathcal{A}_{TTp}(r_{\perp})-\delta_{\lambda,-\bar{\lambda}}\mathcal{A}_{TTf}(r_{\perp})\right)\,,\\\
&\mathcal{A}_{\lambda=0,\bar{\lambda}=0}(\boldsymbol{r}_{\perp},\boldsymbol{\Delta}_{\perp})=eq_{c}\int_{z}\mathrm{e}^{-\mathrm{i}\boldsymbol{\delta}_{\perp}\cdot\boldsymbol{r}_{\perp}}\mathcal{A}_{LL}(r_{\perp})\,,\\\
&\mathcal{A}_{\lambda=\pm
1,\bar{\lambda}=0}(\boldsymbol{r}_{\perp},\boldsymbol{\Delta}_{\perp})=eq_{c}{\lambda\mathrm{e}^{\mathrm{i}\lambda\phi_{r}}}\int_{z}\mathrm{e}^{-\mathrm{i}\boldsymbol{\delta}_{\perp}\cdot\boldsymbol{r}_{\perp}}\mathcal{A}_{TL}(r_{\perp})\,.\\\
\end{split}$ (75)

From (75) we see that for the transverse polarizations of the photon and the
tensor quarkonia both the polarization preserving
($\lambda\to\bar{\lambda}/2$, for the case $\bar{\lambda}=\pm 2$ and
$\lambda\to\bar{\lambda}$ for the case $\bar{\lambda}=\pm 1$) and polarization
flipped ($\lambda\to-\bar{\lambda}/2$, $\lambda\to-\bar{\lambda}$) transitions
are allowed, which explains the notation ($p=$ preserving, $f$ = flipped).

In the final step we compute the amplitudes

$\begin{split}&\left\langle\mathcal{M}_{\lambda=0,\bar{\lambda}=\pm
2}\right\rangle\equiv
q^{-}\widetilde{\mathcal{M}}_{\lambda=0,\bar{\lambda}=\pm 2}\equiv
q^{-}\mathrm{e}^{-\mathrm{i}\bar{\lambda}\phi_{\Delta}}\widetilde{\mathcal{M}}_{LT2}\,,\\\
&\left\langle\mathcal{M}_{\lambda=\pm 1,\bar{\lambda}=\pm
2}\right\rangle\equiv q^{-}\widetilde{\mathcal{M}}_{\lambda=\pm
1,\bar{\lambda}=\pm 2}\equiv
q^{-}\lambda\mathrm{e}^{-\mathrm{i}\lambda\phi_{\Delta}}\delta_{\lambda,\bar{\lambda}/2}\widetilde{\mathcal{M}}_{TT2p}+q^{-}\lambda\mathrm{e}^{3\mathrm{i}\lambda\phi_{\Delta}}\delta_{\lambda,-\bar{\lambda}/2}\widetilde{\mathcal{M}}_{TT2f}\,,\\\
&\left\langle\mathcal{M}_{\lambda=0,\bar{\lambda}=\pm 1}\right\rangle\equiv
q^{-}\widetilde{\mathcal{M}}_{\lambda=0,\bar{\lambda}=\pm 1}\equiv
q^{-}\mathrm{e}^{-\mathrm{i}\bar{\lambda}\phi_{\Delta}}\widetilde{\mathcal{M}}_{LT}\,,\\\
&\left\langle\mathcal{M}_{\lambda=\pm 1,\bar{\lambda}=\pm
1}\right\rangle\equiv q^{-}\widetilde{\mathcal{M}}_{\lambda=\pm
1,\bar{\lambda}=\pm 1}\equiv
q^{-}\delta_{\lambda\bar{\lambda}}\widetilde{\mathcal{M}}_{TTp}+q^{-}\delta_{\lambda,-\bar{\lambda}}\mathrm{e}^{2\mathrm{i}\lambda\phi_{\Delta}}\widetilde{\mathcal{M}}_{TTf}\,,\\\
&\left\langle\mathcal{M}_{\lambda=0,\bar{\lambda}=0}\right\rangle\equiv
q^{-}\widetilde{\mathcal{M}}_{\lambda=0,\bar{\lambda}=0}\equiv
q^{-}\widetilde{\mathcal{M}}_{LL}\,,\\\ &\left\langle\mathcal{M}_{\lambda=\pm
1,\bar{\lambda}=0}\right\rangle\equiv
q^{-}\widetilde{\mathcal{M}}_{\lambda=\pm 1,\bar{\lambda}=0}\equiv
q^{-}(\delta_{\lambda
1}-\delta_{\lambda,-1})\mathrm{e}^{\mathrm{i}\lambda\phi_{\Delta}}\widetilde{\mathcal{M}}_{TL}\,,\end{split}$
(76)

where also the explicit helicity depndence was factored out. The scalar
functions $\widetilde{\mathcal{M}}_{B}$ are given by Eqs. (29).

##  Appendix C The Primakoff contribution in specific kinematic limits

###  C.1 The proof of formula (40)

We begin by considering the amplitude for the decay process
$\mathcal{H}(\bar{\lambda},\Delta)\to\gamma(\lambda,q)\gamma(\lambda^{\prime},\ell)$,
where $\mathcal{H}=\mathcal{S}$, $\mathcal{T}$. In what follows, it will be
useful to separate out the exchanged photon polarization from the amplitude as

$\mathcal{M}_{\lambda\lambda^{\prime}\bar{\lambda}}(\mathcal{H}\to\gamma\gamma)=\epsilon_{\mu}(\lambda^{\prime},\ell)\mathcal{M}^{\mu}_{\lambda\bar{\lambda}}(\mathcal{H}\to\gamma\gamma)\,.$
(77)

Inserting (77) into the standard formula for the decay width [33] we get

$\begin{split}\Gamma(\mathcal{H}\to\gamma\gamma)&=\frac{1}{2J+1}\frac{1}{32\pi
M_{\mathcal{H}}}\sum_{\bar{\lambda}=-J}^{J}\,\sum_{\lambda\lambda^{\prime}}|\mathcal{M}_{\lambda\lambda^{\prime}\bar{\lambda}}(\mathcal{H}\to\gamma\gamma)|^{2}\,,\\\
&=-\frac{1}{2J+1}\frac{1}{32\pi
M_{\mathcal{H}}}\sum_{\bar{\lambda}=-J}^{J}\sum_{\lambda}g_{ii^{\prime}}\mathcal{M}^{i}_{\lambda\bar{\lambda}}(\mathcal{H}\to\gamma\gamma)\mathcal{M}^{i^{\prime}*}_{\lambda\bar{\lambda}}(\mathcal{H}\to\gamma\gamma)\,.\end{split}$
(78)

Next, we consider the amplitude for exclusive photoproduction
$\mathcal{M}_{\lambda\bar{\lambda}}(\gamma p\to\mathcal{H}p)$, see (46). The
cross section reads [33]

$\frac{\mathrm{d}\sigma_{\mathcal{H}}(\gamma
p\to\mathcal{H}p)}{\mathrm{d}|t|}=\frac{1}{2}\sum_{\bar{\lambda}=-J}^{J}\sum_{\lambda}\int\frac{\mathrm{d}\phi_{\ell}}{2\pi}\frac{1}{16\pi
W^{4}}\left|\mathcal{M}_{\lambda\bar{\lambda}}(\gamma
p\to\mathcal{H}p)\right|^{2}\,.$ (79)

where $W^{2}=(q+P)^{2}\simeq 2q\cdot P\gg
M_{\mathcal{H}}^{2},\,P^{2},\,|\ell^{2}|$ in the high energy limit121212The
amplitudes $\mathcal{M}_{\lambda\bar{\lambda}}$ in this section and the color
averaged amplitudes $\widetilde{\mathcal{M}}_{\lambda\bar{\lambda}}$ as
obtained through Eqs. (23), (72) and (76) are related through
$\widetilde{\mathcal{M}}_{\lambda\bar{\lambda}}=\mathcal{M}_{\lambda\bar{\lambda}}/W^{2}$,
compare (79) with (31).. Moreover
$\ell\approx(x_{\mathcal{P}}P^{+},0,\boldsymbol{\ell}_{\perp})$ with
$\phi_{\ell}$ being the respective azimuthal angle. We also have
$x_{\mathcal{P}}\approx M_{\mathcal{H}}^{2}/W^{2}$ (see Eq. (62) below), and
$\ell^{2}=t\approx-\boldsymbol{\ell}_{\perp}^{2}$.

Taking the high energy limit leads to (47). Here we pick up the leading
$F_{1}$ term when $W^{2}\gg|t|$

$\mathcal{M}_{\lambda\bar{\lambda}}(\gamma
p\to\mathcal{H}p)=P_{\mu}\,\mathcal{M}_{\lambda\bar{\lambda}}^{\mu}(\gamma\gamma\to\mathcal{H})\frac{eF_{1}(\ell_{\perp})}{\ell^{2}}\,\frac{1}{q\cdot
P}\bar{u}(P^{\prime})\not{q}u(P)\,,$ (80)

To make the connection (40) we rewrite (80) using QED gauge invariance
$\ell_{\mu}\mathcal{M}_{\lambda\bar{\lambda}}^{\mu}(\gamma\gamma\to\mathcal{H})=0$
in the high-energy limit (the Collins-Ellis trick [89])

$P_{\mu}\mathcal{M}_{\lambda\bar{\lambda}}^{\mu}(\gamma\gamma\to\mathcal{H})\approx-\frac{1}{x_{\mathcal{P}}}\ell_{i}\mathcal{M}_{\lambda\bar{\lambda}}^{i}(\gamma\gamma\to\mathcal{H})\,.$
(81)

so that

$\mathcal{M}_{\lambda\bar{\lambda}}(\gamma
p\to\mathcal{H}p)=-\frac{2eF_{1}(\ell_{\perp})}{t}\frac{W^{2}}{M^{2}_{\mathcal{H}}}\ell_{i}\mathcal{M}_{\lambda\bar{\lambda}}^{i}(\gamma\gamma\to\mathcal{H}).$
(82)

Finally, the result in (79) is isotropic in $\phi_{\ell}$ after the
polarization sum. We can perform the angular average leading to

$\begin{split}\frac{1}{2}\sum_{\bar{\lambda}=-J}^{J}\sum_{\lambda}\int\frac{\mathrm{d}\phi_{\ell}}{2\pi}\left|\ell_{i}\mathcal{M}_{\lambda\bar{\lambda}}^{i}(\gamma\gamma\to\mathcal{H})\right|^{2}&=-\frac{1}{4}\boldsymbol{\ell}_{\perp}^{2}\sum_{\bar{\lambda}=-J}^{J}\sum_{\lambda}g_{ii^{\prime}}\mathcal{M}_{\lambda\bar{\lambda}}^{i}(\gamma\gamma\to\mathcal{H})\mathcal{M}_{\lambda\bar{\lambda}}^{i^{\prime}*}(\gamma\gamma\to\mathcal{H})\\\
&=(2J+1)8\pi
M_{\mathcal{H}}\boldsymbol{\ell}_{\perp}^{2}\Gamma(\mathcal{H}\to\gamma\gamma)\,,\end{split}$
(83)

where in the last line we used (78). Inserting (83) into (79) via (82) gives
Eq. (40), which is the desired result. In the derivation we assumed the
scattering off a proton target, but the formula is valid for any charged
particle.

###  C.2 The $|t|\to 0$ limit for axial vector quarkonia and its connection to
the Landau-Yang theorem

Thanks to the Collins-Ellis trick, the amplitude
$\gamma^{*}\gamma^{*}\to\mathcal{H}$ for scalar and the tensor quarkonia has a
finite $O(\Delta_{\perp})$ contribution. The Collins-Ellis trick is a rather
general and robust procedure in the high-energy limit and so from this
perspective the case of axial quarkonia seems special in that the
$O(\Delta_{\perp})$ contribution vanishes, with the amplitude scaling as
$O(\Delta_{\perp}^{2})$, check (43). In the following we demonstrate that
$O(\Delta_{\perp}^{2})$ scaling follows from a general argument and not
relying to the choice of axial quarkonia wave function. We take Eq. (37) as a
starting point and further factor out the polarization vectors from the
helicity amplitude as

$n_{\nu}\mathcal{M}^{\nu}_{\lambda\bar{\lambda}}(\gamma^{*}\gamma^{*}\to\mathcal{A})=\epsilon_{\mu}(\lambda,q)n_{\nu}E_{\rho}(\bar{\lambda},\Delta_{0})\mathcal{M}^{\mu\nu\rho}(\gamma^{*}\gamma^{*}\to\mathcal{A})\,.$
(84)

The covariant amplitude
$\mathcal{M}^{\mu\nu\rho}(\gamma^{*}\gamma^{*}\to\mathcal{A})$ depends only on
the vectors $q^{\mu}$ and $\ell^{\nu}$. Because $\mathcal{A}$ is an axial
vector quarkonia,
$\mathcal{M}^{\mu\nu\rho}(\gamma^{*}\gamma^{*}\to\mathcal{A})$ has to be
proportional to the $\epsilon$-tensor. In general the covariant decomposition
is as follows

$\begin{split}\mathcal{M}^{\mu\nu\rho}(\gamma^{*}\gamma^{*}\to\mathcal{A})&=\left(q^{\rho}-\ell^{\rho}+\frac{-q^{2}+\ell^{2}}{(q+\ell)^{2}}(q^{\rho}+\ell^{\rho})\right)\epsilon^{\mu\nu
q\ell}M_{\mathcal{A}}F_{TT}(q^{2},\ell^{2})\\\
&+\left(\ell^{\mu}-\frac{q\cdot\ell}{q^{2}}q^{\mu}\right)\epsilon^{\nu\rho
q\ell}\sqrt{-q^{2}}F_{LT}(q^{2},\ell^{2})+\left(q^{\nu}-\frac{q\cdot\ell}{\ell^{2}}\ell^{\nu}\right)\epsilon^{\mu\rho
q\ell}\sqrt{-\ell^{2}}F_{TL}(q^{2},\ell^{2})\,,\end{split}$ (85)

see e. g. [59] and references therein. The decomposition (85) is QED gauge
invariant, that is:
$q_{\mu}\mathcal{M}^{\mu\nu\rho}(\gamma^{*}\gamma^{*}\to\mathcal{A})=0$ and
$\ell_{\nu}\mathcal{M}^{\mu\nu\rho}(\gamma^{*}\gamma^{*}\to\mathcal{A})=0$,
and we also have
$(q^{\rho}+l^{\rho})\mathcal{M}_{\mu\nu\rho}(\gamma^{*}\gamma^{*}\to\mathcal{A})=0$.
Here the notation $F_{TT}$, $F_{TL}$ and $F_{LT}$ refers to the
$\gamma^{*}\gamma^{*}$ polarizations in the center of mass frame. Here only
the form-factor $F_{TT}$ is constrained by the LY theorem: $F_{TT}(0,0)=0$
[75, 76], see also ref. [59]. Due to the Bose symmetry of the full amplitude
(85) we must have $F_{TT}(q^{2},\ell^{2})=-F_{TT}(\ell^{2},q^{2})$ and so
$F_{TT}(q^{2},\ell^{2})\propto q^{2}-\ell^{2}$ for small $q^{2}$ and
$\ell^{2}$. On the other hand, the form-factors $F_{LT}$ and $F_{TL}$ scale as
$F_{LT}(q^{2},\ell^{2})\propto\sqrt{-q^{2}}$,
$F_{TL}(q^{2},\ell^{2})\propto\sqrt{-\ell^{2}}$ in order to avoid kinematic
singularities [90, 59].

Contracting (85) with the longitudinal polarization vector
$E^{\rho}(0,\Delta_{0})$, only the piece proportional to the gauge vector
$n^{\rho}$ from (12) survives. Further contracting with the transverse photon
polarization $\epsilon^{\mu}(\lambda=\pm 1,q)$ and the gauge vector $n^{\nu}$
we obtain

$n_{\nu}\mathcal{M}_{\lambda=\pm
1,\bar{\lambda}=0}^{\nu}(\gamma^{*}\gamma^{*}\to\mathcal{A})=-\frac{M_{0}}{\Delta_{0}\cdot
n}(q\cdot
n)\epsilon^{\epsilon(\lambda,q)nq\ell}\left[\left(1-\frac{q^{2}-\ell^{2}}{\Delta_{0}^{2}}\right)M_{\mathcal{A}}F_{TT}(q^{2},\ell^{2})+\sqrt{-\ell^{2}}F_{TL}(q^{2},\ell^{2})\right]\,,$
(86)

where we have used $\ell\cdot n=0$. The form-factor $F_{LT}(q^{2},\ell^{2})$
decouples after contracting $\epsilon^{\nu\rho ql}$ with
$n_{\nu}E_{\rho}(0,\Delta_{0})$. Taking into account the scalings of the form-
factors the square bracket in (86) is $O(\Delta_{\perp}^{2})$ while the
$\epsilon$-tensor is $O(\Delta_{\perp})$ leading to the overall
$O(\Delta_{\perp}^{3})$ scaling.

For the transverse polarization we find

$\begin{split}n_{\nu}\mathcal{M}_{\lambda=\pm 1,\bar{\lambda}=\pm
1}^{\nu}(\gamma^{*}\gamma^{*}\to\mathcal{A})&=(q-l)\cdot
E(\bar{\lambda},\Delta_{0})\epsilon^{\epsilon(\lambda,q)nq\ell}M_{\mathcal{A}}F_{TT}(q^{2},\ell^{2})+(\epsilon(\lambda,q)\cdot\ell)\epsilon^{nE(\bar{\lambda},\Delta_{0})q\ell}\sqrt{-q^{2}}F_{LT}(q,\ell^{2})\\\
&+(q\cdot
n)\epsilon^{\epsilon(\lambda,q)E(\bar{\lambda},\Delta_{0})q\ell}\sqrt{-\ell^{2}}F_{TL}(q^{2},\ell^{2})\,.\end{split}$
(87)

The first term in (87) is $O(\Delta_{\perp}^{3})$, the second term decouples
in the $q^{2}\to 0$ limit and the third term is $O(\Delta_{\perp}^{2})$ which
is the leading contribution in the $\Delta_{\perp}\to 0$ limit.

###  C.3 The NRQCD limit of the Primakoff cross sections

In this subsection we consider the non-relativisitic QCD (NRQCD) limit of the
Primakoff process by taking the heavy quark limit. We also take $Q^{2}\to 0$
and focus on transverse photon - proton scattering.

In the NRQCD limit all the momenta in the amplitude are considered to be much
smaller than the heavy quark mass. Since the results obtained in the previous
Secs. C.1 and C.2 have already been obtained in the $|t|\to 0$ limit, we take
these as starting point and further expand to second order in $\xi=z-1/2$ and
$\boldsymbol{l}_{\perp}$.

For the scalar $\mathcal{S}$ quarkonia using (42) we obtain

$\widetilde{\mathcal{M}}_{T}\to 8\sqrt{2}\pi\mathrm{i}q_{c}^{2}\alpha
N_{c}\frac{eF_{1}(\Delta_{\perp})}{\Delta_{\perp}}\frac{1}{m_{c}^{3}}\int_{z}\int_{\boldsymbol{l}_{\perp}}(\boldsymbol{l}_{\perp}^{2}+4\xi^{2}m_{c}^{2})\frac{\phi_{\mathcal{S}}(\boldsymbol{l}_{\perp},z)}{z\bar{z}}\,.$
(88)

In the next step we replace the LCWF
$\phi_{\mathcal{S}}(\boldsymbol{l}_{\perp},z)$ by the non-relativistic radial
wave function $u(l)$. Here $\boldsymbol{l}=(\boldsymbol{l}_{\perp},l^{3})$ and
$l$ is its modulus, with $l^{3}=2\xi m_{c}$ as appropriate for the NRQCD
limit. The correspondence is given by [56]

$\int_{z}\int_{\boldsymbol{l}_{\perp}}\frac{\phi_{\mathcal{S}}(\boldsymbol{l}_{\perp},z)}{z\bar{z}}\to\frac{1}{16\pi^{2}\sqrt{m_{c}}}\frac{1}{\sqrt{N_{c}}}\int\mathrm{d}^{3}\boldsymbol{l}\frac{u(l)}{\boldsymbol{l}^{2}}\,.$
(89)

After performing the angular integrations in Eq. (88) we relate $u(l)$ to the
derivative of the radial wave function at the origin, as in Eq. (3.18) of Ref.
[56]:

$\int_{0}^{\infty}\mathrm{d}l\,l^{2}u(l)=3\sqrt{\frac{\pi}{2}}R^{\prime}(0)\,.$
(90)

This leads to

$\widetilde{\mathcal{M}}_{T}\to
12\sqrt{\pi}\mathrm{i}q_{c}^{2}\alpha\sqrt{N_{c}}\frac{eF_{1}(\Delta_{\perp})}{\Delta_{\perp}}\frac{1}{m_{c}^{5/2}}R^{\prime}(0)\,.$
(91)

The cross section is obtained from standard formulas, see Eq. (32). We find

$\frac{\mathrm{d}\sigma(\gamma p\to\mathcal{S}p)}{\mathrm{d}|t|}\to\frac{9\pi
q_{c}^{4}\alpha^{3}N_{c}|R^{\prime}(0)|^{2}F_{1}^{2}(0)}{m_{c}^{7}|t|}\,.$
(92)

The result agrees with Eq. (8b) of Jia et al. [87].

For the axial vector quarkonia there is an additional factor of $\sqrt{3/2}$
on the r.h.s. of Eq. (89). Recalling that the leading contribution in the
$|t|\to 0$ limit comes from the amplitude $\widetilde{\mathcal{M}}_{TT}$ of
(43), and using (33), we find that in the NRQCD limit

$\frac{\mathrm{d}\sigma(\gamma p\to\mathcal{A}p)}{\mathrm{d}|t|}\to\frac{3\pi
q_{c}^{4}\alpha^{3}N_{c}|R^{\prime}(0)|^{2}F_{1}^{2}(0)}{m_{c}^{9}}\,,$ (93)

which also agrees with Eq. (8c) in Ref. [87].

For the tensor quarkonia with $\bar{\lambda}=\pm 1,\pm 2$, the r.h.s. of Eq.
(89) is multiplied by a factor of $\sqrt{3}/M_{\mathcal{T}}$. Using (45) and
(34) we find that the leading contribution comes from the amplitude
$\widetilde{\mathcal{M}}_{TT2,p}$ so that

$\frac{\mathrm{d}\sigma(\gamma p\to\mathcal{T}p)}{\mathrm{d}|t|}\to\frac{12\pi
q_{c}^{4}\alpha^{3}N_{c}|R^{\prime}(0)|^{2}F_{1}^{2}(0)}{m_{c}^{7}|t|}\,.$
(94)

This also agrees with Eq. (8d) in Ref. [87].

##  Appendix D The relative sign of the eikonal photon and Odderon exchange
amplitudes

In this appendix we explain in detail our conventions for covariant
derivatives, field equations, and Wilson lines. Consistent conventions are
important for obtaining the correct relative sign of the Primakoff and Odderon
amplitudes.

We write the covariant derivative in the fundamental representation as

$(D^{\mu})_{ij}=\delta_{ij}\left(\partial^{\mu}+\mathrm{i}eq_{c}A^{\mu}\right)+\mathrm{i}gA^{\mu
a}(t^{a})_{ij}~{}.$ (95)

Here, $i,j=1,\dots,N_{c}$ are fundamental and $a=1,\dots,(N_{c}^{2}-1)$ is an
adjoint color index; $A^{\mu},A^{\mu a}$ represent the electromagnetic and
color fields, respectively, and $t^{a}$ are the traceless generators of color-
SU$(N_{c})$, normalized as $\mathrm{tr}(t^{a}t^{b})=\frac{1}{2}\delta^{ab}$.
$q_{c}$ denotes the fractional electromagnetic charge of the fermion field on
which this covariant derivative acts, $q_{c}=+2/3$ for $c$ quarks.

From the Dirac equation for a massless fermion field, $\psi_{j}$
$\mathrm{i}\gamma_{\mu}(D^{\mu})_{ij}\psi_{j}=0$ in a background with eikonal
(“shockwave”) $A^{\mu}$ and $A^{\mu a}$ fields (where $A^{+}$ is the only non-
vanishing field component and independent of $x^{+}$), $\psi_{j}$ is
proportional to the Wilson line

$V(\boldsymbol{x}_{\perp})={\cal
P}\,\exp\left\\{-\mathrm{i}\int_{-\infty}^{+\infty}\mathrm{d}x^{-}\left[eq_{c}A^{+}(x^{-},\boldsymbol{x}_{\perp})+gA^{+a}(x^{-},\boldsymbol{x}_{\perp})t^{a}\right]\right\\}~{}.$
(96)

The Yang-Mills and Maxwell equations for the shockwave fields in covariant
gauge are

$\displaystyle\partial_{\mu}F^{\mu\nu}_{a}=J^{\nu}_{a}$ $\displaystyle\to$
$\displaystyle-\boldsymbol{\nabla}_{\perp}^{2}A^{+}_{a}=J^{+}_{a}~{},$ (97)
$\displaystyle\partial_{\mu}F^{\mu\nu}=J^{\nu}$ $\displaystyle\to$
$\displaystyle-\boldsymbol{\nabla}_{\perp}^{2}A^{+}=J^{+}~{},$ (98)

where $J^{+a}$ and $J^{+}$ are the plus components of the color and
electromagnetic currents, respectively. Note that here $A^{+}$ denotes the
electromagnetic field sourced by the charges in the proton, from which the
projectile $c\bar{c}$ dipole scatters.

We now define the $S$-matrix for eikonal scattering of the color singlet
dipole, averaged over the colors of the $c,\bar{c}$ quarks

$S(\boldsymbol{x}_{\perp},\boldsymbol{y}_{\perp})=\frac{1}{N_{c}}\,{\rm
tr}\left\langle
V(\boldsymbol{x}_{\perp})V^{\dagger}(\boldsymbol{y}_{\perp})\right\rangle\,.$
(99)

$\langle\cdots\rangle$ denotes the matrix element of the respective operator
between proton states $|P^{+},\boldsymbol{P}_{\perp}\rangle$ and $\langle
P^{+}+\Delta^{+},\boldsymbol{P}_{\perp}+\boldsymbol{\Delta}_{\perp}|$,
stripped of the $\delta$-functions representing conservation of light-cone and
transverse momentum.

We can now define the amplitudes for single photon or three gluon exchange as
follows. Setting the QCD coupling $g=0$ and expanding the Wilson lines to
linear order in $A^{+}$ we have

$V(\boldsymbol{x}_{\perp})=1-\mathrm{i}eq_{c}\int_{-\infty}^{+\infty}\mathrm{d}x^{-}A^{+}(x^{-},\boldsymbol{x}_{\perp})\equiv
1-\mathrm{i}eq_{c}\alpha(\boldsymbol{x}_{\perp})\,.$ (100)

Then,

$\Omega(\boldsymbol{x}_{\perp},\boldsymbol{y}_{\perp})\equiv-\frac{1}{2\mathrm{i}}\left\langle
V(\boldsymbol{x}_{\perp})V^{\dagger}(\boldsymbol{y}_{\perp})-V(\boldsymbol{y}_{\perp})V^{\dagger}(\boldsymbol{x}_{\perp})\right\rangle=eq_{c}\left\langle\alpha(\boldsymbol{x}_{\perp})-\alpha(\boldsymbol{y}_{\perp})\right\rangle\,.$
(101)

corresponds to the scattering amplitude for single photon exchange. With the
field equation (98) we can write
$\langle\alpha(\boldsymbol{x}_{\perp})\rangle=-\left\langle\frac{1}{\boldsymbol{\nabla}^{2}_{\perp}}\rho(\boldsymbol{x}_{\perp})\right\rangle$
where the integrated electric charge density is
$\rho(\boldsymbol{x}_{\perp})=\int_{-\infty}^{+\infty}\mathrm{d}x^{-}J^{+}(x^{-},\boldsymbol{x}_{\perp})$.
Performing a Fourier transform from the transverse coordinate to the
transverse momentum space, the matrix element of this operator is simply the
Dirac electromagnetic form factor of the proton,
$eF_{1}(q_{\perp}^{2})=\langle\rho(\boldsymbol{q}_{\perp})\rangle$. This leads
to

$\Omega(\boldsymbol{r}_{\perp},\boldsymbol{\Delta}_{\perp})=8\pi\mathrm{i}\,\alpha
q_{c}\frac{F_{1}(\Delta_{\perp}^{2})}{\boldsymbol{\Delta}_{\perp}^{2}}\,\sin\left(\frac{\boldsymbol{r}_{\perp}\cdot\boldsymbol{\Delta}_{\perp}}{2}\right)~{}.$
(102)

Here we introduced
$\boldsymbol{r}_{\perp}=\boldsymbol{x}_{\perp}-\boldsymbol{y}_{\perp}$. This
is the expression written in Eq. (36) of the main text. In impact parameter
space,

$\Omega(\boldsymbol{r}_{\perp},\boldsymbol{b}_{\perp})=8\pi\alpha
q_{c}\int\limits_{\boldsymbol{\Delta}_{\perp}}\sin(\boldsymbol{b}_{\perp}\cdot\boldsymbol{\Delta}_{\perp})\,\frac{-F_{1}(\Delta_{\perp}^{2})}{\boldsymbol{\Delta}_{\perp}^{2}}\,\sin\left(\frac{\boldsymbol{r}_{\perp}\cdot\boldsymbol{\Delta}_{\perp}}{2}\right)~{}.$
(103)

Figure 8: The Odderon amplitude
$\mathcal{O}(\boldsymbol{r}_{\perp},\boldsymbol{b}_{\perp})$ obtained from a
light-front quark model [65, 66], when the relative angle between
$\boldsymbol{r}_{\perp}$ and $\boldsymbol{b}_{\perp}$ $\phi_{rb}=0$. At small
$b_{\perp}$ the sign is negative.

We now proceed to the three gluon exchange amplitude by setting $q_{c}=0$ in
the Wilson line (96) followed by an expansion of (5), namely

$\mathcal{O}(\boldsymbol{x}_{\perp},\boldsymbol{y}_{\perp})=-\frac{1}{2\mathrm{i}N_{c}}\,\mathrm{tr}\,\left\langle
V(\boldsymbol{x}_{\perp})V^{\dagger}(\boldsymbol{y}_{\perp})-V(\boldsymbol{y}_{\perp})V^{\dagger}(\boldsymbol{x}_{\perp})\right\rangle\,,$
(104)

to third order in $A^{+a}$:

$\begin{split}\mathcal{O}(\boldsymbol{r}_{\perp},\boldsymbol{b}_{\perp})&=\frac{g^{3}}{4N_{c}}d^{abc}\int_{\boldsymbol{q}_{1\perp}\boldsymbol{q}_{2\perp}\boldsymbol{q}_{3\perp}}\frac{\langle\rho^{a}(\boldsymbol{q}_{1\perp})\rho^{b}(\boldsymbol{q}_{2\perp})\rho^{c}(\boldsymbol{q}_{3\perp})\rangle}{\boldsymbol{q}_{1\perp}^{2}\boldsymbol{q}_{2\perp}^{2}\boldsymbol{q}_{3\perp}^{2}}\sin\left(\boldsymbol{b}_{\perp}\cdot(\boldsymbol{q}_{1\perp}+\boldsymbol{q}_{2\perp}+\boldsymbol{q}_{3\perp})\right)\\\
&\times\left[\sin\left(\frac{1}{2}\boldsymbol{r}_{\perp}\cdot(\boldsymbol{q}_{1\perp}-\boldsymbol{q}_{2\perp}-\boldsymbol{q}_{3\perp})\right)+\frac{1}{3}\sin\left(\frac{1}{2}\boldsymbol{r}_{\perp}\cdot(\boldsymbol{q}_{1\perp}+\boldsymbol{q}_{2\perp}+\boldsymbol{q}_{3\perp})\right)\right]~{}.\end{split}$
(105)

For details see Eq. (77) of Ref. [22], to be multiplied by a factor of
$\mathrm{i}$ (compare our Eq. (104) to their Eq. (76)), and their Appendix B.
Here
$\rho^{a}(\boldsymbol{x}_{\perp})=\int_{-\infty}^{+\infty}\mathrm{d}x^{-}J^{+a}(x^{-},\boldsymbol{x}_{\perp})$
denotes the color charge density integrated along the eikonal path. We now
write the $C$-conjugation odd part of the color charge correlator as

$\left\langle\rho^{a}(\boldsymbol{q}_{1\perp})\rho^{b}(\boldsymbol{q}_{2\perp})\rho^{c}(\boldsymbol{q}_{3\perp})\right\rangle=\frac{1}{4}d^{abc}g^{3}\,G_{3}^{-}(\boldsymbol{q}_{1\perp},\boldsymbol{q}_{2\perp},\boldsymbol{q}_{3\perp})~{}.$
(106)

In a three-quark model this correlator takes the form [22]

$\displaystyle
G_{3}^{-}(\boldsymbol{q}_{1\perp},\boldsymbol{q}_{2\perp},\boldsymbol{q}_{3\perp})$
$\displaystyle=$
$\displaystyle\int[\mathrm{d}x_{i}]\int[\mathrm{d}^{2}p_{i}]\,\psi(x_{1},\boldsymbol{p}_{1\perp};x_{2},\boldsymbol{p}_{2\perp};x_{3},\boldsymbol{p}_{3\perp})$
(107)
$\displaystyle\Bigl{[}\psi^{*}(x_{1},\boldsymbol{p}_{1\perp}+(x_{1}-1)\boldsymbol{q}_{\perp};x_{2},\boldsymbol{p}_{2\perp}+x_{2}\boldsymbol{q}_{\perp};x_{3},\boldsymbol{p}_{3\perp}+x_{3}\boldsymbol{q}_{\perp})$
$\displaystyle-\psi^{*}(x_{1},\boldsymbol{p}_{1\perp}-\boldsymbol{q}_{1\perp}+x_{1}\boldsymbol{q}_{\perp};x_{2},\boldsymbol{p}_{2\perp}-\boldsymbol{q}_{2\perp}-\boldsymbol{q}_{3\perp}+x_{2}\boldsymbol{q}_{\perp};x_{3},\boldsymbol{p}_{3\perp}+x_{3}\boldsymbol{q}_{\perp})$