# Chiral and trace anomalies in Deeply Virtual Compton Scattering II:
QCD factorization and beyond

Shohini Bhattacharya<EMAIL_ADDRESS>RIKEN BNL Research Center, Brookhaven
National Laboratory, Upton, NY 11973, USA  Yoshitaka Hatta<EMAIL_ADDRESS>Physics Department, Brookhaven National Laboratory, Upton, NY 11973, USA
RIKEN BNL Research Center, Brookhaven National Laboratory, Upton, NY 11973,
USA  Werner Vogelsang<EMAIL_ADDRESS>Institute for
Theoretical Physics, Tübingen University, Auf der Morgenstelle 14, 72076
Tübingen, Germany

###### Abstract

We extend the discussion of the recently discovered ‘anomaly poles’ in QCD
Compton scattering. We perform the complete one-loop calculation of the
Compton amplitude using momentum transfer $t$ as the regulator of collinear
divergences. In the gluon channel, we confirm the presence of poles $1/t$ in
both the real and imaginary parts of the amplitude. In the quark channel, we
find unexpected infrared single $1/\epsilon$ and double $1/\epsilon^{2}$
poles. We then perform the one-loop calculation of the leading-twist quark
generalized parton distributions (GPDs) for quark and gluon external states
with the same regulators and find that all these singular terms can be
systematically absorbed into the GPDs, showing that QCD factorization is
restored to this order. Having established this, we discuss the fate of the
$1/t$ poles. We argue that they become the nonperturbative building blocks of
GPDs that encode the chiral and trace anomalies of QCD, in a way consistent
with the known constraints these anomalies impose on the nucleon axial and
gravitational form factors. The scope of research on GPDs can therefore be
expanded to address the manifestation and implications of quantum anomalies in
high-energy exclusive processes.

##  I Introduction

The past several years have witnessed significant progress in the higher-order
calculation of Deeply Virtual Compton Scattering (DVCS). In the flavor-
nonsinglet channel, the three-loop evolution equation for the generalized
parton distributions (GPDs) has been derived Braun et al. (2017) together with
the two-loop coefficient functions Braun et al. (2020). In the flavor-singlet
channel, the two-loop coefficient functions for DVCS have been recently
calculated Braun et al. (2022) and even higher order resummation effects have
been studied Schoenleber (2023). These developments are on a steady path
toward achieving the NNLO accuracy in DVCS that is required for precision GPD
studies at the future Electron-Ion Collider (EIC) Abdul Khalek et al. (2022).

Meanwhile, in a previous paper Bhattacharya et al. (2023), we have explored a
new approach to compute the NLO corrections in DVCS, following an earlier
suggestion in polarized deep inelastic scattering (DIS) Tarasov and
Venugopalan (2020, 2022). The key ingredient is to use momentum transfer
$t=(P_{1}-P_{2})^{2}$ as an infrared cutoff to regulate the collinear
divergence, instead of the usual dimensional regularization. Previously in the
calculation of the Compton amplitude in the Bjorken limit, the variable $t$
had always been neglected when computing partonic amplitudes. Naively, one
would expect that the only new effect of introducing nonzero $t$ would be to
generate higher twist corrections of order $|t|/Q^{2}\ll 1$ where $Q^{2}$ is
the photon virtuality. However, our explicit calculations with nonzero $t$
have revealed ‘anomaly poles’ $1/t$ which had not been detected in the
previous calculations performed at $t=0$ Ji and Osborne (1998); Belitsky and
Mueller (1998); Mankiewicz et al. (1998); Pire et al. (2011); Bertone et al.
(2022), but are consistent with the result in Tarasov and Venugopalan (2020).
Moreover, these poles are accompanied by certain twist-four GPDs but without
an expected suppression factor $1/Q^{2}$. (Rather, $1/Q^{2}$ has been replaced
by $1/t$.) In fact, they can be interpreted as the manifestations of the QCD
chiral Jaffe and Manohar (1989); Tarasov and Venugopalan (2020, 2022) and
trace Bhattacharya et al. (2023) anomalies in high energy scattering. The
finding thus points towards a novel connection between the study of GPDs and
phenomena associated with quantum anomalies such as chiral symmetry breaking
and confinement.

At face value, the emergence of poles is in apparent contradiction with the
QCD factorization theorem Collins and Freund (1999); Ji and Osborne (1998)
which states that the QCD Compton scattering amplitude factorizes into the
perturbatively calculable coefficient functions and the nonperturbative twist-
two GPDs up to higher-twist corrections of order $1/Q^{2}$. However, we have
already suggested in Bhattacharya et al. (2023) that the poles found in the
one-loop calculation may be absorbed into the twist-two GPDs as a part of the
infrared subtraction procedure. The purpose of this paper is to fully
demonstrate that this is indeed the case. We first perform a complete
calculation of the Compton amplitude with nonzero $t$ at one loop, both in the
quark and gluon channels, in both the polarized and unpolarized sectors, and
for the real and imaginary parts of the amplitude. (In Bhattacharya et al.
(2023), we only calculated the imaginary part in the gluon channel.) In the
gluon channel, we find $1/t$ poles also in the real part. Surprisingly, in the
quark channel, we find uncancelled infrared single $1/\epsilon$ and double
$1/\epsilon^{2}$ poles. We next perform the corresponding one-loop calculation
of the unpolarized and polarized quark GPDs for free quark and gluon external
states at finite $t$ and show that all the singular terms can be
systematically absorbed.

Therefore, at least to one loop, the emergence of $1/t$ poles does not
contradict the QCD factorization theorem. The calculation with nonzero $t$ may
be regarded as an alternative factorization scheme. Having established this,
we shift our focus to the fate of the $1/t$ poles absorbed into the twist-two
GPDs. It is well known that the chiral and trace anomalies impose constraints
on the nucleon axial and gravitational form factors, respectively. Since these
form factors are certain moments of the twist-two GPDs, there must be
corresponding constraints directly on GPDs Bhattacharya et al. (2023). A
preliminary discussion of this has been already presented in Bhattacharya et
al. (2023). Our extended treatment here will lend more support to the idea
that this new scheme can uniquely address such profound aspects of QCD in GPD
studies.

##  II Compton scattering

The amplitude for QCD Compton scattering off a proton target,
$\gamma^{*}(q_{1})p(P_{1})\to\gamma^{*}(q_{2})p(P_{2})$, is given by

$\displaystyle T^{\mu\nu}$ $\displaystyle=$ $\displaystyle i\int
d^{4}ye^{iq\cdot y}\langle P_{2}|{\rm
T}\\{j^{\mu}(y/2)j^{\nu}(-y/2)\\}|P_{1}\rangle,$ (1)

where $j^{\mu}=\sum_{q}e_{q}\bar{\psi}_{q}\gamma^{\mu}\psi_{q}$ is the
electromagnetic current and $q^{\mu}=\frac{q^{\mu}_{1}+q^{\mu}_{2}}{2}$ is the
average of the incoming and outgoing photon momenta. The momentum transfer is
denoted as $t=l^{2}$ where
$l^{\mu}=P_{2}^{\mu}-P_{1}^{\mu}=q_{1}^{\mu}-q_{2}^{\mu}$. We introduce the
generalized Bjorken variable $x_{B}$ and the skewness parameter $\xi$,

$\displaystyle x_{B}=\frac{Q^{2}}{2P\cdot
q},\quad\xi=\frac{q_{2}^{2}-q_{1}^{2}}{2P\cdot
q}\approx-\frac{l^{+}}{2P^{+}},$ (2)

where $Q^{2}=-q^{2}$ is the photon virtuality and
$P^{\mu}=\frac{P_{1}^{\mu}+P_{2}^{\mu}}{2}$. In DVCS, $q_{2}^{2}=0$ and
$x_{B}\approx\xi$, but we shall keep general $x_{B}$ and $\xi$ throughout the
paper.

In the generalized Bjorken limit $Q^{2},\,2P\cdot q\to\infty$ with $x_{B}$,
$t$ fixed and $Q^{2}\gg|t|$, the Compton amplitude can be expanded as Diehl
(2003); Belitsky and Radyushkin (2005)

$\displaystyle
T^{\mu\nu}=\frac{g^{\mu\nu}_{\perp}}{2P^{+}}\bar{u}(P_{2})\left[\gamma^{+}{\cal
H}+\frac{i\sigma^{+\nu}l_{\nu}}{2M}{\cal
E}\right]u(P_{2})-i\frac{\epsilon^{\mu\nu}_{\perp}}{2P^{+}}\bar{u}(P_{2})\left[\gamma^{+}\gamma_{5}\tilde{{\cal
H}}+\frac{\gamma^{5}l^{+}}{2M}\tilde{{\cal E}}\right]u(P_{2})+\cdots\,,$ (3)

where $M$ is the proton mass. $g_{\perp}^{\mu\nu}$ and
$\epsilon_{\perp}^{\mu\nu}$ are transverse projectors such that
$g^{ij}_{\perp}=-\delta^{ij}$ and $\epsilon^{ij}_{\perp}=\epsilon^{ij}$ for
transverse indices $i,j=1,2$ and the other components are zero. Our convention
is $\gamma_{5}=i\gamma^{0}\gamma^{1}\gamma^{2}\gamma^{3}$ and
$\epsilon^{0123}=\epsilon^{-+12}=\epsilon^{12}=1$. The ellipsis in (3) stand
for the contributions from the (generalized) longitudinal structure function
and the so-called gluon transversity GPD. As observed in Bhattacharya et al.
(2023), they are not sensitive to anomalies and are thus left for future work.

According to QCD factorization, the Compton form factors ${\cal H}$, ${\cal
E}$, $\tilde{\cal H}$ and $\tilde{\cal E}$ can be written as convolutions of
nonperturbative GPDs and partonic hard-scattering amplitudes. The latter are
commonly calculated in dimensional regularization in $d=4-2\epsilon$
dimensions, with $\epsilon$ regularizing both ultraviolet (UV) and infrared
(IR) divergences. We shall also work in $d$ dimensions, since individual
diagrams contain UV divergences. But we regularize the collinear singularity
by introducing the physical variable $t$ at the partonic level. Such a
calculation is safely justified when $\sqrt{|t|}\gg\Lambda_{\rm QCD}$ (still
keeping $Q\gg\sqrt{|t|}$),111 We consider large but finite $\sqrt{|t|}$, so we
are still in the generalized Bjorken limit $|t|\ll Q^{2}\to\infty$ and the
usual factorization in terms of the twist-two GPDs is expected. This is
different from so-called wide angle Compton scattering (see, e.g., Radyushkin
(1998)) where $Q=0$ and $s\sim-t\sim-u$. In our calculation, $s\sim Q^{2}$ and
we systematically neglect terms of order $t/Q^{2}\sim t/s$.  but we shall
eventually be interested in the region $\sqrt{|t|}\sim\Lambda_{\rm QCD}$. The
result of the one-loop calculation can be summarized in the form

$\displaystyle\begin{pmatrix}{\cal H}(x_{B},\xi,t)\\\ {\cal
E}(x_{B},\xi,t)\end{pmatrix}$ $\displaystyle=$
$\displaystyle\sum_{q}e_{q}^{2}\int_{0}^{1}dx\Biggl{[}\left(C_{0}(x,x_{B})+\frac{\alpha_{s}}{2\pi}C^{q}_{1}(x,x_{B},\xi)\right)\begin{pmatrix}H_{q}(x,\xi,t)-H_{q}(-x,\xi,t)\\\
E_{q}(x,\xi,t)-E_{q}(-x,\xi,t)\end{pmatrix}$ (4)
$\displaystyle\qquad\qquad\qquad+\frac{\alpha_{s}}{2\pi}C^{g}_{1}(x,x_{B},\xi)\begin{pmatrix}H_{g}(x,\xi,t)\\\
E_{g}(x,\xi,t)\end{pmatrix}$
$\displaystyle\qquad\qquad\qquad+\frac{\alpha_{s}}{2\pi}\frac{M^{2}}{t}A(x,x_{B},\xi)\begin{pmatrix}{\cal
F}(x,\xi,t)\\\ -{\cal F}(x,\xi,t)\end{pmatrix}\Biggr{]}+{\cal
O}(1/Q^{2})+{\cal O}(\alpha_{s}^{2}),$
$\displaystyle\begin{pmatrix}\tilde{{\cal H}}(x_{B},\xi,t)\\\ \tilde{{\cal
E}}(x_{B},\xi,t)\end{pmatrix}$ $\displaystyle=$
$\displaystyle\sum_{q}e_{q}^{2}\int_{0}^{1}dx\Biggl{[}\left(\tilde{C}_{0}(x,x_{B})+\frac{\alpha_{s}}{2\pi}\tilde{C}^{q}_{1}(x,x_{B},\xi)\right)\begin{pmatrix}\tilde{H}_{q}(x,\xi,t)+\tilde{H}_{q}(-x,\xi,t)\\\
\tilde{E}_{q}(x,\xi,t)+\tilde{E}_{q}(-x,\xi,t)\end{pmatrix}$ (5)
$\displaystyle\qquad\qquad\qquad+\frac{\alpha_{s}}{2\pi}\tilde{C}_{1}^{g}(x,x_{B},\xi)\begin{pmatrix}\tilde{H}_{g}(x,\xi,t)\\\
\tilde{E}_{g}(x,\xi,t)\end{pmatrix}$
$\displaystyle\qquad\qquad\qquad+\frac{\alpha_{s}}{2\pi}\frac{M^{2}}{t}\tilde{A}(x,x_{B},\xi)\begin{pmatrix}0\\\
\tilde{\cal F}(x,\xi,t)\end{pmatrix}\Biggr{]}+{\cal O}(1/Q^{2})+{\cal
O}(\alpha_{s}^{2}),$

where the notations for the twist-two quark and gluon GPDs $H_{q,g}$,
$E_{q,g}$, $\tilde{H}_{q,g}$, $\tilde{E}_{q,g}$ are standard Diehl (2003);
Belitsky and Radyushkin (2005). (The gluon GPDs are normalized as
$H_{g}(x)=xG(x)$ and $\tilde{H}_{g}(x)=x\Delta G(x)$ in the forward limit,
where $G$, $\Delta G$ are the unpolarized and polarized gluon PDFs.) Note that
(4) and (5) are still in their ‘unsubtracted’ forms in the sense that some of
the coefficients contain divergences in the formal limit $t\to 0$. Their
subtraction is rather unconventional, and to elaborate on this is one of the
main objectives of our paper.

The leading-order coefficient functions are well known:

$\displaystyle C_{0}\left(x,x_{B}\right)$
$\displaystyle=\frac{1}{x-x_{B}+i\epsilon}+\frac{1}{x+x_{B}-i\epsilon},$
$\displaystyle\tilde{C}_{0}\left(x,x_{B}\right)$
$\displaystyle=\frac{1}{x-x_{B}+i\epsilon}-\frac{1}{x+x_{B}-i\epsilon}.$ (6)

The one-loop corrections, $C_{1}^{q}$ etc., have the following generic
structure:

$\displaystyle C^{q}_{1}(x,x_{B},\xi)$
$\displaystyle=\frac{C_{F}}{x}\left(\kappa_{qq}(\hat{x},\hat{\xi})\ln\frac{Q^{2}}{-l^{2}}+\delta
C^{q}_{1}(\hat{x},\hat{\xi})\right),\qquad\tilde{C}^{q}_{1}(x,x_{B},\xi)=\frac{C_{F}}{x}\left(\tilde{\kappa}_{qq}(\hat{x},\hat{\xi})\ln\frac{Q^{2}}{-l^{2}}+\delta\tilde{C}^{q}_{1}(\hat{x},\hat{\xi})\right),$
(7) $\displaystyle C^{g}_{1}(x,x_{B},\xi)$
$\displaystyle=\frac{2T_{R}}{x^{2}}\left(\kappa_{qg}(\hat{x},\hat{\xi})\ln\frac{Q^{2}}{-l^{2}}+\delta
C^{g}_{1}(\hat{x},\hat{\xi})\right),\qquad\tilde{C}^{g}_{1}(x,x_{B},\xi)=\frac{2T_{R}}{x^{2}}\left(\tilde{\kappa}_{qg}(\hat{x},\hat{\xi})\ln\frac{Q^{2}}{-l^{2}}+\delta\tilde{C}^{g}_{1}(\hat{x},\hat{\xi})\right),$
(8)

