University of Warsaw
Faculty of Physics
Oskar Grocholski
Record book number: 382618

Factorization in hard exclusive
processes.
Computation of one-loop corrections to the diphoton photoproduction on proton.
Master’s thesis
in the field of Physics

The thesis was written under the supervision of
Dr hab. Jakub Wagner
National Centre for Nuclear Research
and Prof. Piotr H. Chankowski
University of Warsaw, Faculty of Physics
Institute of Theoretical Physics
Chair of Theory of Hadrons and Leptons

Warsaw, July 2021

Summary

Generalized Parton Distributions (GPDs) carry information on the internal
structure of hadrons such as the angular momentum of quarks and gluons, or
their spacelike distribution. They can be experimentally studied in exclusive
experiments with hadrons, i.e. processes in which all initial and final states
are measured. The tool that creates the necessary bridge between the
theoretical predictions and the experiments is the collinear factorization. It
allows disentangling perturbatively computable parts of the amplitude, which
describes interactions of quark and gluons with the external particles, from
the non-perturbative quantities, which are identified as GPDs.

In this work, I extend the theoretical analysis of the process photoproduction
of photon pairs on a proton to the next-to-leading order in perturbative
Quantum Chromodynamics within the framework of collinear factorization. I give
the proof that all collinear divergences which arise in one-loop computations
cancel at the level of the amplitude. This result enlarges the family of
reactions, which can be studied using the collinear factorization by processes
of the type $2\rightarrow 3$, which have not been previously studied within
this theoretical framework beyond the leading QCD order. Furthermore, I
compute the full form of the amplitude of the discussed process at the one-
loop order. That improves the accuracy of theoretical predictions for this
experiment, which may be used for planning future experiments in JLAB or EIC.

Keywords
Quantum Chromodynamics, Generalised Parton Distributions, Factorization,
Renormalization

The title of the thesis in Polish
Faktoryzacja w twardych procesach ekskluzywnych. Poprawki jednopȩtlowe do
fotoprodukcji par fotonów na protonie.

###### Contents

  1. 1 Introduction
  2. 2 The QCD factorization in hard exclusive processes
    1. 2.1 Leading order factorization  and generalized parton distributions
    2. 2.2 Higher order QCD corrections
      1. 2.2.1 Gauge link
      2. 2.2.2 UV and soft collinear divergences
    3. 2.3 Renormalization of non-local product of operators
    4. 2.4 From evolution of operators in the position space to evolution of GPDs
  3. 3 Photoproduction of photon pairs
    1. 3.1 Kinematics, amplitude, and the differential cross section
    2. 3.2 The LO amplitude
    3. 3.3 Self-energy corrections
      1. 3.3.1 Integral
      2. 3.3.2 Traces
    4. 3.4 Vertex corrections - 3.L and 3.R
    5. 3.5 Vertex corrections - 3.M
    6. 3.6 Box diagrams - 4.L and 4.R
      1. 3.6.1 Localizing divergences inside the integral
      2. 3.6.2 Trace - terms leading to divergences
      3. 3.6.3 The finite part of the amplitude
    7. 3.7 The 5-point diagram
      1. 3.7.1 The divergent part
      2. 3.7.2 The finite part
    8. 3.8 Factorization theorem at the 1-loop order
  4. 4 Conclusions
  5. A Derivation of the finite part of 4- and 5- point loop diagrams.
  6. B The full form of the amplitude

##  Chapter 1 Introduction

As far as we know, hadrons are the smallest composite objects in the Universe.
The first indication that nucleons, unlike electrons, are not pointlike
particles came from the measurement of the magnetic moment of the proton
performed by Stern and collaborators [1] already in 1933. More than twenty
years later, analysis of elastic scattering of several hundreds MeV electrons
off nuclei allowed R. Hofstadter to determine the electric charge
distributions in protons and neutrons [2] (the charge r.m.s. radius of the
proton is now estimated to be approximately $0.87\times 10^{-15}$ m [3]).
Finally, at the end of the sixtieth of the XX century, the experiments, in
some respect analogous to the historic one of Rutherford which revealed the
existence of atomic nuclei, with deeply inelastic scattering (the so-called
DIS) of leptons on nucleons performed at Stanford Linear Accelerator Center
(SLAC) [4] inspired Feynman to formulate the parton model [5, 6, 7, 8] in
which hadrons probed by highly virtual photons can be treated as composed of
almost noninteracting pointlike constituents. These, called by Feynman
partons, were subsequently identified with quarks, first proposed by Zweig and
Gell-Mann [10, 11], and spin 1 bosons (gluons) which are now known to be the
elementary constituents of hadrons.

The proper theory of interactions of quarks and gluons is a non-abelian
quantum gauge field theory known as Quantum Chromodynamics (QCD). The
effective coupling constant of this theory becomes large at long distances
(or, equivalenty, at small transfers of momenta) [12]. Therefore, applications
of the perturbative expansion to QCD are limited to only special classes of
high energy processes which are characterized by the existence of at least one
“hard scale” $Q$, e.g. of a large (compared to a typical hadronic mass scale
set by $\Lambda_{\rm QCD}$ \- the intrinsic energy scale of QCD) four-momentum
transfer $-q^{2}=Q^{2}$ ($\sqrt{|q^{2}|}\gg\Lambda_{\rm QCD}$). This makes the
study of hadrons, their reactions and structure, a great challenge: rather
advanced quantum field theory tools are required in order to extract the
relevant information from QCD. Already to obtain precise theory predictions
for scaling violation in DIS experiments (i.e. to go beyond the simple Feynman
parton picture) one needs to employ the so-called operator product expansion.
It allows in this case to rigorously split the computed cross section of the
inclusive process into the part (Wilson coefficients) which is dominated by
the aforementioned hard scale $Q$ and can, owing to this, be computed
perturbatively, and the other part which can only be phenomenologically
parametrized in terms of the parton distribution functions (PDFs) $f_{p}(x)$
introduced originally at the intuitive level by Feynman in his (so-called
naive) parton model and which are given the interpretation [6, 7, 8, 9] of the
probability densities that inside the hadron a parton $p$ of a given type
carries a fraction $x$ of the total hadron four-momentum. Measurements which
confirmed the scaling violation, i.e. the predicted by QCD changes of the
distributions $f_{p}(x)$ with the changes of the hard scale $Q$, were
important tests of QCD which together with the observation by the DORIS and
Petra experiments in Hamburg of the three-jet events in $e^{-}e^{+}$
collissions essentially established it as a correct theory of strong
interactions.

Amplitudes of more complicated processes, in particular of exclusive processes
in which characteristics of not only the initial but also of the final state
hadron(s) are specified (and measured experimentally), usually cannot be
analysed as rigorously as DIS. Nevertheless, one commonly believes that
amplitudes of the studied process which are characterized by a hard scale can
still be split into a perturbatively computable (within QCD) hard part
(dominated by the scale $Q$) and the rest which is nonperturbative and must be
treated phenomenologically, i.e. parametrized in agreement with the known
symmetries of the underlying theory (QCD), determined experimentally in one
group of experiments and then used to predict results of another group of
experiments. For example, analysis of the data accumulated over years in DIS-
type experiments carried out with the help of the HERA accelerator
(experiments ZEUS and H1) allowed to obtain parametrizations of PDFs of the
parton model and are now routinely being used to predict rates of processes
studied experimentally in the proton-proton inelastic collisions in LHC.

The principle on which the possibility of splitting amplitudes of exclusive
hadronic processes into perturbative and nonperturbative parts relies is the
so-called factorization which operates at the level of amplitudes. It is not
as rigorous as the OPE, the use of which is limited to only a rather narrow
class of directly experimentally measurable processes, but can, at least in
the case of some special processes characterized by a hard scale $Q$, be
proven with a fair degree of rigour (enough to satisfy most physicists but few
pure mathematicians - as my and prof. Chankowski’s hero, A.B. Pippard, would
say). Demonstrating factorization requires, however, a careful analysis of
Feynman diagrams and their singularities in various kinematical regimes, both
ultraviolet (UV), infrared and/or collinear, in order to be able to isolate
them and properly absorb into the nonperturbative part of the amplitude. As
far as exclusive processes with one and the same single hadron (usually
nucleon) in the initial and final state are concerned, the nonperturbative
parts of their amplitudes are parametrized in terms of the Generalized Parton
Distributions (GPDs), which are also known as Off-Forward Parton
Distributions. GPDs which, similarly to PDFs, describe properties of a given
hadron and are, therefore, process-independent, significantly generalize PDFs
(which can be obtained from the GPDs by taking their forward limit) and as
such contain much more detailed information on the hadron structure, e.g. on
angular momenta of partons of a given kind or on their spacelike distribution
inside the hadron [14, 15, 17, 16], than do the PDFs. Still, the strategy of
their use is similar to that of the latter: they should be extracted from the
data on one class of processes and applied to predict rates of other classes
of processes.

Extraction of GPDs from experimental data is complicated, because cross
sections of processes sensitive to parton distributions are rather small and
because what actually can be extracted are only convolutions of GPDs with some
other, process-dependent, distributions [7, 17, 18]. Hence, to reliably
determine generalized distributions, processes of different kinds must be
considered in conjunction. Moreover, since at energy scales of order few GeV
($|q|^{2}\sim 10$ GeV2) relevant for these processes [17, 18, 19, 20], the
strong coupling constant $\alpha_{s}$ is of order $\sim 0.1$, higher order
corrections must necessarily be taken into account. One-loop results
concerning electroproduction of a real photon, of photoproduction of lepton
pair or electroproduction of a meson are already available - they have been
given in the papers [22], [23] and [24], respectively. Recently, the process
of photoproduction of two photons on nucleon has been considered, but only at
the leading order [20, 21].

The main goal of this work is to extend the analysis of factorization of the
amplitude of diphoton photoproduction on proton, i.e. of the process $\gamma
p\rightarrow\gamma^{\prime}\gamma^{\prime\prime}p^{\prime}$, given in [20,
21], to the first nontrivial order. The work consists of two parts. In Chapter
2 I describe in some details the collinear factorization. In particular I
explain how factorization should be handled in higher orders, when the
perturbatively computed part of the amplitude develops in addition to UV also
infrared divergences resulting from diagrams containing loops with propagators
of partons treated as massless quanta. It will be seen that if factorization
is to hold also in higher orders, the infrared (collinear) divergences must
cancel against the ones present in parton distributions defined as (Fourier
transforms - see below - of) matrix elements between hadron states of
appropriately defined renormalized nonlocal composite operators. Parton
distributions defined in this way are frequently called, somewhat
misleadingly, renormalized distributions.111Parton distributions defined in
terms of unrenormalized operators are accordingly and also typically
misleadingly called bare parton distributions. To be sure: whether defined in
the leading order or in higher orders they are extracted from the same
measured quantities. It it thus only their theoretical definition which is
given either in terms of bare or in terms of renormalized operators; the GPDs
are always the same quantities. If such a cancellation does indeed hold, the
amplitude can still be split (i.e. factorized) into the perturbative (hard)
part computable with the help of renormalized perturbative expansion and the
nonperturbative one parametrized now in terms of the renormalized GPDs. Since
this cancellation is not obvious (and there are known situations, in which it
does not obtain [25]), it must be checked for each process independently.
Proofs that factorization holds in the leading order in $|Q|^{-1}$ and up to
an arbitrary order in perturbative QCD can be carried out only for a limited
family of processes. Since renormalization of GPDs (i.e. the construction of
renormalized operators in terms of which GPDs must be defined in higher
orders) is more complicated than, for example, renormalization of conventional
$S$-matrix elements or Green’s functions of elementary or local composite
operators (here the renormalized operators must be defined in a gauge
invariant way, i.e. must be supplemented by the Wilson gauge link operator), I
provide in Section 2.2 a detailed demonstration of the necessary
renormalization procedure at the one-loop level, and derive the resulting
evolution equations (analogous to the ones governing the scale dependence of
the ordinary PDFs) satisfied by the resulting generalized parton densities.

The second part of the work is devoted to the analysis in the next-to-leading
order of photoproduction of photon pairs with large invariant mass (which sets
the hard scale $Q$) on proton. This process is especially interesting from the
point of view of the theoretical analysis, because it is the simplest one of
all $2\rightarrow 3$ processes that can be studied within the framework of
collinear factorization discussed in the first part of the work. The
corresponding amplitude is sensitive only to combinations of partons
distributions, which are odd with respect to the charge conjugation. That
provides the additional source of information on GPDs. The leading-order
analysis [20] shows, that this process can in principle be studied
experimentally at the JLab facilities which provide intense beams of quasi-
real photons and in the planned Electron-Ion Collider (EIC). Moreover, some
experimental data pertaining to this process probably already exists, as
exclusive experiments with protons and highly energetic quasi-real photons
have already been performed to measure another process, the Timelike Compton
Scattering. In Chapter 3 I present a detailed computation of one loop
corrections to the amplitude of the
$\gamma~{}\\!p\rightarrow\gamma^{\prime}\gamma^{\prime\prime}p^{\prime}$
process. I give a proof of the validity of the collinear factorization in it
at NLO, compute the finite part of the amplitude and the differential cross
section of the process. I also introduce some techniques which may, I believe,
prove useful also in theoretical studies of other processes.

Before proceeding with the main part of the work, it is convenient to define
the notation used throughout it. I will introduce here also the definitions of
the Generalized Parton Distributions since, in my opinion, it is better to do
it before explaining how they arise in actual computations of amplitudes of
specific physical processes.

In defining GPDs it is customary to introduce the so-called Sudakov
decomposition of four-vectors $v^{\mu}$ defined in the following way:

$v^{\mu}=v^{+}n^{\mu}_{+}+v^{-}n^{\mu}_{-}+v^{\mu}_{\perp},$ (1.1)

where

$n^{\mu}_{\pm}\equiv\frac{1}{\sqrt{2}}[1,0,0,\pm 1]^{\mu},$ (1.2)

and $v^{\mu}_{\perp}$ denotes the spatial part of the vector projected onto
the $xy$ plane. In the analogous way one introduces the lightcone basis of the
gamma matrices:

$\gamma^{\pm}\equiv\frac{1}{\sqrt{2}}\big{(}\gamma^{0}\pm\gamma^{3}\big{)}.$
(1.3)

Consider now two states $|N(\mathbf{p}_{1},\sigma_{1})\rangle$ and
$|N(\mathbf{p}_{2},\sigma_{2})\rangle$ of the same hadron $N$ with the four-
momenta $p_{1}^{\mu}=(E_{1},\mathbf{p}_{1})$ and
$p_{2}^{\mu}=(E_{2},\mathbf{p}_{2})$
($E_{i}=\sqrt{m^{2}_{N}+\mathbf{p}_{i}^{2}}$, $i=1,2$, where $m_{N}$ is the
mass of the hadron), and $\sigma_{i}$ its spin projection. The ”$+$”
components $p_{i}^{+}=(E_{i}+p_{i}^{z})/\sqrt{2}$, $i=1,2$ of the four-momenta
$p_{1}$ and $p_{2}$ will be parametrized with the help of two parameters:
$\xi$ ($-1<\xi<1$) and $p^{+}\geq 0$ defined by

$p_{1}^{+}=(1+\xi)~{}\\!p^{+},\quad p_{2}^{+}=(1-\xi)~{}\\!p^{+}.$ (1.4)

It follows that $p^{+}=(p_{1}^{+}+p_{2}^{+})/2\geq 0$ and
$\xi=(p_{1}^{+}-p_{2}^{+})/(p_{1}^{+}+p_{2}^{+})$. We assume that both hadrons
move nearly with the speed of light in the direction of the $z$-axis, so that
$p^{+}\gg m_{N}$.

Generalized Parton Distribution can be defined by considering in the limit
$p^{+}\rightarrow\infty$ of the following Fourier transforms of matrix
elements of non-local composite operators:

$\displaystyle\int\\!\frac{dz^{-}}{2\pi}~{}\\!e^{ixp^{+}z^{-}}\left.\\!\\!\langle
N(\mathbf{p}_{2},\sigma_{2})|\bar{\psi}_{q}(-z/2)[-z/2,z/2]\gamma^{+}\psi_{q}(z/2)|N(\mathbf{p}_{1},\sigma_{1})\rangle\right|_{z^{+},z_{\perp}=0},$
(1.5)
$\displaystyle\int\\!\frac{dz^{-}}{2\pi}~{}\\!e^{ixp^{+}z^{-}}\left.\\!\\!\langle
N(\mathbf{p}_{2},\sigma_{2})|\bar{\psi}_{q}(-z/2)[-z/2,z/2]\gamma^{+}\gamma_{5}~{}\\!\psi_{q}(z/2)|N(\mathbf{p}_{1},\sigma_{1})\rangle\right|_{z^{+},z_{\perp}=0},\phantom{aa}$
(1.6)
$\displaystyle\int\\!\frac{dz^{-}}{2\pi}~{}\\!e^{ixp^{+}z^{-}}\left.\\!\\!\langle
N(\mathbf{p}_{2},\sigma_{2})|G^{+\mu}(-z/2)[-z/2,z/2]G^{\phantom{a}+}_{\mu}(z/2)|N(\mathbf{p}_{1},\sigma_{1})\rangle\right|_{z^{+},z_{\perp}=0},$
(1.7)
$\displaystyle\int\\!\frac{dz^{-}}{2\pi}~{}\\!e^{ixp^{+}z^{-}}\left.\\!\\!\langle
N(\mathbf{p}_{2},\sigma_{2})|G^{+\mu}(-z/2)[-z/2,z/2]\tilde{G}^{\phantom{a}+}_{\mu}(z/2)|N(\mathbf{p}_{1},\sigma_{1})\rangle\right|_{z^{+},z_{\perp}=0}.$
(1.8)

$\psi_{q}$ denotes here the bare (i.e. the Heisenberg picture operator the
scale of which is fixed by the canonical kinetic term in the QCD lagrangian
density) quark field operator of flavor $q$, $G^{\mu\nu}$ is the bare field
strength tensor of the gluon field and
$\tilde{G}^{\mu\nu}=\frac{1}{2}\epsilon^{\mu\nu\kappa\lambda}G_{\kappa\lambda}$
the tensor dual to it; $\epsilon^{\mu\nu\kappa\lambda}$ is the completely
antisymmetric tensor (I use $\epsilon^{0123}=+1$). The matrix $\gamma_{5}$ is
defined as $\gamma_{5}=\gamma^{5}=i\gamma^{0}\gamma^{1}\gamma^{2}\gamma^{3}$.
Strictly speaking, in the leading order the operator denoted in (1.5-1.8)
$[-z/2,z/2]$, called the gauge link is replaced by the unit operator. It must
be included in higher orders (together with appropriate operator counterterms
which complete the definitions of the renormalized nonlocal operators) and is
defined in the following way:

$[z_{1},z_{2}]\equiv
P\exp\Big{(}ig\\!\int^{z_{1}}_{z_{2}}\\!d\xi_{\mu}\mathbb{A}^{\mu}(\xi)\Big{)}.$
(1.9)

The symbol $P$ denotes here the path ordering and

$\mathbb{A}^{\mu}\equiv\sum_{a}t^{a}A^{\mu,a},$ (1.10)

where $t^{a}$ are the generators (taken in the fundamental representation) of
the $SU(3)$ gauge group; $a=1,\dots 8$ is the adjoint color index. Presence of
this operator ensures gauge invariance of GPDs. How it arises in actual
calculations is explained in Section 2.2.1. In the literature, see e.g. [29,
30] products of two operators separated by a light-like distance and the gauge
link, as in the matrix elements (1.5-1.8), are often called string operators.

In most applications the hadron $N$ is a spin-$1/2$ particle (proton or
neutron). In such cases the most general forms of the Fourier transforms
(1.5)-(1.8) can be, taking into account Lorentz invariance, parity, time
reversal and charge conjugation symmetries of the underlying theory,
parametrized in terms of the functions $H^{q,g}$, $\tilde{H}^{q,g}$, $E^{q,g}$
and $\tilde{E}^{q,g}$ (their definitions given here agree with those used in
[27]) called generalized parton distributions
($\Delta^{\kappa}\equiv(p_{2}-p_{1})^{\kappa}$):

$\displaystyle\eqref{eq-
gpd-q1}=\frac{1}{p^{+}}\Big{[}H^{q}(x,\xi,t)~{}\\!\bar{u}(\mathbf{p}_{2},\sigma_{2})\gamma^{+}u(\mathbf{p}_{1},\sigma_{1})+E^{q}(x,\xi,t)~{}\\!\bar{u}(\mathbf{p}_{2},\sigma_{2})~{}\\!\frac{i\sigma^{+\kappa}\Delta_{\kappa}}{2m}~{}\\!u(\mathbf{p}_{1},\sigma_{1})\Big{]},$
(1.11) $\displaystyle\eqref{eq-
gpd-q2}=\frac{1}{p^{+}}\Big{[}\tilde{H}^{q}(x,\xi,t)~{}\\!\bar{u}(\mathbf{p}_{2},\sigma_{2})\gamma^{+}\gamma_{5}~{}\\!u(\mathbf{p}_{1},\sigma_{1})+\tilde{E}^{q}(x,\xi,t)\bar{u}(\mathbf{p}_{2},\sigma_{2})~{}\\!\frac{\gamma_{5}\Delta^{+}}{2m}~{}\\!u(\mathbf{p}_{1},\sigma_{1})\Big{]},$
(1.12) $\displaystyle\eqref{eq-
gpd-g1}=\frac{1}{2}\Big{[}H^{g}(x,\xi,t)~{}\\!\bar{u}(\mathbf{p}_{2},\sigma_{2})\gamma^{+}u(\mathbf{p}_{1},\sigma_{1})+E^{g}(x,\xi,t)~{}\\!\bar{u}(\mathbf{p}_{2},\sigma_{2})~{}\\!\frac{i\sigma^{+\kappa}\Delta_{\kappa}}{2m}~{}\\!u(\mathbf{p}_{1},\sigma_{1})\Big{]},$
(1.13) $\displaystyle\eqref{eq-
gpd-g2}=\frac{i}{2}\Big{[}\tilde{H}^{q}(x,\xi,t)~{}\\!\bar{u}(\mathbf{p}_{2},\sigma_{2})\gamma^{+}\gamma_{5}~{}\\!u(\mathbf{p}_{1},\sigma_{1})+\tilde{E}^{q}(x,\xi,t)~{}\\!\bar{u}(\mathbf{p}_{2},\sigma_{2})~{}\\!\frac{\gamma_{5}\Delta^{+}}{2m}~{}\\!u(\mathbf{p}_{1},\sigma_{1})\Big{]}.$
(1.14)