where $C_{F}=\frac{4}{3}$ and $T_{R}=\frac{1}{2}$ are the usual color factors.
We have introduced the partonic variables $\hat{x}=\frac{x_{B}}{x}$ and
$\hat{\xi}=\frac{\xi}{x}$, and set the $\overline{\rm MS}$ renormalization
scale to be $4\pi e^{-\gamma_{E}}\mu^{2}=Q^{2}$. The logarithm
$\ln\frac{Q^{2}}{-l^{2}}$ originates from the collinear singularity and
replaces the $-\tfrac{1}{\epsilon_{\text{\tiny{IR}}}}$ pole in the usual
calculation in dimensional regularization with $t=0$. The coefficients
$\kappa$, $\tilde{\kappa}$ are fixed by the evolution equation of GPDs and
must agree with the known results in the literature. On the other hand, the
coefficient functions $\delta C_{1}^{q}$, $\delta C_{1}^{g}$,
$\delta\tilde{C}_{1}^{q}$ and $\delta\tilde{C}_{1}^{g}$ are potentially scheme
dependent. The results in the $\overline{\rm MS}$ scheme can be found in Ji
and Osborne (1998); Belitsky and Mueller (1998); Mankiewicz et al. (1998);
Belitsky and Radyushkin (2005). Note that, somewhat at variance with the
previous literature, we have used the reflection symmetry in $x$ to restrict
the $x$-integral to the region $0<x<1$. Namely, $\tilde{C}_{0}$,
$\tilde{C}_{1}^{q}$, $C_{1}^{g}$, $H_{g}$, $E_{g}$ are even functions and
$C_{0}$, $C_{1}^{q}$, $\tilde{C}_{1}^{g}$ $\tilde{H}_{g}$, $\tilde{E}_{g}$ are
odd functions, respectively, under $x\to-x$. This is convenient for the
discussion below.

Eqs. (4) and (5) resemble the usual structure dictated by the QCD
factorization theorem except for the ‘anomaly pole’ terms proportional to
$1/t$. These poles are accompanied by the twist-four gluon GPDs ${\cal F}$ and
$\tilde{\cal F}$ defined as Hatta and Zhao (2020); Tarasov and Venugopalan
(2020); Hatta (2020); Radyushkin and Zhao (2021, 2022)

$\displaystyle{\cal F}(x,\xi,t)$
$\displaystyle\equiv\frac{-4xP^{+}}{M}\int\frac{dz^{-}}{2\pi}e^{ixP^{+}z^{-}}\frac{\langle
P_{2}|F^{\mu\nu}(-z^{-}/2)WF_{\mu\nu}(z^{-}/2)|P_{1}\rangle}{\bar{u}(P_{2})u(P_{1})}\,,$
(9) $\displaystyle\tilde{{\cal F}}(x,\xi,t)$
$\displaystyle\equiv\frac{iP^{+}}{M}\int\frac{dz^{-}}{2\pi}e^{ixP^{+}z^{-}}\frac{\langle
P_{2}|F^{\mu\nu}(-z^{-}/2)W\tilde{F}_{\mu\nu}(z^{-}/2)|P_{1}\rangle}{\bar{u}(P_{2})\gamma_{5}u(P_{1})}\,,$
(10)

where $W$ is the straight Wilson line between $[-z^{-}/2,z^{-}/2]$. We have
changed the normalization with respect to Bhattacharya et al. (2023) in order
to make these distributions dimensionless. Despite involving twist-four GPDs,
these terms are not suppressed by $1/Q^{2}$ and apparently cause problems in
the forward limit $t\to 0$. As discussed in Bhattacharya et al. (2023) and
will be further elaborated later, the emergence of poles and their fate may
shed new light on the nonperturbative structure of GPDs, connecting to deep
issues such as chiral symmetry breaking and the origin of hadron masses.

##  III Calculations

Figure 1: Diagrams for the subprocess $\gamma^{*}g\to\gamma^{*}g$ in Compton
scattering. Figure 2: Diagrams for the subprocess $\gamma^{*}q\to\gamma^{*}q$
in Compton scattering. Diagrams with photon lines crossed are not shown.

In this section, we outline our calculation of the perturbative corrections to
Compton scattering at one-loop order. The relevant Feynman diagrams for the
subprocess initiated by the gluons are shown in Fig. 1 and the ones initiated
by the quarks are shown in Fig. 2. For the latter case, we choose to work in
Feynman gauge. (We have also worked in a general covariant gauge and confirmed
that the final results are independent of the gauge as it should be.) As in
Ref. Bhattacharya et al. (2023), we parameterize the incoming and outgoing
momenta as

$\displaystyle q_{1}=q+\frac{l}{2},\quad q_{2}=q-\frac{l}{2},\quad
p_{1}=p-\frac{l}{2},\quad p_{2}=p+\frac{l}{2}.$ (11)

We also define the partonic versions of the Bjorken and skewness variables (2)
as

$\displaystyle\hat{x}=\frac{Q^{2}}{2p\cdot
q}=\frac{x_{B}}{x},\qquad\hat{\xi}=\frac{\xi}{x}=-\hat{x}\frac{q\cdot
l}{Q^{2}}.$ (12)

The incoming and outgoing partons are assumed to be massless,
$p_{1}^{2}=p_{2}^{2}=0$, which leads to the conditions $p^{2}=-l^{2}/4$ and
$p\cdot l=0$. The virtuality of the photons can be written as

$\displaystyle
q_{1}^{2}=-Q^{2}\frac{\hat{x}+\hat{\xi}}{\hat{x}}+\frac{l^{2}}{4},\qquad
q_{2}^{2}=-Q^{2}\frac{\hat{x}-\hat{\xi}}{\hat{x}}+\frac{l^{2}}{4}.$ (13)

We will abbreviate the polarization vectors for the gluons in Fig. 1 as
$\epsilon^{\mu}(p_{1})\equiv\epsilon_{1}^{\mu}$ and
$\epsilon^{*\mu}(p_{2})\equiv\epsilon_{2}^{*\mu}$.

The collinear singularity in the above diagrams will be regularized by
$t=l^{2}$. We emphasize that, in the present ‘handbag’ approximation,
$t=(p_{2}-p_{1})^{2}=(P_{2}-P_{1})^{2}$ is the same at the hadronic and
partonic levels. However, we still have to work in $d=4-2\epsilon$ dimensions
because the individual diagrams will contain UV divergences in the real part.
At the same time, working in $d$-dimensions also helps to check if there are
any leftover IR divergences that are not regularized by nonzero $t$ alone.
This point will be relevant for the quark-channel diagrams in Fig. 2. Our
convention is that, if $\epsilon$ is used for the UV divergences, then
$\epsilon\rightarrow\epsilon_{\text{\tiny{UV}}}>0$, while if it is used for
the IR divergences then $\epsilon\rightarrow\epsilon_{\text{\tiny{IR}}}<0$.

In the following, we shall refer to the two terms in (3) as the symmetric and
antisymmetric parts of the Compton amplitude, respectively. The symmetric part
can be extracted with the help of the projector

$\displaystyle g^{\mu\nu}_{\perp}$
$\displaystyle=g^{\mu\nu}+\frac{1}{q^{2}(1+\gamma^{2})}\left(q^{\mu}-\frac{q^{2}}{p\cdot
q}p^{\mu}\right)\left(q^{\nu}-\frac{q^{2}}{p\cdot
q}p^{\nu}\right)-\frac{q^{\mu}q^{\nu}}{q^{2}}\,,\qquad\gamma^{2}=-\frac{p^{2}q^{2}}{(p\cdot
q)^{2}}=\frac{l^{2}q^{2}}{4(p\cdot q)^{2}},$ (14)

such that

$\displaystyle g^{\mu}_{\perp\mu}=d-2=2(1-\epsilon),\qquad{\cal H},{\cal
E}\sim\frac{1}{2(1-\epsilon)}g_{\mu\nu}^{\perp}T^{\mu\nu}\,.$ (15)

For the antisymmetric part, we use the projector $\epsilon^{\alpha
p\mu\nu}\equiv\epsilon^{\alpha\beta\mu\nu}p_{\beta}$.

We evaluate the above diagrams with the help of the Mathematica package
‘Package-X’ Patel (2015). Below we first discuss the main features of our
results specific to the gluon and quark channels. The complete results will
then be presented in Sec. IV.

###  III.1 Gluon channel

Our results feature (i) a ${1}/{l^{2}}$ pole and (ii) a logarithm
$\ln({Q^{2}/}{-l^{2}})$, both arising from the first and third diagrams in
Fig. 1. For the symmetric case, the UV poles from the first and third diagrams
add up to cancel the one arising from the second diagram. For the
antisymmetric case, the UV poles from the first and third diagram cancel.
There are no leftover ${1}/{\epsilon_{\text{\tiny{IR}}}}$ divergences, meaning
that $l^{2}\neq 0$ functions as a genuine regulator of collinear divergences.

In the symmetric case, the result for the one-loop Compton scattering
amplitude with external gluon polarization vectors $\epsilon_{1}$,
$\epsilon_{2}^{*}$ (Fig. 1) has the following generic structure:

$\displaystyle-\epsilon_{1}\cdot\epsilon^{*}_{2}\left(A\ln\frac{Q^{2}}{-l^{2}}+B\right)+\frac{C}{l^{2}}\epsilon_{1}\cdot
l\epsilon_{2}^{*}\cdot
l=-\epsilon_{1}\cdot\epsilon^{*}_{2}\left(A\ln\frac{Q^{2}}{-l^{2}}+B-\frac{C}{2}\right)+\frac{C}{l^{2}}\left(\epsilon_{1}\cdot
l\epsilon_{2}^{*}\cdot
l-\frac{\epsilon_{1}\cdot\epsilon_{2}^{*}}{2}l^{2}\right),$ (16)

where $A,B,C$ are coefficients that depend on $\hat{x}$ and $\hat{\xi}$. In
the asymmetric case, we find instead

$\displaystyle i\epsilon^{\alpha
p\epsilon^{*}_{2}\epsilon_{1}}\left(\tilde{A}\ln\frac{Q^{2}}{-l^{2}}+\tilde{B}\right)+\tilde{C}\frac{l^{\alpha}}{l^{2}}i\epsilon^{\epsilon_{1}\epsilon_{2}^{*}lp},$
(17)

where
$\epsilon^{\epsilon_{1}\epsilon_{2}^{*}lp}\equiv\epsilon^{\mu\nu\rho\lambda}\epsilon_{1\mu}\epsilon^{*}_{2\nu}l_{\rho}p_{\lambda}$.
The first terms in (16) and (17) can be interpreted as the usual one-loop
corrections to the Compton amplitude $\sim
C_{1}^{g}H_{g},\tilde{C}_{1}^{g}\tilde{H}_{g}$ through the identifications

$\displaystyle-\epsilon_{1}\cdot\epsilon_{2}^{*}\sim\frac{\langle
p_{2}|F^{+\mu}F^{\ +}_{\mu}|p_{1}\rangle}{1-\hat{\xi}^{2}},\qquad
i\epsilon^{+p\epsilon^{*}_{2}\epsilon_{1}}\sim\frac{\langle
p_{2}|iF^{+\mu}\tilde{F}^{\ +}_{\mu}|p_{1}\rangle}{1-\hat{\xi}^{2}}.$ (18)

However, the second terms in (16) and (17) cannot be attributed to twist-two
GPDs. Their structures can only arise from the twist-four operators
$F^{\mu\nu}F_{\mu\nu}$ and $F^{\mu\nu}\tilde{F}_{\mu\nu}$

$\displaystyle\epsilon_{1}\cdot l\epsilon_{2}^{*}\cdot
l-\frac{\epsilon_{1}\cdot\epsilon^{*}_{2}}{2}l^{2}\sim\langle
p_{2}|F^{\mu\nu}F_{\mu\nu}|p_{1}\rangle,\qquad
2i\epsilon^{\epsilon_{1}\epsilon_{2}^{*}lp}\sim\langle
p_{2}|iF^{\mu\nu}\tilde{F}_{\mu\nu}|p_{1}\rangle,$ (19)

and this is how the twist-four GPDs (9), (10) come into play. It should be
mentioned, however, that the present argument only concerns the two-parton
matrix element of the operators $FF$ and $F\tilde{F}$. Further justifications
from other approaches are desirable.

###  III.2 Quark channel

In this case, our results do not contain any terms ${1}/{l^{2}}$. This is
consistent with the expectation that the anomalies, being of purely gluonic
nature, should not affect the quark sector, at least at this order. Quite
unexpectedly though, we find (i) double IR poles
${1}/{\epsilon^{2}_{\text{\tiny{IR}}}}$, (ii) single IR poles
${1}/{\epsilon_{\text{\tiny{IR}}}}$, apart from (iii) a logarithm
$\ln({Q^{2}}/{-l^{2}})$. Besides, we also find (iv) UV poles
${1}/{\epsilon_{\text{\tiny{UV}}}}$. It is interesting to discuss the origin
of these poles. The UV poles arise from all the diagrams in Fig. 2 excluding
the first. To cancel them, we need to include the square root of the self-
energy corrections
$0=\frac{1}{\epsilon_{\text{\tiny{UV}}}}-\frac{1}{\epsilon_{\text{\tiny{IR}}}}$
on the incoming and outgoing (massless) quark lines. This converts the UV
poles into single IR poles that add to the ones in (ii). The double IR poles
arise from the first diagram while the single IR poles arise from all diagrams
except the fourth. It is interesting to note that, in inclusive DIS, the
second and third diagrams also give rise to double poles, canceling the one
from the first diagram, but in the present case they do not because the quark
lines (after re-absorption of gluons) are massive. Another interesting feature
is that, in the usual $\overline{\rm MS}$ scheme, one obtains the evolution
kernel of GPDs as the coefficient of single IR poles for which all the
aforementioned diagrams contribute. In our case, we reproduce the kernel as
the coefficient of the logarithm (iii) for which only the first diagram
contributes.

In the antisymmetric case, we project the result onto
$\bar{u}(p_{2})\gamma^{+}\gamma_{5}u(p_{1})$ in order to extract the twist-two
component. This makes it necessary to specify the treatment of the Dirac
matrix in $d\neq 4$ dimensions. We have used both the ‘naive’ fully
anticommuting $\gamma_{5}$ and the HVBM scheme ’t Hooft and Veltman (1972);
Breitenlohner and Maison (1977). For the latter we have computed the Dirac
traces using the Mathematica package ‘Tracer’ Jamin and Lautenbacher (1993).
The HVBM scheme provides the preferred scheme Weber (1992); Vogelsang (1996)
because, unlike the one with the fully anticommuting $\gamma_{5}$, it is known
to be algebraically consistent. Remarkably, however, in the present case the
result is the same for both schemes. The reason is that all $1/\epsilon$ pole
terms we find enter with a part of the Dirac trace that manifestly gives the
same answer for both treatments of $\gamma_{5}$. All further collinear
singularities are regularized by the logarithm $\ln({Q^{2}}/{-l^{2}})$, so
that this part of the calculation can essentially be carried out in four
dimensions where of course both schemes coincide. The same is true for the
gluonic coefficient function which has no poles in $1/\epsilon$ and can be
obtained in four dimensions. Hence there are no ambiguities related to the
Levi-Civita tensor being an entirely four-dimensional object. The issue of the
scheme nevertheless will show up later when we discuss the one-loop
calculation of the GPDs.

##  IV Results

###  IV.1 One-loop evolution kernels

The coefficients of the logarithm $\ln({Q^{2}/}{-l^{2}})$ in (7), (8) are
dictated by the evolution of the twist-two GPDs and therefore must agree with
the known results in the literature Ji and Osborne (1998); Belitsky and
Mueller (1998); Mankiewicz et al. (1998). We find that this is indeed the
case, meaning that the physical parameter $t$ does the job of regularizing the
collinear singularity associated with GPDs. For completeness, here we
reproduce the results:

$\displaystyle\kappa_{qq}(\hat{x},\hat{\xi})$ $\displaystyle=$
$\displaystyle\frac{3}{2(1-\hat{x})}+\frac{\hat{x}^{2}+1-2\hat{\xi}^{2}}{(1-\hat{\xi}^{2})(1-\hat{x})}\ln\frac{\hat{x}-1}{\hat{x}}+\frac{(\hat{x}-\hat{\xi})(1-\hat{x}^{2}-2\hat{x}\hat{\xi}-2\hat{\xi}^{2})}{(1-\hat{x}^{2})\hat{\xi}(1-\hat{\xi}^{2})}\ln\frac{\hat{x}-\hat{\xi}}{\hat{x}}+(\hat{x}\to-\hat{x}),$
(20) $\displaystyle\tilde{\kappa}_{qq}(\hat{x},\hat{\xi})$ $\displaystyle=$
$\displaystyle\frac{3}{2(1-\hat{x})}+\frac{\hat{x}^{2}+1-2\hat{\xi}^{2}}{(1-\hat{\xi}^{2})(1-\hat{x})}\ln\frac{\hat{x}-1}{\hat{x}}-\frac{(\hat{x}-\hat{\xi})(1+\hat{x}^{2}+2\hat{x}\hat{\xi})}{(1-\hat{x}^{2})(1-\hat{\xi}^{2})}\ln\frac{\hat{x}-\hat{\xi}}{\hat{x}}-(\hat{x}\to-\hat{x}),$
(21) $\displaystyle\kappa_{qg}(\hat{x},\hat{\xi})$ $\displaystyle=$
$\displaystyle\frac{1-2\hat{x}+2\hat{x}^{2}-\hat{\xi}^{2}}{(1-\hat{\xi}^{2})^{2}}\ln\frac{\hat{x}-1}{\hat{x}}+\frac{(\hat{x}-\hat{\xi})(1-2\hat{\xi}\hat{x}-\hat{\xi}^{2})}{\hat{\xi}(1-\hat{\xi}^{2})^{2}}\ln\frac{\hat{x}-\hat{\xi}}{\hat{x}}+(\hat{x}\to-\hat{x}),$
(22) $\displaystyle\tilde{\kappa}_{qg}(\hat{x},\hat{\xi})$ $\displaystyle=$
$\displaystyle\frac{2\hat{x}-1-\hat{\xi}^{2}}{(1-\hat{\xi}^{2})^{2}}\ln\frac{\hat{x}-1}{\hat{x}}-2\frac{\hat{x}-\hat{\xi}}{(1-\hat{\xi}^{2})^{2}}\ln\frac{\hat{x}-\hat{\xi}}{\hat{x}}-(\hat{x}\to-\hat{x}),$
(23)

where as before $\hat{x}=\frac{x_{B}}{x}$, $\hat{\xi}=\frac{\xi}{x}$. Note
that $\hat{x},\hat{\xi}$ are always positive because we restricted to $0<x<1$
in (4) and (5). Also, an infinitesimal, negative imaginary part is understood
in $\hat{x}$, namely, $\hat{x}\to\hat{x}-i\epsilon$ and

$\displaystyle\ln(\hat{x}-1)=\ln(1-\hat{x})-i\pi.$ (24)

###  IV.2 Coefficient functions

The ‘coefficient functions’ in (7), (8) are given as follows:

$\displaystyle\delta C_{1}^{q}(\hat{x},\hat{\xi})$
$\displaystyle=-\dfrac{\bigg{(}\dfrac{Q^{2}}{-l^{2}}\bigg{)}^{\epsilon_{\text{\tiny{IR}}}}}{\epsilon_{\text{\tiny{IR}}}^{2}(1-\hat{x})}-\dfrac{3\bigg{(}\dfrac{Q^{2}}{-l^{2}}\bigg{)}^{\epsilon_{\text{\tiny{IR}}}}}{2\epsilon_{\text{\tiny{IR}}}(1-\hat{x})}+\dfrac{1-2\hat{x}-2\hat{x}^{2}+3\hat{\xi}^{2}}{2(1-\hat{x})(1-\hat{\xi}^{2})}\ln\dfrac{\hat{x}-1}{\hat{x}}+\dfrac{(\hat{x}-\hat{\xi})(-1+\hat{x}^{2}+3\hat{x}\hat{\xi}+3\hat{\xi}^{2})}{(1-\hat{x}^{2})(1-\hat{\xi}^{2})\hat{\xi}}\ln\dfrac{\hat{x}-\hat{\xi}}{\hat{x}}$
$\displaystyle+\dfrac{1+\hat{x}^{2}-2\hat{\xi}^{2}}{2(1-\hat{x})(1-\hat{\xi}^{2})}\ln^{2}\dfrac{\hat{x}-1}{\hat{x}}+\dfrac{\hat{x}}{2(1-\hat{\xi}^{2})\hat{\xi}}\ln^{2}\dfrac{\hat{x}-\hat{\xi}}{\hat{x}}+\dfrac{-1-\hat{x}^{2}+2\hat{\xi}^{2}}{2(1-\hat{x}^{2})(1-\hat{\xi}^{2})}\ln\dfrac{\hat{x}-\hat{\xi}}{\hat{x}}\ln\dfrac{\hat{x}+\hat{\xi}}{\hat{x}}+\dfrac{\pi^{2}-54}{12(1-\hat{x})}$
$\displaystyle+\dfrac{\hat{x}}{(1-\hat{\xi}^{2})\hat{\xi}}\textrm{Li}_{2}\dfrac{-2\hat{\xi}}{\hat{x}-\hat{\xi}}+\dfrac{1+\hat{x}^{2}-2\hat{\xi}^{2}}{(1-\hat{x})(1-\hat{\xi}^{2})}\left(\textrm{Li}_{2}\dfrac{1-\hat{\xi}}{1-\hat{x}}+\textrm{Li}_{2}\dfrac{1+\hat{\xi}}{1-\hat{x}}\right)\,+(\hat{x}\to-\hat{x}),$
(25) $\displaystyle\delta\tilde{C}_{1}^{q}(\hat{x},\hat{\xi})$
$\displaystyle=-\dfrac{\bigg{(}\dfrac{Q^{2}}{-l^{2}}\bigg{)}^{\epsilon_{\text{\tiny{IR}}}}}{\epsilon_{\text{\tiny{IR}}}^{2}(1-\hat{x})}-\dfrac{3\bigg{(}\dfrac{Q^{2}}{-l^{2}}\bigg{)}^{\epsilon_{\text{\tiny{IR}}}}}{2\epsilon_{\text{\tiny{IR}}}(1-\hat{x})}+\dfrac{-1+2\hat{x}-4\hat{x}^{2}+3\hat{\xi}^{2}}{2(1-\hat{x})(1-\hat{\xi}^{2})}\ln\dfrac{\hat{x}-1}{\hat{x}}+\dfrac{(\hat{x}-\hat{\xi})(1+2\hat{x}^{2}+3\hat{x}\hat{\xi})}{(1-\hat{x}^{2})(1-\hat{\xi}^{2})}\ln\dfrac{\hat{x}-\hat{\xi}}{\hat{x}}$
$\displaystyle+\dfrac{1+\hat{x}^{2}-2\hat{\xi}^{2}}{2(1-\hat{x})(1-\hat{\xi}^{2})}\ln^{2}\dfrac{\hat{x}-1}{\hat{x}}+\dfrac{\hat{\xi}}{2(1-\hat{\xi}^{2})}\ln^{2}\dfrac{\hat{x}-\hat{\xi}}{\hat{x}}-\dfrac{\hat{x}(1+\hat{x}^{2}-2\hat{\xi}^{2})}{2(1-\hat{x}^{2})(1-\hat{\xi}^{2})}\ln\dfrac{\hat{x}-\hat{\xi}}{\hat{x}}\ln\dfrac{\hat{x}+\hat{\xi}}{\hat{x}}+\dfrac{\pi^{2}-54}{12(1-\hat{x})}$
$\displaystyle+\dfrac{\hat{\xi}}{1-\hat{\xi}^{2}}\textrm{Li}_{2}\dfrac{-2\hat{\xi}}{\hat{x}-\hat{\xi}}+\dfrac{1+\hat{x}^{2}-2\hat{\xi}^{2}}{(1-\hat{x})(1-\hat{\xi}^{2})}\left(\textrm{Li}_{2}\dfrac{1-\hat{\xi}}{1-\hat{x}}+\textrm{Li}_{2}\dfrac{1+\hat{\xi}}{1-\hat{x}}\right)-(\hat{x}\to-\hat{x}),$
(26) $\displaystyle\delta C_{1}^{g}(\hat{x},\hat{\xi})$
$\displaystyle=-\dfrac{1-2\hat{x}+2\hat{x}^{2}-\hat{\xi}^{2}}{(1-\xi^{2})^{2}}\ln\dfrac{\hat{x}-1}{\hat{x}}+\dfrac{1-2\hat{x}+2\hat{x}^{2}-\hat{\xi}^{2}}{2(1-\hat{\xi}^{2})^{2}}\ln^{2}\dfrac{\hat{x}-1}{\hat{x}}$
$\displaystyle-\dfrac{(\hat{x}-\hat{\xi})(1-2\hat{x}\hat{\xi}-\hat{\xi}^{2})}{\hat{\xi}(1-\hat{\xi}^{2})^{2}}\ln\dfrac{\hat{x}-\hat{\xi}}{\hat{x}}+\dfrac{\hat{x}(1+\hat{\xi}^{2})}{2\xi(1-\hat{\xi}^{2})^{2}}\ln^{2}\dfrac{\hat{x}-\hat{\xi}}{\hat{x}}-\dfrac{1+2\hat{x}^{2}-\hat{\xi}^{2}}{2(1-\hat{\xi}^{2})^{2}}\ln\dfrac{\hat{x}-\hat{\xi}}{\hat{x}}\ln\dfrac{\hat{x}+\hat{\xi}}{\hat{x}}$
$\displaystyle+\dfrac{\hat{x}(1+\hat{\xi}^{2})}{\hat{\xi}(1-\hat{\xi}^{2})^{2}}\textrm{Li}_{2}\dfrac{-2\hat{\xi}}{\hat{x}-\hat{\xi}}+\dfrac{1-2\hat{x}+2\hat{x}^{2}-\hat{\xi}^{2}}{(1-\hat{\xi}^{2})^{2}}\left(\textrm{Li}_{2}\dfrac{1-\hat{\xi}}{1-\hat{x}}+\textrm{Li}_{2}\dfrac{1+\hat{\xi}}{1-\hat{x}}\right)\,+(\hat{x}\to-\hat{x}),$
(27) $\displaystyle\delta\tilde{C}_{1}^{g}(\hat{x},\hat{\xi})$
$\displaystyle=-\dfrac{2\hat{x}-1-\hat{\xi}^{2}}{(1-\hat{\xi}^{2})^{2}}\ln\dfrac{\hat{x}-1}{\hat{x}}+\dfrac{2\hat{x}-1-\hat{\xi}^{2}}{2(1-\hat{\xi}^{2})^{2}}\ln^{2}\dfrac{\hat{x}-1}{\hat{x}}$
$\displaystyle+2\dfrac{\hat{x}-\hat{\xi}}{(1-\hat{\xi}^{2})^{2}}\ln\dfrac{\hat{x}-\hat{\xi}}{\hat{x}}+\dfrac{\hat{\xi}}{(1-\hat{\xi}^{2})^{2}}\ln^{2}\dfrac{\hat{x}-\hat{\xi}}{\hat{x}}-\dfrac{\hat{x}}{(1-\hat{\xi}^{2})^{2}}\ln\dfrac{\hat{x}-\hat{\xi}}{\hat{x}}\ln\dfrac{\hat{x}+\hat{\xi}}{\hat{x}}$
$\displaystyle+\dfrac{2\hat{\xi}}{(1-\hat{\xi}^{2})^{2}}\textrm{Li}_{2}\dfrac{-2\hat{\xi}}{\hat{x}-\hat{\xi}}+\dfrac{2\hat{x}-1-\hat{\xi}^{2}}{(1-\hat{\xi}^{2})^{2}}\left(\textrm{Li}_{2}\dfrac{1-\hat{\xi}}{1-\hat{x}}+\textrm{Li}_{2}\dfrac{1+\hat{\xi}}{1-\hat{x}}\right)\,-(\hat{x}\to-\hat{x}),$
(28)

where ${\rm Li}_{2}$ is the dilogarithm function. As mentioned before, in the
quark channel, we find single $1/\epsilon_{\text{\tiny{IR}}}$ and double
$1/\epsilon_{\text{\tiny{IR}}}^{2}$ infrared poles. Note that, since
$\epsilon_{\text{\tiny{IR}}}<0$,
$({Q^{2}/}{-l^{2}})^{\epsilon_{\text{\tiny{IR}}}}\to 0$ if one takes the
$l^{2}\to 0$ limit first. However, if one keeps $l^{2}$ finite and expands in
$\epsilon_{\text{\tiny{IR}}}$, the first term $\frac{3}{2(1-\hat{x})}$ in (20)
and (21) gets canceled. We discuss below how these problematic terms
eventually disappear. As we also mentioned, the result (26) is independent of
the scheme for $\gamma_{5}$.

###  IV.3 Anomaly pole terms

The coefficients of the ‘anomaly poles’ $1/t$ in (4) and (5) are found to be

$\displaystyle A(x,x_{B},\xi)$ $\displaystyle=$
$\displaystyle\frac{2T_{R}}{x}\left(1+\frac{\hat{x}(1-\hat{x})\ln\dfrac{\hat{x}-1}{\hat{x}}+\hat{x}(\hat{x}-\hat{\xi})\ln\dfrac{\hat{x}-\hat{\xi}}{\hat{x}}+(\hat{x}\to-\hat{x})}{1-\hat{\xi}^{2}}\right),$
(29) $\displaystyle\tilde{A}(x,x_{B},\xi)$ $\displaystyle=$
$\displaystyle\frac{8T_{R}}{x}\frac{(1-\hat{x})\ln\frac{\hat{x}-1}{x}+(\hat{x}-\hat{\xi})\ln\frac{\hat{x}-\hat{\xi}}{\hat{x}}-(\hat{x}\to-\hat{x})}{1-\hat{\xi}^{2}}.$
(30)

The imaginary parts of these expressions agree with our results in
Bhattacharya et al. (2023). Clearly, there are $1/t$ poles also in the real
part of the Compton amplitude.222In Ref. Bhattacharya et al. (2023) we
computed only the imaginary part of the Compton amplitude by directly
extracting the discontinuity across the variables $s=(p+q)^{2}$ and
$q^{2}_{2}$. In this paper, we compute the full amplitude. We have checked
that the imaginary parts of (29), (30) and all the other results in this paper
are consistent with the corresponding results in Bhattacharya et al. (2023).
While (29) and (30) look unfamiliar and complicated, remarkably the
$x$-integrals in (4) and (5) can be exactly rewritten in the following form