Owing to the Lorentz invariance, the distributions $H$, $\tilde{H}$, $E$ and
$\tilde{E}$ depend on only 3 variables: $x$, $\xi$ and $t=\Delta^{2}$;
moreover, as can be shown [28], $-1<x<1$.

As already said, beyond the leading order GPDs (called then renormalized
Generalized Parton Distributions) must be defined in the analogous way but as
matrix elements of renormalized string operators. The form of these
renormalized operators is established in Section 2.3. Because of cancellation
of divergences, which is an essential element of factorization, all
observables must be then expressed in terms of these renormalized GPDs. How
the matrix elements (1.5-1.8), and therefore the generalized parton
distributions, appear in actual computation of amplitudes of physical
processes, and how the generalized parton distributions vary with changes of
the hard scale $Q$ (i.e. what renormalization group they satisfy) will be
shown in Chapter 2.

##  Chapter 2 The QCD factorization in hard exclusive processes

###  2.1 Leading order factorization
and generalized parton distributions

One of the simplest processes allowing to extract GPDs from the experimental
data is the so-called Deeply Virtual Compton Scattering (DVCS) [17] in which a
nucleon $N$ of mass $m_{N}$ absorbs a highly virtual photon carrying the four-
momentum $q$ ($-q^{2}\equiv Q^{2}\gg m^{2}_{N}$) and subsequently emits a real
photon:
$\gamma^{\ast}(q)N(\mathbf{p})\rightarrow\gamma(\mathbf{q}^{\prime})N(\mathbf{p}^{\prime})$.
This process can be experimentally studied as a sub-process shown in Fig. 2.1
of the photon electroproduction that is, in the reaction $e^{-}N\rightarrow
e^{-}\gamma N$.

$N(\mathbf{p})$$N(\mathbf{p}^{\prime})$$e^{-}(\mathbf{k})$$q$$e^{-}(\mathbf{k}^{\prime})$$\gamma(\mathbf{q}^{\prime})$
Figure 2.1: DVCS as a subprocess of the photon electroproduction on nucleon.
$k$ and $k^{\prime}$ are the four-momenta of the initial and final electron.
The hard scale is in this case set by the virtuality of the exchanged photon
$-q^{2}=-(k-k^{\prime})^{2}$. The gray blob is the part of the amplitude which
cannot be computed perturbatively. Another subprocess contributing to the same
photon electroproduction on nucleon is the Bethe-Heitler subprocess in which
nucleon absorbs only the virtual photon and the final state photon is emitted
from one of the electron lines and, therefore, the “hadronic part” of the
amplitude is parametrized in term of well-known elastic nucleon form-factors
[2].

Another process sensitive to the generalized parton distributions [18] is the
Timelike Compton Scattering (TCS),
$\gamma(\mathbf{q})N(\mathbf{p})\rightarrow\gamma^{\ast}(q^{\prime})N(\mathbf{p}^{\prime})$,
which can be regarded as a “reversed” DVCS: the nucleon absorbs a real photon
and emits a virtual photon, which subsequently produces a lepton pair. Both
DVCS and TCS are the limiting cases of a more general process called Double
Deeply Virtual Compton Scattering (DDVCS), which can be represented as
$\gamma^{\ast}(q)N(\mathbf{p})\rightarrow\gamma^{\ast}(q^{\prime})N(\mathbf{p}^{\prime})$,
in which at least one photon is far off-shell, so that $|q^{2}|+|q^{\prime
2}|\gg m^{2}_{N}$. An analysis of DDVCS at the next-to-leading order, together
with the proof of factorization of its amplitude at the 1-loop level, has been
presented in [31].

Before considering the process
$\gamma(q)N(\mathbf{p})\rightarrow\gamma(\mathbf{q}^{\prime})\gamma(\mathbf{q}^{\prime\prime})N(\mathbf{p}^{\prime})$
which is the main goal of this work, let us discuss factorization of the
photon electroproduction amplitude starting from the basic principles of
quantum field theory. In view of the fact that hadrons are not “quanta” of
elementary QCD field operators, the $S$-matrix element

$\displaystyle\langle\left(N(\mathbf{p}^{\prime},\sigma^{\prime})e^{-}(\mathbf{k}^{\prime})\gamma(\mathbf{q}^{\prime})\right)_{\rm
out}|\left(N(\mathbf{p},\sigma)e^{-}(\mathbf{k})\right)_{\rm
in}\rangle~{}\\!,$ (2.1)

(the subscripts “in” and “out” denote the in and out asymptotic eigenstates of
the QCD+QED Hamiltonian) corresponding to the experimentally realizable
process $N(\mathbf{p})e^{-}(\mathbf{k})\rightarrow
N(\mathbf{p}^{\prime})e^{-}(\mathbf{k}^{\prime})\gamma(\mathbf{q}^{\prime})$
must be extracted from the (connected part of the) vacuum Green’s function

$\displaystyle\langle\Omega_{\rm
out}|T\\!\left[O(y^{\prime})O^{\dagger}(y)\psi_{(e)}(w^{\prime})\bar{\psi}_{(e)}(w)A^{\mu}(x^{\prime})\right]\\!|\Omega_{\rm
in}\rangle~{}\\!,$

where $O(y^{\prime})$ is an appropriately constructed composite Heisenberg
picture operator ($\psi_{(e)}$ and $A^{\mu}$ are the electron and photon
Heisenberg picture elementary field operators) using the Lehmann-Symanzik-
Zimmermann reduction formula [33]. If $N$ stands for proton, as $O$ one can
take e.g. the composite operator transforming as the spin $1/2$ representation
of the Lorentz group constructed out of three quark (Heisenberg picture) field
operators
$\varepsilon_{ijk}(C^{-1})^{\beta^{\prime}\beta}\psi^{(u)i}_{\beta^{\prime}}\psi^{(d)j}_{\beta}\psi^{(u)k}_{\alpha}$,
where $C^{-1}$ is the charge conjugation matrix (which makes a Lorentz singlet
out of two quark fields) and $\varepsilon_{ijk}$ makes the operator a colour
singlet. Because the quantum numbers match, the operator $O(y^{\prime})$
certainly has a nonzero matrix element between the vacuum state
$|\Omega\rangle$ and the (in or out) state $|N(\mathbf{p},\sigma)\rangle$ (an
eigenstate of the complete Hamiltonian of QCD+QED) which by symmetry
principles must take the form