$\displaystyle\int_{0}^{1}dxA(x,x_{B},\xi){\cal F}(x,\xi,t)$
$\displaystyle=2T_{R}\int_{0}^{1}dxC_{0}(x,x_{B})\left[\int_{x}^{1}\frac{dx^{\prime}}{x^{\prime}}K\left(\frac{x}{x^{\prime}},\frac{\xi}{x^{\prime}}\right){\cal
F}(x^{\prime},\xi,t)-\theta(\xi-x)\int_{0}^{1}\frac{dx^{\prime}}{x^{\prime}}L\left(\frac{x}{x^{\prime}},\frac{\xi}{x^{\prime}}\right){\cal
F}(x^{\prime},\xi,t)\right]$ $\displaystyle\equiv
2T_{R}\int_{0}^{1}dxC_{0}(x,x_{B})C^{\rm anom}\otimes{\cal F}(x,\xi,t),$ (31)
$\displaystyle\int_{0}^{1}dx\tilde{A}(x,x_{B},\xi)\tilde{\cal F}(x,\xi,t)$
$\displaystyle=2T_{R}\int_{0}^{1}dx\tilde{C}_{0}(x,x_{B})\left[\int_{x}^{1}\frac{dx^{\prime}}{x^{\prime}}\tilde{K}\left(\frac{x}{x^{\prime}},\frac{\xi}{x^{\prime}}\right)\tilde{\cal
F}(x^{\prime},\xi,t)-\theta(\xi-x)\int_{0}^{1}\frac{dx^{\prime}}{x^{\prime}}\tilde{L}\left(\frac{x}{x^{\prime}},\frac{\xi}{x^{\prime}}\right)\tilde{\cal
F}(x^{\prime},\xi,t)\right]$ $\displaystyle\equiv
2T_{R}\int_{0}^{1}dx\tilde{C}_{0}(x,x_{B})\tilde{C}^{\rm
anom}\otimes\tilde{\cal F}(x,\xi,t),$ (32)

where

$\displaystyle K(x,\xi)=\frac{x(1-x)}{1-\xi^{2}}\,,\qquad
L(x,\xi)=\frac{x(\xi-x)}{1-\xi^{2}}\,,$ (33)
$\displaystyle\tilde{K}(x,\xi)=\frac{4(1-x)}{1-\xi^{2}},\qquad\tilde{L}(x,\xi)=\frac{4(\xi-x)}{1-\xi^{2}}.$
(34)

That is, the leading-order kernels $C_{0}$ and $\tilde{C}_{0}$ can be factored
out. The resulting convolution $C^{\rm anom}\otimes{\cal F}$ agrees with what
was anticipated in Bhattacharya et al. (2023) following the general argument
in White (2002) where actually the same integral (31) can be found. We now
have the corresponding result with $\xi\neq 0$ in the polarized sector. As
mentioned already in Bhattacharya et al. (2023), the two terms in $C^{\rm
anom}$ and $\bar{C}^{\rm anom}$ come from the first and third diagrams in Fig.
1. The latter is nonzero only when the outgoing photon becomes timelike,
$q_{2}^{2}>0$, see (13).

The identities (31), (32) guarantee that, if the $1/t$ poles are cancelled in
the imaginary part of the Compton amplitude Bhattacharya et al. (2023), the
same cancellation automatically occurs in the real part as well.

##  V GPD at one loop

We have seen that the Compton scattering amplitudes at one loop contain three
types of singular behaviors: (i) logarithms $\ln(-t)$, (ii) anomaly poles
$1/t$, (iii) single $1/\epsilon$ and double $1/\epsilon^{2}$ infrared poles
(only in the quark channel). The logarithms are as expected, but the other two
are unusual and potentially cause problems with factorization. We now
demonstrate that all these singular structures can be absorbed into the quark
GPDs in the leading-order terms of (4) and (5). Specifically, we compute the
unpolarized and polarized quark GPDs333The variables $x$ and $\xi$ in this
section (and also in the appendix) should better be written as $\hat{x}$ and
$\hat{\xi}$ to be more consistent with the notation in the previous sections.
We however abbreviate $\hat{x},\hat{\xi}\to x,\xi$ for simplicity.

$\displaystyle f_{q}(x,\xi,t)$ $\displaystyle=$
$\displaystyle\int\frac{dz^{-}}{4\pi}e^{ixP^{+}z^{-}}\langle
p_{2}|\bar{q}(-z/2)W\gamma^{+}q(z^{-}/2)|p_{1}\rangle,$ (35)
$\displaystyle\tilde{f}_{q}(x,\xi,t)$ $\displaystyle=$
$\displaystyle\int\frac{dz^{-}}{4\pi}e^{ixP^{+}z^{-}}\langle
p_{2}|\bar{q}(-z/2)W\gamma^{+}\gamma_{5}q(z^{-}/2)|p_{1}\rangle,$ (36)

to one loop for on-shell $p_{1}^{2}=p_{2}^{2}=0$ quark and gluon external
states, keeping $t=(p_{2}-p_{1})^{2}\neq 0$. We need to separately consider
the DGLAP region $0<\xi<x\leq 1$ and the Efremov-Radyushkin-Brodsky-Lepage
(ERBL) region $0<x<\xi$ Diehl (2003). We work in $d=4-2\epsilon$ dimensions to
regularize the UV divergences and any leftover IR divergences. As before, they
are distinguished by $1/\epsilon_{\text{\tiny{UV}}}$ and
$1/\epsilon_{\text{\tiny{IR}}}$, respectively. The $\overline{\rm MS}$ scale
is denoted by $\tilde{\mu}^{2}=4\pi e^{-\gamma_{E}}\mu^{2}$.

###  V.1 Quark matrix element

The divergent part of the quark matrix element in the DGLAP region $\xi<x<1$
is, omitting the common prefactor $\frac{\alpha_{s}C_{F}}{4\pi
p^{+}}\bar{u}(p_{2})\gamma^{+}u(p_{1})$ or $\frac{\alpha_{s}C_{F}}{4\pi
p^{+}}\bar{u}(p_{2})\gamma^{+}\gamma_{5}u(p_{1})$,

$\displaystyle\frac{\left(\frac{\tilde{\mu}^{2}}{-l^{2}}\right)^{\epsilon}}{\epsilon_{\text{\tiny{UV}}}}\left[\frac{1+x^{2}-2\xi^{2}}{(1-\xi^{2})(1-x)_{+}}+\left(\frac{3}{2}-\ln(1-\xi^{2})\right)\delta(1-x)\right]-\delta(1-x)\left(\frac{1}{\epsilon^{2}_{\text{\tiny{IR}}}}+\frac{3}{2}\frac{1}{\epsilon_{\text{\tiny{IR}}}}\right)\left(\frac{\tilde{\mu}^{2}}{-l^{2}}\right)^{\epsilon}.$
(37)

The derivation is outlined in Appendix A. The same result holds for both the
unpolarized and polarized GPDs. Note the double and single infrared poles.
Similar ‘Sudakov’ double poles have been previously encountered in the one-
loop gluon matrix element of the twist-four GPDs (9) and (10) Radyushkin and
Zhao (2021, 2022). After convoluting with $C_{0}$ and $\tilde{C}_{0}$, which
can be done trivially using the delta function $\delta(1-x)$, they exactly
match the double and single poles in (25) and (26). In other words, the poles
in (25) and (26) can be absorbed into the GPDs. The finite terms are given by

$\displaystyle-\frac{1+x^{2}-2\xi^{2}}{1-\xi^{2}}\left(\frac{\ln\frac{(1-x)^{2}}{1-\xi^{2}}}{1-x}\right)_{+}-\frac{1-x}{1-\xi^{2}}-\frac{1}{2}\delta(1-x)\left(\ln^{2}(1-\xi^{2})-\frac{\pi^{2}}{6}\right),$
(38)

in the unpolarized case. In the polarized case, interestingly the result
depends on the scheme for $\gamma_{5}$ in contrast to what we observed with
the coefficient function (26). If one uses the ‘fully anticommuting
$\gamma_{5}$,’ the result is the same as (38). If, on the other hand, one uses
the HVBM scheme, one finds an additional term

$\displaystyle 4\frac{1-x}{1-\xi^{2}}.$ (39)

The $\mathcal{O}(\epsilon)$ difference between the two schemes now leaves a
finite contribution because of the presence of the UV pole
$1/\epsilon_{\text{\tiny{UV}}}$.

In the ERBL region, the divergent terms read

$\displaystyle\frac{\left(\frac{\tilde{\mu}^{2}}{-l^{2}}\right)^{\epsilon}}{\epsilon_{\text{\tiny{UV}}}}\frac{x+\xi}{(1+\xi)}\left(\frac{1}{1-x}+\frac{1}{2\xi}\right),$
(40)

in both the unpolarized and polarized cases. The finite terms are rather
cumbersome. In the unpolarized case, we find

$\displaystyle\frac{1}{2\xi(1-x)(1-\xi^{2})}\Biggl{[}(x-1)(x+\xi)(\xi+1)\ln(\xi^{2}-x^{2})+2\xi
x^{2}\ln\frac{1+\xi}{(x+\xi)(1-x)}+2\xi^{3}\ln\frac{x+\xi}{\xi-x}$
$\displaystyle\qquad+2\xi\ln\frac{(1+\xi)(\xi-x)}{1-x}+4\xi^{3}\ln\frac{1-x}{1+\xi}+2(1-x)(x+\xi^{2})\ln(2\xi)\Biggr{]}-\frac{x+\xi}{2\xi(1+\xi)},$
(41)

and in the polarized case there is an additional term

$\displaystyle 4\frac{x+\xi}{2\xi(1+\xi)},$ (42)

in the HVBM scheme.

The coefficients of the UV pole in (37) and (40) constitute the $q\to q$
evolution kernel of the GPDs. After expanding
$\left(\frac{\tilde{\mu}^{2}}{-l^{2}}\right)^{\epsilon}=1+\epsilon\ln\frac{\tilde{\mu}^{2}}{-l^{2}}$
and convoluting with $C_{0}$ and $\tilde{C}_{0}$, we recover (20) and (21). In
other words, the logarithmic terms in (7) can be absorbed into the GPDs, as
expected. The same comment applies to the $g\to q$ evolution kernel below.

###  V.2 Gluon matrix element, unpolarized

For the gluon matrix elements $\langle
p_{2}\epsilon_{2}|...|p_{1}\epsilon_{1}\rangle$ of the quark GPD (35), (36),
we find it convenient to work in the light-cone gauge
$\epsilon_{1}^{+}=\epsilon_{2}^{+}=0$. The result for the unpolarized GPD in
the DGLAP region is

$\displaystyle f_{q}$ $\displaystyle=$
$\displaystyle\frac{\alpha_{s}T_{R}}{2\pi}\Biggl{[}-(1-\xi^{2})\epsilon_{1}\cdot\epsilon_{2}^{*}\left(\frac{\left(\frac{\tilde{\mu^{2}}}{-l^{2}}\right)^{\epsilon}}{\epsilon_{\rm
UV}}\frac{2x^{2}-2x+1-\xi^{2}}{(1-\xi^{2})^{2}}-\frac{(2x^{2}-2x+1-\xi^{2})\ln\frac{(1-x)^{2}}{1-\xi^{2}}+2x(1-x)}{(1-\xi^{2})^{2}}\right)$
(43)
$\displaystyle\qquad\qquad-\frac{4}{l^{2}}\frac{x(1-x)}{1-\xi^{2}}\left(\epsilon_{1}\cdot
l\epsilon_{2}^{*}\cdot
l-\frac{l^{2}}{2}\epsilon_{1}\cdot\epsilon_{2}^{*}\right)\Biggr{]},$

where we have factored out the structure that represents the twist-two GPD,
see (18). Note the anomaly pole $1/l^{2}$ which was absent in the quark matrix
elements. Its coefficient matches $K$ in (33). (The factor of 4 is from (9).)

In the ERBL region, we find, omitting the prefactor
$\frac{\alpha_{s}T_{R}}{2\pi}$,

$\displaystyle-(1-\xi^{2})\epsilon_{1}\cdot\epsilon_{2}^{*}\frac{\left(\frac{\tilde{\mu}^{2}}{-l^{2}}\right)^{\epsilon}}{\epsilon_{\text{\tiny{UV}}}}\frac{(x+\xi)(1+\xi-2x)}{2\xi(1+\xi)(1-\xi^{2})}-\left(\epsilon_{1}\cdot
l\epsilon_{2}^{*}\cdot
l-\frac{l^{2}}{2}\epsilon_{1}\cdot\epsilon_{2}^{*}\right)\frac{1}{l^{2}}\frac{2x(x+\xi)}{\xi(1+\xi)}-(x\to-x),$
(44)

The first term comes from the same ladder diagram as in (43). The $x\to-x$
term comes from a diagram in which the gluon legs are crossed. The latter
contributes only in the ERBL region. The finite terms read, including the
$x\to-x$ contribution,

$\displaystyle-(1-\xi^{2})\epsilon_{1}\cdot\epsilon_{2}^{*}\frac{1}{\xi(1-\xi^{2})^{2}}\Biggl{[}-x(1+\xi^{2})\ln\frac{\xi^{2}-x^{2}}{4}-\xi(1+2x^{2}-\xi^{2})\ln\frac{(x+\xi)(1-x)}{(\xi-x)(1+x)}$
$\displaystyle+2x\Bigl{(}\xi\ln(1-x^{2})-2\xi\ln(1+\xi)+\xi^{2}-\xi+(1+\xi^{2})\ln\xi\Bigr{)}\Biggr{]}.$
(45)

The coefficient of the anomaly pole is thus

$\displaystyle\frac{1}{l^{2}}\left(\frac{2x(x+\xi)}{\xi(1+\xi)}-\frac{2x(x-\xi)}{\xi(1+\xi)}\right)=\frac{4}{l^{2}}\frac{x}{1+\xi}.$
(46)

This agrees with

$\displaystyle K(x,\xi)-L(x,\xi)=\frac{x}{1+\xi},$ (47)

which is the relevant linear combination in the ERBL region, see (31). These
results, together with the observation (31), show that the $1/t$ pole in (4)
can be absorbed into the unpolarized quark GPDs $H_{q},\,E_{q}$ in the leading
order as a part of the infrared subtraction procedure.

###  V.3 Gluon matrix element, polarized

In the polarized case (36), we find

$\displaystyle\tilde{f}_{q}=\frac{\alpha_{s}T_{R}}{2\pi}\left[(1-\xi^{2})i\epsilon^{+p\epsilon_{2}^{*}\epsilon_{1}}\left(\frac{2x-1-\xi^{2}}{(1-\xi^{2})^{2}}\left(\frac{\left(\frac{\tilde{\mu}^{2}}{-l^{2}}\right)^{\epsilon}}{\epsilon_{\text{\tiny{UV}}}}-\ln\frac{(1-x)^{2}}{1-\xi^{2}}\right)-2\frac{1-x}{(1-\xi^{2})^{2}}\right)+\frac{2il^{+}\epsilon^{\epsilon_{1}\epsilon_{2}^{*}lp}}{l^{2}}\frac{1-x}{1-\xi^{2}}\right],$
(48)

in the DGLAP region. A simplified version of this result (with a different UV
prescription) was already reported in Bhattacharya et al. (2023). The
coefficient of the pole agrees with $\tilde{K}$ in (34). Note that the pole is
proportional to $l^{\mu=+}$, meaning that it contributes to a shift in the GPD
$\tilde{E}_{q}$. In the ERBL region, the singular terms are, omitting the
prefactor $\frac{\alpha_{s}T_{R}}{2\pi}$,

$\displaystyle\frac{x+\xi}{2\xi(1+\xi)}\left((1-\xi^{2})i\epsilon^{+p\epsilon_{2}^{*}\epsilon_{1}}\frac{\left(\frac{\tilde{\mu}^{2}}{-l^{2}}\right)^{\epsilon}}{\epsilon_{\text{\tiny{UV}}}}\frac{-1}{1+\xi}+\frac{2il^{+}\epsilon^{\epsilon_{1}\epsilon_{2}^{*}lp}}{l^{2}}\right)+(x\to-x)$
$\displaystyle\quad=(1-\xi^{2})i\epsilon^{+p\epsilon_{2}^{*}\epsilon_{1}}\frac{\left(\frac{\tilde{\mu}^{2}}{-l^{2}}\right)^{\epsilon}}{\epsilon_{\text{\tiny{UV}}}}\frac{-1}{(1+\xi)^{2}}+\frac{2il^{+}\epsilon^{\epsilon_{1}\epsilon_{2}^{*}lp}}{l^{2}}\frac{1}{1+\xi},$
(49)

where again the $x\to-x$ term comes from the crossed-leg diagram. The finite
terms are, including the $x\to-x$ contribution,

$\displaystyle(1-\xi^{2})i\epsilon^{+p\epsilon_{2}^{*}\epsilon_{1}}\frac{1}{(1-\xi^{2})^{2}}\Biggl{[}-2\xi\ln(\xi^{2}-x^{2})+(1+\xi^{2})\ln(1-x^{2})-2x\ln\frac{(1-x)(x+\xi)}{(1+x)(\xi-x)}$
$\displaystyle\qquad\qquad-2(1+\xi^{2})\ln(1+\xi)+4\xi\ln(2\xi)+2\xi-2\Biggr{]}.$
(50)

The coefficient of the UV pole $\frac{-1}{(1+\xi)^{2}}$ is the correct
evolution kernel in the ERBL region as can be seen by taking the imaginary
part of (23):

$\displaystyle\frac{2x-1-\xi^{2}}{(1-\xi^{2})^{2}}-\frac{2(x-\xi)}{(1-\xi^{2})^{2}}=\frac{-1}{(1+\xi)^{2}}.$
(51)

Again the coefficient of the anomaly pole agrees with

$\displaystyle\tilde{K}(x,\xi)-\tilde{L}(x,\xi)=\frac{4}{1+\xi},$ (52)

from (34). With the help of (32), we can absorb the $1/t$ pole in (5) into the
polarized quark GPD $\tilde{E}_{q}$.

###  V.4 Relation to the $\overline{\rm MS}$ scheme

We have thus shown that all the singular structures $1/t$, $\ln(-t)$,
$1/\epsilon_{\rm IR}$ and $1/\epsilon^{2}_{\rm IR}$ in the ‘unsubtracted’
expressions (4) and (5) can be systematically absorbed into the twist-two GPDs
in the leading order. Since the matrix elements (35) and (36) contain non-
singular terms, one might choose to perform this subtraction also for the
finite terms (25)-(28). An interesting question then arises as to whether,
after such a subtraction, (25)-(28) reduce to the known coefficient functions
in the $\overline{\rm MS}$ scheme.444We thank Vladimir Braun for raising this
question. Here we partially address this question by explicitly performing the
subtraction in the imaginary part of the Compton amplitude.

Let us first consider the DGLAP region $0<\xi<x$. For simplicity, we assume
$x<1$ to avoid the delta function $\delta(1-x)$. The imaginary part of (25) is

$\displaystyle\frac{1-2x-2x^{2}+3\xi^{2}}{2(1-x)(1-\xi^{2})}+\frac{1+x^{2}-2\xi^{2}}{(1-x)(1-\xi^{2})}\ln\frac{1-\xi^{2}}{x(1-x)},$
(53)

where we used

$\displaystyle{\rm Im}\,{\rm
Li}_{2}\frac{1\pm\xi}{1-x+i\epsilon}=-\pi\ln\frac{1\pm\xi}{1-x}.$ (54)

On the other hand, the finite terms in the unpolarized quark GPD are, from
(38),

$\displaystyle-\frac{1+x^{2}-2\xi^{2}}{(1-x)(1-\xi^{2})}\ln\frac{(1-x)^{2}}{1-\xi^{2}}-\frac{1-x}{1-\xi^{2}}.$
(55)

The convolution with the leading-order kernel (6) is trivial for the imaginary
part since ${\rm Im}\,C_{0}\propto\delta(1-x)$. We just need to subtract (55)
from (53) to obtain

$\displaystyle\frac{1+x^{2}-2\xi^{2}}{(1-x)(1-\xi^{2})}\ln\frac{1-x}{x}+\frac{3(1-2x+\xi^{2})}{2(1-x)(1-\xi^{2})}\to
2+x-\frac{3}{2(1-x)}+\frac{1+x^{2}}{1-x}\ln\frac{1-x}{x}+1-x,$ (56)

where we have set $\xi=0$ on the right-hand side. This agrees with the
imaginary part of the $q\to q$ coefficient function in the $\overline{\rm MS}$
scheme Ji and Osborne (1998); Belitsky and Mueller (1998); Mankiewicz et al.
(1998); Pire et al. (2011); Bertone et al. (2022). In particular, the right-
hand side is the familiar $q\to q$ coefficient function for the $F_{1}$
structure function in DIS Altarelli et al. (1979) for $x<1$. Similarly, the
imaginary part of (26) reads

$\displaystyle\frac{-1+2x-4x^{2}+3\xi^{2}}{2(1-x)(1-\xi^{2})}+\frac{1+x^{2}-2\xi^{2}}{(1-x)(1-\xi^{2})}\ln\frac{1-\xi^{2}}{x(1-x)}.$
(57)

The finite terms in the polarized quark PDF depend on the scheme adopted for
the treatment of $\gamma_{5}$. In the HVBM scheme, we find from (38) and (39),

$\displaystyle-\frac{1+x^{2}-2\xi^{2}}{(1-x)(1-\xi^{2})}\ln\frac{(1-x)^{2}}{1-\xi^{2}}+3\frac{1-x}{1-\xi^{2}}.$
(58)

After the subtraction,

$\displaystyle\frac{-7+14x-10x^{2}+3\xi^{2}}{2(1-x)(1-\xi^{2})}+\frac{1+x^{2}-2\xi^{2}}{(1-x)(1-\xi^{2})}\ln\frac{1-x}{x}\to
2+x-\frac{3}{2(1-x)}+\frac{1+x^{2}}{1-x}\ln\frac{1-x}{x}-4(1-x),$ (59)

in agreement with the $q\to q$ coefficient function for the $g_{1}$ structure
function in the HVBM prescription. As was discussed in Refs. Stratmann et al.
(1996); Vogelsang (1996), it is necessary to subtract this term in order to
avoid a conflict with helicity conservation and with the known first-order
correction to the Bjorken sum rule. Incidentally, in the present case, the
result obtained after this finite subtraction coincides with that found for a
fully anticommuting $\gamma_{5}$. Either way, instead of (59) the correct
result becomes

$\displaystyle\frac{1-2x-2x^{2}+3\xi^{2}}{2(1-x)(1-\xi^{2})}+\frac{1+x^{2}-2\xi^{2}}{(1-x)(1-\xi^{2})}\ln\frac{1-x}{x}\to
2+x-\frac{3}{2(1-x)}+\frac{1+x^{2}}{1-x}\ln\frac{1-x}{x}.$ (60)

As for the $g\to q$ coefficients, the imaginary part of (27) is

$\displaystyle\frac{1-2x+2x^{2}-\xi^{2}}{(1-\xi^{2})^{2}}\left(\ln\frac{1-\xi^{2}}{x(1-x)}-1\right).$
(61)

From this, we subtract the finite terms in (43),

$\displaystyle-\frac{1-2x+2x^{2}-\xi^{2}}{(1-\xi^{2})^{2}}\ln\frac{(1-x)^{2}}{1-\xi^{2}}-\frac{2x(1-x)}{(1-\xi^{2})^{2}},$
(62)

to obtain

$\displaystyle\frac{1-2x+2x^{2}-\xi^{2}}{(1-\xi^{2})^{2}}\ln\frac{1-x}{x}+\frac{-1+4x-4x^{2}+\xi^{2}}{(1-\xi^{2})^{2}}\to(1-2x+2x^{2})\left(\ln\frac{1-x}{x}-1\right)+2x(1-x).$
(63)

This agrees with the $g\to q$ coefficient function for the $F_{1}$ structure
function in the $\overline{\rm MS}$ scheme. Finally, the imaginary part of
(28) is

$\displaystyle\frac{2x-1-\xi^{2}}{(1-\xi^{2})^{2}}\left(\ln\frac{1-\xi^{2}}{x(1-x)}-1\right).$
(64)

Subtracting from this the finite terms in (48),

$\displaystyle-\frac{2x-1-\xi^{2}}{(1-\xi^{2})^{2}}\ln\frac{(1-x)^{2}}{1-\xi^{2}}-\frac{2(1-x)}{(1-\xi^{2})^{2}},$
(65)

we find

$\displaystyle\frac{2x-1-\xi^{2}}{(1-\xi^{2})^{2}}\ln\frac{1-x}{x}+\frac{3-4x+\xi^{2}}{(1-\xi^{2})^{2}}\to(2x-1)\left(\ln\frac{1-x}{x}-1\right)+2(1-x),$
(66)

in agreement with the $\overline{\rm MS}$ $g\to q$ coefficient function for
the $g_{1}$ structure function. It is interesting to recall that the last term
$2(1-x)\otimes\Delta G(x)$ in (66) caused a lot of discussion (see, e.g.,
Bodwin and Qiu (1990); Vogelsang (1991)) in the wake of the proton ‘spin
crisis’ in the late 80s. In the standard $\overline{\rm MS}$ calculation in
the forward limit $t=0$, this term arises from the IR region of the loop
diagram, and therefore does not seem to qualify as a part of the ‘hard’
coefficient. In our calculation of the Compton amplitude, this term is
replaced by the pole term $\frac{1}{l^{2}}(1-x)\otimes\tilde{\cal F}(x)$ and
gets absorbed into the GPD $\tilde{E}_{q}$. Nevertheless, the $2(1-x)$ term is
restored after the subtraction because the polarized GPD (36) generates it
from the UV region of the loop momentum. Therefore, even though the final
result (66) is the same, from our perspective the term $2(1-x)$ is
legitimately considered a ‘hard’ contribution.

We have further performed the subtraction of the constant terms (41), (45) and
(50) in the ERBL region $x<\xi$ from the imaginary part of (25)-(28) and
observed consistent agreement with the $\overline{\rm MS}$ coefficient
functions Belitsky and Radyushkin (2005). We have thus partially verified that
the “off-forward” regularization is equivalent to the $\overline{\rm MS}$
scheme after the subtraction of finite terms is made. Extending this to the
real part of the Compton amplitude is left for future work. On the other hand,
since the treatment of finite terms is a matter of scheme choice, one can
choose to subtract only the singular terms. Eqs. (25)-(28), with the single
and double poles omitted, are then the coefficient functions in such a scheme.

##  VI Imprints of anomalies on GPD

Let us discuss the implications of our results. Superficially, it may seem as
if nothing has happened in the end. After absorbing all the singular terms
into the twist-two GPDs, the Compton amplitude will be given by the usual
factorized form with possibly different coefficient functions due to a
different scheme choice. The common attitude is that one does not care about
these singular terms once they have been ‘discarded’ into a parton
distribution, as they will be taken care of by the nonperturbative QCD
dynamics. One can also take the view that the $1/t$ poles should disappear in
the limit $t\to 0$, because nonperturbative effects must intervene when
$\sqrt{|t|}\sim\Lambda_{\rm QCD}$. However, from the point of view of
factorization, technically speaking one is allowed to choose any infrared
regulator that can isolate the collinear divergences, as long as they can be
eventually absorbed into parton distributions when the latter are computed
with the same IR regulator. In this sense, the use of $t$ is no different from
other regulators such as the current quark mass $m_{q}$ and the dimensionality
$d\neq 4$. One may even argue that it is a more physical scheme, since $t\neq
0$ in actual experiments and naturally cuts off collinear divergences.

Technicalities aside, the real reason we are pursuing the off-forward
calculation with nonzero $t$ is that this approach has the potential to
uncover nonperturbative connections between GPDs and QCD anomalies. Indeed,
the very idea that twist-four GPDs are absorbed into twist-two GPDs is quite
non-standard and needs to be investigated further, rather than dismissed as a
routine infrared subtraction procedure. This is all the more so because, as
discussed in Tarasov and Venugopalan (2022); Bhattacharya et al. (2023) and
elaborated further below, the results we shall get are consistent with what we
know about the axial and gravitational form factors which are certain moments
of the twist-two GPDs.

###  VI.1 Axial and gravitational form factors

####  VI.1.1 Isovector axial form factors

In order to motivate our discussion, let us first recall the familiar example
of the isovector axial current
$J_{5a}^{\alpha}=\sum_{q}\bar{q}\gamma^{\alpha}\gamma_{5}\frac{\tau^{a}}{2}q$
where $\tau^{a=1,2,3}$ are the Pauli matrices. Its nucleon matrix element is
parameterized by the axial form factors,

$\displaystyle\langle
P_{2}|J_{5a}^{\alpha}|P_{1}\rangle=\bar{u}(P_{2})\left[\gamma^{\alpha}\gamma_{5}F_{A}(t)+\frac{l^{\alpha}\gamma_{5}}{2M}F_{P}(t)\right]\frac{\tau^{a}}{2}u(P_{1})\,.$
(67)

In QCD with $n_{f}=2$ massless flavors, the current is exactly conserved,
$\partial_{\alpha}J^{\alpha}_{5a}=0$, due to chiral symmetry. This imposes a
constraint among the form factors

$\displaystyle 2MF_{A}(t)+\frac{tF_{P}(t)}{2M}=0.$ (68)

Clearly, $F_{P}(t)$ has a pole at $t=0$:

$\displaystyle F_{P}(t)\approx\frac{-4M^{2}g_{A}^{(3)}}{t}\qquad(t\to 0),$
(69)

where $g^{(3)}_{A}=F_{A}(0)\approx 1.3$ is the isovector axial coupling
constant. The pole is generated by the exchange of the massless pion which is
the Nambu-Goldstone boson of spontaneously broken chiral symmetry. This
requirement leads to the well-known Goldberger-Treiman relation

$\displaystyle g_{A}^{(3)}=\frac{f_{\pi}g_{\pi NN}}{M},$ (70)

where $f_{\pi}$ is the pion decay constant and $g_{\pi NN}$ is the pion-
nucleon coupling. Now recall that $F_{P}(t)$ is the first moment of the
isovector GPD $\tilde{E}$,

$\displaystyle
F_{P}(t)=\int_{-1}^{1}dx\left(\tilde{E}_{u}(x,\xi,t)-\tilde{E}_{d}(x,\xi,t)\right).$
(71)

Barring an unlikely possibility that the pole $1/t$ is generated by the
$x$-integral, the GPDs themselves hence have a massless pole at $t=0$:

$\displaystyle\tilde{E}_{u}(x,\xi,t)-\tilde{E}_{d}(x,\xi,t)\sim\theta(\xi-|x|)\frac{g_{A}^{(3)}}{t}\qquad(t\to
0).$ (72)

Indeed, such a pole has been discussed in the GPD literature (see, e.g.,
Penttinen et al. (2000)) where it has been argued that the pole exists only in
the ERBL region $\xi>x$ where GPDs probe the mesonic degrees of freedom inside
the nucleon. In actual QCD with massive quarks, the pole is shifted to the
physical pion mass, $\frac{1}{t}\to\frac{1}{t-m_{\pi}^{2}}$.

####  VI.1.2 Isoscalar axial form factors

The story is more complicated for the singlet axial current
$J_{5}^{\alpha}=\sum_{q}\bar{q}\gamma^{\alpha}\gamma_{5}q$. The associated
form factors $g_{A}$, $g_{P}$, appearing in the nucleon matrix element via

$\displaystyle\langle
P_{2}|J_{5}^{\alpha}|P_{1}\rangle=\bar{u}(P_{2})\left[\gamma^{\alpha}\gamma_{5}g_{A}(t)+\frac{l^{\alpha}\gamma_{5}}{2M}g_{P}(t)\right]u(P_{1})\,,$
(73)

are related to the flavor-singlet polarized GPDs as

$\displaystyle g_{A}(t)$
$\displaystyle=\sum_{q}\int_{-1}^{1}dx\tilde{H}_{q}(x,\xi,t)=\sum_{q}\int_{0}^{1}dx(\tilde{H}_{q}(x,\xi,t)+\tilde{H}_{q}(-x,\xi,t))\,,$
(74) $\displaystyle g_{P}(t)$
$\displaystyle=\sum_{q}\int_{-1}^{1}dx\tilde{E}_{q}(x,\xi,t)=\sum_{q}\int_{0}^{1}dx(\tilde{E}_{q}(x,\xi,t)+\tilde{E}_{q}(-x,\xi,t))\,.$
(75)

In contrast to the isovector current above, $J_{5}^{\alpha}$ is not conserved
due to the chiral (UA(1)) anomaly,