$\langle\Omega_{\rm
out}|O^{H}_{\alpha}(y^{\prime})|(N(\mathbf{p}^{\prime},\sigma^{\prime})_{\rm
out}\rangle=\mathcal{Z}_{O}^{1/2}e^{-ip^{\prime}\cdot
y^{\prime}}u_{\alpha}(\mathbf{p}^{\prime},\sigma^{\prime}),$ (2.2)

in which $u_{\alpha}(\mathbf{p}^{\prime},\sigma^{\prime})$ is the ordinary
spinor and $\mathcal{Z}_{O}$ is some non-zero constant (which is the property
of the operator $O$). However since the electromagnetic interaction is weak,
it is sufficient to consider the amplitude of the process
$N(\mathbf{p})e^{-}(\mathbf{k})\rightarrow
N(\mathbf{p}^{\prime})e^{-}(\mathbf{k}^{\prime})\gamma(\mathbf{q}^{\prime})$
in the lowest possible order in the coupling $e$. This allows to write it as a
sum of two terms: one the main part of which is the amplitude of the DVCS
subprocess and another one corresponding to the Bethe-Heitler subprocess. The
latter is well-known (it can be written in terms of the elastic nucleon form-
factors) and will not be discussed here. Thus in this approximation111As
usually we consider only connected parts of $S$-matrix elements and of the
corresponding Green’s functions; this will not be indicated explicitly.

$\displaystyle\langle\left(N(\mathbf{p}^{\prime},\sigma^{\prime})e^{-}(\mathbf{k}^{\prime})\gamma(\mathbf{q}^{\prime})\right)_{\rm
out}|\left(N(\mathbf{p},\sigma)e^{-}(\mathbf{k})\right)_{\rm
in}\rangle=(2\pi)^{4}\delta^{4}(q^{\prime}+k^{\prime}+p^{\prime}-k-p)\phantom{aaaaaaaaaaaaaaaa}$
$\displaystyle\times\left\\{\epsilon^{*}_{\nu}(\mathbf{q}^{\prime})~{}\\!(-ie)~{}\\!W^{\mu\nu}(\mathbf{p}^{\prime},\mathbf{p};q^{\prime},q)~{}\\!(-ie)~{}\\!{-i\over
q^{2}}~{}\\!ie~{}\\![\bar{u}(\mathbf{k}^{\prime},s^{\prime})\gamma_{\mu}u(\mathbf{k},s)]-i{\cal
A}_{\rm BeHe}\right\\},$

where ${\cal A}_{\rm BeHe}$ stands for the invariant amplitude of the Bethe-
Heitler subprocess, $\epsilon^{\ast}_{\nu}(\mathbf{q}^{\prime})$ is the
polarization vector of the outgoing photon, $-i/q^{2}$ is the virtual photon
propagator (marked $q$ in Fig. 2.1) and the tensor
$W^{\mu\nu}(\mathbf{p}^{\prime},\mathbf{p};q^{\prime},q)$ is defined by the
double Fourier transform

$\displaystyle(2\pi)^{4}\delta^{(4)}(q^{\prime}+p^{\prime}-q-p)~{}\\!W^{\mu\nu}(\mathbf{p}^{\prime},\mathbf{p};q^{\prime},q)=\int\\!d^{4}x^{\prime}d^{4}x~{}\\!e^{iq^{\prime}\cdot
x^{\prime}}e^{-iq\cdot
x}~{}\\!W^{\mu\nu}(\mathbf{p}^{\prime},\mathbf{p};x^{\prime},x)~{}\\!,$ (2.3)

of the matrix element

$\displaystyle
W^{\mu\nu}(\mathbf{p}^{\prime},\mathbf{p};x^{\prime},x)=\langle(N(\mathbf{p}^{\prime},\sigma^{\prime}))_{\rm
out}|TJ^{\mu}_{H}(x)J^{\nu}_{H}(x^{\prime})|(N(\mathbf{p},\sigma))_{\rm
in}\rangle~{}\\!,$ (2.4)

of the chronological product of two (Heisenberg picture) operators
$J_{H}^{\mu}$ of the hadronic electromagnetic current

$J_{H}^{\mu}(x)=\sum_{f}e_{f}\bar{\psi}^{H}_{(f)}(x)\gamma^{\mu}\psi^{H}_{(f)}(x)~{}\\!,$
(2.5)

where the sum is over all light quark flavours $f$ and $e_{f}$ are the quark
electric charges (in units of $e>0$). The matrix element (2.4) can be
extracted using the LSZ prescription [33] from the appropriate vacuum Green’s
function similarly as the $S$-matrix element (2.1), e.g. as

$\displaystyle
W^{\mu\nu}(\mathbf{p}^{\prime},\mathbf{p};x^{\prime},x)=\left(i\mathcal{Z}_{O}^{-1/2}\right)^{2}\\!\\!\int\\!d^{4}y^{\prime}d^{4}y~{}\\!e^{ip^{\prime}\cdot
y^{\prime}}e^{-ip\cdot
y}~{}\\!\bar{u}(\mathbf{p}^{\prime},\sigma^{\prime})~{}\\!\phantom{aaaaaaaaaaaaaaaaaaaaaaaaaaaaa}$
(2.6)
$\displaystyle(i\overset{\rightarrow}{\not{\partial}}_{y^{\prime}}-m_{N})\langle\Omega_{\rm
out}|TJ^{\mu}_{H}(x)J^{\nu}_{H}(x^{\prime})O_{H}(y^{\prime})\bar{O}_{H}(y)|\Omega_{\rm
in}\rangle(i\overset{\leftarrow}{\not{\partial}}_{y}+m_{N})u(\mathbf{p},\sigma)~{}\\!.$

The application of the Gell-Mann - Low formula222It is interesting to note
that the Gell-Mann - Low formula can be directly applied here owing to nonzero
quark masses which are sources of explicit chiral symmetry breaking which
removes degeneracy of a continuum of QCD vacua which, if quark masses were
vanishing, would exist making impossible adiabatic reaching the full QCD
vacuum from the free one. (see e.g. [33]) allows in principle to compute
vacuum matrix elements like the one in (2.6) perturbatively expressing them
through the chronological products of the interaction picture operators
according to the standard rule

$\displaystyle\langle\Omega_{\rm
out}|TJ^{\mu}_{H}(x)J^{\nu}_{H}(x^{\prime})O_{H}(y^{\prime})\bar{O}_{H}(y)|\Omega_{\rm
in}\rangle\phantom{aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa}$
$\displaystyle=\langle\Omega_{0}|TJ^{\mu}_{I}(x)J^{\nu}_{I}(x^{\prime})O_{I}(y^{\prime})\bar{O}_{I}(y)\exp\\!\left(-i\\!\\!\int\\!d^{4}w~{}\\!\mathcal{H}^{\rm
int}_{I}(w)\\!\right)\\!|\Omega_{0}\rangle~{}\\!,$ (2.7)

in which $\mathcal{H}^{\rm int}_{I}(w)$ is the interaction hamiltonian density
of QCD expressed in terms of the interaction picture quark and gluon field
operators and the operators $O_{I}$ and $J^{\mu}_{I}$ in this picture are
obtained by simply replacing in their Heisenberg counterparts the elementary
Heisenberg picture operators by the interaction picture ones. The Green’s
function can be therefore computed perturbatively by expanding the exponens
and using the Wick theorem. Because of the large value of the QCD coupling
$g_{s}$, it is clear, however, that the expansion obtained in this way cannot
be truncated to its few first terms only and for this reason the entire
prescription for computing
$W^{\mu\nu}(\mathbf{p}^{\prime},\mathbf{p};q^{\prime},q)$ remains rather
formal. To make some progress, one has therefore to reorganize the formal
perturbative expansion obtained from (2.7) in order to pick up this part of it
which through which the large momentum transfer $q$ flows and which, owing to
this, can be reliably computed evaluating only a few subdiagrams (as a small
value of the QCD coupling can be used for this purpose), and the rest which
can be identified as giving the expansion of matrix elements analogous to the
left hand side of (2.7), but with the electromagnetic currents replaced by
some other operators. In other words one wants to represent the right hand
side of the formula (2.3) defining the tensor
$W^{\mu\nu}(\mathbf{p}^{\prime},\mathbf{p};q^{\prime},q)$ by the sum of
convolutions of “hard” and “nonperturbative” parts of the schematic form
(possible contractions of spinor and/or Lorentz indices other than ${\mu\nu}$
between the amplitudes $\mathcal{T}^{\mu\nu}_{W}$ and $W$ are not displayed)

$\displaystyle(2\pi)^{4}\delta^{(4)}(q^{\prime}+p^{\prime}-q-p)~{}\\!W^{\mu\nu}(\mathbf{p}^{\prime},\mathbf{p};q^{\prime},q)\phantom{aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa}$
$\displaystyle=\sum_{W}\\!\int\\!{d^{4}k_{1}\over(2\pi)^{4}}\dots
d^{4}k_{n}\over(2\pi)^{4}~{}\\!(2\pi)^{4}\delta^{(4)}(q^{\prime}-q-k_{n}-...-k_{1})~{}\\!\mathcal{T}^{\mu\nu}_{W}(q^{\prime},q,k_{n},\dots,k_{1})$
$\displaystyle\times(2\pi)^{4}\delta^{(4)}(p^{\prime}+k_{1}+...+k_{n}-p)~{}\\!W(k_{n},\dots,k_{1};\mathbf{p}^{\prime},\mathbf{p})~{}\\!,\phantom{aaaaa}$
(2.8)

where the expressions in the second line are Fourier transforms of amplitudes
defined similarly as (2.4) but with the product of two electromagnetic
currents replaced by various products of different numbers of other operators.
One such amplitude (in fact the most important one in the leading order) is

$\displaystyle
W_{(f)\beta\alpha}(\mathbf{p}^{\prime},\mathbf{p};z)=\langle(N(\mathbf{p}^{\prime},\sigma^{\prime}))_{\rm
out}|T\bar{\psi}^{H}_{(f)\beta}(-z/2)~{}\\!\psi^{H}_{(f)\alpha}(z/2)|(N(\mathbf{p},\sigma))_{\rm
in}\rangle~{}\\!.$ (2.9)

It too, similarly as the amplitude (2.6), can formally be represented by a
perturbative expansion by applying the Gell-Mann - Low theorem to the
appropriate vacuum matrix element. The contribution of (2.9) to
$W^{\mu\nu}(\mathbf{p}^{\prime},\mathbf{p};q^{\prime},q)$ can be written in
the form

$\displaystyle
W^{\mu\nu}(\mathbf{p}^{\prime},\mathbf{p};q^{\prime},q)=\int\\!{d^{4}k\over(2\pi)^{4}}~{}\\!\mathcal{T}^{\mu\nu}_{\beta\alpha}(q^{\prime},q,k)\int\\!d^{4}z~{}\\!e^{ik\cdot
z/2+i(k+q-q^{\prime})\cdot
z/2}~{}\\!W_{(f)\beta\alpha}(\mathbf{p}^{\prime},\mathbf{p};z),$ (2.10)

To see this, let’s write explicitly the corresponding term in the sum (2.8) in
the form (see Fig. 2.2):

$q$$q^{\prime}$$N(\mathbf{p},\sigma_{1})$$N(\mathbf{p}^{\prime},\sigma_{2})$$k_{1}$$k_{2}$Nonperturbative
partHard part Figure 2.2: Graphical representation of the contribution (2.10)
to the amplitude $W^{\mu\nu}(\mathbf{p}^{\prime},\mathbf{p};q^{\prime},q)$: a
quark of flavour $f$ and four-momentum $k_{1}$ emerges from the
nonperturbative part, interacts absorbing a large momentum transfer $q$ (the
“hard” part of the amplitude) and reimmerses in the nonperturbative part with
four-momentum $k_{2}=k_{1}+q-q^{\prime}$.

$\displaystyle\int\\!{d^{4}k_{2}\over(2\pi)^{4}}{d^{4}k_{1}\over(2\pi)^{4}}~{}\\!(2\pi)^{4}\delta^{(4)}(q^{\prime}-q-k_{2}-k_{1})~{}\\!\mathcal{T}^{\mu\nu}_{\beta\alpha}(q^{\prime},q,k)\phantom{aaaaaaaaaaaaaaaaaaaaaaaaa}$
$\displaystyle\int\\!d^{4}z_{2}~{}\\!d^{4}z_{1}~{}\\!e^{ik_{2}\cdot
z_{2}}~{}\\!e^{ik_{1}\cdot
z_{1}}~{}\\!\langle(N(\mathbf{p}^{\prime},\sigma^{\prime}))_{\rm
out}|T\bar{\psi}^{H}_{(f)\beta}(z_{2})~{}\\!\psi^{H}_{(f)\alpha}(z_{1})|(N(\mathbf{p},\sigma))_{\rm
in}\rangle~{}\\!.$

Making the change of variables: $z_{2}={1\over 2}(w-z)$, $z_{1}={1\over
2}(w+z)$ (the Jacobian is $2^{-4}$) and using the Poincaré covariance of the
matrix element

$\displaystyle\langle(N(\mathbf{p}^{\prime},\sigma^{\prime}))_{\rm
out}|T\psi^{H}_{(f)\beta}(z_{2})~{}\\!\bar{\psi}^{H}_{(f)\alpha}(z_{1})|(N(\mathbf{p},\sigma))_{\rm
in}\rangle\phantom{aaaaaaaaaaaaaaaaaaaaaa}$
$\displaystyle=e^{i(p^{\prime}-p)\cdot
w/2}~{}\\!\langle(N(\mathbf{p}^{\prime},\sigma^{\prime}))_{\rm
out}|T\bar{\psi}^{H}_{(f)\beta}(-z/2)~{}\\!\psi^{H}_{(f)\alpha}(z/2)|(N(\mathbf{p},\sigma))_{\rm
in}\rangle~{}\\!,$

enables one to take explicitly the integral over $w$ which gives
$(2\pi)^{4}\delta^{(4)}(q^{\prime}+p^{\prime}-q-p)$ (times $2^{4}$ which
cancels the Jacobian) and comparing the result with the definition (2.3) one
arrives at the contribution (2.10) to the tensor
$W^{\mu\nu}(\mathbf{p}^{\prime},\mathbf{p};q^{\prime},q)$.

The contribution of the amplitude (2.9) to the tensor
$W^{\mu\nu}(\mathbf{p}^{\prime},\mathbf{p};q^{\prime},q)$, taking into account
only the zeroth order (in $g_{s}$) contributions to
$\mathcal{T}^{\mu\nu}_{\beta\alpha}(q^{\prime},q,k)$, is shown in Figs. 2.3a
and 2.3b. Fig. 2.3c shows one of possible contributions to
$W^{\mu\nu}(\mathbf{p}^{\prime},\mathbf{p};q^{\prime},q)$ of the amplitude
analogous to (2.9) but with four quark fields under the chronological product
(and with the lowest order contribution to the associated amplitude
$\mathcal{T}^{\mu\nu}_{\bar{\alpha}\alpha\bar{\beta}\beta}(k_{1},k_{2},q,q^{\prime})$).
This contribution, which can be interpreted as absorption and emission by
nucleon of two different quarks, is subleading in the limit of high $|q^{2}|$
[8], and will therefore be neglected in our analysis.

$q$$q^{\prime}$$(a)$

$q$$q^{\prime}$$(b)$

$q$$q^{\prime}$$(c)$

Figure 2.3: Contributions to the DVCS amplitude (to the tensor
$W^{\mu\nu}(\mathbf{p}^{\prime},\mathbf{p};q^{\prime},q)$) of the
nonperturbative amplitude (2.9) combined with the lowest order contributions
(diagrams a and b) to the associated perturbative amplitude
$\mathcal{T}^{\mu\nu}_{\beta\alpha}(q^{\prime},q,k)$. Diagram c is the
possible contribution of the nonperturbative amplitude with four quark fields.
Such contributions can be shown to be of a lower order in the hard scale [8],
and are, hence, neglected.

The perturbative amplitudes $\mathcal{T}^{\mu\nu}_{W}$ can be computed using
the standard Feynman rules and are simply equivalent to connected Green’s
functions with amputated external lines. The lowest order contribution to
${\cal T}^{\mu\nu}_{(f)\beta\alpha}(q^{\prime},q,k)$ therefore is (cf. Figs
2.3a and 2.3b):

$\displaystyle{\cal
T}^{\mu\nu}_{(f)\beta\alpha}(q^{\prime},q,k)=e_{f}^{2}\left[\gamma^{\nu}~{}\\!{i(\not\\!k~{}\\!+\not\\!q+m_{f})\over(k+q)^{2}-m^{2}_{f}+i0}~{}\\!\gamma^{\mu}+\gamma^{\mu}~{}\\!{i(\not\\!k~{}\\!-\not\\!q^{\prime}+m_{f})\over(k-q^{\prime})^{2}-m^{2}_{f}+i0}~{}\\!\gamma^{\nu}\right]_{\beta\alpha}~{}\\!.$
(2.11)

It can be viewed as the amplitude of the elastic scattering of the (highly)
virtual photon of four-momentum $q$ on an off-shell quark of flavour $f$ and
four-momentum $k$.

In the limit $-q^{2}\rightarrow\infty$ (physically realized when $-q^{2}\gg
m_{N}$) it is convenient to go to the infinite momentum frame, in which the
$"+"$ components of the four-momenta $p$, $p^{\prime}$ and $q$, as well as the
$"-"$ component of $q$ are much larger than $m_{N}$. Recall, that in this
frame the momenta of hadrons are parametrized in the following way:

$p^{\mu}\approx(1+\xi)~{}\\!p^{+}n^{\mu}_{+}~{}\\!,\qquad
p^{\prime\mu}\approx(1-\xi)~{}\\!p^{+}n^{\mu}_{+}~{}\\!.$ (2.12)

At this point one invokes the basic assumption of the parton model, namely
that in the hard amplitude, that is in the lowest order in the formula (2.11),
the four-momentum of the incoming quark can be taken in the form
$k^{\mu}=k^{+}n_{+}^{\mu}$ (see: [7], page 112, Fig. 5.2). Heuristically one
can argue that this quark, being a constituent of a hadron with large ”$+$”
component of its (on-shell) four-momentum $p^{\mu}$, should also have a large
”$+$” component of the four-momentum $k^{\mu}$, while all other components of
$k^{\mu}$ should be much smaller compared to the components of the $q$ and
$q^{\prime}$ four-momenta (note, that $q-q^{\prime}\approx 2\xi p^{+}n_{+}$).
It should, however, be clearly said, that this is a model assumption which
cannot be justified from the first principles. Using this approximation, the
integration over $dk^{-}d^{2}k_{\perp}$ in the formula (2.10) can be
explicitly performed giving the Dirac delta which allows to perform the
integration over $dz^{+}d^{2}z_{\perp}$. Parametrizing the four-momentum of
the quark ’entering’ the hard part using the variable $x$ as
$k^{+}=(x+\xi)p^{+}$, one can represent the contribution (2.10) to
$W^{\mu\nu}(\mathbf{p}^{\prime},\mathbf{p},q^{\prime},q)$ in the form

$\displaystyle
W^{\mu\nu}(\mathbf{p}^{\prime},\mathbf{p},q^{\prime},q)=p^{+}\\!\\!\int\\!dx~{}\\!\mathcal{T}^{\mu\nu}_{\beta\alpha}(q^{\prime},q,k=(x+\xi)p^{+}n_{+})\\!\int\\!{dz^{-}\over
2\pi}~{}\\!e^{ixp^{+}z^{-}}~{}\\!W_{\beta\alpha}(\mathbf{p}^{\prime},\mathbf{p},z=z^{-}n_{-})~{}\\!.\phantom{aa}$
(2.13)

In order to perform the contraction of the spinor indices in (2.13) one can
write the product $\bar{\psi}_{\beta}\psi_{\alpha}$ of field operators in
(2.9) using the Fierz decomposition (see for example [7], p. 33):

$\bar{\psi}_{\beta}\psi_{\alpha}=\frac{1}{4}\big{(}1\big{)}_{\alpha\beta}\bar{\psi}\psi+\frac{1}{4}\big{(}\gamma_{\mu}\big{)}_{\alpha\beta}\bar{\psi}\gamma^{\mu}\psi+\frac{1}{4}\big{(}\gamma_{5}\gamma_{\mu}\big{)}_{\alpha\beta}\bar{\psi}\gamma^{\mu}\gamma_{5}\psi+\frac{1}{4}\big{(}\gamma_{5}\big{)}_{\alpha\beta}\bar{\psi}\gamma_{5}\psi+\frac{1}{8}\big{(}\sigma^{\mu\nu}\big{)}_{\alpha\beta}\bar{\psi}\sigma^{\mu\nu}\psi,$
(2.14)

where $1$ denotes the unit matrix and
$\sigma^{\mu\nu}=\frac{i}{2}\big{(}\gamma^{\mu}\gamma^{\nu}-\gamma^{\nu}\gamma^{\mu}\big{)}$.
If the quark mass $m_{f}$ in the numerator of (2.11) is neglected, only two
structures, the one with $\gamma^{\lambda}$ and that with
$\gamma^{\lambda}\gamma_{5}$, survive after taking the trace. Moreover,
Lorentz invariance implies that the expressions

$\displaystyle\\!\int\\!{dz^{-}\over 2\pi}~{}\\!e^{ixp^{+}z^{-}}\langle
N(\mathbf{p}^{\prime},\sigma^{\prime})_{\rm
out}|T\bar{\psi}^{H}_{(f)}(-z^{-}n_{-}/2)\gamma^{\lambda}~{}\\!\psi^{H}_{(f)}(z^{-}n_{-}/2)|N(\mathbf{p},\sigma)_{\rm
in}\rangle~{}\\!,$ (2.15) $\displaystyle\\!\int\\!{dz^{-}\over
2\pi}~{}\\!e^{ixp^{+}z^{-}}\langle N(\mathbf{p}^{\prime},\sigma^{\prime})_{\rm
out}|T\bar{\psi}^{H}_{(f)}(-z^{-}n_{-}/2)\gamma^{\lambda}\gamma_{5}~{}\\!\psi^{H}_{(f)}(z^{-}n_{-}/2)|N(\mathbf{p},\sigma)_{\rm
in}\rangle~{}\\!,$

can depend on the index $\lambda$ only through $p^{\lambda}$ and
$p^{\prime\lambda}$. Therefore, in the Sudakov frame in which the “tranverse”
(”$\perp$”) and ”$-$” components of $p^{\lambda}$ and $p^{\prime\lambda}$ can
be neglected compared to the ”$+$” one one can approximate (2.13) by

$\displaystyle
W^{\mu\nu}(\mathbf{p}^{\prime},\mathbf{p},q^{\prime},q)\approx\sum_{f}\Bigg{[}{p^{+}\over
4}\\!\int\\!d^{4}x~{}\\!{\rm
tr}\Big{[}\mathcal{T}^{\mu\nu}(q^{\prime},q,(x+\xi)p^{+}n_{+})~{}\\!\gamma^{-}\Big{]}\phantom{aaaaaaaaaaaaaaaaaaaaaaaaaaa}$
(2.16)
$\displaystyle\times\\!\int\\!\frac{dz^{-}}{2\pi}~{}\\!e^{ixp^{+}z^{-}}\\!\langle
N(\mathbf{p}^{\prime},\sigma^{\prime})_{\rm
out}|T\bar{\psi}^{H}_{(f)}(-z^{-}n_{-}/2)\gamma^{+}~{}\\!\psi^{H}_{(f)}(z^{-}n_{-}/2)|N(\mathbf{p},\sigma)_{\rm
in}\rangle$ $\displaystyle+{p^{+}\over 4}\\!\int\\!d^{4}x~{}\\!{\rm
tr}\Big{[}\mathcal{T}^{\mu\nu}(q^{\prime},q,(x+\xi)p^{+}n_{+})~{}\\!\gamma^{-}\gamma_{5}\Big{]}\phantom{aaaaaaaaaaaaaaaaaaaaaaaa}$
$\displaystyle\times\\!\int\\!\frac{dz^{-}}{2\pi}~{}\\!e^{ixp^{+}z^{-}}\\!\langle
N(\mathbf{p}^{\prime},\sigma^{\prime})_{\rm
out}|T\bar{\psi}^{H}_{(f)}(-z^{-}n_{-}/2)\gamma^{+}\gamma_{5}~{}\\!\psi^{H}_{(f)}(z^{-}n_{-}/2)|N(\mathbf{p},\sigma)_{\rm
in}\rangle\Bigg{]}.$

The $1/4$ from the Fierz decomposition is conventionally absorbed in the hard
part. As shown e.g. in [28], in the kinematics of the Sudakov frame in the
matrix elements (2.15) the symbols $T$ of the chronological products can be
removed.

Figure 2.4: The lowest order (in $g_{s}$) diagram which introduces the
dependence of $W^{\nu\mu}(\mathbf{p}^{\prime},\mathbf{p},q^{\prime},q)$ on the
gluon generalized distributions.

Higher order (in $g_{s}$) corrections to the hard amplitude (2.11) will be
discussed in section 2.2. Since gluons can interact with photons through quark
loops, in the sum (2.8) one has to take into account also other terms. In
particular in order $g_{s}^{2}$ there will appear terms corresponding to the
emergence from the nonperturbative part of two gluons. The order $g_{s}^{2}$
contribution to the corresponding hard amplitude
$\mathcal{T}^{\mu\nu}_{\phantom{aa}\rho\kappa}(q^{\prime},q,k)$ is shown in
Fig. 2.4. The additional Lorentz structure of this amplitude (the indices
$\rho\kappa$) can be decomposed into $g_{\rho\kappa}$ and
$\epsilon_{\rho\kappa\sigma\lambda}n_{+}^{\sigma}n_{-}^{\lambda}$ and the
analysis similar to the one discussed above shows that the corresponding
contribution to the $W^{\mu\nu}(\mathbf{p}^{\prime},\mathbf{p},q^{\prime},q)$
amplitude takes the form analogous to (2.16) but with the nonperturbative part
represented by the matrix elements

$\displaystyle\int\\!\frac{dz^{-}}{2\pi}~{}\\!e^{ixp^{+}z^{-}}\\!\bra{N(\mathbf{p}^{\prime},\sigma^{\prime})_{\mathrm{out}}}G^{+\lambda}\big{(}-z^{-}n_{-}/2\big{)}G^{\phantom{a}+}_{\lambda}\big{(}z^{-}n_{-}/2\big{)}\ket{N(\mathbf{p},\sigma)_{\mathrm{in}}},$
(2.17)
$\displaystyle\int\\!\frac{dz^{-}}{2\pi}~{}\\!e^{ixp^{+}z^{-}}\\!\bra{N(\mathbf{p}^{\prime},\sigma^{\prime})_{\mathrm{out}}}G^{+\lambda}\big{(}-z^{-}n_{-}/2\big{)}\tilde{G}^{\phantom{a}+}_{\lambda}\big{(}z^{-}n_{-}/2\big{)}\ket{N(\mathbf{p},\sigma)_{\mathrm{in}}},$
(2.18)

the gluon field strength tensor $G^{\lambda\rho}$ and and its dual
$\tilde{G}^{\lambda\rho}$ are used instead of the gauge potentials – for
reference see Sec. 3.2.6 in [35].

Summarizing, the factorization, sketched above on the example of the DVCS
subprocess, should in general allow to represent the amplitude ${\cal A}$
($W^{\mu\nu}(\mathbf{p}^{\prime},\mathbf{p},q^{\prime},q)$ in the case of
DVCS) of an exclusive process in the schematic form (depicted in Fig. 2.5)

$\mathcal{A}=\sum_{i}\int_{-1}^{1}\\!dx~{}\\!\mathrm{GPD}^{i}(x,\xi,t)~{}\\!\mathcal{T}_{i}(x,\xi,\dots),$
(2.19)

where the summation over $i$ runs over several kinds of quark and gluon
generalized distributions $\mathrm{GPD}^{i}(x,\xi,t)$ ($H^{q,g}$,
$\tilde{H}^{q,g}$, $E^{q,g}$ defined by the decompositions (1.11)-(1.14) in
the case of DVCS), $t=(p^{\prime}-p)^{2}$, and $\mathcal{T}_{i}(x,\xi,\dots)$
are the corresponding hard parts computable perturbatively (the ellipses stand
for all kinematics other kinematical characteristics on which the hard
amplitude may depend).

$=\quad\sum_{i}\int_{-1}^{1}p^{+}dx$$(x+\xi)p^{+}$$(x-\xi)p^{+}$$\mathrm{GPD}^{i}(x,\xi,t)$$\mathcal{T}_{i}$
Figure 2.5: Graphical representation of the factorization of the DVCS
amplitude. $\int_{-1}^{1}p^{+}dx$ corresponds to integration over all possible
values of the ”$+$” components of momentum of the parton entering the hard
part of the amplitude (the upper green blob).

###  2.2 Higher order QCD corrections

In this section we explain, how perturbative corrections in QCD affect the
theoretical predictions of amplitudes of exclusive processes. It turns out,
that whether the collinear factorization can be applied at higher orders of
perturbative calculation is a non-trivial question, since validity of the
factorization requires that infrared and collinear divergences from bare
parton distributions cancel with those present in the hard amplitude. For
example, see [22, 26, 31].

$(a)$

$(b)$

$(c)$

Figure 2.6: Examples of diagrams in the hard-part of the amplitude at the
$\alpha_{S}$ order in QCD. There are two kinds of divergences. The first one –
ultraviolet divergence – is present in diagrams $(a)$ and $(b)$. By power
counting one sees that the diagram $(c)$ is UV safe. However, there are also
so-called collinear divergences present in diagrams $(b)$ and $(c)$, which
result from the integration region, where the momentum of the gluon is nearly
collinear with the one of the on-shell quark. These divergences are discussed
in more detail in Subsec. 2.2.2.

In this work we use the $\overline{MS}$ renormalization scheme together with
the dimensional regularization, where the number of dimensions is
$d=4-\varepsilon$. In general, one should use the renormalized perturbation
theory, so that the hard amplitude is UV-finite, however, at the one-loop
level the renormalization constants cancel, so that in calculations presented
here, for our convenience, we can use the bare fields and unrenormalized
hamiltonian (so that there is no need to include the diagrams with counter-
terms). For reference, see the beginning of Part III in [22]. At higher
orders, that argument does not necessarily work, so that one has to use the
renormalized theory. Computation of regularized loop corrections to the hard
part is a well-known procedure, which, at this point, is not discussed. In
Chapter 3 we present such calculations in the process of photoproduction of
photon pairs. The NLO analysis of DDVCS can be found in [31].

Apart from hard amplitudes, also bare parton distributions develop divergences
at higher orders in strong coupling $\alpha_{S}$. To control them, one needs
to renormalize string operators used in Eqs. (1.5)-(1.8) by finding such
combinations of these operators, that their matrix elements are finite.
Subsequently, one defines renormalized GPDs denoted by $\mathrm{GPD}^{i}_{R}$
using Green’s functions of the renormalized string operators. The relation
between the bare parton distributions used in the formula (2.21) and the
renormalized ones reads:

$\mathrm{GPD}^{i}(x,\xi,t)=\mathrm{GPD}^{i}_{R}(x,\xi,t;\mu_{F})+\frac{2}{\varepsilon}\Big{(}\frac{\mu_{F}^{2}e^{\gamma}}{4\pi\mu_{R}^{2}}\Big{)}^{-\varepsilon/2}\int
dx^{\prime}\>K^{ij}(x,x^{\prime})\mathrm{GPD}^{j}_{R}(x^{\prime},\xi,t;\mu_{F})+\mathcal{O}(\alpha_{S}^{2}),$
(2.20)

where $\gamma$ is the Euler constant, $\mu_{R}$ and $\mu_{F}$ denote the
renormalization scale and factorization scale, respectively. Meaning of these
scales is elaborated in Section 2.3. The matrix $K^{ij}$ is called the
evolution kernel, and there is summation over index $j$. In Section 2.3 it is
explained, how it is derived by considering Green’s functions of string
operators defining GPDs.

Recall the factorization formula (2.19):

$\mathcal{A}=\sum_{i}\int_{-1}^{1}\mathrm{GPD}^{i}(x,\xi,t)\mathcal{T}_{i}(x,\xi,\dots)dx.$
(2.21)

The divergent hard part can be organized in the following way (for simplicity
we omit arguments other than $x$):

$\mathcal{T}^{i}(x)=\mathcal{C}^{i}_{0}(x)+\bigg{[}-\frac{2}{\varepsilon}\Big{(}\frac{|Q^{2}|e^{\gamma}}{4\pi\mu_{R}^{2}}\Big{)}^{-\varepsilon/2}\mathcal{C}^{i}_{coll.}(x)+\mathcal{C}^{i}_{1}(x)\bigg{]},$
(2.22)

where $\mathcal{C}^{i}_{1}$ and $\mathcal{C}^{i}_{coll.}$ are of order
$\alpha_{S}$. Let us write the divergent part of the amplitude at the order
$\alpha_{S}$:

$\displaystyle\mathcal{A}_{div}=$
$\displaystyle-\frac{2}{\varepsilon}\sum_{j}\int_{-1}^{1}dx\>\mathrm{GPD}^{j}_{R}(x)\mathcal{C}^{j}_{coll.}(x)+$
(2.23) $\displaystyle+$
$\displaystyle\frac{2}{\varepsilon}\sum_{i,j}\int_{-1}^{1}dx\>dx^{\prime}\>K^{ij}(x,x^{\prime})\mathrm{GPD}^{j}_{R}(x^{\prime})\mathcal{C}^{i}_{0}(x^{\prime}).$

The cancellation of divergences occurs if

$\mathcal{C}^{j}_{coll.}(x)=\sum_{i}\int_{-1}^{1}K^{ij}(y,x)\mathcal{C}^{i}_{0}(y)dy.$
(2.24)

If that is the case, the the amplitude can be expressed using the finite
quantities:

$\sum_{i}\int_{-1}^{1}\mathrm{GPD}^{i}_{R}(x,\xi,t;\mu_{F})\bigg{[}\mathcal{C}^{i}_{0}(x,\xi,\dots)+\mathcal{C}^{i}_{1}(x,\xi,\dots)+\log\Big{(}\frac{|Q^{2}|}{\mu_{F}^{2}}\Big{)}\mathcal{C}^{i}_{coll.}(x,\xi,\dots)\bigg{]}.$
(2.25)

Dependence on the scale $\mu_{R}^{2}$ vanished due to the simple relation

$\log\Big{(}\frac{|Q^{2}|}{\mu_{R}^{2}}\Big{)}-\log\Big{(}\frac{\mu_{F}^{2}}{\mu_{R}^{2}}\Big{)}=\log\Big{(}\frac{|Q^{2}|}{\mu_{F}^{2}}\Big{)}.$

Generalization to higher QCD orders is straightforward, and can be described
in the following points:

  1. 1.

Using the renormalized fields and couplings, compute all diagrams in the hard
part of the amplitude within a given regularization scheme. Renormalization
makes them only UV-finite, but the still possess infrared/collinear
divergences.

  2. 2.

Using the same regularization, express the bare parton densities in terms of
renormalized ones. In particular, the divergent part of it must be expressed
as a product of an infinite constant and a term linear in renormalized GPDs.

  3. 3.

Verify cancellation of infrared/collinear divergences and write the amplitude
in terms of finite quantities.

####  2.2.1 Gauge link

It turns out, that at higher orders in perturbative expansion the leading (in
the hard scale) contributions stem not only from terms, in which there are
just two field operators between the hadronic states. Equally important are
graphs with an arbitrary number of gluons leaving the hard diagram and
entering the nonperturbative part, presented graphically in Fig. (2.7). In
this part we argue, that these diagrams result in presence of the gauge link
in operators defining GPDs.

$\dots$$q$$q^{\prime}$$N(\mathbf{p},\sigma_{1})$$N(\mathbf{p}^{\prime},\sigma_{2})$$k_{1}$$k_{2}$Nonperturbative
partHard part Figure 2.7: At higher orders one needs to include contributions,
in which apart from two “main” partons (denoted by straight lines connecting
the hard and nonperturbative parts) there is an arbitrary number of gluons
leaving the hard part and entering the nonperturbative matrix element.

To simplify these contributions, let us work on the matrix element from Eq.
(2.7) in the coordinate space. Let $\mathcal{T}^{(i)}(z_{1},z_{2})$ denote a
hard diagram with one parton $i$ incoming at the position $z_{1}$ and one (of
the same kind) outgoing at $z_{2}$ (we do not explicitly write free spinor or
vector indices that get contracted with their counterparts in the
nonperturbative part). By
$\mathcal{T}_{n}^{(i)\mu_{1}\dots\mu_{n}}(z_{1},z_{2};y_{1},\dots,y_{n})$ we
denote the sum of all diagrams with $n$ outgoing gluons obtained by attaching
gluon fields $\mathbb{A}_{\mu_{i}}(y_{i})$ at positions $y_{1},\dots,y_{n}$ on
propagators lines inside of the diagram $\mathcal{T}^{(i)}(z_{1},z_{2})$. To
make this statement more clear, let us present an example: let
$\mathcal{T}(z_{1},z_{2})$ correspond to a diagram presented in Fig. 2.8.

Figure 2.8: The exemplary diagram $\mathcal{T}(z_{1},z_{2})$.

Then, $\mathcal{T}_{1}^{\mu}(z_{1},z_{2};y_{1})$ is a sum of diagrams shown in
Fig. 2.9.

$+$$\mathbb{A}_{\mu}(y_{1})$

$\mathbb{A}_{\mu}(y_{1})$$+\quad\mathrm{other}\>\mathrm{attachments.}$

Figure 2.9: The sum $\mathcal{T}_{1}^{\mu}(z_{1},z_{2};y)$ defined by
$\mathcal{T}(z_{1},z_{2})$ in Fig. 2.8.

For a given kind of parton $i$, summation over all possible hard part diagrams
involving this parton and having an arbitrary number of outgoing gluon lines
is equivalent to summation over all possible diagrams $\mathcal{T}^{(i)}$ and
the corresponding $\mathcal{T}^{(i)}_{n}$. Inside of the time ordered product
of operators in the matrix element (2.7) one finds the expression

$\sum_{\mathcal{T}^{(i)}}\bigg{(}\mathcal{T}^{(i)}(z_{1},z_{2})+\sum_{n=1}^{\infty}\Big{(}\prod_{j=1}^{n}\int
d^{4}y_{j}\mathbb{A}_{\mu_{j}}(y_{j})\Big{)}\mathcal{T}_{n}^{(i)\mu_{1}\dots\mu_{n}}(z_{1},z_{2};y_{1},\dots,y_{n})\bigg{)}.$
(2.26)

Note, that in the analysis presented in Section 2.1 only the first term from
(2.26) was present. The reasoning shown in the work of Efremov and Radyushkin
[32] allows to prove, that for any diagram $\mathcal{T}^{(i)}$, the sum in
(2.26) can be written as

$\eqref{eq-sum-
T-y}=\mathcal{T}^{(i)}(z_{1},z_{2})P\>exp\Big{(}ig\int_{z_{2}}^{z_{1}}dz_{\mu}\mathbb{A}^{\mu}(z)\Big{)}+\quad\mathrm{remainder},$
(2.27)

where the remainder consists of operators which matrix elements yield sub-
leading terms.

To show it, one has to analyze how attachments of gluon vertices modify
propagators of the particles. Consider the sum of quark propagators with
additional attachments of gluons, as it is shown in Fig. 2.2.1.

$+$

$+$

$+\>\dots$

Sum of such diagrams can be written as

$\mathcal{S}(x_{1}-x_{2})=S(x_{1}-x_{2})+g\int
d^{4}y_{1}S(x_{1}-y_{1})\not{\mathbb{A}}(y_{1})S(y_{1}-x_{2})+\dots$ (2.28)

Observe that for a massless quark the free propagator $S$ is the solution of
the equation

$i\not{\partial}S(x_{1}-x_{2})=-\delta^{4}(x_{1}-x_{2}),$ (2.29)

while $\mathcal{S}$ is the propagator of the quark in the background gauge
field, and

$i\not{D}\mathcal{S}(x_{1}-x_{2})=-\delta^{4}(x_{1}-x_{2}),$ (2.30)

where $D_{\mu}=\partial_{\mu}+ig\mathbb{A}_{\mu}$ is the covariant derivative.
It can be shown [32] that

$\mathcal{S}(x_{1}-x_{2})=Pexp\Big{(}ig\int_{x_{2}}^{x_{1}}dx_{\mu}\mathbb{A}^{\mu}(x)\Big{)}\big{(}S(x_{1}-x_{2})+R(x_{1},x_{2})\big{)},$
(2.31)

where $P\>exp\Big{(}ig\int_{x_{2}}^{x_{1}}dx_{\mu}\mathbb{A}^{\mu}(x)\Big{)}$
is the path-ordered Wilson line (also know as the gauge link). The “remainder”
$R(x_{1},x_{2})$ depends on the gauge field by the field strength tensor only,
which gives a sub-leading terms when inserted between the hadron states. The
analogous reasoning applied to gluon and ghost propagators yields the same
result. Hence, summation over all possible attachments of external gluon lines
is (to the leading order in the hard scale) equivalent to substituting of
propagators inside of the diagram $\mathcal{T}^{(i)}(z_{1},z_{2})$ with ones
multiplied by the gauge link. That, on the other hand can be shown to give the
overall factor

$P\>exp\Big{(}ig\int_{z_{2}}^{z_{1}}dz_{\mu}\mathbb{A}^{\mu}(z)\Big{)}$

standing by the amplitude $T^{(i)}(z_{1},z_{2})$. It is crucial that it does
not depend on the parton type $i$ or the exact form of the diagram
$\mathcal{T}^{(i)}$, and hence can included inside the nonperturbative part
without any further complication. Therefore, after using the collinear
approximation and applying the LSZ reduction formula to the vacuum Green’s
function of the product of parton fields, hadron operators $O$ and the gauge
link one finally obtains the full definition of Generalised Parton
Distributions shown in Formulas (1.5)-(1.8).

####  2.2.2 UV and soft collinear divergences

Here we briefly describe, how divergences arise in all the considered
amplitudes and how to classify them. It allows to clarify, why collinear
divergences are present in our calculations, and say in which diagram they
occur.

A general form of a momentum integral encountered in a diagram with $L$ loops,
which contain $N$ propagators (all of massless particles), can be written in
the following form:

$G\big{(}\\{p_{i}\\}\big{)}=\prod_{l=1}^{L}\int\frac{d^{d}k_{l}}{(2\pi)^{d}}\frac{P\big{(}\\{k_{i}\\}\big{)}}{\prod_{n=1}^{N}\Big{(}Q_{n}^{2}\big{(}\\{k_{i}\\}\big{)}+i0\Big{)}}~{}\\!,$
(2.32)

where $\\{p_{i}\\}$ denotes the set of external momenta entering the loops,
$Q_{n}^{2}\big{(}\\{k_{i}\\}\big{)}$ is the momentum squared in $n$-th
propagator, which depends on momenta over which we integrate, and
$P\big{(}\\{k_{i}\\}\big{)}$ denotes momenta-dependent polynomial. There are
two classes of regions, which may lead to divergences. The first one,
corresponding to the limit $k^{2}\rightarrow\infty$ is the well-known
ultraviolet divergence, which is fixed by the hamiltonian renormalization. The
other source of divergence are poles of the integrand corresponding to
$Q^{2}_{n}=0$ (particle on the mass shell). They can be studied using the
method developed by Landau [34], which we describe now. Using the Feynman
parametrization, one can write formula (2.32) as

$G\big{(}\\{p_{i}\\}\big{)}=\prod_{j=1}^{N}\int_{0}^{1}d\alpha_{j}\delta\Big{(}1-\sum_{i=1}^{N}\alpha_{i}\Big{)}\prod_{l=1}^{L}\frac{d^{d}k_{l}}{(2\pi)^{d}}\frac{P\big{(}\\{k_{i}\\},\\{p_{i}\\}\big{)}}{\Big{[}\sum_{n=1}^{N}\alpha_{n}Q_{n}^{2}\big{(}\\{k_{i}\\}\big{)}+i0\Big{]}^{N}}.$
(2.33)

Let us denote the denominator as $D\big{(}\\{k_{i}\\}\big{)}+i0$. The
integrand is an analytical function of $\\{k_{i}\\}$ and $\\{\alpha_{i}\\}$,
so that one can modify the contour of the integration in a way that it avoids
singular points. For example, assume that, for some fixed momenta, $D$ becomes
$0$ and there is some $\alpha_{j}$ such that at this problematic point
$Q_{j}^{2}\neq 0$. Then, if $\alpha_{j}\neq 0$ (it is not at the endpoint),
one can shift the contour to avoid the singularity, see Fig. 2.10

$\mathrm{Re}\>\alpha_{j}$ Figure 2.10: An example of how one can change the
contour avoiding integration in a neighborhood of the pole, which is marked by
the cross slightly above the real axis.

Modification of the integration contour avoiding the pole is impossible, if
that occurs at one of the endpoints ($0$ or $1$) or the term $Q_{j}^{2}$,
which multiplies $\alpha_{j}$ is zero. In the case $\alpha_{j}=1$, all other
Feynman parameters are equal to $0$. It means that, if $D=0$, then
$Q_{j}^{2}=0$. Hence, modification of $\alpha_{j}$ does not allow to avoid the
singular point if

$\alpha_{j}=0\quad\mathrm{or}\quad
Q_{j}^{2}=0\quad\mathrm{for}\>\mathrm{all}\>j.$ (2.34)

If that is the case, one can try to avoid the singularity by modifying the
contour of integration over some component of the loop momentum $k_{l}^{\mu}$.
Since the domain of the integration is the whole real axis, then there are no
problems with the endpoints. However, there still can be a situation, when it
is impossible to avoid the singularity – because $D$ is a quadratic function
of momenta $\\{k_{i}\\}$, then for each considered component of momentum there
are 2 corresponding roots of $D$. If their real values are equal (and their
corresponding imaginary parts are separated only by an infinitesimally small
value $i0$), then the line of integration lays between these two poles and one
cannot change the contour, see Fig. 2.11. One can say, that these two poles
’pinch’ the line of the integration – hence the name ’pinch singularities’.

$\mathrm{Re}\>k_{l}^{\mu}$

$\mathrm{Re}\>k_{l}^{\mu}$

Figure 2.11: The illustration of a pinch singularity: in the case presented on
the left hand side, one can modify the contour avoiding poles, while in the
case on the right hand side it is impossible, because the initial contour is
’pinched’ between two poles.

The pinch singularity corresponds to the situation, when all components of
considered loop momenta $k^{\mu}_{l}$ are double roots of $D$, which
translates into the condition

$\frac{\partial}{\partial k_{l}^{\mu}}D=0.$ (2.35)

If we assume that all loop momenta come with the positive sign into propagator
momenta $Q$, we obtain that Eq. (2.35) is equivalent to

$\sum_{n=0}^{N}\alpha_{n}Q_{n}^{\mu}=0.$ (2.36)

Equations (2.34) and (2.36) are known as Landau equations. We use them to
analyze, what kind of divergences are present in the 1-loop hard diagrams in
DVCS.

$p$$k$

Let us start from the simplest case, namely the self-energy correction, where
the momentum of the quark $p$ is off-shell. In dimension $d\geq 4$ this
diagram naturally has the UV divergence, but we want to focus only on the on-
shell singularities. The corresponding denominator can be written using the
Feynman parameters as

$D=\alpha_{1}k^{2}+\alpha_{2}(p+k)^{2}.$ (2.37)

Using the Landau equations we find, that there are two singular points, which
cannot be avoided:

$\displaystyle\alpha_{1}=0,\quad p^{\mu}+k^{\mu}=0,$ (2.38)
$\displaystyle\alpha_{2}=0,\quad k^{\mu}=0.$

However, note that this is a singularity at a single point, therefore, for a
number of dimensions $d$ greater than $2$ it is integrable and there are no
resulting infrared divergences. Hence, the considered self-energy diagram has
only UV divergence for $d\geq 4$, is finite for $d=3$ and develops the
infrared divergence for $d=2$.

Now let us consider a vertex diagram with the incoming quark line on-shell.
Assume that $p=p^{+}n_{+}^{\mu}$.

$p$$k$$q$

Landau equation again tell us that there are pinch singularities when $k$,
$p+k$ or $p+q+k$ are zero, but in dimension $4$ it does not lead to IR
divergences. However, there is another solution which corresponds to the
situation, when the gluon’s momentum is collinear with that of the incoming
quark: $k^{\mu}=k^{+}n_{+}^{\mu}$. In such situation the pinch singularity
occurs at

$\alpha_{1}k^{+}+\alpha_{2}(p^{+}+k^{+})=0,\quad\alpha_{3}=0,\quad\iff\quad\alpha_{2}=-\frac{k^{+}}{p^{+}}.$
(2.39)

To obtain the second equation we used the fact that
$\alpha_{1}+\alpha_{2}+\alpha_{3}=1$. That leads to so-called collinear
divergences, which correspond to the situation when two or more momenta of
massless particles become collinear and on-shell. Note, that if the momentum
of the incoming quark was off-shell, then the vertex diagram would have only
the UV divergence. Collinear divergence is present also in the case of the box
diagram (i.e. the graph $(c)$ in Fig. 2.6). As it was shown in the previous
part, the key point of factorization theorem is proving that collinear
divergences from the hard part can be absorbed in the divergences present in
the non-renormalized parton distributions.

###  2.3 Renormalization of non-local product of operators

As it was discussed at the beginning of Section 2.2, to gain control over
divergences in parton distributions, one need to find renormalized string
operators used in the definition of GPDs. In this Section we describe the
procedure of construction of these operators. For further convenience, let use
denote string operators according to the following formula:

$\displaystyle\mathcal{O}^{(f)}(z_{1},z_{2}):=\bar{\psi}_{(f)}(z_{1})\gamma^{+}[z_{1},z_{2}]\psi_{(f)}(z_{2}),$
(2.40)
$\displaystyle\tilde{\mathcal{O}}^{(f)}(z_{1},z_{2}):=\bar{\psi}_{(f)}(z_{1})\gamma^{+}\gamma_{5}[z_{1},z_{2}]\psi_{(f)}(z_{2}),$
(2.41)
$\displaystyle\mathcal{O}^{gg}(z_{1},z_{2}):=G^{+\mu}(z_{1})[z_{1},z_{2}]G^{\phantom{a}+}_{\mu}(z_{2}),$
(2.42)
$\displaystyle\tilde{\mathcal{O}}^{gg}(z_{1},z_{2}):=G^{+\mu}(z_{1})[z_{1},z_{2}]\tilde{G}^{\phantom{a}+}_{\mu}(z_{2}).$
(2.43)

Our goal is to express them in terms of renormalized operators minus counter-
terms. To do so, we follow the reasoning shown in Appendix G in [35]. It
allows to write bare GPDs in terms of renormalized parton distributions plus
divergent terms.

Let us start with description of methods used in renormalization of local
operators. For example will serve the operator $\phi^{3}(0)$ in the scalar
theory $\phi^{4}$ in dimension 4. We demand that all vacuum Green’s function
of this operator be finite to a given order in perturbation theory. At the
first order in perturbation theory, the following two Green’s functions become
divergent:

$\displaystyle G^{1}(p):=\int
d^{4}xe^{-ipx}\bra{\Omega_{\mathrm{out}}}T\phi^{3}_{H}(0)\phi_{H}(x)\ket{\Omega_{\mathrm{in}}},$
(2.44) $\displaystyle G^{3}(p_{1},p_{2},p_{3}):=\int
d^{4}x_{1}\>d^{4}x_{2}\>d^{4}x_{3}\>e^{-ip_{1}x_{1}-ip_{2}x_{2}-ip_{3}x_{3}}\bra{\Omega_{\mathrm{out}}}T\phi^{3}_{H}(0)\phi_{H}(x_{1})\phi_{H}(x_{2})\phi_{H}(x_{3})\ket{\Omega_{\mathrm{in}}},$
(2.45)

Using the Gell-Mann and Low theorem [33] one can perform computation using
standard Feynman rules. The corresponding diagrams leading to divergences are
shown in Fig. 2.12.

Figure 2.12: Diagrams contributing to Green’s function. The gray dots denote
insertion of the operator $\phi^{3}$. The first diagram from the left
represents $G^{1}(p)$ and two other: $G^{3}(p_{1},p_{2},p_{3})$ – note that
for each of them there are 3 possible permutations of insertions of momenta.
The divergence from the self-energy correction in the third diagram is
cancelled by the field-strength renormalization constant. Divergences present
in the rest of diagrams are not cured by the lagrangian renormalization, so
that additional counter-terms are required.

The renormalized operator $\phi^{3}(0)$ at the arbitrary order needs to be
corrected by counter-terms proportional to renormalized operators $\phi_{R}$
and $\phi^{3}_{R}$:

$\big{[}\phi^{3}(0)\big{]}_{R}=\phi^{3}_{R}(0)+c_{1}\phi_{R}(0)+c_{3}\phi^{3}_{R}(0),$
(2.46)

where $c_{1/3}$ are the renormalization constants. Equivalently, one can write
it in terms of the bare operators $\phi_{0}$, using
$\phi_{R}=\mathcal{Z}_{\phi}^{-1/2}\phi_{0}$, where $\mathcal{Z}_{\phi}$ is
the renormalization constant.

In the more general case of two operators separated by a lightlike vector
$z_{1}-z_{2}$, one has to use a more complicated form of counter-terms, which
become an product of operators integrated over positions:

$\mathrm{counterterm}\sim\int_{0}^{1}d\alpha_{1}d\alpha_{2}K(\alpha_{1},\alpha_{2})\mathcal{O}\big{(}\alpha_{1}z_{1}+(1-\alpha_{1})z_{2},\alpha_{2}z_{2}+(1-\alpha_{2})z_{1}\big{)}.$
(2.47)

Distributions $K(\alpha_{1},\alpha_{2})$ are called the evolution kernels.
It turns out, that the Green’s functions of fields of different kinds do not
vanish, hence, for example, the renormalized quark operator consists of gluon
fields. Because of that, it is convenient to consider combinations of
operators defined in Eqs (2.40)-(2.43) of definite charge-parity. The non-
singlet operators, which are $\mathcal{C}$-odd and do not mix under
renormalization, are the following:

$\displaystyle\mathcal{O}^{(f)}_{NS}(z_{1},z_{2})=\mathcal{O}^{(f)}(z_{1},z_{2})+\mathcal{O}^{(f)}(z_{2},z_{1}),$
(2.48)
$\displaystyle\tilde{\mathcal{O}}^{(f)}_{NS}(z_{1},z_{2})=\tilde{\mathcal{O}}^{(f)}(z_{1},z_{2})-\tilde{\mathcal{O}}^{(f)}(z_{2},z_{1}).$
(2.49)

$\mathcal{C}$-even (singlet) combinations are:

$\displaystyle\mathcal{O}_{S}^{(f)}(z_{1},z_{2})=\mathcal{O}^{(f)}(z_{1},z_{2})-\mathcal{O}^{(f)}(z_{2},z_{1})$
(2.50)
$\displaystyle\mathcal{O}_{S}^{gg}(z_{1},z_{2})=\mathcal{O}^{gg}(z_{1},z_{2}),$
(2.51)
$\displaystyle\tilde{\mathcal{O}}_{S}^{(f)}(z_{1},z_{2})=\tilde{\mathcal{O}}^{(f)}(z_{1},z_{2})+\tilde{\mathcal{O}}^{(f)}(z_{2},z_{1})$
(2.52)
$\displaystyle\tilde{\mathcal{O}}_{S}^{gg}(z_{1},z_{2})=\tilde{\mathcal{O}}^{gg}(z_{1},z_{2}),$
(2.53)

(it is worth to stress, that gluon operators are always charge-even).

Similarly as in the simpler case of local operator, one needs to find counter-
terms necessary to make the vacuum Green’s functions of the time-ordered
product of bare fields and operators
$\bar{\psi}(z_{1})[z_{1},z_{2}]\psi(z_{2})$,
$G^{+\mu}(z_{1})[z_{1},z_{2}]G^{\phantom{a}+}_{\mu}(z_{2})$ finite. In the
case of the quark string operator one considers two kinds of such Green’s
functions:

$\displaystyle G^{ff^{\prime}}(p_{1},p_{2},z_{1},z_{2}):=\int d^{4}x_{1}\int
d^{4}x_{2}e^{-ip_{1}x_{1}-ip_{2}x_{2}}\bra{\Omega_{\mathrm{out}}}T\mathcal{O}^{(f)}(z_{1},z_{2})\bar{\psi}_{(f^{\prime})}(x_{1})\psi_{(f^{\prime})}(x_{2})\big{\\}}\ket{\Omega},$
(2.54) $\displaystyle G^{fg}(p_{1},p_{2},z_{1},z_{2}):=\int d^{4}x_{1}\int
d^{4}x_{2}e^{-ip_{1}x_{1}-ip_{2}x_{2}}\bra{\Omega_{\mathrm{out}}}T\mathcal{O}^{(f)}(z_{1},z_{2})G^{+\mu}(z_{1})G^{\phantom{a}+}_{\mu}(x_{2})\ket{\Omega},$