$\displaystyle\partial_{\alpha}J_{5}^{\alpha}=-\frac{n_{f}\alpha_{s}}{4\pi}F^{\mu\nu}\tilde{F}_{\mu\nu}.$
(76)

This leads to the following exact relation:

$\displaystyle 2Mg_{A}(t)+\frac{tg_{P}(t)}{2M}=i\frac{\langle
P_{2}|\frac{n_{f}\alpha_{s}}{4\pi}F\tilde{F}|P_{1}\rangle}{\bar{u}(P_{2})\gamma_{5}u(P_{1})}\,.$
(77)

We see that, in the absence of the anomaly (i.e., if the right-hand side were
zero), $g_{P}(t)$ would have a pole at $t=0$,

$\displaystyle\frac{g_{P}(t)}{2M}\approx-\frac{2M\Delta\Sigma}{t}\qquad(t\to
0),$ (78)

where $\Delta\Sigma=g_{A}(0)$ is the quark helicity contribution to the
nucleon spin. Such a pole can be interpreted as due to the exchange of the
massless ninth Nambu-Goldstone boson, the ‘primordial’ $\eta_{0}$ meson.
Moreover, (75) suggests that already the flavor-singlet GPD
$\sum_{q}\tilde{E}_{q}$ would have a pole $1/t$, just like (72).

In reality, however, the UA(1) axial symmetry is explicitly broken due to the
anomaly, and $g_{P}(t)$ exhibits a pole at the physical $\eta^{\prime}$ meson
mass $t=m^{2}_{\eta^{\prime}}$. The exact mechanism behind this scenario was a
great debate in the late 70s culminating in the works of Witten Witten (1979)
and Veneziano Veneziano (1979). In a nutshell, $\eta_{0}$ acquires mass via a
resummation Veneziano (1979)

$\displaystyle\frac{1}{t}+\frac{m_{\eta^{\prime}}^{2}}{t^{2}}+\frac{m_{\eta^{\prime}}^{4}}{t^{3}}+\cdots=\frac{1}{t-m_{\eta^{\prime}}^{2}}=-\left(\frac{1}{t}\frac{m_{\eta^{\prime}}^{2}}{m_{\eta^{\prime}}^{2}-t}-\frac{1}{t}\right),$
(79)

due to its coupling with the gluonic topological fluctuations
$m_{\eta^{\prime}}^{2}\propto\langle(F\tilde{F})^{2}\rangle$. On the right-
hand side, we have deliberately expressed the resulting propagator as the
difference of two poles at $t=0$. Now let us compare this with (77) which can
be identically rewritten in the form

$\displaystyle\frac{g_{P}(t)}{2M}$ $\displaystyle=$
$\displaystyle\frac{1}{t}\left(i\frac{\langle
P_{2}|\frac{n_{f}\alpha_{s}}{4\pi}F\tilde{F}|P_{1}\rangle}{\bar{u}(P_{2})\gamma_{5}u(P_{1})}-2Mg_{A}(t)\right)$
(80) $\displaystyle=$ $\displaystyle\frac{1}{t}\left(i\frac{\langle
P_{2}|\frac{n_{f}\alpha_{s}}{4\pi}F\tilde{F}|P_{1}\rangle}{\bar{u}(P_{2})\gamma_{5}u(P_{1})}-i\left.\frac{\langle
P_{2}|\frac{n_{f}\alpha_{s}}{4\pi}F\tilde{F}|P_{1}\rangle}{\bar{u}(P_{2})\gamma_{5}u(P_{1})}\right|_{t=0}\right)+2M\frac{g_{A}(0)-g_{A}(t)}{t}\,.$

We neglect the last term assuming $g_{A}(t)\approx g_{A}(0)=\Delta\Sigma$ to
be varying only slowly with $t$.555A partial justification of this comes from
the large-$N_{c}$ approximation where $m_{\eta^{\prime}}\sim{\cal
O}(1/\sqrt{N_{c}})$ is considered as small, at least parametrically, compared
to the singlet axial vector meson masses $m_{A}\sim{\cal O}(N_{c}^{0})$. Thus,
as long as one is interested in the region $|t|\sim m_{\eta^{\prime}}^{2}$,
the variation of $g_{A}(t)\sim 1/(t-m_{A}^{2})$ can be neglected. In practice,
however, the $\eta^{\prime}(957)$ is only slightly lighter than the
$f_{1}(1285)$.  The right-hand side can then be interpreted as a cancellation
of two poles at $t=0$, just like (79), between the ‘anomaly pole’ (first term)
and the naive pole (78) from the massless $\eta_{0}$ meson exchange (second
term). Eqs. (79) and (80) are actually identical in the single-pole
approximation where (80) is saturated by

$\displaystyle\frac{g_{P}(t)}{2M}\approx\frac{2M\Delta\Sigma}{m_{\eta^{\prime}}^{2}-t}\qquad
i\frac{\langle
P_{2}|\frac{n_{f}\alpha_{s}}{4\pi}F\tilde{F}|P_{1}\rangle}{\bar{u}(P_{2})\gamma_{5}u(P_{1})}\approx
2M\Delta\Sigma\frac{m_{\eta^{\prime}}^{2}}{m_{\eta^{\prime}}^{2}-t}.$ (81)

In the context of polarized DIS, the cancellation of poles just described had
been originally envisaged in Jaffe and Manohar (1989) and further elaborated
in Tarasov and Venugopalan (2022) to resolve issues with the $g_{1}$ structure
function. Compton scattering and GPDs offer a more general setup to explore
the physics of the anomaly to its full extent.

####  VI.1.3 Gravitational form factors

We now point out that one can repeat the same story for the QCD energy
momentum tensor $\Theta^{\alpha\beta}$ and its nucleon matrix element that
defines the gravitational form factors,

$\displaystyle\langle P_{2}|\Theta^{\alpha\beta}|P_{1}\rangle$
$\displaystyle=$
$\displaystyle\bar{u}(P_{2})\left[A(t)\frac{P^{\alpha}P^{\beta}}{M}+(A(t)+B(t))\frac{P^{(\alpha}i\sigma^{\beta)\lambda}l_{\lambda}}{2M}+D(t)\frac{l^{\alpha}l^{\beta}-g^{\alpha\beta}t}{4M}\right]u(P_{1}),$
(82)

where $a^{(\mu}b^{\nu)}=\frac{1}{2}(a^{\mu}b^{\nu}+a^{\nu}b^{\mu})$. Taking
the trace, we find an exact constraint among the form factors:

$\displaystyle\langle
P_{2}|(\Theta)^{\alpha}_{\alpha}|P_{1}\rangle=M\left(A(t)+\frac{B(t)}{4M^{2}}t-\frac{3D(t)}{4M^{2}}t\right)\bar{u}(P_{2})u(P_{1})=\langle
P_{2}|\frac{\beta(g)}{2g}F^{\mu\nu}F_{\mu\nu}|P_{1}\rangle.$ (83)

The right-hand side, on which $\beta(g)$ is the QCD beta function, is the
trace anomaly which signifies the explicit breaking of conformal symmetry. If
one naively neglects it, one finds a massless pole in $D(t)$ at $t=0$:

$\displaystyle\frac{3D(t)}{4}\approx\frac{M^{2}}{t}\qquad(t\to 0),$ (84)

where the conditions $A(0)=1$ and $B(0)=0$ have been used (so that one can
omit $tB(t)$ as $t\to 0$). By analogy with the massless $\eta_{0}$ pole in
(78), one might interpret the pole in (84) as due to the exchange of spin-0
glueballs which would couple to the operator $\Theta^{\alpha\beta}$ and which
would have been massless in the absence of the trace anomaly Bhattacharya et
al. (2023). In reality, however, the anomaly modifies (84) to

$\displaystyle\frac{3D(t)}{4}$ $\displaystyle\approx$
$\displaystyle\frac{M^{2}}{t}\left(A(t)-\frac{\langle
P_{2}|\frac{\beta(g)}{2g}F^{2}|P_{1}\rangle}{M\bar{u}(P_{2})u(P_{1})}\right)$
(85) $\displaystyle=$ $\displaystyle-\frac{M}{t}\left(\frac{\langle
P_{2}|\frac{\beta(g)}{2g}F^{2}|P_{1}\rangle}{\bar{u}(P_{2})u(P_{1})}-\left.\frac{\langle
P_{2}|\frac{\beta(g)}{2g}F^{2}|P_{1}\rangle}{\bar{u}(P_{2})u(P_{1})}\right|_{t=0}\right)+M^{2}\frac{A(t)-A(0)}{t}.$

Note the similarity to (80). The $D(t)$-form factor can be interpreted as the
difference of two poles at $t=0$, between the ‘anomaly pole’ (first term in
the brackets) and the naive glueball pole (84) (second term in the brackets).
As a result of this cancellation, the pole in $D(t)$ is shifted from $t=0$ to
physical glueball masses $t=m_{G}^{2}$ presumably in a way similar to (79).
However, unlike the situation in (80), in the present case the last term
$A(t)-A(0)$ of (85), which is related to spin-2 glueballs Bhattacharya et al.
(2023), is likely important at least from the large-$N_{c}$ perspective. Since
the trace anomaly cannot be turned off in the large-$N_{c}$ limit, the masses
of $2^{++}$ and $0^{++}$ glueballs are both ${\cal O}(N_{c}^{0})$. Moreover,
the analysis in Fujita et al. (2022) suggests that the single pole
approximation (cf., (81)) may not be a good approximation. The $D(t)$-form
factor thus exhibits ‘glueball dominance’

$\displaystyle
D(t)=\sum_{i}^{0^{++}}\frac{a_{i}}{m_{G_{i}}^{2}-t}+\sum_{j}^{2^{++}}\frac{b_{j}}{m_{G_{j}}^{2}-t},$
(86)

where the two contributions come from the $\langle F^{2}\rangle$ and $A(t)$
terms in (85), respectively. Incidentally, by taking the $t\to 0$ limit of
(85), one finds Cebulla et al. (2007)

$\displaystyle\frac{3D(0)}{4}=-M\left.\frac{d}{dt}\frac{\langle
P_{2}|\frac{\beta(g)}{2g}F^{2}|P_{1}\rangle}{\bar{u}(P_{2})u(P_{1})}\right|_{t=0}+\left.M^{2}\frac{dA(t)}{dt}\right|_{t=0}.$
(87)

The slope of a form factor at $t=0$ defines a ‘radius’ of the hadron. Eq. (87)
shows that the D-term $D(t=0)$ is related to the difference between two radii,
one defined by the scalar form factor $\langle F^{2}\rangle$ (related to the
$0^{++}$ glueball masses) and the other by the $A$-form factor (related to the
$2^{++}$ glueball masses), see recent discussions in Ji (2021); Mamo and Zahed
(2021); Kharzeev (2021); Fujita et al. (2022).

As we have seen in the above three examples, the existence or not of a
massless pole in form factors teaches us fundamental insights into the
nonperturbative dynamics of QCD. However, despite the known connections
between form factors and GPDs, the corresponding discussion at the GPD level
has been limited to the isovector sector (72) in the literature (see, however,
Bass (2002)). Our main purpose is to extend this argument to the singlet
sector.

###  VI.2 Anomaly poles in GPDs

Let us now return to our context. We have argued in Section V that the pole
$1/t$ in the one-loop Compton amplitude (5) should be absorbed into
$\tilde{E}_{q}$. This means that $\tilde{E}_{q}$ acquires a component related
to the twist-four GPD $\tilde{\cal F}$:

$\displaystyle\sum_{q}(\tilde{E}_{q}(x,\xi,t)+\tilde{E}_{q}(-x,\xi,t))=\frac{T_{R}n_{f}\alpha_{s}}{\pi}\frac{M^{2}}{t}\tilde{C}^{\rm
anom}\otimes\tilde{\cal F}(x,\xi,t)+\cdots,$ (88)

where $\tilde{C}^{\rm anom}$ is defined in (32). Integrating over $x$, we
exactly reproduce the first term of (80). Moreover, (80) suggests that there
is another, ‘primordial’ pole in $\tilde{E}_{q}$ which exactly cancels the
$1/t$ pole to make $\tilde{E}_{q}$ finite for all values of $x$ and $\xi$ in
the limit $t\to 0$. A simple, yet ad-hoc fix consistent with (80) is to add a
‘counterterm’

$\displaystyle\sum_{q}(\tilde{E}_{q}(x,\xi,t)+\tilde{E}_{q}(-x,\xi,t))\approx\frac{T_{R}n_{f}\alpha_{s}}{\pi}\frac{M^{2}}{t}\tilde{C}^{\rm
anom}\otimes(\tilde{\cal F}(x,\xi,t)-\tilde{\cal F}(x,\xi,0)).$ (89)

This may be viewed as the non-local version of the local relation (80). The
second, added term is an analog of (72), but interestingly, in the present
case the pole is not limited to the ERBL region $x<\xi$. We postulate (89) as
a nonperturbative relation between the twist-two and twist-four GPDs mediated
by the chiral anomaly.

The fate of the $1/t$ pole in the unpolarized sector and its connection to the
trace anomaly are more involved. This is partly because the QCD energy
momentum tensor consists of a quark and a gluon part,
$\Theta^{\alpha\beta}=\sum_{q}\Theta_{q}^{\alpha\beta}+\Theta_{g}^{\alpha\beta}$,
in contrast to $J_{5}^{\alpha}$ which is purely a quark operator. Accordingly,
one can define gravitational form factors separately for quarks and gluons Ji
(1997):

$\displaystyle\langle P_{2}|\Theta_{q,g}^{\alpha\beta}|P_{1}\rangle$
$\displaystyle=$
$\displaystyle\bar{u}(P_{2})\left[A_{q,g}(t)\frac{P^{\alpha}P^{\beta}}{M}+(A_{q,g}(t)+B_{q,g}(t))\frac{P^{(\alpha}i\sigma^{\beta)\lambda}l_{\lambda}}{2M}+D_{q,g}(t)\frac{l^{\alpha}l^{\beta}-g^{\alpha\beta}t}{4M}+\bar{C}_{q,g}(t)Mg^{\alpha\beta}\right]u(P_{1}).$

They are related to the second moments of the unpolarized quark GPDs,

$\displaystyle\int_{-1}^{1}dxxH_{q}(x,\xi,t)$
$\displaystyle=\int_{0}^{1}dxx(H_{q}(x,\xi,t)-H_{q}(-x,\xi,t))=A_{q}(t)+\xi^{2}D_{q}(t)\,,$
(91) $\displaystyle\int_{-1}^{1}dxxE_{q}(x,\xi,t)$
$\displaystyle=\int_{0}^{1}dxx(E_{q}(x,\xi,t)-E_{q}(-x,\xi,t))=B_{q}(t)-\xi^{2}D_{q}(t)\,,$
(92)

and similarly for the gluon GPDs. Taking the trace of (82), we find

$\displaystyle\langle
P_{2}|\sum_{q}(\Theta_{q})^{\alpha}_{\alpha}|P_{1}\rangle$ $\displaystyle=$
$\displaystyle\sum_{q}M\left(A_{q}(t)+4\bar{C}_{q}(t)+\frac{B_{q}(t)}{4M^{2}}t-\frac{3D_{q}(t)}{4M^{2}}t\right)\bar{u}(P_{2})u(P_{1})$
(93) $\displaystyle=$ $\displaystyle\langle
P_{2}|\frac{\beta_{q}(g)}{2g}F^{2}|P_{1}\rangle\approx\langle
P_{2}|\frac{T_{R}n_{f}\alpha_{s}}{6\pi}F^{2}|P_{1}\rangle,$

where $\frac{\beta_{q}}{2g}$ is the quark part of the trace anomaly that can
be systematically calculated in perturbation theory Hatta et al. (2018);
Tanaka (2019); Ahmed et al. (2023); Tanaka (2023). To lowest order, it is
simply the $n_{f}$ term of the beta function:

$\displaystyle(\Theta_{q})^{\alpha}_{\alpha}+(\Theta_{g})^{\alpha}_{\alpha}=\frac{\beta(g)}{2g}F^{\mu\nu}F_{\mu\nu}=-\frac{\alpha_{s}}{8\pi}\left(\frac{11N_{c}}{3}-\frac{4T_{R}n_{f}}{3}\right)F^{2}+\cdots.$
(94)

Clearly, (93) is not as constraining as (83) because of the new form factors
$B_{q}(t),\bar{C}_{q}(t)$. (Note that $B_{q}(0),\bar{C}_{q}(0)\neq 0$ although
$B_{q}(0)+B_{g}(0)=\bar{C}_{q}(t)+\bar{C}_{g}(t)=0$.) Nevertheless we may try
to rewrite it in a way similar to (85)

$\displaystyle\sum_{q}\frac{3D_{q}(t)-B_{q}(t)}{4}=-\frac{M}{t}\left(\frac{\langle
P_{2}|\frac{\beta_{q}}{2g}F^{2}|P_{1}\rangle}{\bar{u}(P_{2})u(P_{1})}-\left.\frac{\langle
P_{2}|\frac{\beta_{q}}{2g}F^{2}|P_{1}\rangle}{\bar{u}(P_{2})u(P_{1})}\right|_{t=0}\right)+\frac{M^{2}}{t}\sum_{q}\Bigl{(}A_{q}(t)+4\bar{C}_{q}(t)-A_{q}(0)-4\bar{C}_{q}(0)\Bigr{)}.$
(95)

Let us now discuss how the constraint (95) from the trace anomaly is encoded
in the GPDs. We have argued that the anomaly poles in (4) should be absorbed
into the unpolarized GPDs,

$\displaystyle\sum_{q}(H_{q}(x,\xi,t)-H_{q}(-x,\xi,t))$
$\displaystyle=\frac{T_{R}n_{f}\alpha_{s}}{\pi}\frac{M^{2}}{t}C^{\rm
anom}\otimes{\cal F}(x,\xi,t)+\cdots,$ (96)
$\displaystyle\sum_{q}(E_{q}(x,\xi,t)-E_{q}(-x,\xi,t))$
$\displaystyle=-\frac{T_{R}n_{f}\alpha_{s}}{\pi}\frac{M^{2}}{t}C^{\rm
anom}\otimes{\cal F}(x,\xi,t)+\cdots,$ (97)

where $C^{\rm anom}$ is defined in (31). Taking the second moment and
comparing with (91), we find

$\displaystyle\left.\sum_{q}A_{q}(t)\right|_{\rm
pole}=-\left.\sum_{q}B_{q}(t)\right|_{\rm pole}=-\frac{M}{t}\frac{\langle
P_{2}|\frac{T_{R}n_{f}\alpha_{s}}{6\pi}F^{\mu\nu}(i\overleftrightarrow{D}^{+})^{2}F_{\mu\nu}|P_{1}\rangle}{(P^{+})^{2}\bar{u}(P_{2})u(P_{1})},$
(98) $\displaystyle\left.\sum_{q}D_{q}(t)\right|_{\rm
pole}=-\frac{M}{t}\frac{\langle
P_{2}|\frac{T_{R}n_{f}\alpha_{s}}{6\pi}F^{2}|P_{1}\rangle}{\bar{u}(P_{2})u(P_{1})}.$
(99)

Eq. (99) seems to reproduce the first term of (95) after $\beta_{q}$ is
expanded to lowest order. However, apparently there is a factor $\frac{3}{4}$
mismatch in the normalization. Besides, the previous argument around (85) did
not hint at the possible existence of an anomaly pole in the $A_{q},B_{q}$
form factors.

In order to understand these differences, we quote the one-loop result for the
energy momentum tensor matrix element between on-shell gluon (not nucleon)
states Giannotti and Mottola (2009):

$\displaystyle\langle
p_{2}|\Theta_{q}^{\alpha\beta}|p_{1}\rangle=-\frac{T_{R}\alpha_{s}}{6\pi}\left(\frac{p^{\alpha}p^{\beta}}{t}+\frac{l^{\alpha}l^{\beta}-tg^{\alpha\beta}}{4t}\right)\langle
p_{2}|F^{\mu\nu}F_{\mu\nu}|p_{1}\rangle+\cdots\,.$ (100)

A superficial comparison with (VI.2) suggests that poles of equal magnitude
are induced in the $A_{q},B_{q},D_{q}$ form factors

$\displaystyle A_{q}(t)\approx-B_{q}(t)\approx
D_{q}(t)\sim\frac{\langle\alpha_{s}F^{2}\rangle}{t},$ (101)

and the issue of the factor $\frac{3}{4}$ goes away because
$\frac{3D_{q}(t)-B_{q}(t)}{4}\approx D_{q}(t)$ on the left hand side of (95).
Taking the trace of (100), we find

$\displaystyle\langle
p_{2}|(\Theta_{q})^{\alpha}_{\alpha}|p_{1}\rangle=\langle
p_{2}|\frac{T_{R}\alpha_{s}}{6\pi}F^{2}|p_{1}\rangle,$ (102)

which is the correct trace anomaly relation to this order. To obtain (102), it
is important to use the on-shell condition $p^{2}=-t/4$ of the external
states, so that the two terms in (100) contribute $\frac{1}{4}$ and
$\frac{3}{4}$ of the total anomaly, respectively. Going from gluon to nucleon
targets, we see that the way the trace anomaly relation (93) is fulfilled
among various form factors is highly nontrivial. A different, spin-2 operator
$F(D^{+})^{2}F$ is involved in the $A_{q},B_{q}$ form factor (98) due to the
convolution integral in $x$. Moreover, a naive identification
$p^{\alpha}p^{\beta}\to P^{\alpha}P^{\beta}$ is precarious because the nucleon
is massive $P^{2}=M^{2}-t/4$. While the difference is negligible when
$\sqrt{|t|}\gg M$, this obscures the fate of the poles in (98) as $t$ gets
smaller.

On the other hand, the tensor $l^{\alpha}l^{\beta}$ is formally identical for
both the nucleon and gluon targets,
$l^{\alpha}=p_{2}^{\alpha}-p_{1}^{\alpha}=P_{2}^{\alpha}-P_{1}^{\alpha}$. We
may therefore expect that the anomaly relation at the partonic level is better
reflected in the $D_{q}(t)$ form factor even at the hadronic level, just like
the $g_{P}(t)$ form factor which is the coefficient of $l^{\alpha}$. Indeed,
the opposite signs in (97) suggests that the pole terms mainly feed into the
so-called Polyakov-Weiss D-term Polyakov and Weiss (1999) of the unpolarized
GPDs,

$\displaystyle H^{\rm PW}_{q}(x,\xi,t)=-E^{\rm
PW}_{q}(x,\xi,t)=\theta(\xi-|x|)D_{q}(x/\xi,t).$ (103)

The distribution $D_{q}(z,t)$ is odd in $z$ and is solely responsible for the
highest power of $\xi$ in the moments of GPDs. In order to extract it, we take
the $n$-th moment of (97) with odd integers $n$,

$\displaystyle\sum_{q}\int_{-1}^{1}dxx^{n}H_{q}(x,\xi,t)$
$\displaystyle\approx$
$\displaystyle\frac{T_{R}n_{f}\alpha_{s}}{\pi}\frac{M^{2}}{t}\int_{0}^{1}dx\frac{x^{n}}{(n+2)(n+3)}\frac{1-\left(\frac{\xi}{x}\right)^{n+3}}{1-\frac{\xi^{2}}{x^{2}}}\left({\cal
F}(x,\xi,t)-{\cal F}(x,\xi,0)\right)$ (104) $\displaystyle\equiv$
$\displaystyle\sum_{i=0}^{n+1}\sum_{q}h_{qn}^{i}\xi^{i},$

where we have minimally subtracted the pole at $t=0$ as in (89). The highest
power $h_{qn}^{n+1}$ is related to $D_{q}(z,t)$ as

$\displaystyle\sum_{q}\int_{-1}^{1}dzz^{n}D_{q}(z,t)$ $\displaystyle\approx$
$\displaystyle\sum_{q}h^{n+1}_{qn}(t)$ (105) $\displaystyle=$
$\displaystyle\frac{T_{R}n_{f}\alpha_{s}}{\pi}\frac{M^{2}}{t}\frac{1}{(n+2)(n+3)}\int_{0}^{1}\frac{dx}{x}\left({\cal
F}(x,\xi,t)-{\cal F}(x,\xi,0)\right)$ $\displaystyle=$
$\displaystyle-2\frac{T_{R}n_{f}\alpha_{s}}{\pi}\frac{M}{t}\frac{1}{(n+2)(n+3)}\left(\frac{\langle
P_{2}|F^{2}|P_{1}\rangle}{\bar{u}(P_{2})u(P_{1})}-\left.\frac{\langle
P_{2}|F^{2}|P_{1}\rangle}{\bar{u}(P_{2})u(P_{1})}\right|_{t=0}\right).$

By definition, the $n=1$ moment is the gravitational form factor
$\int_{-1}^{1}dzzD_{q}(z,t)=D_{q}(t)$. Inverting the Mellin transform (105)
and noting that $D_{q}(z,t)$ is an odd function of $z$, we obtain

$\displaystyle\sum_{q}D_{q}(z,t)\approx-\frac{T_{R}n_{f}\alpha_{s}}{\pi}z(1-|z|)\frac{M}{t}\left(\frac{\langle
P_{2}|F^{2}|P_{1}\rangle}{\bar{u}(P_{2})u(P_{1})}-\left.\frac{\langle
P_{2}|F^{2}|P_{1}\rangle}{\bar{u}(P_{2})u(P_{1})}\right|_{t=0}\right),$ (106)

and in particular,

$\displaystyle\sum_{q}D_{q}(t)\approx-\frac{M}{t}\left(\frac{\langle
P_{2}|\frac{T_{R}n_{f}\alpha_{s}}{6\pi}F^{2}|P_{1}\rangle}{\bar{u}(P_{2})u(P_{1})}-\left.\frac{\langle
P_{2}|\frac{T_{R}n_{f}\alpha_{s}}{6\pi}F^{2}|P_{1}\rangle}{\bar{u}(P_{2})u(P_{1})}\right|_{t=0}\right).$
(107)

Since $\langle P|F^{2}|P\rangle<0$ in QCD and the form factor $\langle
P_{2}|F^{2}|P_{1}\rangle$ is a decreasing function of $|t|$, the right-hand
side of (106) is positive, whereas $D_{q}(t)$ is usually believed to be
negative. While we expect that eventually the leading-order coefficient
$\frac{T_{R}n_{f}\alpha_{s}}{6\pi}$ will be replaced by $\frac{\beta_{q}}{2g}$
after including higher-order corrections, according to the three-loop analyses
in Tanaka (2019); Metz et al. (2020); Ahmed et al. (2023); Tanaka (2023), the
sign does not flip $\frac{\beta_{q}}{2g}>0$. This suggests that the other
terms in (95) that were neglected in the above minimal subtraction procedure
may be numerically important as we already suspected in the argument below
(85). Note also that the sign does flip if one includes the gluon contribution
to recover the full beta function of QCD
$\frac{\beta_{q}}{2g}\to\frac{\beta}{2g}<0$.

##  VII Conclusions

In this work, we have performed a complete one-loop calculation of the Compton
scattering amplitude using momentum transfer $t$ as the regulator of the
collinear singularity. Our approach differs from all the previous calculations
in the GPD literature where one typically uses dimensional regularization to
isolate the collinear singularity and sets $t=0$ right from the start,
assuming that nonzero $t$ only generates higher-twist corrections of order
$t/Q^{2}$. In practice, the introduction of an additional variable $t$ makes
the calculation more cumbersome and brings in unusual features. In the gluon
initiated channel, we have found anomaly poles $1/t$ (29), (30) accompanied by
twist-four GPDs (9), (10) in both the real and imaginary parts of the Compton
amplitude, confirming and extending our previous finding Bhattacharya et al.
(2023). In the quark initiated channel, we have unexpectedly found uncancelled
single and double IR poles in the ‘coefficient functions’ (25), (26). Each of
these features potentially implies the violation of factorization. However, we
have also performed the one-loop calculation of GPDs for quark and gluon
states with the same set of regulators and showed how all these poles can be
systematically absorbed into the GPDs themselves. This shows that QCD
factorization is restored at least to this order.

This is however not the end of the story. We have also explored connections
between GPDs and anomalies, as a natural and necessary consequence of the
known connections between form factors and anomalies. We have argued that once
the poles $1/t$ have been absorbed into GPDs, they become a part of the GPDs.
In other words, anomalies nonperturbatively relate twist-two and twist-four
GPDs. Such relations, once integrated over $x$, are expected to reproduce the
constraints among the corresponding form factors. This scenario seems to be
working for the polarized GPD $\tilde{E}_{q}$ and its connection to the chiral
anomaly. Relation (89), partly supported by the large-$N_{c}$ argument, can be
viewed as the $x$-dependent generalization of the form factor relation (77).
The situation is more complicated (and more interesting) for the unpolarized
GPDs $H_{q},E_{q}$ and their relation to the trace anomaly. We have argued
that the anomaly mostly constrains the $D_{q}(t)$ form factor and its GPD
analog, the Polyakov-Weiss D-term. The results we have arrived at (106) (107)
are roughly consistent with the anomaly relation (95), but they differ in
detail. Further investigation in this direction is necessary.

In conclusion, we have proposed finite-$t$ regularization as an alternative
factorization scheme that elucidates the physics of anomalies. This is a
scheme where we are able to unravel novel connections between twist-two and
twist-four GPDs mediated by the anomalies of QCD. Admittedly, the calculation
is more cumbersome than the standard dimensional regularization with $t=0$.
Still, the chiral and trace anomalies are among the most fascinating phenomena
of QCD with far-reaching consequences, and we believe that research on GPDs is
enriched by incorporating such fundamental problems. There are a number of
directions along which the current work can be refined or extended, in
addition to the aforementioned tension between (95) and (107). First, we
strongly suspect that anomaly poles are present in higher order perturbation
theory. Especially in the symmetric case, we expect that each additional loop
provides the corresponding term in the expansion of the (quark part of the)
beta function. A related question is whether there are anomaly poles in the
gluon GPDs that complement the quark ones to restore the full beta function
$\beta=\beta_{q}+\beta_{g}$ Hatta et al. (2018). Another important question
which has not been addressed at all in this paper is how to understand the new
relations from a renormalization group point of view. In the present scheme,
the mixing between the twist-two and twist-four GPDs occurs as a result of a
finite subtraction rather than the DGLAP evolution of GPDs. The UV properties
of the twist-four GPDs (9) and (10) have been studied in Hatta and Zhao
(2020); Hatta (2020); Radyushkin and Zhao (2021, 2022), but more work is
certainly needed. Furthermore, it is well known that at twist-3 accuracy, the
amplitude for DVCS off the nucleon contains twist-3 GPDs apart from the usual
twist-2 GPDs. It is also interesting to pursue whether or not there are
imprints of anomalies on twist-3 GPDs and related observables. Finally,
constraints from anomalies should be implemented in the modeling of GPDs. In
particular, the specific functional form given in (106) might be helpful to
model this poorly constrained distribution.

###### Acknowledgements.

We are very grateful to Vladimir Braun and Anatoly Radyushkin for many useful
discussions. We also thank Kornelija Passek-Kumeri$\check{\rm c}$ki, Swagato
Mukherjee, Kazuhiro Tanaka, Raju Venugopalan and Christian Weiss for
discussions. S. B. and Y. H. are supported by the U.S. Department of Energy
under Contract No. DE-SC0012704, and Laboratory Directed Research and
Development (LDRD) funds from Brookhaven Science Associates. Y. H. is also
supported by the framework of the Saturated Glue (SURGE) Topical Theory
Collaboration. W. V. has been supported by Deutsche Forschungsgemeinschaft
(DFG) through the Research Unit FOR 2926 (project 409651613).

##  Appendix A Derivation of Eq. (37)