in the limit when $p_{1,2}$ have only the ”$+$” component of momentum (that is
due to the fact, that we assume partons taking part in the hard process to be
massless and nearly collinear with the hadron). The function $G^{fg}$
corresponds to the mixing of quark $f$ and gluon operators and is present only
in the case of the singlet combination. At the $\alpha_{S}$ order, Green’s
functions with quark operators do not vanish only as $f=f^{\prime}$, however,
at higher orders there are also contributions responsible for mixing of quark
flavors in the singlet sector, as it is presented in Fig. 2.13.

Figure 2.13: Diagrams appearing in computation of Green’s functions of string
operators at higher orders, which contribute to the mixing of the operators.
The diagram in the middle corresponds to mixing of different quark flavors,
while the other two result in mixing of quarks and gluons.

For simplicity, we consider amputated Green’s functions, which are ordinary
Green’s functions, but with removed external propagators and spinors (or
polarization vectors). For example, in the leading order in the strong
coupling constant, the only non-vanishing Green’s function can be represented
by the following diagram:

$\leavevmode\hbox to91.72pt{\vbox
to22.67pt{\pgfpicture\makeatletter\hbox{\hskip
131.21968pt\lower-39.78577pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\feynman \vertex(z1) at (0, 0)
{$z_{1}$}; \vertex(z2) at (3, 0) {$z_{2}$}; \vertex(p1) at (0, -2);
\vertex(p2) at (3, -2); \diagram* { (p1) -- [fermion, insertion=1.] (z1); (p2)
-- [anti fermion, insertion = 1.] (z2); (p1) -- [fermion, momentum'=$p_{1}$]
(z1); (p2) -- [anti fermion, momentum'=$p_{2}$] (z2); };
{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-127.88667pt}{-32.45276pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{$G^{ff}_{amp.}(p_{1},p_{2},z_{1},z_{2})\Big{|}_{\alpha_{S}^{0}}=$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys<EMAIL_ADDRESS>(2.55)

Crosses labeled $z_{1}$ and $z_{2}$ denote insertions of operators
$\bar{\psi}(z_{1})$ and $\psi(z_{2})$. It is easy to find, that

$\eqref{eq-green-lo}=\gamma^{+}e^{-iz_{1}^{-}p_{1}^{+}-iz_{2}^{-}p_{2}^{+}}.$
(2.56)

At the 1-loop order, when computing $G^{ff}_{amp.}(p_{1},p_{2},z_{1},z_{2})$,
one has to consider diagrams presented in Fig. 2.14.

$z_{1}$$z_{2}$$\xi$$p_{1}$$k$$p_{2}$$(a)$

$z_{1}$$z_{2}$$\xi$$k$$p_{1}$$p_{2}$$(b)$

$z_{1}$$z_{2}$$k$$p_{1}$$p_{2}$$(c)$

$z_{1}$$z_{2}$$p_{1}$$p_{2}$$(d)$

$z_{1}$$z_{2}$$p_{1}$$p_{2}$$(e)$

Figure 2.14: Diagrams present at the 1-loop computation of
$G^{ff}_{amp.}(p_{1},p_{2},z_{1},z_{2})$. Diagrams $(a)$ and $(b)$ contain
insertion of a gluon field coming from expanding the gauge link up to the
order $g$ (it is denoted by the cross labelled by $\xi$). Diagram $(c)$
corresponds to gluon exchange between two on-shell quarks. Diagrams $(d)$ and
$(e)$ are the standard self-energy corrections. As it was said at the
beginning of Section 2.2, at the one-loop level one can use the unrenormalized
perturbation theory, so that here we do not include diagrams with the self-
energy counter-terms.

Let us present a detailed computation of the diagram $(a)$. It allows to show
methods used in computation of considered Green’s functions and explain, how
one obtains the necessary counter-terms. In the Sudakov basis, which is the
most convenient for analysis of operators on the lightcone, usually one has to
treat ”$\pm$” and ”$\perp$” components independently. It turns out, that the
simplest way to achieve it, is to use the Schwinger parametrization of
propagator:

$\frac{i}{k^{2}+i0}=\int_{0}^{\infty}e^{iy(k^{2}+i0)}\>dy.$ (2.57)

When more propagators are present, it is useful to perform the following
substitution in the integration:

$\prod_{j=1}^{N}\int_{0}^{\infty}dy_{j}\rightarrow\int_{0}^{\infty}d\rho\rho^{N-1}\prod_{j=1}^{N}\int_{0}^{1}d\alpha_{j}\delta\Big{(}\sum_{j=1}^{N}\alpha_{j}-1\Big{)}.$
(2.58)

Because of divergences, we use the dimensional regularization. For gluon field
we use the Feynman gauge. Averaging over colors we get that the contribution
to the Green’s function from the considered diagram, which we denote by
$\mathcal{M}^{ff}_{a}$, is

$\mathcal{M}^{ff}_{a}=g^{2}C_{F}\int_{z_{2}^{-}}^{z_{1}^{-}}d\xi^{-}\mu_{R}^{4-d}\int\frac{d^{d}k}{(2\pi)^{d}}\gamma^{+}e^{-iz_{1}^{-}(p_{1}^{+}-k^{+})-iz_{2}^{-}p_{2}^{+}}e^{-i\xi^{-}k^{+}}\frac{\not{p}_{1}-\not{k}}{(p_{1}-k)^{2}+i0}\gamma^{+}\frac{1}{k^{2}+i0}~{}\\!,$
(2.59)

where $C_{F}=\frac{4}{3}$ is the color factor and $\mu_{R}$ is the
renormalization scale added to keep the proper dimension of the amplitude.
Phases are the result of insertion of fields with given momenta at points
$z_{1}$, $z_{2}$ and $\xi$ – the parameter $\xi$ is the position, at which one
inserts the gluon field from the gauge link. Note that, since
$(\gamma^{+})^{2}=0$, we can write

$\gamma^{+}\big{(}\not{p}_{1}-\not{k}\big{)}\gamma^{+}=2\big{(}p_{1}^{+}-k^{+}\big{)}\gamma^{+}~{}\\!,$
(2.60)

so that the gamma matrix structure of (2.59) is the same as that of (2.55).
Performing the integral over $d\xi^{-}$ we get

$\eqref{eq-nlo-
qq-b}=2ig^{2}C_{F}\gamma^{+}\mu_{R}^{4-d}\int\frac{d^{d}k}{(2\pi)^{d}}\frac{p_{1}^{+}-k^{+}}{k^{2}(p_{1}-k)^{2}}e^{-iz_{1}^{-}(p_{1}^{+}-k^{+})-iz_{2}^{-}p_{2}^{+}}\frac{1}{k^{+}}\Big{(}e^{-iz_{1}^{-}k^{+}}-e^{-iz_{2}^{-}k^{+}}\Big{)}.$
(2.61)

Let us denote $\Delta z=z_{2}^{-}-z_{1}^{-}$ and use the following identity:

$\frac{1}{k^{+}}\Big{(}1-e^{-i\Delta zk^{+}}\Big{)}=i\Delta z\int_{0}^{1}d\tau
e^{-i\Delta z^{-}k^{+}\tau}.$ (2.62)

The resulting integrals can be done using the previously discussed Schwinger
parametrization:

$\displaystyle\int\frac{d^{d}k}{(2\pi)^{d}}\int_{0}^{1}d\tau e^{-i\Delta
zk^{+}\tau}\frac{p_{1}^{+}-k^{+}}{k^{2}(p_{1}-k)^{2}}=$ (2.63)
$\displaystyle-\int\frac{d^{d}k}{(2\pi)^{d}}\int_{0}^{\infty}d\rho\int_{0}^{1}d\tau\int_{0}^{1}d\alpha\big{(}(1-\alpha)p_{1}^{+}-k^{+}\big{)}e^{i\rho
k^{2}}e^{-i\tau\Delta z(k^{+}+\alpha p_{1}^{+})\tau}.$

Going to the second line we used the substitution $k\rightarrow k-\alpha
p_{1}$. Integration over $dk^{-}$ produces a delta distribution:

$\int\frac{dk^{-}}{2\pi}e^{2i\rho
k^{+}k^{-}}=\frac{1}{2}\rho^{-1}\delta(k^{+})~{}\\!,$ (2.64)

which allows to easily perform the integral over $dk^{+}$. Keep in mind that
$k^{2}=2k^{+}k^{-}-k_{\perp}^{2}$. Integration over perpendicular momenta is
just a standard Gaussian integral:

$\int\frac{d^{d-2}k_{\perp}}{(2\pi)^{d-2}}e^{-i\rho
k^{2}_{\perp}}=\frac{1}{(4\pi)^{d/2-1}}\big{(}-i\rho\big{)}^{-d/2+1}.$ (2.65)

The resulting integral over $d\rho$ is ill-defined:

$\frac{1}{(4\pi)^{d/2-1}}\int_{0}^{\infty}d\rho\>(-i\rho)^{-d/2+1}.$ (2.66)

To solve this problem, let us multiply the integrand by
$e^{-i\rho\mu_{F}^{2}}$, where $\mu_{F}$ is called the factorization scale. It
is equivalent with the assumption that all partons in the diagram have masses
$\mu_{F}$ [35]. This way, instead of the integral in formula (2.66), we have

$\frac{1}{(4\pi)^{d/2-1}}\int_{0}^{\infty}d\rho\>(-i\rho)^{-d/2+1}e^{-i\mu_{F}^{2}\rho}=\frac{i}{4\pi}\Big{(}\frac{\mu_{F}^{2}}{4\pi}\Big{)}^{d/2-2}\Gamma\Big{(}2-\frac{d}{2}\Big{)},$
(2.67)

where $\Gamma(z)$ is the gamma function. It is clear, that in the limit of
$d\rightarrow 4$, this expression becomes divergent. Now we perform the
integration over the parameter $d\tau$:

$p_{1}^{+}(1-\alpha)\int_{0}^{1}d\tau e^{-i\alpha p_{1}^{+}\Delta
z\tau}=\frac{i}{\Delta z}\frac{1-\alpha}{\alpha}\big{(}e^{-i\alpha
p_{1}^{+}\Delta z}-1\big{)}.$ (2.68)

All these integrals yield the following result:

$\eqref{eq-ren-
momentum}=-\frac{i}{(4\pi)^{2}}\Big{(}\frac{\mu_{F}^{2}}{4\pi}\Big{)}^{d/2-2}\Gamma\Big{(}2-\frac{d}{2}\Big{)}\frac{i}{\Delta
z}\frac{1-x}{x}\big{(}e^{-ixp_{1}^{+}\Delta z}-1\big{)}.$ (2.69)

Hence, after multiplying by the remaining constants, we obtain

$\eqref{eq-nlo-
qq-b}=\frac{1}{8\pi^{2}}C_{F}g^{2}\Big{(}\frac{\mu_{F}^{2}}{4\pi\mu_{R}^{2}}\Big{)}^{d/2-2}\Gamma\Big{(}2-\frac{d}{2}\Big{)}\gamma^{+}\int_{0}^{1}d\alpha\frac{1-\alpha}{\alpha}\Big{(}e^{-i(z_{2}^{-}-z_{1}^{-})\alpha
p_{1}^{+}}-1\Big{)}e^{-iz_{1}^{-}p_{1}^{+}+iz_{2}^{-}p_{2}^{+}}.$ (2.70)

It can be conveniently written using so-called ’plus prescription’:

$\int_{0}^{1}d\alpha\>\Big{[}\frac{1-\alpha}{\alpha}\Big{]}_{+}f(\alpha):=\int_{0}^{1}d\alpha\frac{1-\alpha}{\alpha}\big{(}f(\alpha)-f(0)\big{)}.$
(2.71)

Moreover, using the Laurent expansion of the gamma function
($d=4-\varepsilon$),

$\Gamma\Big{(}2-\frac{d}{2}\Big{)}=\Gamma\Big{(}\frac{\varepsilon}{2}\Big{)}=\frac{2}{\varepsilon}-\gamma+\mathcal{O}(\epsilon),$
(2.72)

we can write

$\Big{(}\frac{\mu_{F}^{2}}{4\pi\mu_{R}^{2}}\Big{)}^{d/2-2}\Gamma\Big{(}2-\frac{d}{2}\Big{)}=\frac{2}{\varepsilon}\Big{(}\frac{e^{\gamma}\mu_{F}^{2}}{4\pi\mu_{R}^{2}}\Big{)}^{-\varepsilon/2}.$
(2.73)

Finally, using $\alpha_{S}=g^{2}/4\pi$, we can re-express the result in a
compact form

$\eqref{eq-nlo-qq-b-pre-
final}=\frac{2}{\varepsilon}\frac{\alpha_{S}}{2\pi}C_{F}\Big{(}\frac{e^{\gamma}\mu_{F}^{2}}{4\pi\mu_{R}^{2}}\Big{)}^{-\varepsilon/2}\gamma^{+}\int_{0}^{1}d\alpha\>\Big{[}\frac{1-\alpha}{\alpha}\Big{]}_{+}e^{-ip_{1}^{+}\big{(}(1-\alpha)z_{1}^{-}+\alpha
z_{2}^{-}\big{)}-iz_{2}^{-}p_{2}^{+}}.$ (2.74)

At the first sight, this expression looks complicated and finding the counter-
term may appear a cumbersome task. However, let us notice, that Eq. (2.74) can
be written in terms of the LO amputated Green’s function, see Eq. (2.56):

$\eqref{eq-Ma-
fin}=\frac{2}{\varepsilon}\frac{\alpha_{S}}{2\pi}C_{F}\Big{(}\frac{e^{\gamma}\mu_{F}^{2}}{4\pi\mu_{R}^{2}}\Big{)}^{-\varepsilon/2}\int_{0}^{1}d\alpha\>\Big{[}\frac{1-\alpha}{\alpha}\Big{]}_{+}G^{ff}_{amp.}(p_{1},p_{2},(1-\alpha)z_{1}^{-}+\alpha
z_{2}^{-},z_{2})\Big{|}_{\alpha_{S}^{0}}.$ (2.75)

Using this result, one can easily notice, that the counter-term necessary to
cancel divergence from the diagram $(a)$ at the order $\alpha_{S}^{1}$ is

$-\frac{2}{\varepsilon}\Big{(}\frac{e^{\gamma}\mu_{F}^{2}}{4\pi\mu_{R}^{2}}\Big{)}^{-\varepsilon/2}\int_{0}^{1}d\alpha\>\frac{\alpha_{S}}{2\pi}C_{F}\Big{[}\frac{1-\alpha}{\alpha}\Big{]}_{+}\mathcal{O}^{ff}\big{(}(1-\alpha)z_{1}^{-}+\alpha
z_{2}^{-},z_{2}\big{)}.$ (2.76)

After considering all diagrams presented in Fig. 2.14 one obtains the
renormalized non-singlet operator:

$\displaystyle\mathcal{O}^{(f)}_{NS,R}(z_{1},z_{2};\mu_{F})=$
$\displaystyle\>\mathcal{O}^{(f)}_{NS}(z_{1},z_{2})\>+$ (2.77)
$\displaystyle-\frac{2}{\varepsilon}\Big{(}\frac{e^{\gamma}\mu_{F}^{2}}{4\pi\mu_{R}^{2}}\Big{)}^{-\varepsilon/2}\int_{0}^{1}d\alpha_{1}\>d\alpha_{2}\>K^{ff}_{NS}(\alpha_{1},\alpha_{2})\>\mathcal{O}^{(f)}_{NS}\big{(}\bar{\alpha}_{1}z_{1}+\alpha_{1}z_{2},\bar{\alpha}_{2}z_{2}+\alpha_{2}z_{1}\big{)}~{}\\!,$

where we denote $\bar{\alpha}_{i}=1-\alpha_{i}$ and the evolution kernel is
given by

$K^{ff}_{NS}(\alpha_{1},\alpha_{2})=\frac{\alpha_{S}}{2\pi}C_{F}\bigg{(}1+\delta(\alpha_{1})\Big{[}\frac{\bar{\alpha}_{2}}{\alpha_{2}}\Big{]}_{+}+\delta(\alpha_{2})\Big{[}\frac{\bar{\alpha}_{1}}{\alpha_{1}}\Big{]}_{+}+\frac{3}{2}\delta(\alpha_{1})\delta(\alpha_{2})\bigg{)}.$
(2.78)

Other evolution kernels in the position space can be found in [27, 35]. The
relation between renormalized and bare operators can be inverted:

$\displaystyle\mathcal{O}^{(f)}_{NS}(z_{1},z_{2})=$
$\displaystyle\>\mathcal{O}^{(f)}_{NS,R}(z_{1},z_{2};\mu_{F})\>+$ (2.79)
$\displaystyle+\frac{2}{\varepsilon}\Big{(}\frac{e^{\gamma}\mu_{F}^{2}}{4\pi\mu_{R}^{2}}\Big{)}^{-\varepsilon/2}\int_{0}^{1}d\alpha_{1}\>d\alpha_{2}\>K^{ff}_{NS}(\alpha_{1},\alpha_{2})\>\mathcal{O}^{(f)}_{NS,R}\big{(}\bar{\alpha}_{1}z_{1}+\alpha_{1}z_{2},\bar{\alpha}_{2}z_{2}+\alpha_{2}z_{1};\mu_{F}\big{)}~{}\\!,$

where we used the fact that, at the order $\alpha_{S}$, in the counter-term we
can use either renormalized operator, or the bare one (the difference is of
order $\alpha_{S}^{2}$). This result allows to write bare GPDs in terms of
renormalized parton distributions and their convolutions with the
corresponding evolution kernels.

###  2.4 From evolution of operators in the position space to evolution of
GPDs

In this section we finally obtain the relation between bare and renormalized
parton distributions and formulate their evolution equations. We define
renormalized GPDs in the same way as bare one, but using renormalized string
operator. For example, the renormalized non-singlet vector quark distribution
is

$F^{f}_{R}(x,\xi,t;\mu_{F}):=\int\frac{dz^{-}}{2\pi}e^{ixp^{+}z^{-}}\bra{N(\mathbf{p}_{2},\sigma_{2})_{\mathrm{out}}}\mathcal{O}^{(f),H}_{NS,R}(-z^{-}/2,z^{-}/2;\mu_{F})\ket{N(\mathbf{p}_{1},\sigma_{1})_{\mathrm{in}}}.$
(2.80)

From the Formula (2.79) we obtain the difference between bare and renormalized
GPD:

$\displaystyle F^{f}_{NS,R}(x,\xi,t;\mu_{F})-F^{f}_{NS}(x,\xi,t)=$ (2.81)
$\displaystyle=$
$\displaystyle-\frac{2}{\varepsilon}\Big{(}\frac{e^{\gamma}\mu_{F}^{2}}{4\pi\mu_{R}^{2}}\Big{)}^{-\varepsilon/2}\int\frac{dz^{-}}{2\pi}e^{ixp^{+}z^{-}}\int_{0}^{1}d\alpha_{1}\>d\alpha_{2}\>K^{ff}_{NS}(\alpha_{1},\alpha_{2})\times$
$\displaystyle\bra{N(\mathbf{p}_{2},\sigma_{2})_{\mathrm{out}}}\mathcal{O}^{(f),H}_{NS,R}\big{(}-z^{-}(1/2-\alpha_{1}),z^{-}(1/2-\alpha_{2});\mu_{F}\big{)}\ket{N(\mathbf{p}_{1},\sigma_{1})_{\mathrm{in}}}.$

To write it as a parton distribution convoluted with some kernel, let us first
use the Poincaré covariance of the matrix element:

$\displaystyle\bra{N(\mathbf{p}_{2},\sigma_{2})_{\mathrm{out}}}\mathcal{O}^{(f)}_{NS,R}\big{(}-z^{-}(1/2-\alpha_{1}),z^{-}(1/2-\alpha_{2})\big{)}\ket{N(\mathbf{p}_{1},\sigma_{1})_{\mathrm{in}}}=$
(2.82) $\displaystyle=$ $\displaystyle
e^{-ip^{+}z^{-}\xi(\alpha_{1}-\alpha_{2})}\bra{N(\mathbf{p}_{2},\sigma_{2})_{\mathrm{out}}}\mathcal{O}^{(f)}_{NS,R}\big{(}-z^{-}(1-\alpha_{1}-\alpha_{2})/2,z^{-}(1-\alpha_{1}-\alpha_{2})/2;\mu_{F}\big{)}\ket{N(\mathbf{p}_{1},\sigma_{1})_{\mathrm{in}}}.$

Using the substitution $z^{\prime}=z(1-\alpha_{1}-\alpha_{2})$ we rewrite the
integrals from the Eq. (2.81) in the following form:

$\int\frac{dz^{\prime-}}{2\pi}e^{ip^{+}z^{\prime-}\frac{x-\xi(\alpha_{1}-\alpha_{2})}{1-\alpha_{1}-\alpha_{2}}}\int_{0}^{1}d\alpha_{1}\>d\alpha_{2}\>K^{ff}_{NS}(\alpha_{1},\alpha_{2})\bra{N(\mathbf{p}_{2},\sigma_{2})_{\mathrm{out}}}\mathcal{O}^{(f)}_{NS,R}\big{(}-z^{\prime-}/2),z^{\prime-}/2;\mu_{F}\big{)}\ket{N(\mathbf{p}_{1},\sigma_{1})_{\mathrm{in}}}.$
(2.83)

It is easy to verify that it can be written as

$\displaystyle\int_{-1}^{1}dx^{\prime}\>\int_{0}^{1}d\alpha_{1}d\alpha_{2}\>\delta\Big{(}x-x^{\prime}(1-\alpha_{1}-\alpha_{2})-\xi(\alpha_{1}-\alpha_{2})\Big{)}K^{ff}_{NS}(\alpha_{1},\alpha_{2})F^{f}_{NS,R}(x^{\prime},\xi,t)=$
(2.84) $\displaystyle=$
$\displaystyle\int_{-1}^{1}dx^{\prime}\>\frac{1}{|\xi|}V_{NS}\Big{(}\frac{x}{\xi},\frac{x^{\prime}}{\xi}\Big{)}F_{NS,R}^{f}(x^{\prime},\xi,t;\mu_{F}),$

where the evolution kernel for the non-singlet GPD is defined as

$\frac{1}{|\xi|}V_{NS}\Big{(}\frac{x}{\xi},\frac{x^{\prime}}{\xi}\Big{)}=\int_{0}^{1}d\alpha_{1}d\alpha_{2}\>\delta\Big{(}x-x^{\prime}(1-\alpha_{1}-\alpha_{2})-\xi(\alpha_{1}-\alpha_{2})\Big{)}K^{ff}_{NS}(\alpha_{1},\alpha_{2}).$
(2.85)

The resulting formula takes the same form as the one in Eq. (2.20):

$F_{NS}^{f}(x,\xi,t)=F_{NS,R}^{f}(x,\xi,t;\mu_{F})+\frac{2}{\varepsilon}\Big{(}\frac{e^{\gamma}\mu_{F}^{2}}{4\pi\mu_{R}^{2}}\Big{)}^{-\varepsilon/2}\int_{-1}^{1}dx^{\prime}\>\frac{1}{|\xi|}V_{NS}\Big{(}\frac{x}{\xi},\frac{x^{\prime}}{\xi}\Big{)}F_{NS,R}^{f}(x^{\prime},\xi,t;\mu_{F}).$
(2.86)

Demanding that the bare operators are independent of the factorization scale,
one obtains the evolution equation:

$\mu_{F}^{2}\frac{d}{d\mu_{F}^{2}}F_{NS,R}^{f}(x,\xi,t;\mu_{F})=\int_{-1}^{1}dx^{\prime}\>\frac{1}{|\xi|}V_{NS}\Big{(}\frac{x}{\xi},\frac{x^{\prime}}{\xi}\Big{)}F_{NS,R}^{f}(x^{\prime},\xi,t;\mu_{F}).$
(2.87)

In the case of singlet GPDs all these relations take the same form, but with
the evolution kernels in the form of a matrix with indices corresponding to
quark and gluon distributions. All kernels can be found in [27, 35].

##  Chapter 3 Photoproduction of photon pairs

After introducing the collinear factorization, generalized parton
distributions and methods used in the computation of higher-order QCD
corrections, we present a detailed analysis of photoproduction of photon pair
with large invariant mass on a nucleon, i.e. in the process

$\gamma(\mathbf{q},\epsilon)\>N(\mathbf{p}_{1},\sigma_{1})\longrightarrow\gamma^{\prime}(\mathbf{q}_{1},\epsilon_{1}(\mathbf{q}_{1}))\>\gamma^{\prime\prime}(\mathbf{q}_{2},\epsilon(\mathbf{q}_{2}))\>N(\mathbf{p}_{2},\sigma_{2}).$

$\epsilon$ denotes polarization vectors of photons. It is the simplest
reaction, which allows to study the factorization in a $2\rightarrow 3$
process – and hence a good starting point in studies of this general class of
processes. Because of the charge conjugation symmetry, it is sensitive only to
non-singlet combinations of GPDs (in contrary to e.g. DVCS, TCS or DVMP), so
it might be useful in extraction of GPDs from experimental data. In Section
3.1 we describe kinematics of this process. In Section 3.2 we show the
leading-order results [20, 21]. The remaining part of the work are my original
results. Throughout Sections 3.3 \- 3.7 we present a detailed computation of
one loop QCD corrections to the hard-part. The full result can be found in
Appendix B. The proof of validity of the collinear factorization at NLO is
shown in Section 3.8.

###  3.1 Kinematics, amplitude, and the differential cross section

We parametrize particles’ momenta and polarizations using the following light-
cone four-vectors:

$p^{\mu}=\frac{\sqrt{s}}{2}(1,0,0,1),\quad\quad
n^{\mu}=\frac{\sqrt{s}}{2}(1,0,0,-1),$ (3.1)

the parameter $s$ will be described in more detail later. We use the Sudakov
decomposition:

$\displaystyle
p_{1}^{\mu}=(1+\xi)p^{\mu}+\frac{m_{N}^{2}}{s(1+\xi)}n^{\mu},\qquad
p_{2}^{\mu}=(1-\xi)p^{\mu}+\frac{m_{N}^{2}+\vec{\Delta}_{t}^{2}}{s(1-\xi)}n^{\mu}+\Delta_{\perp}^{\mu},\qquad
q^{\mu}=n^{\mu},$ (3.2) $\displaystyle q_{1}^{\mu}=\alpha
n^{\mu}+\frac{(\vec{p}_{t}-\frac{1}{2}\vec{\Delta}_{t})^{2}}{\alpha
s}p^{\mu}+p_{\perp}^{\mu}-\frac{1}{2}\Delta_{\perp}^{\mu},$ $\displaystyle
q_{2}^{\mu}=\bar{\alpha}n^{\mu}+\frac{(\vec{p}_{t}+\frac{1}{2}\vec{\Delta}_{t})^{2}}{\bar{\alpha}s}p^{\mu}-p_{\perp}^{\mu}-\frac{1}{2}\Delta_{\perp}^{\mu}.$

$\Delta^{\mu}=p_{2}^{\mu}-p_{1}^{\mu}$, $m_{N}$ is the mass of the nucleon,
and the parameter $\xi$ is called skewedness. We denote the Euclidean product
as $\vec{p}_{t}^{\phantom{a}2}=-p_{\perp}^{2}\geq 0$. We use the gauge
$A^{-}=0$ for the electromagnetic field. It implies that for each photon’s
polarization vector

$\epsilon\cdot p=0.$

Hence, the polarization vectors are the following:

$\displaystyle\epsilon^{\mu}(\mathbf{q}_{1})=\epsilon_{\perp}^{\mu}(\mathbf{q}_{1})+\frac{2\vec{\epsilon}_{t}(\mathbf{q}_{1})\big{(}\vec{p}_{t}-\frac{1}{2}\vec{\Delta}_{t}\big{)}}{\alpha
s}p^{\mu},$ (3.3)
$\displaystyle\epsilon^{\mu}(\mathbf{q}_{2})=\epsilon_{\perp}^{\mu}(\mathbf{q}_{1})-\frac{2\vec{\epsilon}_{t}(\mathbf{q}_{2})\big{(}\vec{p}_{t}+\frac{1}{2}\vec{\Delta}_{t}\big{)}}{\bar{\alpha}s}p^{\mu},$
$\displaystyle\epsilon^{\mu}(\mathbf{q})=\epsilon_{\perp}^{\mu}(\mathbf{q}).$

The process is described in terms of 4 Mandelstam variables:

$S_{\gamma N}=(p_{1}+q)^{2},\quad t=(p_{1}-p_{2})^{2},\quad
M_{\gamma\gamma}^{2}=(q_{1}+q_{2})^{2},\quad u^{\prime}=(q_{2}-q)^{2}.$ (3.4)

For further convenience we relate them to other invariants:

$\displaystyle S_{\gamma N}=(1+\xi)s+m_{N}^{2}\>\implies\>s=\frac{S_{\gamma
N}-m_{N}^{2}}{1+\xi},$ (3.5) $\displaystyle
t^{\prime}=(q_{1}-q)^{2}=t-M_{\gamma\gamma}^{2}-u^{\prime}.$

Like in the work [20], we consider the situation when $M_{\gamma\gamma}^{2}$
is large and $t$ is small (with respect to the mass squared of the nucleon).
It implies that transverse momenta of outgoing photons are large and
approximately opposite. Let us discuss the case $\Delta_{\perp}^{\mu}=0$. It
is justified, since in the considered kinematics one has
$|p_{\perp}^{2}|\gg|\Delta_{\perp}^{2}|$. The simplified kinematical relations
used in the calculation of the hard part are:

$\displaystyle\xi=\frac{\tau}{2-\tau},\qquad\tau=\frac{M_{\gamma\gamma}^{2}-t}{S_{\gamma
N}-m_{N}^{2}},$ (3.6)
$\displaystyle-t^{\prime}=\frac{\vec{p}_{t}^{\phantom{a}2}}{\alpha},\qquad-u^{\prime}=\frac{\vec{p}_{t}^{\phantom{a}2}}{\bar{\alpha}},\qquad\alpha+\bar{\alpha}=1.$

The second line of Eq. (LABEL:eq-kin-simplified) yields

$\vec{p}_{t}^{\phantom{a}2}=-\frac{u^{\prime}t^{\prime}}{u^{\prime}+t^{\prime}},\qquad\alpha=\frac{u^{\prime}}{u^{\prime}+t^{\prime}},\qquad\bar{\alpha}=\frac{t^{\prime}}{u^{\prime}+t^{\prime}}.$
(3.7)

Because of the charge conjugation symmetry, the amplitude of the
photoproduction of two real photons on a proton is sensitive only to non-
singlet (odd with respect to the charge conjugation) parton distributions. To
avoid proliferation of indices, from now on by

$H^{q}(x,\xi,t;\mu_{F})\quad\mathrm{and}\quad E^{q}(x,\xi,t;\mu_{F})$

we denote the renormalized, non-singlet vector GPDs corresponding to quark of
flavour $q$, see Eq. (1.11). We neglect the axial GPDs, since the LO analysis
[20] has shown, that their contribution in this process is negligible. The
scattering amplitude of the full process $\mathcal{M}$ can be written in the
following way:

$\mathcal{M}=\frac{1}{2s}\sum_{q}\Big{[}\mathcal{H}^{q}\bar{u}(\mathbf{p}_{2},\sigma_{2})\gamma^{+}u(\mathbf{p}_{1},\sigma_{1})+\mathcal{E}^{q}\bar{u}(\mathbf{p}_{2},\sigma_{2})\frac{i\sigma^{\mu+}\Delta_{\mu}}{2m_{N}}u(\mathbf{p}_{1},\sigma_{1})\Big{]},$
(3.8)

where $u(\mathbf{p},\sigma)$ denotes the bispinor of the nucleon, and

$\displaystyle\mathcal{H}^{q}=\int_{-1}^{1}dx\>H^{q}(x,\xi,t;\mu_{F})\mathcal{T}^{q}_{R}(x,\dots;\mu_{F}),$
(3.9)
$\displaystyle\mathcal{E}^{q}=\int_{-1}^{1}dx\>E^{q}(x,\xi,t;\mu_{F})\mathcal{T}^{q}_{R}(x,\dots;\mu_{F}),$

where $\mathcal{T}^{q}_{R}(x,\dots;\mu_{F})$ denotes the “renormalized”
amplitude of the hard sub-process:

$\mathcal{T}_{R}^{i}(x,\xi,\dots;\mu_{F})\equiv\mathcal{C}^{i}_{0}(x,\xi,\dots)+\mathcal{C}^{i}_{1}(x,\xi,\dots)+\log\Big{(}\frac{s}{\mu_{F}^{2}}\Big{)}\mathcal{C}^{i}_{coll.}(x,\xi,\dots),$
(3.10)

as it is described by the Formula (2.25). Ellipsis denotes all parameters
other than $x$, on which depends the hard amplitude. Here, the hard scale is
provided by $s$ instead of $|Q^{2}|$ used in Chapter 2.

The amplitude squared, summed (averaged) over the final (initial) nucleon
spins $\sigma_{1}$, $\sigma_{2}$ reads:

$\frac{1}{2}\sum_{\sigma_{1},\sigma_{2}}\big{|}\mathcal{T}\big{|}^{2}=\frac{1}{4}\Big{[}(1-\xi^{2})\big{|}\sum_{q}\mathcal{H}^{q}\big{|}^{2}-2\xi^{2}\mathrm{Re}\sum_{q}\mathcal{H}^{q}\sum_{q^{\prime}}\mathcal{E}^{q^{\prime}\>*}+\frac{\xi^{4}}{1-\xi^{2}}\big{|}\sum_{q}\mathcal{E}^{q}\big{|}^{2}\Big{]}.$
(3.11)

The differential cross section, averaged over the initial photon polarization
and summed over the final polarizations, has the following form [20]:

$\displaystyle\frac{d\sigma}{dM_{\gamma\gamma}^{2}dtd(-u^{\prime})}=\frac{1}{2}\frac{1}{32(2\pi)^{3}S_{\gamma
N}^{2}M_{\gamma\gamma}^{2}}\sum_{\mathrm{pol.}}\sum_{\sigma_{1},\sigma_{2}}\frac{1}{4}\big{|}\mathcal{T}\big{|}^{2}=$
(3.12) $\displaystyle=\frac{1}{(16\pi)^{3}S_{\gamma
N}^{2}M_{\gamma\gamma}^{2}}\sum_{\mathrm{pol.}}\Big{[}(1-\xi^{2})\big{|}\sum_{q}\mathcal{H}^{q}\big{|}^{2}-2\xi^{2}\mathrm{Re}\sum_{q}\mathcal{H}^{q}\sum_{q^{\prime}}\mathcal{E}^{q\>*}+\frac{\xi^{4}}{1-\xi^{2}}\big{|}\sum_{q}\mathcal{E}^{q}\big{|}^{2}\Big{]}.$

In Section 3.2 we recall results of [20] and [21] concerning the LO amplitude,
in particular the term $\mathcal{C}^{q}_{0}(x,\xi,\dots)$. In the subsequent
part we present the computation of the NLO amplitude of the hard part, i.e.
$\mathcal{C}^{q}_{1}(x,\xi,\dots)$ and $\mathcal{C}^{q}_{coll.}(x,\xi,\dots)$.

To reduce the number of considered NLO diagrams (originally 48), we consider a
process with 3 incoming photons with momenta $(k_{1},k_{2},k_{3})$ and
polarisation vectors $(\epsilon_{1},\epsilon_{2},\epsilon_{3})$, and then we
sum the amplitudes over all combinations

$\big{(}(k_{1},\epsilon_{1}),(k_{2},\epsilon_{2}),(k_{3},\epsilon_{3})\big{)}$

corresponding to all permutations of

$\big{(}(-q_{1},\epsilon^{*}(\mathbf{q}_{1})),(-q_{2},\epsilon^{*}(\mathbf{q}_{2})),(q,\epsilon(\mathbf{q}))\big{)}.$

We introduce the following parametrization of these invariants:

$2pk_{i}:=\beta_{i}s,\quad 2k_{1}k_{2}:=\kappa_{3}s,\quad
2k_{2}k_{3}:=\kappa_{1}s,\quad 2k_{1}k_{3}:=\kappa_{2}s.$ (3.13)

It is easy to check that

$\kappa_{i}=-2\xi\beta_{i}.$

Scalar products of the corresponding momenta are:

$\displaystyle 2pq_{1}=\alpha s,\qquad 2pq_{2}=\bar{\alpha}s,\qquad 2pq=s;$
(3.14) $\displaystyle
2q_{1}q_{2}=M_{\gamma\gamma}^{2},\qquad-2qq_{2}=u^{\prime},\qquad-2qq_{1}=t^{\prime}.$

Hence, parameters $\beta_{i}$ and $\kappa_{i}$ are the following permutations:

$\displaystyle\\{\beta_{i}\\}=\>\mathrm{permutations}\>\mathrm{of}\>\Big{\\{}-\alpha,-\bar{\alpha},1\Big{\\}},$
(3.15)
$\displaystyle\\{\kappa_{i}\\}=\>\mathrm{permutations}\>\mathrm{of}\>\Big{\\{}\frac{M_{\gamma\gamma}^{2}}{s},\frac{u^{\prime}}{s},\frac{t^{\prime}}{s}\Big{\\}}.$

We often use the following abbreviations: $\bar{p}:=(x+\xi)p$,
$\underline{p}:=(x-\xi)p$.

###  3.2 The LO amplitude

At the leading order there are 6 tree diagrams, each corresponding to a
different permutation of photons attached to the fermion line, as it is shown
in Fig. 3.1.

$k_{1}$$k_{2}$$k_{3}$$(x+\xi)p$$(x-\xi)p$ Figure 3.1: The general form of the
LO diagram. As it was described in Sec. 3.1, to simplify the computation we
first consider a diagram with 3 incoming photons, and at the end of the
computations we use appropriate substitution and sum over permutations.

The hard-part amplitude corresponding to a given permutation reads

$Tr\big{[}i\mathcal{M}^{0}_{1,2,3}~{}\not{p}\big{]}=ie_{q}^{3}\frac{1}{2(x+\xi)pk_{1}+i0}\frac{1}{-2(x-\xi)pk_{3}+i0}Tr\Big{(}\not{\epsilon}_{3}(\underline{\not{p}}-\not{k}_{3})\not{\epsilon}_{2}(\bar{\not{p}}+\not{k}_{1})\not{\epsilon}_{1}\not{p}\Big{)}.$
(3.16)

Let us denote $i0_{i}=i0\times sgn(\beta_{i})$. Using definitions from Sec.
3.1 and $\not{p}\not{p}=0$, we rewrite (3.16) in a compact form:

$Tr\big{[}i\mathcal{M}^{0}_{1,2,3}~{}\not{p}\big{]}=ie_{q}^{3}s^{-2}\frac{1}{\beta_{1}\beta_{3}}\frac{1}{x+\xi+i0_{1}}\frac{1}{x-\xi-i0_{3}}Tr\Big{(}\not{\epsilon}_{3}\not{k}_{3}\not{\epsilon}_{2}\not{k}_{1}\not{\epsilon}_{1}\not{p}\Big{)}.$
(3.17)

The ”$-$” sign standing by $\not{k}_{3}$ in the trace got absorbed by the
”$-$” in the propagator $\frac{1}{-(x-\xi)+i0_{3}}$. It turns out that the
trace structure present here often appears in loop calculations, and because
of that we write

$4\mathcal{A}_{1,2,3}:=-Tr\Big{(}\not{\epsilon}_{3}\not{k}_{3}\not{\epsilon}_{2}\not{k}_{1}\not{\epsilon}_{1}\not{p}\Big{)}.$
(3.18)

Let us observe that the interchange of indices $1\leftrightarrow 3$ in Eq.
(3.17) does not change that expression if $sgn(\beta_{1})=sgn(\beta_{3})$, and
results in complex conjugation of the term

$\frac{1}{x+\xi+i0_{1}}\frac{1}{x-\xi-i0_{3}}$

if $sgn(\beta_{1})=-sgn(\beta_{3})$. The full amplitude of the hard sub-
process is obtained by summing over all permutations of photons $k_{i}$, so
that the only terms proportional to the imaginary part of propagators which
survive, are those from diagrams in which the incoming photon is the middle
one (so that $k_{2}=q$). It turns out, that terms proportional to the real
part of propagators cancel as well, since in the considered kinematics

$\sum_{perm.}\frac{1}{\beta_{1}\beta_{3}}Tr\Big{(}\not{\epsilon}_{3}\not{k}_{3}\not{\epsilon}_{2}\not{k}_{1}\not{\epsilon}_{1}\not{p}\Big{)}=\frac{1}{\beta_{1}\beta_{2}\beta_{3}}\sum_{perm.}\beta_{2}Tr\Big{(}\not{\epsilon}_{3}\not{k}_{3}\not{\epsilon}_{2}\not{k}_{1}\not{\epsilon}_{1}\not{p}\Big{)}=0.$
(3.19)

That can be checked by a straightforward computation, for reference see [21]
(note that the original formula in [20] contained a mistake). Hence, after
summing over all the possible photons’ permutations, one arrives at the
following hard part:

$\displaystyle\mathcal{C}^{q}_{0}(x,\xi,\dots)=\sum_{\mathrm{perm.}}Tr\big{[}\mathcal{M}^{0}_{1,2,3}~{}\not{p}\big{]}=\frac{e_{q}^{3}}{s^{2}\alpha\bar{\alpha}}\bigg{(}\frac{1}{x+\xi-i0}\frac{1}{x-\xi+i0}-c.c.\bigg{)}Tr\Big{(}\not{\epsilon}^{*}(\mathbf{q}_{2})\not{q}_{2}\not{\epsilon}(\mathbf{q})\not{q}_{1}\not{\epsilon}^{*}(\mathbf{q}_{1})\not{p}\Big{)}$
(3.20) $\displaystyle\phantom{\mathcal{C}^{q}_{0}(x,\xi,\dots)}=-i\frac{2\pi
e_{q}^{3}}{s\alpha\bar{\alpha}}\frac{1}{\xi}\Big{(}\delta(x-\xi)+\delta(x+\xi)\Big{)}\times$
$\displaystyle\qquad\bigg{[}(\alpha-\bar{\alpha})\big{(}\vec{\epsilon^{*}}_{t}(\mathbf{q}_{1})\vec{\epsilon^{*}}_{t}(\mathbf{q}_{2})\big{)}\big{(}\vec{p}_{t}\vec{\epsilon}_{t}(\mathbf{q})\big{)}-\big{(}\vec{p}_{t}\vec{\epsilon^{*}}_{t}(\mathbf{q}_{1})\big{)}\big{(}\vec{\epsilon}_{t}(\mathbf{q})\vec{\epsilon^{*}}_{t}(\mathbf{q}_{2})\big{)}+\big{(}\vec{p}_{t}\vec{\epsilon^{*}}_{t}(\mathbf{q}_{2})\big{)}\big{(}\vec{\epsilon}_{t}(\mathbf{q})\vec{\epsilon^{*}}_{t}(\mathbf{q}_{1})\big{)}\bigg{]}.$

###  3.3 Self-energy corrections

Before presenting the computation of NLO diagrams, let us recapitulate our
conventions. We use the dimensional regularization with the number of
dimensions $d=4-\varepsilon$. It results in the modification of the strong
coupling $g\rightarrow g\mu^{2-d/2}$, which ensures the correct dimensionality
of the amplitude. For gamma matrices in $d$ dimensions one has

$\displaystyle\gamma^{\mu}\gamma^{\nu}\gamma_{\mu}=(2-d)\gamma^{\nu}=(-2+\varepsilon)\gamma^{\nu},$
(3.21)
$\displaystyle\gamma^{\mu}\gamma^{\alpha}\gamma^{\beta}\gamma^{\delta}\gamma_{\mu}=-2\gamma^{\delta}\gamma^{\beta}\gamma^{\alpha}+\varepsilon\gamma^{\alpha}\gamma^{\beta}\gamma^{\delta}.$

We use the Feynman formula for combining denominators

$\frac{1}{a_{1}\dots
a_{n}}=\Gamma(n)\int_{0}^{1}du_{1}\dots\int_{0}^{1}du_{n}\delta\big{(}1-u_{1}-\dots-
u_{n}\big{)}\frac{1}{\big{(}u_{1}a_{1}+\dots u_{n}a_{n}\big{)}^{n}}.$ (3.22)

Finally, we perform momentum integration according to the following formulas:

$\displaystyle\int\frac{d^{d}k}{(2\pi)^{d}}\frac{1}{\Big{(}k^{2}-\Delta+i0\Big{)}^{n}}=\frac{i(-1)^{n}}{(4\pi)^{d/2}}\Big{(}\Delta-i0\Big{)}^{d/2-n}\frac{\Gamma\big{(}n-d/2\big{)}}{\Gamma(n)},$
(3.23)
$\displaystyle\int\frac{d^{d}k}{(2\pi)^{d}}\frac{k^{\mu}k^{\nu}}{(k^{2}+\Delta+i0)^{n}}=-\frac{i}{2}g^{\mu\nu}\frac{(-1)^{n}}{(4\pi)^{d/2}}\frac{\Gamma(n-1-d/2)}{\Gamma(n)}\Big{(}-\Delta-i0\Big{)}^{d/2+1-n}.$

Let us emphasise that, since the hard amplitude is integrated over the
parameter $x$ and may contain singularities, it is important to keep the $i0$
terms through all the calculation.
We use the following naming convention: a diagram with a loop at
left/right/middle involving $n$ propagators will be denoted by superscript
$n.L$/ $n.R$/ $n.M$. Let us recall one more time, that we often use the short-
hand notation $\bar{p}:=(x+\xi)p$, $\underline{p}:=(x-\xi)p$.

After this discussion, we can proceed with the computation of loop
corrections. The simplest ones, self-energy diagrams, are shown in Fig. 3.2.

$k_{1}$$k_{2}$$k_{3}$$(x+\xi)p$$(x-\xi)p$$k$$k_{1}$$k_{2}$$k_{3}$$(x+\xi)p$$(x-\xi)p$$k$

Figure 3.2: Diagrams $2.L$ and $2.R$.

Amplitude corresponding to $2.L$ reads

$\displaystyle
Tr\Big{[}i\mathcal{M}^{2.L}_{1,2,3}~{}\not{p}\Big{]}=-C_{F}(ie_{q})^{3}(ig)^{2}\Big{(}\frac{i}{\big{(}(x+\xi)p+k_{1}\big{)}^{2}+i0}\Big{)}^{2}\frac{i}{\big{(}(x-\xi)p-k_{3}\big{)}^{2}+i0}\times$
(3.24)
$\displaystyle\phantom{Tr\Big{[}i\mathcal{M}^{2.L}_{1,2,3}~{}\not{p}\Big{]}=}Tr\Big{\\{}\not{\epsilon}_{3}\big{[}(x-\xi)\not{p}-\not{k}_{3}\big{]}\not{\epsilon}_{2}\big{[}(x+\xi)\not{p}+\not{k}_{1}\big{]}\times$
$\displaystyle\phantom{Tr\Big{[}i\mathcal{M}^{2.L}_{1,2,3}~{}\not{p}\Big{]}=}\mu^{4-d}\int\frac{i^{2}d^{d}k}{(2\pi)^{d}}\frac{\gamma^{\mu}[(x+\xi)\not{p}+\not{k}_{1}+\not{k}]\gamma_{\mu}}{(k^{2}+i0)(((x+\xi)p+k_{1}+k)^{2}+i0)}\big{[}(x+\xi)\not{p}+\not{k}_{1}\big{]}\not{\epsilon}_{1}\not{p}\Big{\\}}.$

####  3.3.1 Integral

Using the discussed integration formula (3.23), we first employ the Feynman
parametrization and shift momentum to complete the square in denominator:

$I=\mu^{4-d}\int\frac{d^{d}k}{(2\pi)^{d}}\frac{\gamma^{\mu}(\not{P}+\not{k})\gamma_{\mu}}{(k^{2}+i0)((P+k)^{2}+i0)}=\mu^{4-d}\int\frac{d^{d}k}{(2\pi)^{d}}\int_{0}^{1}dx\frac{(2-d)(1-x)\not{P}}{[k^{2}+x(1-x)P^{2}+i0]^{2}},$
(3.25)

and then integrate over the loop momentum $k$:

$I=\mu^{4-d}(2-d)\int_{0}^{1}dx\frac{i\not{P}(1-x)}{(4\pi)^{d/2}}\big{(}-x(1-x)P^{2}-i0\big{)}^{d/2-2}\Gamma(2-d/2).$
(3.26)

In the next step we use the fact that $x(1-x)\geq 0$, so that it can be taken
outside the parenthesis without changing the sign of $i0$

$I=(-2+\varepsilon)\frac{i\not{P}}{(4\pi)^{2}}\int_{0}^{1}dx\big{(}-P^{2}/4\pi\mu^{2}-i0\big{)}^{-\varepsilon/2}x^{-\varepsilon/2}(1-x)^{1-\varepsilon/2}\Gamma(\varepsilon/2),$
(3.27)

and finally we perform the integral over $dx$ obtaining the result expressed
in terms of Euler gamma functions:

$I=(-2+\varepsilon)\frac{i\not{P}}{(4\pi)^{2}}\big{(}-P^{2}/4\pi\mu^{2}-i0\big{)}^{-\varepsilon/2}\Gamma(\varepsilon/2)\frac{\Gamma(1-\varepsilon/2)\Gamma(2-\varepsilon/2)}{\Gamma(3-\varepsilon)}.$
(3.28)

Taking the limit $\varepsilon\rightarrow 0$ we get

$(-2+\varepsilon)\Gamma(\varepsilon/2)\frac{\Gamma(1-\varepsilon/2)\Gamma(2-\varepsilon/2)}{\Gamma(3-\varepsilon)}=-\frac{2}{\varepsilon}-1+\gamma+\mathcal{O}(\varepsilon).$
(3.29)

Finally, we have

$I=\frac{i\not{p}}{(4\pi)^{2}}\big{(}-\frac{2}{\varepsilon}+log(-P^{2}/4\pi\mu^{2}-i0)-1+\gamma\big{)}.$
(3.30)

####  3.3.2 Traces

In the case of the self-energy contributions we have the following traces:

$\displaystyle
Tr\Big{\\{}\not{\epsilon}_{3}\big{[}(x-\xi)\not{p}-\not{k}_{3}\big{]}\not{\epsilon}_{2}\big{[}(x+\xi)\not{p}+\not{k}_{1}\big{]}^{3}\not{\epsilon}_{1}\not{p}\Big{\\}},$
(3.31) $\displaystyle
Tr\Big{\\{}\not{\epsilon}_{3}\big{[}(x-\xi)\not{p}-\not{k}_{3}\big{]}^{3}\not{\epsilon}_{2}\big{[}(x+\xi)\not{p}+\not{k}_{1}\big{]}\not{\epsilon}_{1}\not{p}\Big{\\}}.$

They are simplified using $\not{P}\not{P}=P^{2}$, and the problem reduces to
computation of

$Tr\Big{\\{}\not{\epsilon}_{3}\big{[}(x-\xi)\not{p}-\not{k}_{3}\big{]}\not{\epsilon}_{2}\big{[}(x+\xi)\not{p}+\not{k}_{1}\big{]}\not{\epsilon}_{1}\not{p}\Big{\\}}.$
(3.32)

Using $\not{p}\not{p}=0$ we obtain

$-Tr\Big{\\{}\not{\epsilon}_{3}\not{k}_{3}\not{\epsilon}_{2}\not{k}_{1}\not{\epsilon}_{1}\not{p}\Big{\\}}=4\mathcal{A}_{1,2,3}.$
(3.33)

Putting together the results of the integration and the trace yields

$\displaystyle Tr\big{[}i\mathcal{M}^{2.L}_{1,2,3}~{}\not{p}\big{]}=$
$\displaystyle-
C_{F}e_{q}^{3}g^{2}\frac{i}{(4\pi)^{2}}\frac{1}{(x+\xi)pk_{1}+i0}\frac{1}{(x-\xi)pk_{3}-i0}\mathcal{A}_{1,2,3}\times$
$\displaystyle\Big{(}-\frac{2}{\varepsilon}+log(-(x+\xi)(pk_{1})/2\pi\mu^{2}-i0)-1+\gamma\Big{)}.$

It can be further simplified using $pk_{1}/2\pi=s\beta_{1}/4\pi$ and

$-\frac{2}{\varepsilon}+log(-(x+\xi)s\beta_{1}/4\pi\mu^{2}-i0)+\gamma=\Big{(}\frac{se^{\gamma}}{4\pi\mu^{2}}\Big{)}^{-\frac{\varepsilon}{2}}\bigg{(}-\frac{2}{\varepsilon}+log\Big{(}-(x+\xi)\beta_{1}-i0\Big{)}\bigg{)}+\mathcal{O}(\varepsilon),$
(3.34)

so that

$\displaystyle
Tr\big{[}i\mathcal{M}^{2.L}_{1,2,3}~{}\not{p}\big{]}=-C_{F}e_{q}^{3}g^{2}\frac{i}{(4\pi)^{2}}\frac{1}{x+\xi+i0_{1}}\frac{1}{x-\xi-i0_{3}}\frac{4}{s^{2}\beta_{1}\beta_{3}}\mathcal{A}_{1,2,3}\times$
(3.35)
$\displaystyle\phantom{Tr\big{[}i\mathcal{M}^{2.L}_{1,2,3}~{}\not{p}\big{]}=}\Big{(}\frac{se^{\gamma}}{4\pi\mu^{2}}\Big{)}^{-\frac{\varepsilon}{2}}\bigg{(}-\frac{2}{\varepsilon}+log\Big{(}-(x+\xi)\beta_{1}-i0\Big{)}-1\bigg{)}$
$\displaystyle\phantom{Tr\big{[}i\mathcal{M}^{2.L}_{1,2,3}~{}\not{p}\big{]}}=Tr[i\mathcal{M}^{0}_{1,2,3}~{}\not{p}]\cdot\frac{g^{2}C_{F}}{(4\pi)^{2}}\Big{(}\frac{se^{\gamma}}{4\pi\mu^{2}}\Big{)}^{-\frac{\varepsilon}{2}}\bigg{(}-\frac{2}{\varepsilon}+log\Big{(}-(x+\xi)\beta_{1}-i0\Big{)}-1\bigg{)}$

Using this results one can write the second amplitude in a straightforward
way:

$\displaystyle
Tr\big{[}i\mathcal{M}^{2.R}_{1,2,3}~{}\not{p}\big{]}=-C_{F}(ie_{q})^{3}(ig)^{2}\frac{i}{\big{(}(x+\xi)p+k_{1}\big{)}^{2}+i0}\bigg{(}\frac{i}{\big{(}(x-\xi)p-k_{3}\big{)}^{2}+i0}\bigg{)}^{2}\times$
(3.36)
$\displaystyle\phantom{Tr\big{[}i\mathcal{M}^{2.R}_{1,2,3}~{}\not{p}\big{]}=}Tr\Big{\\{}\not{\epsilon}_{3}\big{[}(x-\xi)\not{p}-\not{k}_{3}\big{]}\mu^{4-d}\int\frac{i^{2}d^{d}k}{(2\pi)^{d}}\frac{\gamma^{\mu}[(x-\xi)\not{p}-\not{k}_{3}+\not{k}]\gamma_{\mu}}{(k^{2}+i0)(((x-\xi)p-k_{3}+k)^{2}+i0)}\times$
$\displaystyle\phantom{Tr\big{[}i\mathcal{M}^{2.R}_{1,2,3}~{}\not{p}\big{]}=}\big{[}(x-\xi)\not{p}-\not{k}_{3}\big{]}\not{\epsilon}_{2}\big{[}(x+\xi)\not{p}+\not{k}_{1}\big{]}\not{\epsilon}_{1}\not{p}\Big{\\}}=$
$\displaystyle=-C_{F}e_{q}^{3}g^{2}\frac{i}{(4\pi)^{2}}\frac{1}{x+\xi+i0_{1}}\frac{1}{x-\xi-i0_{3}}\frac{4}{s^{2}\beta_{1}\beta_{3}}\mathcal{A}_{1,2,3}\Big{(}-\frac{2}{\varepsilon}+log\big{(}(x-\xi)(pk_{3})/2\pi\mu^{2}-i0\big{)}-1+\gamma\Big{)}$
$\displaystyle=\phantom{-}Tr[i\mathcal{M}^{0}_{1,2,3}~{}\not{p}]\cdot\frac{g^{2}C_{F}}{(4\pi)^{2}}\Big{(}\frac{se^{\gamma}}{4\pi\mu^{2}}\Big{)}^{-\frac{\varepsilon}{2}}\bigg{(}-\frac{2}{\varepsilon}+log\Big{(}(x-\xi)\beta_{3}-i0\Big{)}-4\bigg{)}.$

###  3.4 Vertex corrections - 3.L and 3.R

The vertex corrections with the gluon line attached to the on-shell quark
develop both the UV and collinear divergence, as it was explained in Subsec.
2.2.2. The considered diagrams are shown in Fig. 3.3

$k_{1}$$k_{2}$$k_{3}$$(x+\xi)p$$(x-\xi)p$$k$$k_{1}$$k_{2}$$k_{3}$$(x+\xi)p$$(x-\xi)p$$k$

Figure 3.3: Diagrams 3.L and 3.R.

The full form of the amplitude corresponding to $3.L$ is

$\displaystyle
Tr\big{[}i\mathcal{M}^{3.L}_{1,2,3}~{}\not{p}\big{]}=-C_{F}(ie_{q})^{3}(ig)^{2}\frac{i}{\big{(}(x+\xi)p+k_{1}\big{)}^{2}+i0}\frac{i}{\big{(}(x-\xi)p-k_{3}\big{)}^{2}+i0}Tr\Big{\\{}\not{\epsilon}_{3}\big{[}(x-\xi)\not{p}-\not{k}_{3}\big{]}\not{\epsilon}_{2}\times$
(3.37)
$\displaystyle\big{[}(x+\xi)\not{p}+\not{k}_{1}\big{]}\mu^{4-d}\int\frac{d^{d}k}{(2\pi)^{d}}\frac{i^{3}\gamma_{\mu}\big{(}(x+\xi)\not{p}+\not{k}_{1}+\not{k}\big{)}\not{\epsilon}_{1}\big{(}(x+\xi)\not{p}+\not{k}\big{)}\gamma^{\mu}}{\big{(}((x+\xi)p+k)^{2}+i0\big{)}\big{(}((x+\xi)p+k_{1}+k)^{2}+i0\big{)}\big{(}k^{2}+i0\big{)}}\not{p}\Big{\\}}.$

First, let us first simplify the traces present in Eq. (3.37). Using the
formula

$\gamma_{\mu}\not{A}\not{B}\not{C}\gamma^{\mu}=-2\not{C}\not{B}\not{A}+\varepsilon\not{A}\not{B}\not{C},$

for arbitrary numbers $a,b,c,d$ we obtain

$\displaystyle
Tr\Big{\\{}\not{\epsilon}_{3}\big{[}(x-\xi)\not{p}-\not{k}_{3}\big{]}\not{\epsilon}_{2}\big{[}(x+\xi)\not{p}+\not{k}_{1}\big{]}\big{[}a\not{p}+b\not{k}_{1}\big{]}\not{\epsilon}_{1}\big{[}c\not{p}+d\not{k}_{1}\big{]}\not{p}\Big{\\}}=-adTr\Big{\\{}\not{\epsilon}_{3}\not{k}_{3}\not{\epsilon}_{2}\not{k}_{1}\not{p}\not{\epsilon}_{1}\not{k}_{1}\not{p}\Big{\\}}=$
(3.38)
$\displaystyle=-2ad(pk_{1})Tr\big{(}\not{\epsilon}_{3}\not{k}_{3}\not{\epsilon}_{2}\not{k}_{1}\not{p}\not{\epsilon}_{1}\big{)}=-8ad(pk_{1})\mathcal{A}_{1,2,3}.$

We proceed with the computation by applying the Feynman parametrization to the
integral:

$\displaystyle\mu^{4-d}\int\frac{d^{d}k}{(2\pi)^{d}}\frac{\gamma_{\mu}\big{(}(x+\xi)\not{p}+\not{k}_{1}+\not{k}\big{)}\not{\epsilon}_{1}\big{(}(x+\xi)\not{p}+\not{k}\big{)}\gamma^{\mu}}{\big{(}((x+\xi)p+k)^{2}+i0\big{)}\big{(}((x+\xi)p+k_{1}+k)^{2}+i0\big{)}\big{(}k^{2}+i0\big{)}}$
(3.39)
$\displaystyle=2\mu^{4-d}\int_{0}^{1}dy\int_{0}^{1-y}dz\int\frac{d^{d}k}{(2\pi)^{d}}\frac{\gamma_{\mu}\big{(}(x+\xi)\not{p}+\not{k}_{1}+\not{k}\big{)}\not{\epsilon}_{1}\big{(}(x+\xi)\not{p}+\not{k}\big{)}\gamma^{\mu}}{\Big{(}\big{(}k+(y+z)(x+\xi)p+yk_{1}\big{)}^{2}+2(x+\xi)(pk_{1})(1-y-z)y+i0\Big{)}^{3}}.$

Substituting $k^{\prime}=k+(y+z)(x+\xi)p+yk_{1}$ and using the trace formula
(3.38) we reduce the problem of computation of the integral inside the trace
to the following one:

$I=2\mu^{4-d}\int_{0}^{1}dy\int_{0}^{1-y}dz\int\frac{d^{d}k}{(2\pi)^{d}}\frac{-2(1-y-z)(x+\xi)(1+(\varepsilon/2-1)y)\not{p}\not{\epsilon}_{1}\not{k}_{1}+(-2+\varepsilon)k^{\mu}k^{\nu}\gamma_{\mu}\not{\epsilon}_{1}\gamma_{\nu}}{\Big{(}k^{2}+2(x+\xi)(pk_{1})(1-y-z)y+i0\Big{)}^{3}}.$
(3.40)

It breaks into 2 term which are proportional to different gamma structures:

$I\equiv(x+\xi)\not{p}\not{\epsilon}_{1}\not{k}_{1}I_{1}+I_{2}\not{\epsilon}_{1}.$
(3.41)

After performing the momentum integration, $I_{1}$ is found to be

$I_{1}=\frac{2i\mu^{4-d}}{(4\pi)^{d/2}}\Gamma\big{(}3-\frac{d}{2}\big{)}\int_{0}^{1}dy\int_{0}^{1-y}dz\frac{(1-y-z)(1+(\varepsilon/2-1)y)}{\big{(}-2(x+\xi)(pk_{1})(1-y-z)y-i0\big{)}^{3-d/2}}.$
(3.42)

The integration over Feynman parameters $y$ and $z$ can be simplified by the
substitution:

$y=\bar{z}(1-\bar{y}),\quad\quad
z=\bar{y}\bar{z},\quad\quad\int_{0}^{1}dy\int_{0}^{1-y}dz\rightarrow\int_{0}^{1}d\bar{y}\int_{0}^{1}\bar{z}d\bar{z}.$
(3.43)

$I_{1}$ has now the following form (for simplicity we omit the lines over the
reparametrized variables $\bar{y},\bar{z}$):

$I_{1}=\Gamma\big{(}3-\frac{d}{2}\big{)}\frac{2i\mu^{4-d}}{(4\pi)^{d/2}}\int_{0}^{1}dy\int_{0}^{1}zdz\frac{(1-z)\big{(}1+(\varepsilon/2-1)z(1-y)\big{)}}{\big{(}-2(x+\xi)(pk_{1})(1-z)(1-y)z-i0\big{)}^{3-d/2}}.$

The integration over $dx$ and $dz$ can be now factorized into a product of 2
integrals:

$\displaystyle I_{1}=$
$\displaystyle\frac{2i\mu^{4-d}}{(4\pi)^{d/2}}\big{(}-2(x+\xi)(pk_{1})-i0\big{)}^{d/2-3}\Gamma\big{(}3-\frac{d}{2}\big{)}\Big{(}\int_{0}^{1}(1-y)^{d/2-3}dy\int_{0}^{1}z^{d/2-2}(1-z)^{d/2-2}$
$\displaystyle+(\varepsilon/2-1)\int_{0}^{1}(1-y)^{d/2-2}\int_{0}^{1}z^{d/2-1}(1-z)^{d/2-2}dz\Big{)}.$

These integral are easily performed and one obtains

$\displaystyle I_{1}=$
$\displaystyle\frac{2i\mu^{4-d}}{(4\pi)^{d/2}}\big{(}-2(x+\xi)(pk_{1})-i0\big{)}^{d/2-3}\times$
(3.44)
$\displaystyle\Big{(}\frac{\Gamma\big{(}3-\frac{d}{2}\big{)}}{d/2-2}\frac{\Gamma(d/2-1)\Gamma(d/2-1)}{\Gamma(d-2)}+(\varepsilon/2-1)\frac{\Gamma\big{(}3-\frac{d}{2}\big{)}}{d/2-1}\frac{\Gamma(d/2)\Gamma(d/2-1)}{\Gamma(d-1)}\Big{)}.$

For $d=4-\varepsilon$, in the limit $\varepsilon\rightarrow 0$ one has

$\displaystyle\frac{\Gamma\big{(}3-\frac{d}{2}\big{)}}{d/2-2}\frac{\Gamma(d/2-1)\Gamma(d/2-1)}{\Gamma(d-2)}\approx-\frac{2}{\varepsilon}-2+\gamma,$
(3.45)
$\displaystyle(\varepsilon/2-1)\frac{\Gamma\big{(}3-\frac{d}{2}\big{)}}{d/2-1}\frac{\Gamma(d/2)\Gamma(d/2-1)}{\Gamma(d-1)}\approx-\frac{1}{2}.$

In this way, we obtain the final form of $I_{1}$:

$\displaystyle I_{1}$
$\displaystyle=\frac{2i}{(4\pi)^{2}}\big{(}-2(x+\xi)(pk_{1})-i0\big{)}^{-1}\Big{(}1-\frac{\varepsilon}{2}log\big{(}-(x+\xi)(pk_{1})/2\pi\mu^{2}-i0\big{)}\Big{)}\Big{(}-\frac{2}{\varepsilon}+\gamma-\frac{5}{2}\Big{)}$
(3.46)
$\displaystyle=\frac{2i}{(4\pi)^{2}}\big{(}-2(x+\xi)(pk_{1})-i0\big{)}^{-1}\Big{(}-\frac{2}{\varepsilon}+log\big{(}-(x+\xi)(pk_{1})/2\pi\mu^{2}-i0\big{)}+\gamma-\frac{5}{2}\Big{)}.$

Let us note, that the loop momentum integral inside of $I_{1}$ was UV safe,
from the argument of power counting, and hence the divergence present therein
is the result of the region of integration, where $k$ and $(x+\xi)p+k$ become
collinear and lightlike (i.e. the collinear divergence). In such case, the
divergence manifested in the integration over the Feynman parameters – one can
easily check, that if one assumed a non-zero mass of partons, the relevant
integral would be finite.

The second integral, after using
$\gamma^{\mu}\not{\epsilon}_{1}\gamma_{\mu}=(-2+\varepsilon)\not{\epsilon}_{1}$
and integrating over momenta becomes

$I_{2}\not{\epsilon}_{1}=\frac{i\mu^{4-d}}{2(4\pi)^{d/2}}\Gamma\big{(}2-d/2\big{)}(2-\varepsilon)^{2}\not{\epsilon}_{1}\int_{0}^{1}dy\int_{0}^{1-y}dz\frac{1}{\big{(}-2(x+\xi)(pk_{1})(1-y-z)y-i0\big{)}^{2-d/2}}.$

We see that this integral is not UV-safe and there is a divergent term
$\Gamma\big{(}2-d/2\big{)}$ resulting from the integration over the loop
momentum. Now we use the same substitution as in Eq. (3.43):

$I_{2}\not{\epsilon}_{1}=\frac{i\mu^{4-d}}{2(4\pi)^{d/2}}\Gamma\big{(}2-d/2\big{)}(2-\varepsilon)^{2}\not{\epsilon}_{1}\int_{0}^{1}dy\int_{0}^{1}zdz\big{(}-2(x+\xi)(pk_{1})(1-z)(1-y)z-i0\big{)}^{d/2-2},$

because $(1-z)(1-y)z\geq 0$, we can bring that term outside the parenthesis
without changing the sign of $i0$, and as a result we get

$I_{2}\not{\epsilon}_{1}=\frac{i\mu^{4-d}}{2(4\pi)^{d/2}}\Gamma\big{(}2-d/2\big{)}(2-\varepsilon)^{2}\not{\epsilon}_{1}\big{(}-2(x+\xi)(pk_{1})-i0\big{)}^{d/2-2}\frac{\Gamma(d/2)\Gamma(d/2-1)}{\Gamma(d-1)}\frac{1}{d/2-1}.$

We expand with respect to $\varepsilon$:

$\Gamma\big{(}2-d/2\big{)}(2-\varepsilon)^{2}\frac{\Gamma(d/2)\Gamma(d/2-1)}{\Gamma(d-1)}\frac{1}{d/2-1}\approx
2\big{(}\frac{2}{\varepsilon}+1-\gamma\big{)}.$

In this case, the integrals over Feynman parameters are finite, since the
$k^{2}$ from the denominator cancels with the $k^{2}$ from the trace, and the
integrand does not have a collinear singularity. The final expression for the
second integral is then

$\displaystyle I_{2}\not{\epsilon}_{1}$
$\displaystyle=\frac{i\mu^{4-d}}{2(4\pi)^{d/2}}\not{\epsilon}_{1}\big{(}-2(x+\xi)(pk_{1})-i0\big{)}^{d/2-2}2\big{(}\frac{2}{\varepsilon}+1-\gamma\big{)}=$
(3.47)
$\displaystyle=\frac{i}{(4\pi)^{2}}\not{\epsilon}_{1}\Big{(}\frac{2}{\varepsilon}+1-\gamma-
log\big{(}-(x+\xi)pk_{1}/2\pi\mu^{2}-i0\big{)}\Big{)}.$

Let us recall that

$I=(x+\xi)\not{p}\not{\epsilon}_{1}\not{k}_{1}I_{1}+I_{2}\not{\epsilon}_{2}.$

Using the previous results concerning traces we obtain

$\displaystyle Tr\big{[}i\mathcal{M}^{3.L}_{1,2,3}~{}\not{p}\big{]}$
$\displaystyle=-C_{F}e_{q}^{3}g^{2}\frac{i}{(4\pi)^{2}}\frac{1}{\bar{p}k_{1}+i0}\frac{1}{\underline{p}k_{3}-i0}\mathcal{A}_{1,2,3}\Big{(}-\frac{2}{\varepsilon}+log\big{(}-(x+\xi)(pk_{1})/2\pi\mu^{2}-i0\big{)}+\gamma-4\Big{)}$
(3.48)
$\displaystyle=Tr[i\mathcal{M}^{0}_{1,2,3}~{}\not{p}]\cdot\frac{g^{2}C_{F}}{(4\pi)^{2}}\Big{(}\frac{se^{\gamma}}{4\pi\mu^{2}}\Big{)}^{-\frac{\varepsilon}{2}}\bigg{(}-\frac{2}{\varepsilon}+log\Big{(}-(x+\xi)\beta_{1}-i0\Big{)}-4\bigg{)}$

The computation of the diagram $3.R$ is done in the same way, and it yields

$\displaystyle Tr\big{[}i\mathcal{M}^{3.R}_{1,2,3}~{}\not{p}\big{]}$
$\displaystyle=-C_{F}e_{q}^{3}g^{2}\frac{i}{(4\pi)^{2}}\frac{1}{\bar{p}k_{1}+i0}\frac{1}{\underline{p}k_{3}-i0}\mathcal{A}_{1,2,3}\Big{(}-\frac{2}{\varepsilon}+log\big{(}(x-\xi)(pk_{3})/2\pi\mu^{2}-i0\big{)}+\gamma-4\Big{)}$
(3.49)
$\displaystyle=Tr[i\mathcal{M}^{0}_{1,2,3}~{}\not{p}]\cdot\frac{g^{2}C_{F}}{(4\pi)^{2}}\Big{(}\frac{se^{\gamma}}{4\pi\mu^{2}}\Big{)}^{-\frac{\varepsilon}{2}}\bigg{(}-\frac{2}{\varepsilon}+log\Big{(}(x-\xi)\beta_{3}-i0\Big{)}-4\bigg{)}.$

###  3.5 Vertex corrections - 3.M

The considered diagram is shown in fig. 3.4

$k_{1}$$k_{2}$$k_{3}$$(x+\xi)p$$(x-\xi)p$$k$ Figure 3.4: Diagram 3.M.

The fact that the gluon line in the diagram $3.M$ is connected to quarks that
are off-shell implies that in this term there is no collinear divergence here.
The integration over the Feynman parameters becomes more complicated, but it
turns out that an appropriate substitution allows to simplify it
significantly.
The considered amplitude reads

$\displaystyle
Tr\big{[}i\mathcal{M}_{3.M}~{}\not{p}\big{]}=-C_{F}(ie_{q})^{3}(ig)^{2}\frac{i}{\big{(}(x+\xi)p+k_{1}\big{)}^{2}+i0}\frac{i}{\big{(}(x-\xi)p-k_{3}\big{)}^{2}+i0}Tr\Big{\\{}\not{\epsilon}_{3}\big{[}(x-\xi)\not{p}-\not{k}_{3}\big{]}\times$
(3.50)
$\displaystyle\mu^{4-d}\int\frac{d^{d}k}{(2\pi)^{d}}\frac{i^{3}\gamma_{\mu}\big{[}(x+\xi)\not{p}+\not{k}_{1}+\not{k}_{2}+\not{k}\big{]}\not{\epsilon}_{2}\big{[}(x+\xi)\not{p}+\not{k}_{1}+\not{k}\big{]}\gamma^{\mu}}{\big{(}((x+\xi)p+k_{1}+k_{2}+k)^{2}+i0\big{)}\big{(}((x+\xi)p+k_{1}+k)^{2}+i0\big{)}\big{(}k^{2}+i0\big{)}}\big{[}(x+\xi)\not{p}+\not{k}_{1}\big{]}\not{\epsilon}_{1}\not{p}\Big{\\}}.$

First, we simplify the trace:

$\displaystyle
Tr\Big{\\{}\not{\epsilon}_{3}(-\not{k}_{3})\big{[}a(x+\xi)\not{p}+a\not{k}_{1}+b\not{k}_{2}\big{]}\not{\epsilon}_{2}\big{[}a(x+\xi)\not{p}+a\not{k}_{1}+c\not{k}_{2}\big{]}\not{k}_{1}\not{\epsilon}_{1}\not{p}\Big{\\}}=$
(3.51)
$\displaystyle=2pk_{1}Tr\Big{\\{}\not{\epsilon}_{3}\not{k}_{3}\big{[}a(x+\xi)\not{p}+a\not{k}_{1}+b\not{k}_{2}\big{]}\not{\epsilon}_{2}\big{[}a(x+\xi)\not{p}+c\not{k}_{2}\big{]}\not{\epsilon}_{1}\Big{\\}}+$
$\displaystyle-
Tr\Big{\\{}\not{\epsilon}_{3}\not{k}_{3}\big{[}a(x+\xi)\not{p}+a\not{k}_{1}+b\not{k}_{2}\big{]}\not{\epsilon}_{2}c\not{k}_{2}\not{p}\not{k}_{1}\not{\epsilon}_{1}\Big{\\}}.$

Using identities $(x+\xi)p+k_{1}=(x-\xi)p-k_{2}-k_{3}$ and
$\not{k}_{i}\not{k}_{i}=0$ we obtain

$\displaystyle\eqref{eq-
tr-3m}=2pk_{1}Tr\Big{\\{}\not{\epsilon}_{3}\not{k}_{3}\big{[}a(x-\xi)\not{p}+(b-a)\not{k}_{2}\big{]}\not{\epsilon}_{2}\big{[}a(x+\xi)\not{p}+c\not{k}_{2}\big{]}\not{\epsilon}_{1}\Big{\\}}-Tr\Big{\\{}\not{\epsilon}_{3}\not{k}_{3}a(x-\xi)\not{p}\not{\epsilon}_{2}c\not{k}_{2}\not{p}\not{k}_{1}\not{\epsilon}_{1}\Big{\\}}$
(3.52)
$\displaystyle=2a(x-\xi)c(pk_{1})Tr\Big{\\{}\not{\epsilon}_{3}\not{k}_{3}\not{p}\not{\epsilon}_{2}\not{k}_{2}\not{\epsilon}_{1}\Big{\\}}+2a(x+\xi)(b-a)(pk_{1})Tr\Big{\\{}\not{\epsilon}_{3}\not{k}_{3}\not{k}_{2}\not{\epsilon}_{2}\not{p}\not{\epsilon}_{1}\Big{\\}}+$
$\displaystyle-2ac(x-\xi)(pk_{2})Tr\Big{\\{}\not{\epsilon}_{3}\not{k}_{3}\not{p}\not{\epsilon}_{2}\not{k}_{1}\not{\epsilon}_{1}\Big{\\}}\equiv
a\big{(}\mathcal{B}(b-a)+\mathcal{C}c\big{)}.$

$\mathcal{B}$ and $\mathcal{C}$ are auxiliary parameters, which do not appear
in the final expression for the amplitude. Let us further simplify
$\mathcal{C}$. First, apply $pk_{2}=-p(k_{1}+k_{3})$ to the second trace
inside the parenthesis:

$\displaystyle\mathcal{C}$
$\displaystyle=2(x-\xi)\bigg{(}(pk_{1})Tr\Big{\\{}\not{\epsilon}_{3}\not{k}_{3}\not{p}\not{\epsilon}_{2}\not{k}_{2}\not{\epsilon}_{1}\Big{\\}}-(pk_{2})Tr\Big{\\{}\not{\epsilon}_{3}\not{k}_{3}\not{p}\not{\epsilon}_{2}\not{k}_{1}\not{\epsilon}_{1}\Big{\\}}\bigg{)}=$
$\displaystyle=2(x-\xi)\bigg{(}(pk_{1})Tr\Big{\\{}\not{\epsilon}_{3}\not{k}_{3}\not{p}\not{\epsilon}_{2}(\not{k}_{1}+\not{k}_{2})\not{\epsilon}_{1}\Big{\\}}+(pk_{3})Tr\Big{\\{}\not{\epsilon}_{3}\not{k}_{3}\not{p}\not{\epsilon}_{2}\not{k}_{1}\not{\epsilon}_{1}\Big{\\}}\bigg{)}.$

Then, from the relation
$\not{p}\not{\epsilon}_{2}(\not{k}_{1}+\not{k}_{2})=-\not{p}\not{\epsilon}_{2}\not{k}_{3}$
follows that

$\mathcal{C}=2(x-\xi)\bigg{(}-(pk_{1})Tr\Big{\\{}\not{\epsilon}_{3}\not{k}_{3}\not{p}\not{\epsilon}_{2}\not{k}_{3}\not{\epsilon}_{1}\Big{\\}}+(pk_{3})Tr\Big{\\{}\not{\epsilon}_{3}\not{k}_{3}\not{p}\not{\epsilon}_{2}\not{k}_{1}\not{\epsilon}_{1}\Big{\\}}\bigg{)}.$

Finally, using a straightforward trace algebra one verifies that

$\displaystyle\mathcal{C}$
$\displaystyle=2(x-\xi)\bigg{(}-8(pk_{1})(pk_{3})(k_{3}\epsilon_{1})(\epsilon_{2}\epsilon_{3})+(pk_{3})Tr\Big{\\{}\not{\epsilon}_{3}\not{k}_{3}\not{p}\not{\epsilon}_{2}\not{k}_{1}\not{\epsilon}_{1}\Big{\\}}\bigg{)}=$
(3.53)
$\displaystyle=2(x-\xi)pk_{3}Tr\Big{\\{}\not{\epsilon}_{3}\not{p}\not{\epsilon}_{2}\not{k}_{1}\not{k}_{3}\not{\epsilon}_{1}\Big{\\}}.$

To perform the integration, let us denote $P=(x+\xi)p+k_{1}$, $Q=P+k_{2}$.
Obviously $(P-Q)^{2}=0$, so that $P^{2}+Q^{2}=2PQ$. We take the integral
(3.50), employ the Feynman parametrization and shift the loop momentum to
complete the square. The part of the resulting integral without the $k^{2}$
term in numerator, after computing the trace according to the previous
results, reads

$I_{1}=-2i\mu^{4-d}\int\frac{d^{d}k}{(2\pi)^{d}}\int_{0}^{1}dy\int_{0}^{1-y}dz\frac{-2(1-y-z)\big{(}(\mathcal{C}-\mathcal{B})+(1-\varepsilon/2)\mathcal{B}y-(1-\varepsilon/2)\mathcal{C}z\big{)}}{\Big{(}k^{2}+(1-y-z)(yP^{2}+zQ^{2})+i0\Big{)}^{3}}.$
(3.54)

Performing the momenta integration and making the substitution (3.43) one
obtains

$\displaystyle I_{1}$
$\displaystyle=\frac{2\mu^{4-d}}{(4\pi)^{d/2}}\int_{0}^{1}dy\int_{0}^{1}zdz\Big{(}-(1-z)z\big{(}(1-y)P^{2}+yQ^{2}\big{)}-i0\Big{)}^{d/2-3}\Gamma(3-d/2)\times$
(3.55)
$\displaystyle(1-z)\Big{(}(\mathcal{C}-\mathcal{B})+z\big{(}1-\frac{\varepsilon}{2}\big{)}\big{(}\mathcal{B}(1-y)-y\mathcal{C}\big{)}\Big{)}.$

The integrals over $dz$ are finite:

$\displaystyle I_{1}$
$\displaystyle=\frac{2\mu^{4-d}}{(4\pi)^{d/2}}\Gamma(3-d/2)\bigg{(}\frac{\Gamma(d/2-1)\Gamma(d/2-1)}{\Gamma(d-2)}\int_{0}^{1}dy\Big{(}y(P^{2}-Q^{2})-(P^{2}+i0)\Big{)}^{d/2-3}(\mathcal{C}-\mathcal{B})+$
(3.56)
$\displaystyle+\frac{\Gamma(d/2)\Gamma(d/2-1)}{\Gamma(d-1)}(1-\frac{\varepsilon}{2})\Big{[}\mathcal{B}\int_{0}^{1}dy\Big{(}y(P^{2}-Q^{2})-(P^{2}+i0)\Big{)}^{d/2-3}+$
$\displaystyle-(\mathcal{C}+\mathcal{B})\int_{0}^{1}dy\Big{(}y(P^{2}-Q^{2})-(P^{2}+i0)\Big{)}^{d/2-3}y\Big{]}\bigg{)}.$

Now let us calculate integrals over $dy$. The first one (second line of Eq.
(3.56)) is

$\displaystyle\frac{\mu^{4-d}}{{(4\pi)^{d/2}}}\int_{0}^{1}dy\Big{(}y(P^{2}-Q^{2})-(P^{2}+i0)\Big{)}^{d/2-3}=$
$\displaystyle=\frac{\mu^{4-d}}{{(4\pi)^{d/2}}}\frac{1}{d/2-2}\frac{1}{P^{2}-Q^{2}}\Big{(}\big{(}-Q^{2}-i0\big{)}^{d/2-2}-\big{(}-P^{2}-i0\big{)}^{d/2-2}\Big{)}.$

In the limit $d\rightarrow 4$ this expression is finite and equal to

$\frac{1}{(4\pi)^{2}}\frac{1}{P^{2}-Q^{2}}\Big{(}log\big{(}-Q^{2}/4\pi\mu^{2}-i0\big{)}-log\big{(}-P^{2}/4\pi\mu^{2}-i0\big{)}^{d/2-2}\Big{)}.$
(3.57)

The second integral (in the last line of Eq. (3.56)) yields

$\displaystyle\frac{\mu^{4-d}}{{(4\pi)^{d/2}}}\int_{0}^{1}ydy\Big{(}y(P^{2}-Q^{2})-(P^{2}+i0)\Big{)}^{d/2-3}=$
(3.58)
$\displaystyle=\frac{\mu^{4-d}}{{(4\pi)^{d/2}}}\frac{1}{d/2-2}\frac{1}{P^{2}-Q^{2}}\Big{(}(-Q^{2}-i0)^{d/2-2}-\frac{1}{d/2-1}\frac{1}{P^{2}-Q^{2}}\big{(}(-Q^{2}-i0)^{d/2-1}-(-P^{2}-i0)^{d/2-1}\big{)}\Big{)}$
$\displaystyle=\frac{1}{(4\pi)^{2}}\frac{1}{(P^{2}-Q^{2})^{2}}\Big{(}P^{2}\log\frac{Q^{2}+i0}{P^{2}+i0}+P^{2}-Q^{2}\Big{)}.$

Notice that both values are finite. Finally, we get:

$\displaystyle I_{1}=$
$\displaystyle\frac{4}{(4\pi)^{2}}\bigg{(}2(\mathcal{C}-\frac{1}{2}\mathcal{B})\frac{log(-Q^{2}/4\pi\mu^{2}-i0)-log(-P^{2}/4\pi\mu^{2}-i0)}{P^{2}-Q^{2}}+$
(3.59)
$\displaystyle-(\mathcal{B}+\mathcal{C})\frac{1}{(P^{2}-Q^{2})^{2}}\Big{(}P^{2}\log\frac{Q^{2}+i0}{P^{2}+i0}+P^{2}-Q^{2}\Big{)}\bigg{)}.$

Now let us compute the second part of the loop momentum integral with the term
$k^{2}$ in the numerator.

$\displaystyle
I_{2}=2\mu^{4-d}\int\frac{d^{d}k}{(2\pi)^{d}}\int_{0}^{1}dy\int_{0}^{1-y}dz\frac{\gamma_{\mu}\not{k}\not{\epsilon}_{2}\not{k}\gamma^{\mu}}{\big{(}k^{2}+(1-y-z)(yP^{2}+zQ^{2})+i0\big{)}^{3}}$
(3.60)
$\displaystyle=\not{\epsilon}_{2}\frac{i\mu^{4-d}}{(4\pi)^{d/2}}\frac{\Gamma(2-d/2)}{2}(-2+\varepsilon)^{2}\int_{0}^{1}dy\int_{0}^{1}zdz\not{\epsilon}_{2}\big{(}-(1-z)z((1-y)P^{2}+yQ^{2})-i0\big{)}^{d/2-2}$
$\displaystyle=\not{\epsilon}_{2}\frac{i\mu^{4-d}}{(4\pi)^{d/2}}\frac{\Gamma(2-d/2)}{2}(-2+\varepsilon)^{2}\frac{\Gamma(d/2)\Gamma(d/2-1)}{\Gamma(d-1)}\frac{1}{d/2-1}\frac{(-Q^{2}-i0)^{d/2-1}-(-P^{2}-i0)^{d/2-1}}{P^{2}-Q^{2}}$
$\displaystyle=\not{\epsilon}_{2}\frac{i}{(4\pi)^{2}}\Big{(}\frac{2}{\varepsilon}+1-\gamma+\frac{Q^{2}log(-Q^{2}/4\pi\mu^{2}-i0)-P^{2}log(-P^{2}/4\pi\mu^{2}-i0)}{P^{2}-Q^{2}}\Big{)}.$

After taking the traces one obtains

$\displaystyle
Tr\big{[}i\mathcal{M}^{3.M}_{1,2,3}~{}\not{p}\big{]}=-C_{F}e_{q}^{3}g^{2}\frac{i}{(4\pi)^{2}}\frac{1}{\bar{p}k_{1}+i0}\frac{1}{\underline{p}k_{3}-i0}\times$
(3.61)