In this appendix we give an outline of the derivation of (37) and (38). The
other results in Section V can be derived similarly. For the quark matrix
elements, we work in the Feynman gauge. The ladder diagram reads, up to a
prefactor,

$\displaystyle\mu^{2\epsilon}\int\frac{dk^{-}d^{2-2\epsilon}k_{\perp}}{(2\pi)^{3-2\epsilon}}\bar{u}(p+l/2)\frac{\gamma_{\mu}({\ooalign{\hfil/\hfil\crcr$k$}}+{\ooalign{\hfil/\hfil\crcr$l$}}/2)\gamma^{+}({\ooalign{\hfil/\hfil\crcr$k$}}-{\ooalign{\hfil/\hfil\crcr$l$}}/2)\gamma^{\mu}}{(p-k)^{2}(k-l/2)^{2}(k+l/2)^{2}}u(p-l/2),$
(116)

where

$\displaystyle k^{+}=xp^{+},\quad l^{+}=-2\xi p^{+},\quad l^{-}=-\frac{\xi
l^{2}}{4p^{+}},\quad\vec{l}_{\perp}^{2}=(\xi^{2}-1)l^{2},\quad
p^{-}=-\frac{l^{2}}{8p^{+}}.$ (117)

In the DGLAP region $\xi<x<1$, the $k^{-}$ integral can be done by picking up
the pole of $(p-k)^{2}+i\epsilon=0$ at

$\displaystyle k^{-}=-\frac{k_{\perp}^{2}+\frac{1-x}{4}l^{2}}{2(1-x)p^{+}},$
(118)

in the upper half plane. The remaining propagators can be combined as

$\displaystyle\frac{1}{1-x}\int_{0}^{1}da\frac{A(k_{\perp})}{(k^{2}+(1-2a)k\cdot
l+\frac{l^{2}}{4})^{2}}=\frac{1-x}{(1-\xi(1-2a))^{2}}\int_{0}^{1}da\frac{A(k^{\prime}_{\perp}-\frac{(1-2a)(1-x)}{2(1-\xi(1-2a))}l_{\perp})}{\left(k^{\prime
2}_{\perp}-\frac{(1-a)a(1-x)^{2}l^{2}}{(1-\xi(1-2a))^{2}}\right)^{2}},$ (119)

where in the denominator we have shifted momentum $k_{\perp}\to
k^{\prime}_{\perp}$ to complete the square. In the numerator we have projected
onto the twist-two component,

$\displaystyle\bar{u}(p+l/2)\cdots u(p-l/2)\to\left[(1-\epsilon)k^{\prime
2}_{\perp}+(B-\epsilon C)l^{2}\right]\bar{u}\gamma^{+}u\equiv
A\bar{u}\gamma^{+}u,$ (120)

with

$\displaystyle
B=\frac{(1-\xi+a(x+2\xi-1))(x-\xi-a(x-2\xi-1))}{(1-\xi(1-2a))^{2}},\quad
C=\frac{(1-a)a(1-x)^{2}}{(1-\xi(1-2a))^{2}}.$ (121)

The terms linear in $k^{\prime}_{\perp}$ have been dropped in (120) since they
vanish after the $k^{\prime}_{\perp}$ integral:

$\displaystyle\mu^{2\epsilon}\int_{0}^{1}da\frac{1-x}{(1-\xi(1-2a))^{2}}\int\frac{d^{2-2\epsilon}k^{\prime}_{\perp}}{(2\pi)^{2-2\epsilon}}\frac{(1-\epsilon)k^{\prime
2}_{\perp}+(B-\epsilon C)l^{2}}{\left(k^{\prime 2}_{\perp}-Cl^{2}\right)^{2}}$
$\displaystyle\approx\left(\frac{\mu^{2}}{-l^{2}}\right)^{\epsilon}\int_{0}^{1}da\frac{1-x}{(1-\xi(1-2a))^{2}}\frac{1}{(4\pi)^{1-\epsilon}}\left(\frac{(1-2\epsilon)\Gamma(\epsilon)}{C^{\epsilon}}-\frac{B\Gamma(1+\epsilon)}{C^{1+\epsilon}}\right).$
(122)

The first integral gives a UV pole:

$\displaystyle\frac{\Gamma(\epsilon_{\text{\tiny{UV}}})}{4\pi}\left(\frac{4\pi\mu^{2}}{-l^{2}}\right)^{\epsilon}\int_{0}^{1}da\frac{(1-x)(1-2\epsilon)}{(1-\xi(1-2a))^{2}C^{\epsilon}}$
$\displaystyle=\frac{1}{4\pi}\left(\frac{\tilde{\mu}^{2}}{-l^{2}}\right)^{\epsilon}\left(\frac{1}{\epsilon_{\text{\tiny{UV}}}}\frac{1-x}{1-\xi^{2}}-\frac{(1-x)(2\ln(1-x)-\ln(1-\xi^{2}))}{1-\xi^{2}}\right),$
(123)

while the second integral gives double and single IR poles:

$\displaystyle\left(\frac{4\pi\mu^{2}}{-l^{2}}\right)^{\epsilon}\int_{0}^{1}da\frac{1-x}{(1-\xi(1-2a))^{2}}\frac{-B\Gamma(1+\epsilon)}{C^{1+\epsilon}}$
(124) $\displaystyle=$
$\displaystyle\left(\frac{\tilde{\mu}^{2}}{-l^{2}}\right)^{\epsilon}\frac{1}{(1-x)^{1+2\epsilon}}\left(\frac{2(x-\xi^{2})}{1-\xi^{2}}\frac{1}{\epsilon_{\text{\tiny{IR}}}}-\frac{(1-x)^{2}-2(x-\xi^{2})\ln(1-\xi^{2})}{1-\xi^{2}}+\epsilon
f(x)\right)$ $\displaystyle=$
$\displaystyle\left(\frac{\tilde{\mu}^{2}}{-l^{2}}\right)^{\epsilon}\Biggl{[}-\delta(1-x)\left(\frac{1}{\epsilon_{\text{\tiny{IR}}}^{2}}+\frac{\ln(1-\xi^{2})}{\epsilon_{\text{\tiny{IR}}}}\right)+\frac{2(x-\xi^{2})}{(1-x)_{+}(1-\xi^{2})}\frac{1}{\epsilon_{\rm
IR}}-\frac{4(x-\xi^{2})}{1-\xi^{2}}\left(\frac{\ln(1-x)}{1-x}\right)_{+}$
$\displaystyle-\frac{(1-x)^{2}-2(x-\xi^{2})\ln(1-\xi^{2})}{(1-\xi^{2})(1-x)_{+}}-\frac{1}{2}\delta(1-x)\left(\ln^{2}(1-\xi^{2})-\frac{\pi^{2}}{6}\right)\Biggr{]}.$

Here $f(x)$ is a certain function whose value at $x=1$ is the only thing we
need.

In Feynman gauge, there are two other diagrams, giving

$\displaystyle\int\frac{d^{d}k}{(2\pi)^{d-1}}\bar{u}(p+l/2)\frac{\gamma^{+}{\ooalign{\hfil/\hfil\crcr$k$}}\gamma^{+}\bigl{(}\delta(k^{+}-(x+\xi)p^{+})-\delta((1-x)p^{+})\bigr{)}}{k^{2}(p-l/2-k)^{2}((1+\xi)p^{+}-k^{+})}u(p-l/2),$
(127)

and the corresponding contribution for its mirror diagram. They can be
similarly evaluated. The result is

$\displaystyle\left(\frac{2(x-\xi^{2})}{(1-\xi^{2})(1-x)_{+}}+\delta(1-x)\left(2-\ln(1-\xi^{2})\right)\right)\frac{\left(\frac{\tilde{\mu}^{2}}{-l^{2}}\right)^{\epsilon}}{4\pi}\left(\frac{1}{\epsilon_{\text{\tiny{UV}}}}-\frac{1}{\epsilon_{\text{\tiny{IR}}}}\right).$
(128)

Adding also the quark self energy diagrams on the external legs, we arrive at
(37) and (38). We note that the calculation can also be performed in Landau
gauge which has the advantage that the self-energy diagrams vanish
identically.

In the ERBL region $\xi>x$, the pole of $(k+l/2)^{2}$ in (116) moves to the
upper half plane. We thus pick up the pole of $(k-l/2)^{2}$ at

$\displaystyle
k^{-}_{c}=\frac{k_{\perp}^{2}-(x\xi+1)\frac{l^{2}}{4}-\vec{k}_{\perp}\cdot\vec{l}_{\perp}}{2(x+\xi)p^{+}},$
(129)

in the lower half plane. The first term in (127) also contributes (but not its
mirror diagram).

## References

  * Braun et al. (2017) V. M. Braun, A. N. Manashov, S. Moch, and M. Strohmaier, JHEP 06, 037 (2017), eprint 1703.09532.
  * Braun et al. (2020) V. M. Braun, A. N. Manashov, S. Moch, and J. Schoenleber, JHEP 09, 117 (2020), [Erratum: JHEP 02, 115 (2022)], eprint 2007.06348.
  * Braun et al. (2022) V. M. Braun, Y. Ji, and J. Schoenleber, Phys. Rev. Lett. 129, 172001 (2022), eprint 2207.06818.
  * Schoenleber (2023) J. Schoenleber, JHEP 02, 207 (2023), eprint 2209.09015.
  * Abdul Khalek et al. (2022) R. Abdul Khalek et al., Nucl. Phys. A 1026, 122447 (2022), eprint 2103.05419.
  * Bhattacharya et al. (2023) S. Bhattacharya, Y. Hatta, and W. Vogelsang, Phys. Rev. D 107, 014026 (2023), eprint 2210.13419.
  * Tarasov and Venugopalan (2020) A. Tarasov and R. Venugopalan, Phys. Rev. D 102, 114022 (2020), eprint 2008.08104.
  * Tarasov and Venugopalan (2022) A. Tarasov and R. Venugopalan, Phys. Rev. D 105, 014020 (2022), eprint 2109.10370.
  * Ji and Osborne (1998) X.-D. Ji and J. Osborne, Phys. Rev. D 57, 1337 (1998), eprint hep-ph/9707254.
  * Belitsky and Mueller (1998) A. V. Belitsky and D. Mueller, Phys. Lett. B 417, 129 (1998), eprint hep-ph/9709379.
  * Mankiewicz et al. (1998) L. Mankiewicz, G. Piller, E. Stein, M. Vanttinen, and T. Weigl, Phys. Lett. B 425, 186 (1998), [Erratum: Phys.Lett.B 461, 423–423 (1999)], eprint hep-ph/9712251.
  * Pire et al. (2011) B. Pire, L. Szymanowski, and J. Wagner, Phys. Rev. D 83, 034009 (2011), eprint 1101.0555.
  * Bertone et al. (2022) V. Bertone, H. Dutrieux, C. Mezrag, J. M. Morgado, and H. Moutarde, Eur. Phys. J. C 82, 888 (2022), eprint 2206.01412.
  * Jaffe and Manohar (1989) R. L. Jaffe and A. Manohar, Phys. Lett. B 223, 218 (1989).
  * Collins and Freund (1999) J. C. Collins and A. Freund, Phys. Rev. D 59, 074009 (1999), eprint hep-ph/9801262.
  * Diehl (2003) M. Diehl, Phys. Rept. 388, 41 (2003), eprint hep-ph/0307382.
  * Belitsky and Radyushkin (2005) A. V. Belitsky and A. V. Radyushkin, Phys. Rept. 418, 1 (2005), eprint hep-ph/0504030.
  * Radyushkin (1998) A. V. Radyushkin, Phys. Rev. D 58, 114008 (1998), eprint hep-ph/9803316.
  * Hatta and Zhao (2020) Y. Hatta and Y. Zhao, Phys. Rev. D 102, 034004 (2020), eprint 2006.02798.
  * Hatta (2020) Y. Hatta, Phys. Rev. D 102, 094004 (2020), eprint 2009.03657.
  * Radyushkin and Zhao (2021) A. Radyushkin and S. Zhao, JHEP 12, 010 (2021), eprint 2111.00887.
  * Radyushkin and Zhao (2022) A. Radyushkin and S. Zhao, JHEP 02, 163 (2022), eprint 2201.02181.
  * Patel (2015) H. H. Patel, Comput. Phys. Commun. 197, 276 (2015), eprint 1503.01469.
  * ’t Hooft and Veltman (1972) G. ’t Hooft and M. J. G. Veltman, Nucl. Phys. B 44, 189 (1972).
  * Breitenlohner and Maison (1977) P. Breitenlohner and D. Maison, Commun. Math. Phys. 52, 11 (1977).
  * Jamin and Lautenbacher (1993) M. Jamin and M. E. Lautenbacher, Comput. Phys. Commun. 74, 265 (1993).
  * Weber (1992) A. Weber, Nucl. Phys. B 382, 63 (1992).
  * Vogelsang (1996) W. Vogelsang, Nucl. Phys. B 475, 47 (1996), eprint hep-ph/9603366.
  * White (2002) B. E. White, J. Phys. G 28, 203 (2002), eprint hep-ph/0102121.
  * Altarelli et al. (1979) G. Altarelli, R. K. Ellis, and G. Martinelli, Nucl. Phys. B 157, 461 (1979).
  * Stratmann et al. (1996) M. Stratmann, A. Weber, and W. Vogelsang, Phys. Rev. D 53, 138 (1996), eprint hep-ph/9509236.
  * Bodwin and Qiu (1990) G. T. Bodwin and J.-W. Qiu, Phys. Rev. D 41, 2755 (1990).
  * Vogelsang (1991) W. Vogelsang, Z. Phys. C 50, 275 (1991).
  * Penttinen et al. (2000) M. Penttinen, M. V. Polyakov, and K. Goeke, Phys. Rev. D 62, 014024 (2000), eprint hep-ph/9909489.
  * Witten (1979) E. Witten, Nucl. Phys. B 156, 269 (1979).
  * Veneziano (1979) G. Veneziano, Nucl. Phys. B 159, 213 (1979).
  * Fujita et al. (2022) M. Fujita, Y. Hatta, S. Sugimoto, and T. Ueda, PTEP 2022, 093 (2022), eprint 2206.06578.
  * Cebulla et al. (2007) C. Cebulla, K. Goeke, J. Ossmann, and P. Schweitzer, Nucl. Phys. A 794, 87 (2007), eprint hep-ph/0703025.
  * Ji (2021) X. Ji, Front. Phys. (Beijing) 16, 64601 (2021), eprint 2102.07830.
  * Mamo and Zahed (2021) K. A. Mamo and I. Zahed, Phys. Rev. D 103, 094010 (2021), eprint 2103.03186.
  * Kharzeev (2021) D. E. Kharzeev, Phys. Rev. D 104, 054015 (2021), eprint 2102.00110.
  * Bass (2002) S. D. Bass, Phys. Rev. D 65, 074025 (2002), eprint hep-ph/0102036.
  * Ji (1997) X.-D. Ji, Phys. Rev. Lett. 78, 610 (1997), eprint hep-ph/9603249.
  * Hatta et al. (2018) Y. Hatta, A. Rajan, and K. Tanaka, JHEP 12, 008 (2018), eprint 1810.05116.
  * Tanaka (2019) K. Tanaka, JHEP 01, 120 (2019), eprint 1811.07879.
  * Ahmed et al. (2023) T. Ahmed, L. Chen, and M. Czakon, JHEP 01, 077 (2023), eprint 2208.01441.
  * Tanaka (2023) K. Tanaka, JHEP 03, 013 (2023), eprint 2212.09417.
  * Giannotti and Mottola (2009) M. Giannotti and E. Mottola, Phys. Rev. D 79, 045014 (2009), eprint 0812.0351.
  * Polyakov and Weiss (1999) M. V. Polyakov and C. Weiss, Phys. Rev. D 60, 114017 (1999), eprint hep-ph/9902451.
  * Metz et al. (2020) A. Metz, B. Pasquini, and S. Rodini, Phys. Rev. D 102, 114042 (2020), eprint 2006.11171.