# Probing Hadron Structure in Proton-Nucleus Collisions

David Zaslavsky

(2014)

_set_boolean_T:n outer_mult [datatype=bibtex] [fieldsource=howpublished,
match=.+Stack, final=true] [typesource=misc, typetarget=stackexchange,
final=true] stackexchange [fieldsource=author, match=.+*$https?://.+$]
[fieldset=nameaddon, origfieldval] [fieldsource=author,
match=(.+)*$https?://.+$, replace=$1$] [fieldsource=nameaddon,
match=.+*$(https?://.+)$, replace=$1$] [fieldsource=note,
match=.+$version:\s\d{4}-\d{2}-\d{2}$] [fieldset=date, origfieldval]
[fieldset=note, null] [fieldsource=date,
match=.+$version:\s(\d{4}-\d{2}-\d{2})$, replace=$1$] [fieldset=eprint, null]
Physics December 2014 Dissertation The Graduate School Eberly College of
Science 4 [Dissertation Advisor, Chair of Committee]Anna M. StaśtoAssociate
Professor of Physics []Irina MocioiuAssociate Professor of Physics []Mark
StrikmanDistinguished Professor of Physics []James BrannickProfessor of
Mathematics [Graduate Program Chair]Richard RobinettProfessor of Physics

abstract.tex

acknowledgements.tex  This document is distributed under the Creative Commons
Attribution 4.0 License.

## Introduction

Ever since quantum chromodynamics, the theory of strong interactions, was
developed in the 1960’s, high-energy physicists have endeavored to understand
the internal structure of strongly bound particles. This effort was driven, in
large part, by experimental results from increasingly powerful particle
accelerators. But as accelerator technology has become more and more complex,
so has the depth of structure displayed by the particles being collided.
Instead of having a fixed set of constituents, as one might classically
expect, protons and atomic nuclei act as though they contain ever-increasing
numbers of smaller and smaller particles as we probe them at higher and higher
energies.

Physicists traditionally characterize this apparently changing composition
using a set of parton distributions, functions which characterize how the
momentum is distributed among the constituents (partons) of the strongly bound
particle. Although these functions emerge from fundamentally nonperturbative
physics, it is possible to describe how they vary with the parameters of the
collision using the methods of perturbative QCD. At high energies, which
typically correlate to small momentum fractions $x$, a renormalization group
equation derived by Balitsky, Fadin, Kuraev, and Lipatov [3, 4, 5, 6, 7]
governs the evolution of the parton distributions as $x$ changes. Physically,
this equation, called the BFKL evolution equation, represents the process of
parton branching, in which partons effectively divide into other partons with
smaller fractions of the original hadron’s momentum. The BFKL equation
predicts a rapid rise in the parton distributions, especially the gluon
distribution, as $x$ becomes small.

However, this growth in the gluon distribution is not sustainable. At very low
values of $x$ it leads to violations of the Froissart bound [8, 9, 10], a
fundamental limit on the growth of cross sections with energy. Violation of
the Froissart bound is equivalent to breaking unitarity — that is, to
nonconservation of probability. In the 1990’s, Balitsky and Kovchegov [11, 12,
13, 14, 15, 16] investigated the dynamics of the gluon distribution as it
approaches these very low values of $x$, finding an important contribution
from gluon recombination, the inverse of the branching process that underlies
BFKL evolution. They derived a correction term to the BFKL equation that
represents this recombination, a term which restricts the growth of the gluon
distribution at small $x$, in a phenomenon known as gluon saturation.
Accelerator technology has only recently advanced to the point where we may be
able to experimentally detect the effects of gluon saturation. Naturally, one
of the exciting topics in the world of modern particle physics is the search
for these effects.

Detecting experimental signatures of saturation generally proceeds along these
lines:

  1. 1.

Develop a theoretical model for the gluon distribution which exhibits
saturation

  2. 2.

Insert that model into a formula which relates the gluon distribution to some
measurable quantity, and calculate results for various parameters

  3. 3.

Identify features in the results which correspond to characteristics of the
gluon distribution that identify saturation

  4. 4.

Compare the results to corresponding experimental data points to see if those
features are present

For the process to be successful, the formula relating the gluon distribution
to a measurable quantity needs to be sensitive to small changes in the gluon
distribution in the region of interest, and it also needs to be precise in the
sense that it does not introduce significant theoretical uncertainties. There
is an ongoing effort in the saturation physics community to find measurements
which satisfy these criteria.

In this dissertation, I’ll describe two such measurements. First, I’ll examine
the derivation of, and predictions for, the azimuthal angular correlation in
Drell-Yan lepton pair-hadron production. This quantity is essentially the
cross section for the process AhX, expressed as a function of the angular
separation between two outgoing momenta: the total momentum of the lepton pair
, and the momentum of the hadron h. We’ll see that the correlation exhibits
characteristic features which are sensitive to details of the gluon
distribution in the saturation regime, including a double peak for back-to-
back emission. Furthermore, examining the angular correlation in the back-to-
back emission region allows us to nearly isolate the effect of the low-
momentum part of the gluon distribution, and the correlation in the near-
parallel emission region gives us a window on the high-momentum gluon
distribution.

Afterwards, I’ll describe the calculation of the inclusive hadron production
cross section, AhX, in the saturation formalism, with the complete next-to-
leading order (NLO) corrections (except next-to-leading corrections to BK
evolution). Unlike the angular correlation calculation, which only
incorporates leading order terms corresponding to tree-level diagrams, the NLO
inclusive hadron production calculation involves loop diagrams, which come
with their own set of theoretical challenges. In particular, the formulas
contain divergences in rapidity and momentum which need to be canceled by
subtractions, and then analytical singularities need to be removed to produce
formulas which can be evaluated numerically. Although the inclusive hadron
production cross section does not separate the low- and high-momentum regions
of the gluon distribution as well as the correlation does, it does provide a
very precise test since the NLO corrections act to reduce the theoretical
uncertainty in the formula.

However, the LO+NLO calculation for the inclusive hadron cross section is only
applicable in a limited kinematic regime. We’ll see that it produces negative
results at high transverse momenta, so evidently the formula is not useful as
a probe of the gluon distribution in those conditions. This implies that
further corrections to the formula are necessary to use the high-momentum
cross section to accurately probe the gluon distribution. I’ll finish by
discussing some of the first work toward that goal, which examines two types
of higher-order corrections: a simple resummation of the most prominent
contributions of higher-order loop diagrams, and an exact kinematic formula
that matches the collinear factorization result which is known to be positive.
The exact kinematic formula, in particular, is quite successful at describing
the experimental data at high transverse momentum where the simple LO+NLO
calculation fails.

At each step, I’ll present the theoretical analysis along with comparisons to
experimental data. Currently, the state of the art in experimental data for
small-$x$ physics comes from the BRAHMS and STAR detectors at the Relativistic
Heavy Ion Collider (RHIC). These detectors share several common
characteristics that make them suitable for the studies I will describe in
this dissertation: they are capable of measuring particles emitted at forward
angles close to the beamline, equivalent to high pseudorapidity $\eta\sim
3-4$, and they have good resolution when measuring the momentum of emitted
pions and other light hadrons at such forward pseudorapidities. There are also
detectors at the Large Hadron Collider (LHC) with these capabilities, namely
ALICE, ATLAS, CMS, and LHCf, although most of the relevant experimental
results from these detectors have not yet been published, so when comparing to
data I will concentrate on the BRAHMS and STAR results.

The body of the dissertation is structured as follows. Chapter Hadron
Structure From QCD lays out the basic principles of small-$x$ and saturation
physics, together with its underlying background in quantum chromodynamics,
and then chapter The Color Glass Condensate addresses specific concepts of
saturation physics in more detail, in particular the color glass condensate
model which describes target nucleons and nuclei under saturation conditions.
The remaining chapters describe the calculations previously mentioned. Chapter
Drell-Yan Correlations shows the Drell-Yan lepton pair-hadron angular
correlation, which illustrates the basic principles of calculations in
small-$x$ physics. Then, chapter NLO Inclusive Pion Production presents the
cross section for single inclusive pion production incorporating next-to-
leading order corrections, and describes some of the difficulties inherent in
treating the higher-order contributions. Chapter Computational Results and
Analysis shows the numerical results of this calculation, and then chapter
Extensions Beyond NLO discusses how additional contributions might be used to
remedy some of the shortcomings of the LO+NLO formula. Finally, chapter
Summary and Conclusion summarizes the results and presents an outlook toward
further development of these calculations.

## Hadron Structure From QCD

Quantum chromodynamics (QCD), the theory that describes the behavior of
strongly interacting particles, was first developed in the 1970s by physicists
exploring generalizations of quantum electrodynamics (the quantum theory of
electromagnetism). Within a few short years, it became clear that constituents
of the proton with small fractions $x$ of its momentum would play an
increasingly important role in the physics as collision energy increased. This
chapter lays the groundwork for the discussion of small-$x$ physics in chapter
The Color Glass Condensate by motivating the idea of a target hadron or
nucleus as a dressed bound state of quarks, and introducing the integrated and
unintegrated parton distributions used to describe it.

### Quark and Gluon Fields

QCD is based on the $\mathrm{SU}(3)$ transformation group, whose fundamental
representation is the set of all matrices $U$ that satisfy three properties:

  * •

Special

$\det U=1$ (1)

  * •

Unitary

$\herm{U}U=1$ (2)

  * •

Three-dimensional

$U=\begin{pmatrix}U_{11}&U_{12}&U_{13}\\\ U_{21}&U_{22}&U_{23}\\\
U_{31}&U_{32}&U_{33}\end{pmatrix}$ (3)

where $U_{ij}$ (with two subscripts) are complex numbers.

To construct a physical theory from this group, we first need a set of quantum
fields which $\mathrm{SU}(3)$ transformations can act upon. Since the simplest
nontrivial representation of $\mathrm{SU}(3)$ is the three-dimensional
fundamental representation, the simplest nontrivial set of fields we can use
will have three members: $q_{r}$, $q_{g}$, and $q_{b}$.

$\displaystyle\begin{pmatrix}q^{\prime}_{r}\\\ q^{\prime}_{g}\\\
q^{\prime}_{b}\end{pmatrix}$
$\displaystyle=\begin{pmatrix}U_{rr}&U_{rg}&U_{rb}\\\ U_{gr}&U_{gg}&U_{gb}\\\
U_{br}&U_{bg}&U_{bb}\end{pmatrix}\begin{pmatrix}q_{r}\\\ q_{g}\\\
q_{b}\end{pmatrix}$ (4) $\displaystyle q^{\prime}_{k}$
$\displaystyle=U_{kj}q_{j}$ (5)

where I use the usual convention of summing over repeated indices. These are
the quark fields.

We then put these fields into a Lagrangian in combinations which are
$\mathrm{SU}(3)$ invariant: when performing an $\mathrm{SU}(3)$ transformation
on the fields, the Lagrangian should remain unchanged. For a Dirac field like
the quarks, the Lagrangian will be

$\mathcal{L}=\bar{q}i(\gamma^{\mu}\mathcal{D}_{\mu}-m)q$ (6)

where I’ve omitted the $\mathrm{SU}(3)$ index for notational simplicity, but
we should keep in mind that $q$ still refers to the $\mathrm{SU}(3)$ vector of
three quark fields.

While the mass term $\bar{q}imq$ is trivially $\mathrm{SU}(3)$ invariant (any
$\mathrm{SU}(3)$ transformation made to the quark fields $q$ cancels out with
the inverse transformation made to $\bar{q}$), the kinetic term
$\bar{q}i(\gamma^{\mu}\mathcal{D}_{\mu})q$ is not, at least not if we use the
simple partial derivative $\mathcal{D}_{\mu}=\partial_{\mu}$. Applying an
$\mathrm{SU}(3)$ transformation (4) to that term, and expanding
$\partial_{\mu}Uq$ using the product rule, produces

$i\bar{q}\herm{U}\gamma^{\mu}\partial_{\mu}Uq=i\bar{q}\gamma^{\mu}\partial_{\mu}q+i\bar{q}\herm{U}\gamma^{\mu}(\partial_{\mu}U)q$
(7)

The first term is the same one, involving the partial derivative, that we had
before the transformation, but this residual contribution
$i\bar{q}\herm{U}\gamma^{\mu}(\partial_{\mu}U)q$ spoils the gauge invariance.

To restore gauge invariance, we’ll need to add to $\mathcal{D}_{\mu}$ an extra
term $\Delta_{\mu}$ which transforms in a way that will cancel out that
residual contribution. Specifically, the product $\bar{q}\Delta_{\mu}q$ needs
to change by $-\bar{q}\herm{U}\partial_{\mu}Uq$ under an $\mathrm{SU}(3)$
transformation:

$\bar{q}\Delta_{\mu}q\to\bar{q}\herm{U}\Delta_{\mu}Uq-\bar{q}\herm{U}(\partial_{\mu}U)q$
(8)

Now, an arbitrary $\mathrm{SU}(3)$ transformation $U$ can be written as
$e^{i\alpha_{c}T^{c}}$, where $\alpha_{c}$ is an eight-component vector
(corresponding to the eight generators $T^{c}$ of $\mathrm{SU}(3)$). From
this, we can explicitly write
$\partial_{\mu}U=i(\partial_{\mu}\alpha_{c})T^{c}U$, which means the
transformation of $\Delta_{\mu}$ has to satisfy

$\bar{q}\Delta_{\mu}q\to\bar{q}\herm{U}\Delta_{\mu}Uq-\bar{q}\herm{U}i(\partial_{\mu}\alpha_{c})T^{c}Uq$
(9)

or

$\Delta_{\mu}\to\herm{U}[\Delta_{\mu}-i(\partial_{\mu}\alpha_{c})T^{c}]U$ (10)

Given this form, it makes sense to write
$\Delta_{\mu}=ig_{\text{s}}A_{c\mu}T^{c}$, and specify that the gauge field,
or gluon field for QCD specifically, transforms as

$A_{c\mu}\to A_{c\mu}-\frac{1}{g_{\text{s}}}\partial_{\mu}\alpha_{c}$ (11)

The covariant derivative then becomes

$\mathcal{D}_{\mu}=\partial_{\mu}+ig_{\text{s}}T^{c}A_{c\mu}$ (12)

where $g_{\text{s}}$ is the coupling, $T^{c}$ is an $\mathfrak{su}(3)$
generator, and $A_{c\mu}$ is a new field we will call the gluon field (or the
gauge field). In this way, the requirement of local $\mathrm{SU}(3)$ gauge
invariance, where the transformation matrix $U$ can be a function of spacetime
position while still preserving the symmetry, leads to the existence of eight
gluon fields.

As explained in standard references [17, 18, 19], we can derive a kinetic term
for the gluon fields by defining $F^{i}_{\mu\nu}$ from

$\commut{\mathcal{D}_{\mu}}{\mathcal{D}_{\nu}}=-ig_{\text{s}}F^{c}_{\mu\nu}\frac{\sigma^{c}}{2}$
(13)

or, explicitly,

$F^{c}_{\mu\nu}=\partial_{\mu}A^{c}_{\nu}-\partial_{\nu}A^{c}_{\mu}+g_{\text{s}}\epsilon^{cjk}A^{j}_{\mu}A^{k}_{\nu}$
(14)

The kinetic term, which must be a gauge-invariant function of $F$, works out
to be $-\frac{1}{4}F^{2}$. Putting this all together, we have the Lagrangian
of QCD

$\mathcal{L}=\bar{q}i(\gamma^{\mu}\mathcal{D}_{\mu}-m)q-\frac{1}{4}\bigl{(}F^{c}_{\mu\nu}\bigr{)}^{2}$
(15)

with three quark fields and eight gluon fields. Expanding the terms shows the
following types of interactions:

$\displaystyle-m\bar{q}q$ quark mass (16a) $\displaystyle-
g_{\text{s}}\bar{q}\gamma^{\mu}T^{c}A_{c\mu}q$ quark-quark-gluon vertex (16b)
$\displaystyle
g_{\text{s}}(\partial_{\mu}A^{c}_{\nu}-\partial_{\nu}A^{c}_{\mu})\epsilon^{cjk}A^{j\mu}A^{k\nu}$
triple gluon vertex (16c) $\displaystyle
g_{\text{s}}^{2}\epsilon^{cjk}A^{j\mu}A^{k\nu}\epsilon^{cj^{\prime}k^{\prime}}A^{j^{\prime}}_{\mu}A^{k^{\prime}}_{\nu}$
quadruple gluon vertex (16d)

The gluon self-interactions (16c) and (16d) are a unique feature of a
nonabelian theory; they don’t occur in quantum electrodynamics, for example.

#### Running Coupling

One feature of QCD (or any renormalizable quantum field theory) that will play
an important role in calculations is the running of the strong coupling. So
far we have been treating $g_{\text{s}}$, or equivalently
$\alpha_{\text{s}}=g_{\text{s}}^{2}/4\pi$, as a constant, but in fact it has
to vary with energy. Here is the reasoning: suppose we have some observable
quantity $O$ which depends on a single energy scale $Q$. $O$ may be
proportional to some power of $Q$, determined by its units, but any more
complex energy dependence will have to occur in the form $Q/\mu$, where $\mu$
is the renormalization scale of the theory (the scale which separates high-
energy behavior we can ignore from low-energy behavior which we cannot).
However, $O$ should be independent of $\mu$. So we need its expression to
satisfy

$\biggl{[}\mu^{2}\pd{}{\mu^{2}}+\mu^{2}\pd{\alpha_{\text{s}}}{\mu^{2}}\pd{}{\alpha_{\text{s}}}\biggr{]}O=0$
(17)

This can be rewritten as a differential equation for the coupling, which takes
the form

$\pd{\alpha_{\text{s}}}{\ln(Q^{2}/\mu^{2})}=\mu^{2}\pd{\alpha_{\text{s}}}{\mu^{2}}=-b\alpha_{\text{s}}^{2}(1+b^{\prime}\alpha_{\text{s}}+b^{\prime\prime}\alpha_{\text{s}}^{2}+\cdots)$
(18)

With this in hand, as shown in e.g. reference [19], we can take the observable
$O$, expressed in terms of the coupling at an energy scale $Q$, and write it
as a perturbation series in the coupling constant at the renormalization
scale, $\alpha_{\text{s}}(\mu)$:

$O(\alpha_{\text{s}}(Q^{2}))=O_{1}\bigl{[}1-c_{1}\alpha_{\text{s}}(\mu^{2})+c_{2}(\alpha_{\text{s}}(\mu^{2}))^{2}+\cdots\bigr{]}$
(19)

with $O=O_{1}\alpha_{\text{s}}+\cdots$. This becomes relevant when we discuss
the BFKL and BK equations, where we will interpret the running of the coupling
as providing next-to-leading order corrections to the operator that implements
the BFKL/BK evolution.

For calculations up to next-to-leading order, we can safely take the first
term on the right from equation (18) and solve it to obtain

$\alpha_{\text{s}}=\frac{1}{b\ln\frac{Q^{2}}{\mu^{2}}}$ (20)

A conventional choice for $\mu$ is the parameter $\Lambda_{\text{QCD}}$ which
has a value around $\SI{200}{MeV}$, although as a renormalization scale
parameter, one can make different, equally valid choices for $\mu$,
corresponding to a theoretical uncertainty in calculations that use the
running coupling.

#### Dressed Particles

At this point we have to make a distinction between bare particles, which are
the particles corresponding to the fields in the Lagrangian, and dressed
particles, which are the physical particles that propagate through space and
time.

Physically, we can understand the difference as follows. The vertex in
equation (16b) tells us that a bare quark, for example, will occasionally emit
gluons as it propagates. Those gluons in turn emit other gluons, and some of
these emitted gluons will turn into quark-antiquark pairs. In the end, the
quark will be surrounded by a “cloud” of copropagating particles. But the
particles from this cloud can interact with each other and be reabsorbed by
the original particle. If the relative amounts of quarks, antiquarks, and
gluons are balanced, all these emissions and reabsorptions reach an
equilibrium, and we have a whole cloud of particles propagating together, not
just a single disturbance in a quark field. This propagating cloud constitutes
a dressed quark.

Mathematically, this situation arises because excitations in the bare fields
of the QCD Lagrangian are not eigenstates of the Hamiltonian, not when
interactions are included. The eigenstates are linear combinations of bare
fields, and these linear combinations correspond to dressed particles.

When a dressed quark interacts with another particle, there is some
probability that the other particle actually interacts with a bare gluon, or a
pair of bare gluons, or a bare quark and antiquark, or even a “subcloud” that
corresponds to a dressed gluon inside the top-level dressed quark.
Accordingly, we’ll say that there is some probability of finding, say, a gluon
“inside” a quark, which will be given by the splitting wavefunction $\psi$
later in this chapter.

At a higher level, we can actually think of a proton or any other hadron as a
dressed bound state of three bare quarks. The overarching goal of hadron
structure research, as described in the following sections, is understanding
the nature of this dressing.

### Structure of Protons and Nuclei

Modern hadron structure research has its origin in 1969, with the development
of the parton model by Feynman and Björken [20, 21, 22]. This model treats a
target proton or nucleus as a collection of pointlike constituents which
incoherently interact with the projectile.

As motivation for the parton model, consider the severe Lorentz contraction
that affects protons and nuclei in a high-energy collider. The proton or
nucleus is contracted into an ellipsoid with longitudinal extent
$\frac{1}{\gamma}$ of its radius, and due to time dilation, its internal
dynamics — the branching and recombining of bare particles in a dressed proton
— are slowed by a factor of $\gamma$. For example, a target proton or nucleus
at RHIC, in the reference frame of the projectile, has
$\gamma=\frac{E}{m_{p}}=\frac{\SI{200}{GeV}}{\SI{938}{MeV}}=213$. In
$\mathrm{p}A$ collisions at the LHC, the corresponding value is
$\gamma=\frac{\SI{2.76}{TeV}}{\SI{938}{MeV}}=2940$. These values of $\gamma$
are large enough that we can take the $\gamma\to\infty$ limit, which
corresponds to analyzing the collision in the infinite momentum frame. In this
frame, the target acts as a static configuration of quantum fields on a flat
disk, and as a result, we can consider the scattering event, the hard process,
to be effectively instantaneous, equivalent to scattering off a pointlike
constituent.

Figure 1: A proton or nucleus behaves like a flat disk in a collision

#### The Parton Model and Integrated Parton Distributions

As the originators of hadron structure research did, let’s consider deep
inelastic scattering (DIS), a process which features an _electron_ projectile
impacting a target proton or nucleus and breaking it up: AX. Under the parton
model, the behavior of the target hadron in a DIS event is determined by two
pieces: the hard scattering event involving the target parton, and the
distribution of those partons within the target hadron.111This division
between hard scattering and the parton distribution is formalized in quantum
field theory by factorization theorems [23, 24, 25] which have nontrivial
proofs. However, we simply need to know that factorization works; most of the
details are beyond the scope of this dissertation, and those that are not
(e.g. the role of the factorization scale $\mu$) will be introduced later as
needed. See e.g. reference [26] for a more involved discussion of
factorization.

To understand the roles of these pieces, consider the simplified situation of
scattering off a fixed electric charge distribution $\rho(\vec{r})$. The hard
scattering is represented by the cross section for scattering off a point
elementary charge, called the Mott cross section,

$\biggl{(}\udd{\sigma}{\Omega}\biggr{)}_{\text{Mott}}=\frac{\alpha_{\text{em}}^{2}E^{2}}{4k^{4}\sin^{4}(\theta/2)}\biggl{(}1-\frac{k^{2}}{E^{2}}\sin^{2}\frac{\theta}{2}\biggr{)}$
(21)

where $E$ is the incoming particle’s energy and $k$ is the magnitude of its
momentum, and $\theta$ is the angle by which it is scattered. Now, we _could_
in principle represent the distribution of partons within the target with the
charge distribution $\rho(\vec{r})$, and then the scattering amplitude would
be a convolution of the Mott amplitude with the charge distribution:

$\mathcal{A}(\vec{r})=\iiint\mathcal{A}_{\text{Mott}}(\vec{r}-\vec{r}^{\prime})\rho(\vec{r^{\prime}})\udddc\vec{r}^{\prime}$
(22)

with $\udd{\sigma}{\Omega}\sim\mathcal{A}$. However, we already have the Mott
cross section in terms of momentum, and we would like the end result in
momentum space as well, so it’s easier to just rewrite that whole relation in
momentum space, by which the convolution becomes a simple product.

$\mathcal{A}(\vec{q})=\mathcal{A}_{\text{Mott}}(\vec{q})G_{E}(\vec{q})$ (23)

The function $G_{E}(\vec{q})$, called the electric form factor, is the Fourier
transform of the charge distribution
$\iiint\rho(\vec{r})\exp(-i\vec{q}\cdot\vec{r})\udddc\vec{r}$. in this limit
where the target’s recoil can be neglected. This simplified model gives the
overall cross section as

$\udd{\sigma}{\Omega}=\biggl{(}\udd{\sigma}{\Omega}\biggr{)}_{\text{Mott}}\abs{G_{E}(\vec{q})}^{2}$
(24)

The first factor represents the hard scattering event, and the second factor
represents the distribution of partons. It’s conventional to denote that
second factor as $F(\vec{q})\coloneqq\abs{G_{E}(\vec{q})}^{2}$, called the
structure function.

For realistic collisions, this model fails for several reasons: the partons
have magnetic moments in addition to their electric charges, spin interactions
between the projectile and target partons mean we can no longer use the Mott
cross section as the “reference” cross section for a single parton, and the
target will recoil. But the general framework, of expressing a cross section
in terms of structure functions, still works. These complications can be
accomodated by generalizing the formula for the differential cross section of
DIS to include two structure functions, $F_{1}$ and $F_{2}$.

$\frac{\uddc\sigma}{\udc x_{g}\udc
Q^{2}}=\frac{4\pi\alpha_{\text{s}}^{2}}{Q^{4}}\biggl{[}[1+(1-y)^{2}]F_{1}(x_{g},Q^{2})+\frac{1-y}{x_{g}}\bigl{[}F_{2}(x_{g},Q^{2})-2x_{g}F_{1}(x_{g},Q^{2})\bigr{]}\biggr{]}$
(25)

where $y=\frac{q\cdot p}{k\cdot p}$ in the language of figure 2. However, the
specific form of the coefficients is not as important as the fact that the
cross section can be expressed as multiples of two structure functions in this
way.

##### Collision Kinematics

At this stage I need to take a brief interlude to explain the meaning of the
kinematic variables $x_{g}$ and $Q^{2}$ which I’ve just introduced.

Figure 2: Illustration of a deep inelastic scattering event with the momentum
variables used to describe it. $k^{\mu}$ and $k^{\prime\mu}$ are the momenta
of the incoming and outgoing electron, and similarly with $p^{\mu}$ and
$p^{\prime\mu}$

When analyzing a DIS collision, we have four momenta to work with, as shown in
figure 2: the incoming and outgoing momenta of the electron, $k^{\mu}$ and
$k^{\prime\mu}$ respectively, the incoming momentum of the proton $p^{\mu}$,
and the outgoing momentum $p^{\prime\mu}$ of the assorted hadrons produced by
the breakup of the proton. Between the conservation law
($k^{\mu}+p^{\mu}=k^{\prime\mu}+p^{\prime\mu}$), the on-shell relations for
incoming and outgoing momenta ($k^{2}=k^{\prime 2}=m_{e}^{2}$,
$p^{2}=m_{p}^{2}$), and the freedom to choose a suitable reference frame, we
can eliminate all but three degrees of freedom from these
momenta.222Unfortunately, there is almost no consistency in the pedagogical
literature when it comes to parametrizing those two degrees of freedom. [27,
ch. 8] Following a derivation in any of the standard textbooks can be rather
difficult. (For elastic scattering, we can use $p^{\prime
2}=m_{\text{proton}}^{2}$ to get down to two.)

Figure 3: Snapshots of the structure of a target proton or nucleus as $x_{g}$
and $Q^{2}$ vary. The upper left corner is the saturation regime, which we
will encounter in the next chapter.

In the analysis of deep inelastic scattering, the parametrization of choice is
the photon virtuality $Q^{2}$ and the Björken variable $x_{\text{Bj}}$,
defined in terms of the momentum transfer $q^{\mu}\coloneqq
k^{\mu}-k^{\prime\mu}$:

$\displaystyle Q^{2}$ $\displaystyle\coloneqq-q^{2}$ $\displaystyle
x_{\text{Bj}}$ $\displaystyle\coloneqq\frac{Q^{2}}{2p\cdot q}$ (26)

These variables are particularly useful in combination with the parton model
because both of them have straightforward physical interpretations, at least
in leading order calculations. Suppose the photon interacts with a parton
carrying a fraction $x_{g}$ of the proton’s momentum,
$p_{q}^{\mu}=x_{g}p^{\mu}$, and let us neglect the parton’s mass. We can write
the squared four-momentum of the parton _after_ the interaction as

$\displaystyle 0\approx p_{q}^{\prime 2}=(p_{q}+q)^{2}\approx
q^{2}+2p_{q}\cdot q=-2p\cdot qx_{\text{Bj}}+2x_{g}p\cdot q=2p\cdot
q(x_{g}-x_{\text{Bj}})$ (27)

which implies that $x_{\text{Bj}}=x_{g}$: the Björken variable represents the
momentum fraction of the struck parton.

Knowing $x_{g}$ and $Q^{2}$ of a collision allows us to characterize the
behavior of the target in that collision, as shown in figure 3.

##### Integrated Parton Distributions

We now have the behavior of a target proton or nucleus in a high-energy
collision parametrized in terms of two structure functions,
$F_{1}(x_{g},Q^{2})$ and $F_{2}(x_{g},Q^{2})$. If the target consisted only of
quarks with a fixed fraction $x_{g}$ of the target’s momentum, we would find
[19, ch. 4] that

$F_{2}(x_{g},Q^{2})=2x_{g}F_{1}(x_{g},Q^{2})=\delta(1-x_{g})$ (28)

Experimentally, we know this is not the case; structure functions can be
measured and we find that they form a broad distribution. It will be
convenient to separate the contributions from each of the different flavors of
partons to the structure function, which we do by expressing it as a sum of
parton distributions,

$F_{2}(x_{g})=\sum_{\text{flavors }i}e_{i}^{2}x_{g}f_{i/h}(x_{g})$ (29)

These functions $x_{g}f_{i/h}(x_{g})$ are the integrated parton distributions
(PDFs).

PDFs cannot be completely calculated within perturbative QCD. Although the
$Q^{2}$ dependence of parton distributions is given by perturbation theory, in
the form of the DGLAP evolution equations [28, 29, 30, 31] (see also [17, 18,
19, 32, 33] for pedagogical descriptions), the initial conditions for these
equations emerge from nonperturbative physics acting at low $Q^{2}$.
Thereform, the parton distributions used in calculations are developed by
postulating a functional form, which satisfies the DGLAP equations, with
several free parameters, and fitting that function to a wide variety of
collider data. Different research groups have different methods of performing
the fit.333The Durham PDF server at http://hepdata.cedar.ac.uk/pdfs/ is a good
starting point for information about PDF fitting groups and their results. A
representative set are shown in figure 4.

Figure 4: Parton distributions $x_{g}f_{i/h}(x_{g},Q^{2})$ derived from fits
to measurements at the HERA collider. These are not the independent parton
distributions for each parton flavor, but the valence quark distributions,
$x_{g}u_{v}=x_{g}(u-\bar{u})$ and $x_{g}d_{v}=x_{g}(d-\bar{d})$ (differences
between quark and antiquark distributions of the same flavor), and the sea
quark distribution, $\sum_{i}(f_{i/h}+\bar{f}_{i/h})/2$ (average of all quark
and antiquark distributions). The gluon PDF is kept separate, and it is clear
that it is by far the largest one at small $x$. These plots were originally
published in reference [34] and are reused under the terms of the Creative
Commons Attribution 4.0 license.

At small values of $x_{g}\lesssim\num{1e-2}$, we can see from figure 4 that
the gluon distribution exceeds all others by several orders of magnitude.
Accordingly, small-$x$ physics focuses on determining the behavior of the
gluon distribution.

#### Unintegrated Parton Distributions

Factorizing cross sections into hard factors and integrated parton
distributions works well for inclusive measurements, where we measure the
total cross section, but the procedure has proved unable to explain features
of certain exclusive measurements, where the cross section is broken down by
the kinematic properties of one final-state particle. It was proposed in the
1990s [35, 36] that parton distributions which depend on the transverse
momentum transfer could provide the additional flexibility needed. Since that
time, a large body of work has been developed focusing on these transverse
momentum-dependent (TMD) distributions and their associated factorization
procedure.

These unintegrated gluon distributions can be built up from two fundamental
building blocks [37, 38], the Weizsäcker-Williams distribution [39, 40], which
measures the number density of gluons, and the dipole distribution. Both
distributions have the same asymptotic behavior at large $k_{\perp}$ [41], but
they contribute to different physical processes. This dissertation focuses on
the dipole gluon distribution specifically.

Probably the simplest expression of the unintegrated dipole gluon distribution
is the derivative of the integrated gluon distribution from the previous
section (29), [42]

$Q^{2}\mathcal{F}_{x_{g}}(x_{g},Q^{2})=\frac{2\pi^{2}\alpha_{\text{s}}}{N_{c}}\pd{}{Q^{2}}[x_{g}g(x_{g},Q^{2})]$
(30)

There are a variety of other notations used in the literature, among them
$\phi$, which is related to $\mathcal{F}_{x_{g}}$ as

$\mathcal{F}_{x_{g}}(\vec{}k_{\perp})=\frac{1}{8\pi}\lapl\phi(x_{g},\vec{}k_{\perp})+\delta^{(2)}(\vec{}k_{\perp})$
(31)

One will also occasionally see $\varphi=\frac{1}{2}\phi$.

It’s quite common to use the position space version of the distribution,
$S_{x_{g}}^{(2)}$, which is expressed as a Fourier integral

$\mathcal{F}_{x_{g}}(\vec{k}_{\perp})=\int\frac{\uddc\vec{x}_{\perp}\uddc\vec{y}_{\perp}}{(2\pi)^{2}}e^{-i\vec{k}_{\perp}\cdot(\vec{x}_{\perp}-\vec{y}_{\perp})}S_{x_{g}}^{(2)}(\vec{x}_{\perp},\vec{y}_{\perp})$
(32)

as well as the dipole amplitude

$\mathcal{N}(\vec{}x_{\perp},\vec{}y_{\perp},x_{g})=1-S_{x_{g}}^{(2)}(\vec{}x_{\perp},\vec{}y_{\perp})$
(33)

##### Analytic Models

We’ve seen that parton distributions cannot be calculated using perturbative
QCD, but we still need some concrete form for the PDF if we’re going to use it
in the calculation of a cross section. Developing and refining these concrete
forms is an active area of research. The general procedure is to postulate
some analytic expression with several free parameters that can then be
determined from fits to collider data. In this way, various researchers have
developed models for the unintegrated gluon distribution, with varying levels
of complexity.

Figure 5: Plots of the GBW and MV gluon distributions in momentum space. The
key feature of the MV distribution is the power tail at large $k_{\perp}$, as
opposed to the exponential falloff of the GBW model. These distributions use
the same saturation scale of $Q_{s}^{2}=cA^{1/3}(x_{0}/x_{g})^{\lambda}$ with
$c=0.85$, $A=208$, $x_{0}=0.000304$, and $\lambda=0.288$.

###### GBW Model

One model which is particularly common due to its simplicity was developed by
Golec-Biernat and Wüsthoff [43, 44]. The model is defined by two components.
First, an exponential form for the gluon distribution

$S_{x_{g}}^{(2)}(\vec{}r_{\perp})\coloneqq\exp\biggl{(}\frac{r_{\perp}^{2}Q_{s}^{2}}{4}\biggr{)}$
(34)

where $Q_{s}$ is the saturation scale, a momentum scale that characterizes the
transition between $S_{x_{g}}^{(2)}\approx 1$ and $S_{x_{g}}^{(2)}\ll 1$.
We’ll see the significance of the saturation scale in the next chapter. The
other component is the dependence of the saturation momentum on $x_{g}$,

$Q_{s}^{2}=Q_{0}^{2}\biggl{(}\frac{x_{0}}{x_{g}}\biggr{)}^{\lambda}$ (35)

$Q_{0}$ is the saturation momentum at some reference momentum fraction
$x_{0}$, so one of $x_{0}$ and $Q_{0}$ can be chosen arbitrarily. The choice
made by Golec-Biernat and Wüsthoff, based on a fit to DIS data collected at
HERA, is to set $Q_{0}=\SI{1}{GeV}$ for a proton, for which the corresponding
value of $x_{0}$ from the fit is $x_{0}=0.000304$.

For a nucleus, the formula for $Q_{s}^{2}$ is modified by adding a factor of
$cA^{1/3}$, where $A$ is the atomic mass number and $c$ is the centrality
coefficient, $c=\sqrt{1-R_{\perp}^{2}/R^{2}}$ with $R_{\perp}$ being the
impact parameter and $R$ being the radius of the nucleus. Effectively,
$A^{1/3}$ is proportional to the number of nucleons that fit across a diameter
of the nucleus, and thus $cA^{1/3}$ can be thought of as a measure of the
number of nucleons the projectile will encounter on its way through the
nucleus.444A more precise argument for why the saturation scale of a nucleus
scales as $A^{1/3}$ can be found in reference [32] and many papers on the
subject.

The momentum space distribution for the GBW model,

$\mathcal{F}_{x_{g}}(\vec{}k_{\perp})=\frac{1}{\pi
Q_{s}^{2}}\exp\biggl{(}-\frac{k_{\perp}^{2}}{Q_{s}^{2}}\biggr{)}$ (36)

is plotted in figure 5. It exhibits a characteristic sharp exponential falloff
at high $k_{\perp}$.

###### MV Model

The GBW model has one significant shortcoming: it falls off exponentially at
high momenta $k_{\perp}\gtrsim Q_{s}$, which conflicts with what we would
expect from rather general arguments based on perturbative QCD. To remedy
this, we can use the MV model [39, 40, 45, 46], which modifies the exponent to
produce a power law (not exponential) decay at high momentum:

$S_{x_{g}}^{(2)}(\vec{}r_{\perp})\coloneqq\exp\biggl{[}\frac{r_{\perp}^{2}Q_{s}^{2}}{4}\ln\biggl{(}e+\frac{1}{r_{\perp}\Lambda_{\text{MV}}}\biggr{)}\biggr{]}$
(37)

With this model as well, the $x_{g}$ dependence of the saturation scale
$Q_{s}$ is unconstrained. Conventionally we take it to have the same
functional dependence as in the GBW model (35).

There is no symbolic expression for the momentum space distribution
$\mathcal{F}_{x_{g}}$ in the MV model, but it can be computed numerically. The
result is plotted in figure 5. Unlike the GBW model, this gluon distribution
has a softer drop at high $k_{\perp}$, not exponential but proportional to a
power law $k_{\perp}^{-a}$ with $a\approx 4.5$.

###### AAMQS Model

Recent work by Albacete, Armesto, Milhano, Quiroga, and Salgado [47] has
generalized the preceding two models by adding an anomalous dimension,
$\gamma$.

$S_{x_{g}}^{(2)}(\vec{}r_{\perp})\coloneqq\exp\biggl{[}\frac{(r_{\perp}^{2}Q_{s}^{2})^{\gamma}}{4}\ln\biggl{(}e+\frac{1}{r_{\perp}\Lambda_{\text{MV}}}\biggr{)}\biggr{]}$
(38)

The MV model can be considered a special case of this with $\gamma=1$, and the
GBW model a further specialization with $\Lambda_{\text{MV}}\to\infty$.

In practice, the AAMQS model is usually used as an initial condition to the BK
evolution equation (to be introduced later), rather than as an independent
model of the gluon distribution in its own right. The best fit with this model
yields $\gamma$ between $\approx 0.9$ and $\approx 1.4$ depending on the
values of other parameters, suggesting that the additional degree of freedom
in the model may play a useful role in describing data. I mention it here
because it appears frequently in the literature, although I will not use it
for the fits in the remainder of this dissertation.

## The Color Glass Condensate

The previous chapter presented a very high-level overview of how the picture
of hadron structure emerges from QCD. In this chapter, I’m going to revisit
those concepts, showing the details of how they emerge from the underlying
quantum field theory. These derivations will also serve as a demonstration of
how the cross sections in later chapters can be calculated.

I’ll begin with a discussion of light cone perturbation theory, a method of
approaching quantum field theory that is commonly used in the field. The focus
of this chapter is on the color glass condensate, which is a model developed
to describe the behavior of colliding nucleons and nuclei when their gluon
distributions are large. Having discussed the CGC prepares us to analyze the
small-$x$ evolution of the unintegrated gluon distribution, given by the BK
equation, which will play a key role in the results presented in the remainder
of the dissertation.

### Light Cone Perturbation Theory

Calculations in saturation physics are often performed in light cone
perturbation theory,555For an accessible introduction to light cone
perturbation theory, including a set of rules analogous to those used for
Feynman diagrams, see reference [32], which is based on the more detailed
treatment in references [48, 49]. which is similar to, but slightly different
from, the “traditional” perturbation theory that uses Feynman diagrams.

The basis of light cone perturbation theory is the light cone coordinate
system, in which any four-vector $a^{\mu}$ is expressed as the components
$(a^{+},a^{-},\vec{a}_{\perp})$:

$\displaystyle a^{+}$ $\displaystyle=a_{t}+a_{z}$ $\displaystyle a^{-}$
$\displaystyle=a_{t}-a_{z}$ $\displaystyle\vec{a}_{\perp}$
$\displaystyle=(a_{x},a_{y})$ (39)

Light cone coordinates are described in more detail in appendix 615.

This choice of coordinates is particularly convenient to describe two
particles moving at the speed of light in opposite directions along the $z$
axis. For one particle, which we’ll take to be the projectile, the minus
component of position is constant, and the plus component can be used to
parametrize the path, filling the role of time. And of course for the target,
the reverse is true: $x^{+}$ is constant and $x^{-}$ parametrizes the path.

Light cone perturbation theory comes with its own rules for representing
physical processes by Feynman-like diagrams. For comparison, consider a
Feynman diagram in plain old $\phi^{3}$ field theory.

This diagram represents the total contribution of _all_ possible processes in
which the particle lines are connected this way, no matter where in spacetime
the vertices occur. In particular, even processes where vertex 2 is
chronologically _before_ vertex 1 ($t_{2}<t_{1}$) are included in the
diagram’s contribution. On the other hand, there is a different procedure,
time ordered perturbation theory, in which the vertices have to be ordered by
time, so we use two diagrams to represent the same process: one for all the
contributions where $t_{2}>t_{1}$, and one for all the contributions where
$t_{2}<t_{1}$.

Light cone perturbation theory is quite similar, except that instead of being
ordered by time, the vertices are ordered by the light cone time $x^{+}$. The
light cone version of the diagram on the left represents all processes where
$x_{2}^{+}>x_{1}^{+}$, and the one on the right, all processes where
$x_{2}^{+}<x_{1}^{+}$.

Coupled with this choice of coordinates and diagrams is the use of light cone
gauge for the gluon field. In general, the light cone gauge is defined as

$\eta^{\mu}A_{\mu}^{a}=0$ (40)

for some four-vector $\eta$ which is lightlike, satisfying
$\eta^{\mu}\eta_{\mu}=0$. A common choice, which I’ll be using in this work,
is

$\eta=(0,2,\vec{0})\qquad\implies\qquad A^{+}=0$ (41)

#### The Light Cone Wavefunction

A quantum state in light cone perturbation theory is a superposition of Fock
states, each of which corresponds to a definite number of particles of
specific types with definite quantum numbers. For example, a state with one
gluon with plus momentum $k^{+}$, transverse momentum $\vec{}k_{\perp}$,
polarization $\lambda$, and color indices $\vec{a}$ might be represented
$\ket{k^{+},\vec{}k_{\perp},\lambda,\vec{a}}$. More generally, let’s denote a
state with $n_{g}$ gluons and $n_{q}$ quarks, where the individual momenta are
irrelevant, as $\ket{n_{g},n_{q}}$.

All possible Fock states form a complete basis, so if we sum
$\ket{n_{g},n_{q}}\bra{n_{g},n_{q}}$ over all possible particle numbers, and
integrate over all possible momenta, polarizations, and color indices
(represented by $\udc\mathcal{PS}$, for “phase space”) at each particle
number, the resulting operator is the identity operator

$1=\sum_{n_{g},n_{q}}\int\udc\mathcal{PS}_{n_{g},n_{q}}\ket{n_{g},n_{q}}\bra{n_{g},n_{q}}$
(42)

So, if $\ket{p}$ is the state that represents, say, a proton, we can write

$\ket{p}=\sum_{n_{g},n_{q}}\int\udc\mathcal{PS}_{n_{g},n_{q}}\ket{n_{g},n_{q}}\braket{n_{g},n_{q}|p}$
(43)

This is reminiscent of the way we define a position-space wavefunction in
basic quantum mechanics:

$\ket{\psi}=\int\udddc\vec{r}\ket{\vec{r}}\braket{\vec{r}|\psi}$ (44)

In this case we would define the complex-valued function
$\psi(\vec{r})\coloneqq\braket{\vec{r}|\psi}$ as the wavefunction.
Analogously, from equation (43), we define the complex-valued function

$\Psi_{p}(n_{g},n_{q},\text{momenta})\coloneqq\braket{n_{g},n_{q}|p}$ (45)

as the light cone wavefunction. Just as the normal wavefunction represents the
amplitude of a state at a particular position, the light-cone wavefunction
represents the amplitude of a state to have a particular particle content with
specified momenta.

#### Compound Diagrams

$\displaystyle\leavevmode\hbox{\lower
41.64282pt\hbox{\includegraphics{generated-
figures/compounddiagram1}}}\quad\times\quad\leavevmode\hbox{\lower
42.16016pt\hbox{\includegraphics{generated-figures/compounddiagram2}}}\quad$
$\displaystyle=\quad\leavevmode\hbox{\lower
42.16016pt\hbox{\includegraphics{generated-figures/compounddiagram3}}}$
$\displaystyle\leavevmode\hbox{\lower
41.64282pt\hbox{\includegraphics{generated-
figures/compounddiagram4}}}\quad\times\quad\leavevmode\hbox{\lower
42.16016pt\hbox{\includegraphics{generated-figures/compounddiagram5}}}\quad$
$\displaystyle=\quad\leavevmode\hbox{\lower
42.16016pt\hbox{\includegraphics{generated-figures/compounddiagram6}}}$ Figure
6: Two examples of compound diagrams for the process

Each individual diagram in light cone perturbation theory (as in Feynman
perturbation theory) represents a contribution $\mathcal{A}_{i}$ to the
overall scattering amplitude $\mathcal{A}$.

$\mathcal{A}=\mathcal{A}_{1}+\mathcal{A}_{2}+\cdots$ (46)

When we calculate the cross section, we’re going to multiply the amplitude by
its complex conjugate, which produces a series of terms of the form
$(\text{conjugate amplitude})\times(\text{amplitude})$.

$\abs{\mathcal{A}}^{2}=(\conj{\mathcal{A}_{1}}+\conj{\mathcal{A}_{2}}+\cdots)(\mathcal{A}_{1}+\mathcal{A}_{2}+\cdots)=\conj{\mathcal{A}_{1}}\mathcal{A}_{1}+\conj{\mathcal{A}_{2}}\mathcal{A}_{2}+\cdots+\conj{\mathcal{A}_{1}}\mathcal{A}_{2}+\conj{\mathcal{A}_{2}}\mathcal{A}_{1}+\cdots$
(47)

If we adopt the convention that the mirror image of a diagram represents the
complex conjugate of the corresponding term, then each of these
$(\text{conjugate amplitude})\times(\text{amplitude})$ terms can be
represented using a “compound diagram” that is constructed from one normal and
one mirrored diagram. Figure 6 shows a couple of examples. In these compound
diagrams, a dashed line, the cut, indicates where the normal diagram ends and
the mirrored one begins, or equivalently, the intermediate state of the
compound diagram which corresponds to the final state of the normal diagram.

We can apply the rules of perturbation theory to construct the mathematical
expression corresponding to the compound diagram, and we’ll get precisely the
product of the normal and mirrored diagrams. In this way, we can use the rules
of our field theory to obtain the terms of $\abs{\mathcal{A}}^{2}$ directly,
which is often easier than finding $\mathcal{A}$ and multiplying it by its
complex conjugate.

There is one caveat: when constructing the expression that corresponds to a
compound diagram, we would ordinarily allow the momenta of the intermediate
particles to be off-shell (i.e. to not satisfy $p^{\mu}p_{\mu}=m^{2}$).
However, in reality, these “intermediate” particles are actually final-state
particles, which will necessarily have on-shell momenta. So we’ll have to
manually insert the constraint that the momentum is on-shell for each particle
line that crosses the cut.

### Calculations in the Color Glass Condensate

In roughly the past decade or so, the color glass condensate (CGC) model has
emerged as a promising candidate for describing the behavior of protons and
nuclei in high-energy collisions where the gluon distributions are large. The
model derives its name from its treatment of the target nucleon or nucleus as
a densely packed gluon fluid, or “color glass.” The appeal of the CGC is its
ability to efficiently handle situations where multiple gluon interactions
play a large role in determining the physics.

In this section I’ll introduce some of the underlying principles of CGC
calculations and demonstrate their use in calculations.

#### Wilson Lines

Much of the formalism of the CGC is based on Wilson lines. A Wilson line is
defined as the path-ordered exponential of a field along a path:

$U(\vec{}x_{\perp})=\operatorname{\mathcal{P}}\exp\biggl{(}ig_{S}\int\udc
s\ud{x^{\mu}}{s}T^{c}A_{c\mu}(x^{\nu})\biggr{)}$ (48)

Physically, this represents the accumulated effect of the field $A_{c\mu}$ on
a particle passing through it along the path parametrized by $s$, where the
fundamental interaction between the particle and the field contributes the
factor of $T^{c}$. For example, the Wilson line with $T^{c}$ in the
fundamental representation of $\mathrm{SU}(3)$ represents the superposition of
all possible multiple gluon interactions between a propagating quark and the
target nucleon or nucleus.

(49)

Alternatively, a Wilson line where the generator is in the adjoint
representation, denoted $\tilde{T}^{c}$, represents the superposition of
interactions between a propagating _gluon_ and the field of the nucleus. In
lieu of equation (48), this gives us

$W(\vec{}x_{\perp})=\operatorname{\mathcal{P}}\exp\biggl{(}ig_{S}\int\udc
s\ud{x^{\mu}}{s}\tilde{T}^{c}A_{c\mu}(x^{\nu})\biggr{)}$ (50)

Adjoint Wilson lines can be expressed in terms of fundamental Wilson lines
using the identity (609).

All the degrees of freedom in the gluon field can also be represented as
Wilson lines — in other words, any function of $A_{c\mu}$ can also be written
in terms of $U$ (or its equivalent in a different representation). But Wilson
lines become particularly appealing when we write them as Wilson loops, traces
of products of Wilson lines which form a closed path. The fact that the path
is closed, together with the trace operation, means that Wilson loops are
automatically gauge invariant, and in fact _any_ gauge-invariant quantity in
the theory can be expressed in terms of Wilson loops. [17, sec. 15.3]

#### Gluon Distributions as Wilson Line Correlators

Figure 7: Coordinates for color dipoles. This represents the splitting of a
single color dipole $\vec{}r_{\perp}$ into two, $\vec{}s_{\perp}$ and
$\vec{}t_{\perp}$.

Unintegrated parton distributions are, of course, gauge invariant quantities,
so we can express them in terms of Wilson loops. In this framework, the dipole
gluon distribution is a correlator of two Wilson lines,

$S_{x_{g}}^{(2)}(\vec{}x_{\perp},\vec{}y_{\perp})=\frac{1}{N_{c}}\braket{\trace[U(\vec{}x_{\perp})\herm
U(\vec{}y_{\perp})]}_{x_{g}}$ (51)

where $\braket{\cdots}_{x_{g}}$ denotes a “target averaging operation” [50,
sec. 2.5] which works like this: from the gluon field in the nucleus, extract
the components corresponding to gluons having longitudinal momentum fraction
$x_{g}$. Then average over all possible configurations of these components of
the field. Here $\vec{}x_{\perp}$ and $\vec{}y_{\perp}$ can be thought of as
representing the coordinates of the endpoints of a color dipole, such as a
quark-antiquark pair.

We’re also going to encounter the quadrupole gluon distribution, which is a
correlator of four Wilson lines,

$S_{x_{g}}^{(4)}(\vec{}x_{\perp},\vec{}b_{\perp},\vec{}y_{\perp})=\frac{1}{N_{c}^{2}}\braket{\trace[U(\vec{}x_{\perp})\herm
U(\vec{}b_{\perp})]\trace[U(\vec{}b_{\perp})\herm
U(\vec{}y_{\perp})]}_{x_{g}}$ (52)

As before, $\vec{}x_{\perp}$ and $\vec{}y_{\perp}$ represent the coordinates
of the endpoints of a parent color dipole, which has emitted a new dipole at
transverse coordinate $\vec{}b_{\perp}$, as shown in figure 7. This dipole
splitting process will be discussed in more detail in section Dynamics of
Small-$x$ Gluons.

As with the dipole distribution, this can be written in momentum space,

$\mathcal{G}_{x_{g}}(\vec{}k_{\perp},\vec{}k_{\perp}^{\prime})=\int\frac{\uddc\vec{}x_{\perp}\uddc\vec{}y_{\perp}\uddc\vec{}b_{\perp}}{(2\pi)^{4}}e^{-i\vec{}k_{\perp}\cdot\vec{(}x_{\perp}-\vec{}b_{\perp})}e^{-i\vec{}k_{\perp}^{\prime}\cdot(\vec{}b_{\perp}-\vec{}y_{\perp})}S_{x_{g}}^{(4)}(\vec{}x_{\perp},\vec{}b_{\perp},\vec{}y_{\perp})$
(53)

#### Constructing Diagrams in the CGC

Calculations using the CGC start with the corresponding diagrams in light cone
perturbation theory. In most cases, they involve at least some combination of
the following three components:

  * •

the splitting wavefunction $\psi(p^{+},k^{+},\vec{}r_{\perp})$, which is a
light cone wavefunction giving the amplitude for a dressed quark with plus-
component momentum $p_{+}$ to contain a gluon with plus-component momentum
$k^{+}\coloneqq(1-\xi)p^{+}$ and a dressed quark with plus-component momentum
$\xi p^{+}$ with transverse separation $\vec{}r_{\perp}$. Effectively, this
represents the emission of the gluon from the quark.

  * •

the fundamental Wilson line $U(\vec{}x_{\perp})$, from equation (48), which
represents the interaction between a quark at transverse position
$\vec{}x_{\perp}$ and the chromodynamic field of the nucleus. Again, this can
be modeled as the superposition of all possible numbers of gluons interacting
with the quark.

  * •

the adjoint Wilson line $W(\vec{}y_{\perp})$, from equation (50), which is the
same thing with a gluon as the propagating particle.

As an example of how to put these pieces together, consider the process AA, in
which a quark from the projectile interacts with the target nucleus and emits
a gluon. These two events could happen in either order. If the interaction
with the gluon field of the nucleus comes first, the process is represented by
the diagram

This will be represented by the factor

$\psi(p^{+},k^{+},\vec{}r_{\perp})U(\vec{}x_{\perp})$ (54)

Alternatively, the interaction could take place after the emission of the
gluon, in which case the chromodynamic field interacts with _both_ the quark
and the emitted gluon.

We can think of this as the construction shown in figure 8: a superposition of
all different numbers of nuclear gluons that could interact with the
propagating partons. Mathematically, the superposition in figure 8 corresponds
to the expression of the form

$U_{0}W_{0}+U_{1}W_{0}+U_{0}W_{1}+U_{2}W_{0}+U_{1}W_{1}+U_{0}W_{2}+\cdots=UW$
(55)

So the interaction of the nuclear chromodynamic field with a quark and gluon
propagating together is the product of the fundamental and adjoint Wilson
lines. Accordingly, the expression corresponding to this diagram is

$\psi(p^{+},k^{+},\vec{}x_{\perp}-\vec{}b_{\perp})U(\vec{}x_{\perp})W(\vec{}b_{\perp})$
(56)

where $\vec{}x_{\perp}-\vec{}b_{\perp}$ appears in the splitting wavefunction
because it is the transverse separation between the fundamental and adjoint
Wilson lines.

Figure 8: Decomposition of the interaction of a quark and gluon with the
chromodynamic field of a target

#### Impact Parameter Independence

A significant simplifying assumption, made in many studies of gluon
distributions, is that the gluon field of the nucleus has an infinite extent
in the transverse plane, and is homogeneous except for local fluctuations. For
these gluon distributions, this corresponds to the statement that the Wilson
line correlators depend only on
$\vec{}r_{\perp}\coloneqq\vec{}x_{\perp}-\vec{}y_{\perp}$ but are independent
of the impact parameter
$\vec{}R_{\perp}\coloneqq\frac{1}{2}(\vec{}x_{\perp}+\vec{}y_{\perp})$.666References
[51, 52, 53, 54, 55] numerically investigate the consequences of dropping this
assumption. If we make this assumption, the dipole distribution can be
expressed as a function of $\vec{}r_{\perp}$ only, and its Fourier transform
(32) can be simplified to

$\mathcal{F}_{x_{g}}(\vec{}k_{\perp})=\int\frac{\uddc\vec{}R_{\perp}}{4}\int\frac{\uddc\vec{}r_{\perp}}{(2\pi)^{2}}e^{-i\vec{}k_{\perp}\cdot\vec{}r_{\perp}}S_{x_{g}}^{(2)}(\vec{}r_{\perp})=S_{\perp}\int\frac{\uddc\vec{}r_{\perp}}{(2\pi)^{2}}e^{-i\vec{}k_{\perp}\cdot\vec{}r_{\perp}}S_{x_{g}}^{(2)}(\vec{}r_{\perp})$
(57)

where $S_{\perp}$ is the transverse area of the target over which the
interaction takes place. We can express the inverse relation as

$S_{x_{g}}^{(2)}(\vec{}r_{\perp})=\frac{1}{S_{\perp}}\int\uddc\vec{}k_{\perp}e^{i\vec{}k_{\perp}\cdot\vec{}r_{\perp}}\mathcal{F}_{x_{g}}(\vec{}k_{\perp})$
(58)

For the quadrupole distribution, we can do the same thing with
$\vec{}s_{\perp}\coloneqq\vec{}x_{\perp}-\vec{}b_{\perp}$ and
$\vec{}t_{\perp}\coloneqq\vec{}y_{\perp}-\vec{}b_{\perp}$, where
$\vec{}b_{\perp}$ is the transverse coordinate of the emitted color pole). The
resulting inverse transform is

$S_{x_{g}}^{(4)}(\vec{}s_{\perp},\vec{}t_{\perp})=\frac{1}{S_{\perp}}\int\uddc\vec{}k_{\perp}\uddc\vec{}k_{\perp}^{\prime}e^{i\vec{}k_{\perp}\cdot\vec{}s_{\perp}}e^{i\vec{}k_{\perp}^{\prime}\cdot\vec{}t_{\perp}}\mathcal{G}_{x_{g}}(\vec{}k_{\perp},\vec{}k_{\perp}^{\prime})$
(59)

Assuming impact parameter independence also allows us to use the mean field
approximation, which we’ll encounter later (83), to factor the momentum space
quadrupole distribution as follows:

$\mathcal{G}_{x_{g}}(\vec{}k_{\perp},\vec{}k_{\perp}^{\prime})=\mathcal{F}_{x_{g}}(\vec{}k_{\perp})\mathcal{F}_{x_{g}}(\vec{}k_{\perp}^{\prime})$
(60)

### Dynamics of Small-$x$ Gluons

Figure 9: A “ladder diagram” — a compound diagram showing the interaction of
gluons between two quarks. This can also be interpreted as in figure 10, as
emission of a gluon from one of the quarks, where the gluon branches into
other gluons with progressively smaller momentum fractions on its way to
interacting with the other quark.

Even though parton distributions are fundamentally nonperturbative, it is
possible to determine their behavior in high-energy collisions using
perturbative QCD. The effort to do so begins with the work of Balitsky, Fadin,
Kuraev, and Lipatov in the late 1970s. Their original derivations [3, 4, 5, 6,
7] focused on the behavior of cross sections, not specifically the
unintegrated parton distributions, but the underlying physics is the same in
both cases. What they discovered is that the high-energy behavior of hadron
scattering cross sections is primarily a result of ladder diagrams, of the
type shown in figure 9. These diagrams give contributions proportional to
$(\alpha_{\text{s}}\ln x_{g})^{n}$. Even though $\alpha_{\text{s}}$ is small
in these high-energy collisions, for sufficiently small values of $x_{g}$ the
combination $\alpha_{\text{s}}\ln x_{g}$ can be large, so that a leading order
($n=0$) or next-to-leading order ($n=0,1$) approximation to the series won’t
be sufficient. We need a resummation of these logarithmic contributions — that
is, a single expression giving the sum of the entire infinite series of
diagrams — to produce a good approximation to the cross section.

Figure 10: A color dipole cascade, which is another way of looking at one side
of a ladder diagram (figure 9). As the gluon (the color dipole) branches,
moving from left to right, the lowest gluon has progressively smaller
fractions $x_{g}$ of the target’s longitudinal momentum. Each circle at the
bottom shows a view, in the transverse plane, of the color dipole cascade at
the corresponding value of $x_{g}$ indicated by the dashed line through the
top part of the diagram.

This resummation is provided, in a sense, by the BFKL equation, which is often
written in terms of the dipole scattering amplitude (33) as

$\pd{\mathcal{N}(\vec{}r_{\perp},x_{g})}{\ln(1/x_{g})}=\frac{\alpha_{\text{s}}N_{c}}{2\pi^{2}}\int\uddc\vec{}b_{\perp}\frac{r_{\perp}^{2}}{s_{\perp}^{2}t_{\perp}^{2}}\bigl{[}\mathcal{N}(\vec{}s_{\perp},x_{g})+\mathcal{N}(\vec{}t_{\perp},x_{g})-\mathcal{N}(\vec{}r_{\perp},x_{g})\bigr{]}$
(61)

(this is, properly speaking, only the leading logarithmic, angle-independent
approximation to the full BFKL equation). Intuitively, this can be interpreted
as the color dipole splitting process pictured in figure 10. As we step toward
smaller values of the gluon momentum fraction $x_{g}$, an existing color
dipole of size $\vec{}r_{\perp}$ (corresponding to a gluon, or a quark-
antiquark pair) disappears and two new ones of sizes $\vec{}s_{\perp}$ and
$\vec{}t_{\perp}$ take its place.

#### Unitarity Violation

The BFKL equation describes the behavior of parton distributions at small
$x_{g}$, corresponding to high collision energy $s$, under a range of
kinematic conditions, but it’s not difficult to see that it cannot be a
_complete_ description. It’s possible to solve the equation for $\mathcal{N}$
in terms of the eigenfunctions of the integral kernel, which are simply powers
of the dipole size, i.e.

$\int\uddc\vec{}b_{\perp}\frac{r_{\perp}^{2}}{s_{\perp}^{2}t_{\perp}^{2}}\bigl{[}s_{\perp}^{\gamma}+t_{\perp}^{\gamma}-r_{\perp}^{\gamma}\bigr{]}=2\pi\bigl{[}\underbrace{2\psi(1)-\psi(\gamma)-\psi(1-\gamma)}_{\chi(\gamma)}\bigr{]}r_{\perp}^{\gamma}$
(62)

where $\psi(x)=\ud{\ln\Gamma(x)}{x}$. A general solution of (61) is a linear
combination of these functions,

$\mathcal{N}(\vec{}r_{\perp},x_{g})=\int_{-\infty}^{\infty}\udc\nu
C(\nu,x_{g})r_{\perp}^{\gamma}$ (63)

and plugging this into (61) gives us the $x_{g}$ evolution of the
coefficients:

$\pd{C(\gamma,x_{g})}{\ln(1/x_{g})}=\frac{\alpha_{\text{s}}N_{c}}{\pi}C(\gamma,x_{g})\chi(\gamma)\qquad\implies\qquad
C(\gamma,x_{g})=\biggl{(}\frac{1}{x_{g}}\biggr{)}^{\alpha_{\text{s}}N_{c}\chi(\gamma)/\pi}$
(64)

So as $x_{g}$ gets smaller, $C(\gamma,x_{g})$ grows exponentially. The largest
contribution to this exponential growth comes from $\gamma=\frac{1}{2}$, where
$\chi(\gamma)=4\ln 2$. When we run through the mathematical machinery to
convert $\mathcal{N}$ into a cross section, we can expect that, in the limit
of large collision energy $s$, the cross section will grow as $\sigma\sim
s^{4\ln 2\alpha_{\text{s}}N_{c}/\pi-1}$. [32]

This growth conflicts with an elementary result of quantum field theory called
the Froissart-Martin bound [8, 9, 10], which says that the unitarity of the
S-matrix (equivalent to conservation of probability) implies that cross
sections cannot grow any faster than $(\ln s)^{2}$ in the high-energy limit.
So there must be at least another term in the equation that acts to limit the
growth of the solution as $\mathcal{N}$ becomes large.

One can solve this problem by incorporating a negative nonlinear term which
can be interpreted as representing multiple scatterings of the probe off
target partons, or as interactions between the partons in the target prior to
the scattering. This solution was first proposed by Gribov, Levin, and Ryskin
[56] (reviewed in [57]), and eventually developed into the Balitsky-Kovchegov
(BK) equation [11, 12, 13, 14, 15, 16]:

$\pd{\mathcal{N}(\vec{}r_{\perp},x_{g})}{\ln(1/x_{g})}=\frac{\alpha_{\text{s}}N_{c}}{2\pi^{2}}\int\uddc\vec{}b_{\perp}\frac{r_{\perp}^{2}}{s_{\perp}^{2}t_{\perp}^{2}}\bigl{[}\mathcal{N}(\vec{}s_{\perp},x_{g})+\mathcal{N}(\vec{}t_{\perp},x_{g})-\mathcal{N}(\vec{}r_{\perp},x_{g})-\mathcal{N}(\vec{}s_{\perp},x_{g})\mathcal{N}(\vec{}t_{\perp},x_{g})\bigr{]}$
(65)

The nonlinear term in this equation limits the growth of $\mathcal{N}$ as
$x_{g}$ becomes small. Because of that term,
$\mathcal{N}(\vec{}r_{\perp},x_{g})=1$ is a fixed point of this equation, so
the nonlinear evolution keeps $\mathcal{N}<1$ which is sufficient to preserve
the Froissart bound.

#### Saturation

The BK equation exhibits an interesting feature called geometric scaling,
which we can see most easily by rewriting the equation in momentum space. In
terms of the momentum space distribution $\phi$, the leading logarithmic BK
equation takes the form [16]

$\pd{\phi(k_{\perp},x_{g})}{\ln(1/x_{g})}=\frac{\alpha_{\text{s}}N_{c}}{\pi}\biggl{[}\chi\biggl{(}-\pd{}{\ln
k_{\perp}^{2}}\biggr{)}\phi(k_{\perp})-\phi(k_{\perp})^{2}\biggr{]}$ (66)

We have seen that the high-energy behavior of the BK solution is dominated by
$\chi(\gamma)$ near $\gamma=\frac{1}{2}$, so we can substitute in the
expansion

$\chi(\gamma)=4\ln
2+\frac{\chi_{2}}{2}\biggl{(}\gamma-\frac{1}{2}\biggr{)}^{2}+\orderof{\gamma^{4}}$
(67)

where $\chi_{2}\coloneqq\chi\bigl{(}\frac{1}{2}\bigr{)}$, to write (66) as
[58, 59, 60]

$\pd{\phi(k_{\perp},x_{g})}{\ln(1/x_{g})}=\frac{\alpha_{\text{s}}N_{c}}{\pi}\biggl{[}(4\ln
2)\phi(k_{\perp},x_{g})+\frac{\chi_{2}}{2}\biggl{(}\pd{}{\ln
k_{\perp}^{2}}+\frac{1}{2}\biggr{)}^{2}\phi(k_{\perp},x_{g})-\phi(k_{\perp},x_{g})^{2}\biggr{]}$
(68)

This general form is known as the Fisher-Komolgorov-Petrov-Piscounov (FKPP)
equation [61, 62]. Its asymptotic solutions, in the limit of large
$\ln\frac{1}{x_{g}}$, can be written as a one-variable function
$k_{\perp}^{2}/Q_{s}^{2}$, where the saturation scale $Q_{s}$ marks the
transition between the high-momentum region where $\phi$ is small, and the
low-momentum region where it is large. For a proton, the dependence of this
variable on $x_{g}$ works out to [58]

$Q_{s}^{2}\biggl{(}\ln\frac{1}{x_{g}}\biggr{)}=Q_{0}^{2}\biggl{(}\ln\frac{1}{x_{g}}\biggr{)}^{\frac{3}{2(1-\gamma_{0})}}\exp\biggl{(}-\frac{\alpha_{\text{s}}N_{c}}{\pi}\frac{\chi(\gamma_{0})}{1-\gamma_{0}}\ln\frac{1}{x_{g}}\biggr{)}$
(69)

with $\gamma_{0}=1-\frac{1}{2}\sqrt{1+32\ln 2/\chi_{2}}$ and $Q_{0}$ some
constant. Geometric scaling is the observation that $\phi$ depends on
$k_{\perp}$ and $x_{g}$ only through the combination
$k_{\perp}^{2}/Q_{s}^{2}$. This scaling behavior is also exhibited by deep
inelastic scattering data from HERA [63], supporting the idea that the BK
equation is physically relevant to these processes.

For a nucleus, recall that $Q_{s}^{2}$ is enhanced by a factor of $cA^{1/3}$,
where $c$ is the centrality parameter of the collision and $A$ is the mass
number of the nucleus.

#### Dipole Derivation

Having explored the properties of BFKL and BK evolution, let’s examine its
derivation in more detail. A commonly referenced derivation of the BFKL
equation777In most introductory sources, such as [33, 19, 32], the BFKL
equation is presented as a generalization of the DGLAP equation. However, this
dissertation only concerns BFKL dynamics, and the derivation I present stands
on its own, so I’ve omitted any discussion of DGLAP evolution. comes from the
color dipole model, proposed by Mueller in reference [64]. In this section
I’ll present a somewhat simplified argument from reference [50], which is
nevertheless equivalent to the dipole model derivation and illustrates its
essential features.

Figure 11: Leading order diagram for interaction of a photon with a nucleus

Consider the process in which a photon fluctuates into a quark-antiquark pair
which interacts with the gluon field of the nucleus and then turns back into a
photon. The leading order diagram for this process is shown in figure 11.
Following the rules of light cone perturbation theory, we construct the
corresponding amplitude by starting at an arbitrary point on the quark loop
and following it around, picking up a factor for each interaction we encounter
as we go:

  * •

The splitting of the photon into the quark-antiquark pair contributes a factor
of the photon splitting wavefunction
$\psi_{li}(p^{+},k^{+},\vec{}x_{\perp}-\vec{}y_{\perp})$, where $k^{+}$ is the
longitudinal momentum of the quark (or we could choose it to be the antiquark)

  * •

The interaction of the upper quark line with the gluon field of the nucleus
contributes the Wilson line $U_{ij}(\vec{}x_{\perp})$

  * •

The reforming of the photon from the quark-antiquark pair contributes the
conjugated photon splitting wavefunction
$\conj{\psi_{jk}}(p^{+},k^{+},\vec{}x_{\perp}-\vec{}y_{\perp})$

  * •

The interaction of the lower quark line with the gluon field contributes the
conjugate Wilson line $\herm U_{kl}(\vec{}y_{\perp})$

The color indices on each factor come from the two lines that connect to it.
We don’t include propagators for the photons because they are external lines.
Putting everything together, we get

$\psi_{li}(p^{+},k^{+},\vec{}r_{\perp})U_{ij}(\vec{}x_{\perp})\conj{\psi_{jk}}(p^{+},k^{+},\vec{}r_{\perp})\herm
U_{kl}(\vec{}y_{\perp})$ (70)

Since photons are color-neutral, their splitting wavefunctions are diagonal in
the color indices, $\psi_{ij}\propto\delta_{ij}$, so we can simplify this to

$\abs{\psi(p^{+},k^{+},\vec{}r_{\perp})}^{2}\trace[U(\vec{}x_{\perp})\herm
U(\vec{}y_{\perp})]$ (71)

Figure 12: Next-to-leading order diagrams contributing to the interaction of a
color dipole with the gluon field of a nucleus. In the top set of diagrams,
the emitted gluon does not interact with the nucleus, whereas in the bottom
set of diagrams it does, adding an adjoint Wilson line to the amplitude. In
the small diagrams on the right, the nucleus itself is not shown to save
space.

The next-to-leading order corrections to this are represented by diagrams with
a gluon line, shown in figure 12. Each emission or absorption of a gluon
brings in a vertex factor of
$2g_{\text{s}}T^{a}\frac{\vec{\epsilon}_{\lambda}\cdot\vec{}q_{\perp}}{q_{\perp}^{2}}$,
where $\vec{}q_{\perp}$ is the transverse momentum difference between the
quark and gluon; however, to be consistent with our other expressions we
should exchange transverse momentum for transverse position. Performing the
Fourier transform gets us to the splitting factor [50]

$\int\frac{\uddc\vec{}q_{\perp}}{(2\pi)^{2}}e^{i\vec{}q_{\perp}\cdot\vec{}r_{\perp}}2g_{\text{s}}T^{a}\frac{\vec{\epsilon}_{\lambda}\cdot\vec{}q_{\perp}}{q_{\perp}^{2}}=\frac{2ig_{\text{s}}}{2\pi}T^{a}\frac{\vec{\epsilon}_{\lambda}\cdot\vec{}r_{\perp}}{r_{\perp}^{2}}$
(72)

It’s also necessary to sum over gluon polarizations using equation (658).

Putting all this together, the top left diagram in figure 12 gives the
expression

$\int^{k^{+}}\frac{\udc
q^{+}}{q^{+}}\int\uddc\vec{}b_{\perp}\psi_{ni}(p^{+},k^{+},\vec{}r_{\perp})\biggl{(}\frac{2ig_{\text{s}}}{2\pi}T_{ij}^{a}\frac{\vec{\epsilon}_{\lambda}\cdot\vec{}s_{\perp}}{s_{\perp}^{2}}\biggr{)}\\\
\times\biggl{(}-\frac{2ig_{\text{s}}}{2\pi}T_{jk}^{a}\frac{\vec{\epsilon}_{\lambda}\cdot\vec{}s_{\perp}}{s_{\perp}^{2}}\biggr{)}U_{kl}(\vec{}x_{\perp})\conj{\psi_{lm}}(p^{+},k^{+},\vec{}r_{\perp})\herm{U_{mn}}(\vec{}y_{\perp})$
(73)

where all repeated indices are summed over. The upper limit on the $q^{+}$
integral comes from the constraint that the gluon momentum can’t be greater
than that of the quark it was emitted from.

Performing the integral over $q^{+}$ gives $\ln k^{+}+\text{constant}$. Now,
the quark momentum $k^{+}$ is related to the momentum of the target gluons it
interacts with: high-momentum quarks will interact with high-momentum gluons,
and same for low momentum. We can therefore conclude that $k^{+}\propto
x_{g}$, and rewrite $\ln k^{+}+\text{constant}$ as
$-\ln(1/x_{g})+\text{constant}$. A suitable shift in reference frame then
allows us to eliminate the constant.

Taking into account that $\psi_{ij}\propto\delta_{ij}$, and that the
polarizations are orthogonal to the gluon momentum so that the last term in
(658) vanishes, this becomes

$-\frac{4\alpha_{\text{s}}}{\pi}\ln\frac{1}{x_{g}}\abs{\psi(\vec{}r_{\perp})}^{2}\trace[T^{a}T^{a}U(\vec{}x_{\perp})\herm{U}(\vec{}y_{\perp})]\int\frac{\uddc\vec{}b_{\perp}}{(2\pi)^{2}}\frac{1}{s_{\perp}^{2}}$
(74)

where $\alpha_{\text{s}}=\frac{g_{\text{s}}^{2}}{4\pi}$ and
$\vec{}s_{\perp}=\vec{}x_{\perp}-\vec{}b_{\perp}$. By a similar procedure, the
bottom left diagram gives

$\int^{k^{+}}\frac{\udc
q^{+}}{q^{+}}\int\uddc\vec{}b_{\perp}\psi_{ni}(p^{+},k^{+},\vec{}r_{\perp})\biggl{(}\frac{2ig_{\text{s}}}{2\pi}T_{ij}^{a}\frac{\vec{\epsilon}_{\lambda}\cdot\vec{}s_{\perp}}{s_{\perp}^{2}}\biggr{)}\\\
\times
U_{jk}(\vec{}x_{\perp})\conj{\psi_{kl}}(p^{+},k^{+},\vec{}r_{\perp})\biggl{(}-\frac{2ig_{\text{s}}}{2\pi}T_{lm}^{b}\frac{\vec{\epsilon}_{\lambda}\cdot\vec{}t_{\perp}}{t_{\perp}^{2}}\biggr{)}\herm{U_{mn}}(\vec{}y_{\perp})W^{ab}(\vec{}b_{\perp})$
(75)

which reduces to

$-\frac{4\alpha_{\text{s}}}{\pi}\ln\frac{1}{x_{g}}\abs{\psi(\vec{}r_{\perp})}^{2}\trace[T^{a}U(\vec{}x_{\perp})T^{b}\herm{U}(\vec{}y_{\perp})]\int\frac{\uddc\vec{}b_{\perp}}{(2\pi)^{2}}\frac{\vec{}s_{\perp}\cdot\vec{}t_{\perp}}{s_{\perp}^{2}t_{\perp}^{2}}W^{ab}(\vec{}b_{\perp})$
(76)

When all seven diagrams in figure 12 are taken into account, the total NLO
contribution, from explicit gluon emission, works out to

$-\frac{\alpha_{\text{s}}}{2\pi^{2}}\ln\frac{1}{x_{g}}\abs{\psi(\vec{}r_{\perp})}^{2}\int\uddc\vec{}b_{\perp}\frac{r_{\perp}^{2}}{s_{\perp}^{2}t_{\perp}^{2}}\Bigl{\\{}N_{c}\trace[U(\vec{}x_{\perp})\herm
U(\vec{}y_{\perp})]\\\ -\trace[U(\vec{}x_{\perp})\herm
U(\vec{}b_{\perp})]\trace[U(\vec{}b_{\perp})\herm
U(\vec{}y_{\perp})]\Bigr{\\}}$ (77)

Recall that the leading order contribution (71) was

$\abs{\psi(p^{+},k^{+},\vec{}r_{\perp})}^{2}\trace[U(\vec{}x_{\perp})\herm
U(\vec{}y_{\perp})]$

This expression includes the contribution from emission of gluons with
momentum fractions _less_ than $x_{g}$, which are resummed by the Wilson
lines. Therefore, the incremental contribution from emission of gluons with
momentum fraction _equal_ to $x_{g}$, within a range
$\Delta\bigl{[}\ln\frac{1}{x_{g}}\bigr{]}$, will be

$\abs{\psi(p^{+},k^{+},\vec{}r_{\perp})}^{2}\pd{\trace[U(\vec{}x_{\perp})\herm
U(\vec{}y_{\perp})]}{\ln(1/x_{g})}\Delta\biggl{[}\ln\frac{1}{x_{g}}\biggr{]}$
(78)

However, this is precisely the same quantity represented by the single gluon
emission contributions (77), if we take the difference between two closely
spaced values of $x_{g}$:

$-\frac{\alpha_{\text{s}}}{2\pi^{2}}\Delta\biggl{[}\ln\frac{1}{x_{g}}\biggr{]}\abs{\psi(\vec{}r_{\perp})}^{2}\int\uddc\vec{}b_{\perp}\frac{r_{\perp}^{2}}{s_{\perp}^{2}t_{\perp}^{2}}\Bigl{\\{}N_{c}\trace[U(\vec{}x_{\perp})\herm
U(\vec{}y_{\perp})]\\\ -\trace[U(\vec{}x_{\perp})\herm
U(\vec{}b_{\perp})]\trace[U(\vec{}b_{\perp})\herm
U(\vec{}y_{\perp})]\Bigr{\\}}$ (79)

Equating expressions (78) and (79), we find the following evolution equation
for the Wilson loop:

$\pd{\trace[U(\vec{}x_{\perp})\herm
U(\vec{}y_{\perp})]}{\ln(1/x_{g})}=-\frac{\alpha_{\text{s}}}{2\pi^{2}}\int\uddc\vec{}b_{\perp}\frac{r_{\perp}^{2}}{s_{\perp}^{2}t_{\perp}^{2}}\Bigl{\\{}N_{c}\trace[U(\vec{}x_{\perp})\herm
U(\vec{}y_{\perp})]\\\ -\trace[U(\vec{}x_{\perp})\herm
U(\vec{}b_{\perp})]\trace[U(\vec{}b_{\perp})\herm
U(\vec{}y_{\perp})]\Bigr{\\}}$ (80)

##### CGC Target Averaging

Technically, this equation only applies to one specific quantum state of the
gluon field of the nucleus. But we can’t prepare a nucleus in a specific
state; in a collider we get an average over all possible states. So to get an
equation that describes physically measurable results, we need to use the
target averaging procedure $\braket{\cdots}_{x_{g}}$ to get

$\pd{\braket{\trace[U(\vec{}x_{\perp})\herm
U(\vec{}y_{\perp})]}_{x_{g}}}{\ln(1/x_{g})}=-\frac{\alpha_{\text{s}}}{2\pi^{2}}\int\uddc\vec{}b_{\perp}\frac{r_{\perp}^{2}}{s_{\perp}^{2}t_{\perp}^{2}}\Bigl{\\{}N_{c}\braket{\trace[U(\vec{}x_{\perp})\herm
U(\vec{}y_{\perp})]}_{x_{g}}\\\ -\braket{\trace[U(\vec{}x_{\perp})\herm
U(\vec{}b_{\perp})]\trace[U(\vec{}b_{\perp})\herm
U(\vec{}y_{\perp})]}_{x_{g}}\Bigr{\\}}$ (81)

or, in terms of the dipole (51) and quadrupole (52) gluon distributions,

$\pd{S_{x_{g}}^{(2)}(\vec{}r_{\perp})}{\ln(1/x_{g})}=-\frac{\alpha_{\text{s}}N_{c}}{2\pi^{2}}\int\uddc\vec{}b_{\perp}\frac{r_{\perp}^{2}}{s_{\perp}^{2}t_{\perp}^{2}}\Bigl{[}S_{x_{g}}^{(2)}(\vec{}r_{\perp})-S_{x_{g}}^{(4)}(\vec{}s_{\perp},\vec{}t_{\perp})\Bigr{]}$
(82)

This differential equation will allow us to compute the evolution of
$S_{x_{g}}^{(2)}$ in $x_{g}$, provided that we can compute $S_{x_{g}}^{(4)}$.
So we might be inspired to look for a corresponding equation for
$S_{x_{g}}^{(4)}$. This equation does exist, but it expresses the evolution of
$S_{x_{g}}^{(4)}$ in terms of $S_{x_{g}}^{(6)}$, a correlator of two
fundamental and _two_ adjoint Wilson lines, $\trace[U\herm{U}W\herm{W}]$. And
the evolution equation for $S_{x_{g}}^{(6)}$ would be expressed in terms of
$\trace[U\herm{U}WW\herm{W}]$, and so on. This series of equations involving
progressively larger Wilson loops is called the Balitsky-JIMWLK hierarchy [11,
65, 66, 67, 68, 69].

In practice, to avoid the full complexity of the Balitsky-JIMWLK hierarchy, we
can use the mean field approximation, which is simply the assumption that the
two Wilson loops in the quadrupole gluon distribution are uncorrelated, and
thus the target average of the product,
$\braket{\trace[U\herm{U}]\trace[U\herm{U}]}_{x_{g}}$ factors into the product
of two averages,
$\braket{\trace[U\herm{U}]}_{x_{g}}\braket{\trace[U\herm{U}]}_{x_{g}}$. The
validity of this assumption in the limit of large $N_{c}$ is supported by
numerical studies [70, 71, 72]. This allows us to reduce the quadrupole gluon
distribution as follows:

$S_{x_{g}}^{(4)}(\vec{}x_{\perp},\vec{}b_{\perp},\vec{}y_{\perp})=S_{x_{g}}^{(2)}(\vec{}x_{\perp},\vec{}b_{\perp})S_{x_{g}}^{(2)}(\vec{}b_{\perp},\vec{}y_{\perp})$
(83)

so that the evolution equation is a self-contained differential equation for
the dipole distribution

$\pd{S_{x_{g}}^{(2)}(\vec{}r_{\perp})}{\ln(1/x_{g})}=-\frac{\alpha_{\text{s}}N_{c}}{\pi}\int\uddc\vec{}b_{\perp}\frac{r_{\perp}^{2}}{s_{\perp}^{2}t_{\perp}^{2}}\Bigl{[}S_{x_{g}}^{(2)}(\vec{}r_{\perp})-S_{x_{g}}^{(2)}(\vec{}s_{\perp})S_{x_{g}}^{(2)}(\vec{}t_{\perp})\Bigr{]}$
(84)

This is called the Balitsky-Kovchegov (BK) equation [11, 12, 13, 14, 15, 16]
in the leading logarithmic (LL) approximation.

The designation as LL reflects the fact that there is a single factor of
$\alpha_{\text{s}}$ for each factor of $\ln\frac{1}{x_{g}}$. One can go beyond
this and incorporate subleading terms which are proportional to
$\alpha_{\text{s}}\biggl{(}\alpha_{\text{s}}\ln\frac{1}{x_{g}}\biggr{)}^{n}$,
producing the next-to-leading logarithmic (NLL) BK equation [73, 74, 75],
though the numerical implementation of the full NLL BK evolution is quite
complex. Some progress has been made toward such an implementation but it is
not yet ready for use in phenomenology. [76]

It turns out that incorporating selected NLL corrections, specifically those
which result from the running of the QCD coupling, is expected to give a
reasonable approximation to the full NLL evolution. Many numerical studies,
including those to be discussed in the remainder of this dissertation, use
this rcBK (BK with running coupling) gluon distribution. A widely used sample
implementation is available in e.g. reference [77].

#### Numerical Analysis of BK Evolution

Figure 13: The gluon distribution from the GBW model (top), and the output of
the numerical BK evolution with fixed coupling (middle) and running coupling
(bottom), for central ($c=0.85$) $\mathrm{pPb}$ ($A=208$) collisions. Each
curve shows $\mathcal{F}_{x_{g}}(k_{\perp})$ for
$x_{g}=x_{\text{init}}\exp(-n\delta Y)$, where $n$ is the number that labels
the curve and $\delta Y=\frac{1}{400}\ln(10^{8})\approx 0.04605$ is the step
size in rapidity used by the numerical integrator. The curve labeled $0$ is
the initial condition. Notable features include the transient “spike” at the
low $k_{\perp}^{2}$ region of the solution to the BK evolution (which is a
numerical artifact), and the significant enhancement at large $k_{\perp}^{2}$
relative to the GBW model. This figure is adapted from reference [78].

Figure 14: Analogous to figure 13, these plots show the dipole gluon
distribution in momentum space, $\mathcal{F}_{x_{g}}$, computed by solving the
equation in position space (65) and numerically Fourier transforming it into
momentum space. The top plot shows the analytic result from the GBW model, the
middle one shows the result from the BK equation with fixed coupling
$\alpha_{\text{s}}=0.1$ (65), and the bottom one shows the result with running
coupling corrections incorporated. Each curve in this plot shows
$\mathcal{F}_{x_{g}}(k_{\perp})$ for $x_{g}=x_{\text{init}}\exp(-n)$, where
$n$ is the number that labels the curve, so the numbers labeling the curves
are not directly comparable to those in figure 13. In addition, these curves
are for minimum bias ($c=0.56$) $\mathrm{pPb}$ ($A=208$) collisions, so the
normalization is not the same as in figure 13, though the trend of the
small-$x$ evolution is clear. The most notable feature of these plots is the
numerical instability at high $k_{\perp}$ which results from the Fourier
transform.

Given the importance of small-$x$ evolution in calculating cross sections, it
should come as no surprise that there are many, many numerical studies of the
solutions to both the BFKL and BK equation (as well as the Balitsky-JIMWLK
hierarchy, but that goes beyond the scope of this dissertation). One can
either solve the equation in position space, using the form (65) directly, or
convert it to momentum space, in terms of $\phi$, and solve that. The latter
approach is taken in reference [79], where the momentum space BK equation is
written as

$\pd{\phi(k_{\perp},x_{g})}{\ln(1/x_{g})}=\frac{\alpha_{\text{s}}N_{c}}{\pi}\int_{0}^{\infty}\frac{\udc
k_{\perp}^{\prime 2}}{k_{\perp}^{\prime 2}}\biggl{[}\frac{k_{\perp}^{\prime
2}\phi(k_{\perp}^{\prime},x_{g})-k_{\perp}^{2}\phi(k_{\perp},x_{g})}{\abs{k_{\perp}^{2}-k_{\perp}^{\prime
2}}}+\frac{k_{\perp}^{2}\phi(k_{\perp},x_{g})}{\sqrt{4k_{\perp}^{\prime
4}+k_{\perp}^{4}}}\biggr{]}-\frac{\alpha_{\text{s}}N_{c}}{\pi}\phi(k_{\perp})^{2}$
(85)

Then a numerical Fourier transform allows us to convert either solution to the
other.

Figures 13 and 14 show two computations of $\mathcal{F}_{x_{g}}$: in figure
13, the solution is computed from the momentum space BK equation (85) and then
$\mathcal{F}_{x_{g}}$ is calculated using equation (31), whereas the plots in
figure 14 come from the solution to the position space BK equation (65) which
is then Fourier transformed. These two specific sets of plots have different
normalization, so they are not directly comparable, but looking at the
mid-$k_{\perp}$ region we can see some slight differences in how they depart
from the initial condition as $x_{g}$ gets smaller.

One important difference between the two solutions is the present of numerical
instabilities at large $k_{\perp}$ in the Fourier transform from the position
space solution (figure 14). The instability results from floating-point
roundoff error in the Fourier transform. This limits the validity of this
solution to moderate values of $k_{\perp}$, though this will not be a problem
for the calculations we will use this solution for in chapter NLO Inclusive
Pion Production. Conversely, in the direct momentum-space solution (figure
13), we see a spike at low $k_{\perp}$, which is a numerical artifact that
owes its existence to the lower boundary condition on $k_{\perp}$. In this
case, the solution can only be used above a certain momentum, though again,
this is not a problem for the calculations in chapter Drell-Yan Correlations
because it only samples the gluon distribution in the moderate $k_{\perp}$
range.

We can characterize the speed of the BK evolution by finding the saturation
scale $Q_{s}$ as a function of $\ln\frac{1}{x_{g}}$. In the GBW model, of
course, this is purely an exponential dependence,

$Q_{s}^{2}\biggl{(}\ln\frac{1}{x_{g}}\biggr{)}=Q_{0}^{2}x_{0}^{\lambda}\exp\biggl{(}\lambda\ln\frac{1}{x_{g}}\biggr{)}$
(86)

The analytic solution to the BK equation predicts [58, 59, 80] a more
complicated form,

$Q_{s}^{2}\biggl{(}\ln\frac{1}{x_{g}}\biggr{)}=Q_{0}^{2}\biggl{(}\ln\frac{1}{x_{g}}\biggr{)}^{\frac{3}{2(1-\gamma_{0})}}\exp\biggl{(}-\frac{\alpha_{\text{s}}N_{c}}{\pi}\frac{\chi(\gamma_{0})}{1-\gamma_{0}}\ln\frac{1}{x_{g}}\biggr{)}$
(87)

Figure 15: Saturation scales for the GBW model and the solution to the BK
equation with running coupling, solved in position space and Fourier
transformed into momentum space, extracted by finding $k_{\perp}$ at which
$\mathcal{F}_{x_{g}}(k_{\perp})=0.5\mathcal{F}_{x_{g}}(0)$. The flat portion
at the left end of the plot is where the initial condition of the BK evolution
is still strongly influencing the result. We can see that the evolution is
slower in the BK solution.

Figure 16: This plot, reproduced from reference [78], shows how adjustment of
the parameters of the BK evolution alters the $x_{g}$ dependence of the
saturation scale. Each curve shows the peak of the momentum distribution,
$k_{\text{max}}=\max k\phi(k^{2},Y)$, computed from the analytic formula for
the GBW model, from the fixed coupling BK evolution for selected values of
$\alpha_{s}$, and from the running coupling BK evolution with selected values
of $\Lambda_{\text{QCD}}^{2}$. For the BK evolution curves, the slope in the
upper range of rapidities decreases as $\alpha_{s}$ or
$\Lambda_{\text{QCD}}^{2}$ decreases, and the closest match to the slope of
the GBW model curve at $Y>15$ is achieved with $\alpha_{s}=0.062$ for fixed
coupling or $\Lambda_{\text{QCD}}^{2}=0.001$ for running coupling. The jagged
“steps” in the curves reflect the finite spacing of the momentum grid used in
the evolution.

To do this calculation we need a way to extract the saturation scale from the
solution to the BK equation. There are a couple of complementary prescriptions
we can use, depending on which form of the equation we’re working from. A
simple method is to define the saturation scale as the momentum at which the
solution crosses some particular threshold; specifically, in position space,

$S_{x_{g}}^{(2)}(1/Q_{s})=S_{0}$ (88)

or in momentum space,

$\mathcal{F}_{x_{g}}(Q_{s})=F_{0}\mathcal{F}_{x_{g}}(0)$ (89)

Alternatively, one can define the saturation scale as the momentum which
maximizes the product $k_{\perp}\phi(k_{\perp})$, as in reference [79]. The
results from the various prescriptions differ slightly in their overall
normalization, but they all predict essentially the same trend.

## Drell-Yan Correlations

In the previous chapters, we saw how the color glass condensate provides a
theoretical model of the effective structure of protons and nuclei. Of course,
the model is useless without the ability to connect it to measurable
quantities, namely cross sections. With that goal in mind, this chapter shows
the basic process of computing the cross section for the Drell-Yan process
AhX, from the unintegrated gluon distribution.

Drell-Yan processes, first predicted in 1970 [81], are those in which the
decay of a virtual photon or Z boson produces a lepton and antilepton. Lepton
pair production, as it’s called, is particularly appealing as a probe of
hadron structure for a few reasons:

  * •

Electrons or muons and their antiparticles make very distinctive signatures in
the detector, so they’re easy to identify

  * •

Not much else produces lepton-antilepton pairs, so background contamination is
low

  * •

Photons and Z bosons are uncharged, so they are not subject to strong
interactions on their way out of the collision and therefore final state
interactions with the hadron remnants are absent

  * •

The decay of the photon or Z into the lepton-antilepton pair is represented by
a factor which can be exactly calculated from electroweak theory

Accordingly, lepton pair production gives us a fairly direct window on the
unintegrated dipole gluon distribution. The correlation, as opposed to a
simple cross section, is particularly sensitive to the transverse momentum
dependence of the distribution, making it especially suitable as a probe of
saturation.

### Kinematics of $\mathrm{p}A$ collisions

Figure 17: Feynman diagram for AX at leading order

In a proton-nucleus collision, we assume the proton comes in from the $-z$
direction and the nucleon (which is part of a nucleus) comes in from the $+z$
direction, with purely longitudinal momenta

$\displaystyle p_{\textmd{\Pproton}}^{\mu}$
$\displaystyle=(E_{p},0,0,p_{\textmd{\Pproton}})$ $\displaystyle
p_{\textmd{\HepGenParticle{A}{}{}\xspace}}^{\mu}$
$\displaystyle=(E_{A},0,0,-p_{\textmd{\HepGenParticle{A}{}{}\xspace}})$ (90)

respectively. Here $p_{\textmd{\Pproton}}$ and
$p_{\textmd{\HepGenParticle{A}{}{}\xspace}}$ represent the magnitudes of the
proton’s 3-momentum $\vec{p}_{\textmd{\Pproton}}$ and the nucleon’s 3-momentum
$\vec{p}_{\textmd{\HepGenParticle{A}{}{}\xspace}}$, respectively — note that
$p_{\textmd{\HepGenParticle{A}{}{}\xspace}}$ refers to the momentum of _one_
nucleon, not the entire nucleus.

It’s easiest to analyze the collision in the longitudinal center of mass
frame, in which the center of mass of the proton and nucleon (not nucleus) is
at rest:
$p_{\textmd{\Pproton}z}+p_{\textmd{\HepGenParticle{A}{}{}\xspace}z}=0$. Since
the incoming proton and nucleon have zero transverse momentum and are
traveling at near light speed, we can neglect the masses of the proton and
nucleon and _approximately_ say

$\displaystyle E_{p}$
$\displaystyle=p_{\textmd{\Pproton}}=p_{\textmd{\Pproton}z}$ $\displaystyle
E_{A}$
$\displaystyle=p_{\textmd{\HepGenParticle{A}{}{}\xspace}}=-p_{\textmd{\HepGenParticle{A}{}{}\xspace}z}$
(91a) $\displaystyle p_{\textmd{\Pproton}}^{+}$
$\displaystyle=2E_{p}=2p_{\textmd{\Pproton}z}$ $\displaystyle
p_{\textmd{\HepGenParticle{A}{}{}\xspace}}^{+}$ $\displaystyle=0$ (91b)
$\displaystyle p_{\textmd{\Pproton}}^{-}$ $\displaystyle=0$ $\displaystyle
p_{\textmd{\HepGenParticle{A}{}{}\xspace}}^{-}$
$\displaystyle=2E_{A}=-2p_{\textmd{\HepGenParticle{A}{}{}\xspace}z}$ (91c)

according to the light-cone coordinate definitions in appendix 615. Given our
assumption that $\vec{p}_{p\perp}=\vec{p}_{A\perp}=0$, this implies that
$E_{p}=E_{A}=p_{\textmd{\Pproton}}=p_{\textmd{\HepGenParticle{A}{}{}\xspace}}$.
This enables us to write the Mandelstam variable $s$ for the proton and
nucleon (also denoted $s_{NN}$ in many sources) as

$s=(p_{\textmd{\Pproton}}^{\mu}+p_{\textmd{\HepGenParticle{A}{}{}\xspace}}^{\mu})^{2}=p_{\textmd{\Pproton}}^{+}p_{\textmd{\HepGenParticle{A}{}{}\xspace}}^{-}$
(92)

which in combination with equations (91) means that

$p_{\textmd{\Pproton}}^{+}=p_{\textmd{\HepGenParticle{A}{}{}\xspace}}^{-}=\sqrt{s}$
(93)

Our eventual goal is to determine the momentum fractions $x_{p}$, $x_{g}$, and
$z_{2}$, which we need for the parton distribution and gluon distribution, in
terms of observable quantities. For this reaction, those are $p_{\Ppi\perp}$,
$\sqrt{s}$, $\eta_{\Ppi}$, $\eta_{\Pphoton}$, and $M$, where $\eta$ is the
pseudorapidity, defined as

$\eta=\frac{1}{2}\ln\tan\frac{\theta}{2}$ (94)

Pseudorapidity and rapidity are related by the formula

$Y=\ln\frac{\sqrt{m^{2}+p_{\perp}^{2}\cosh^{2}\eta}+p_{\perp}\sinh\eta}{\sqrt{m^{2}+p_{\perp}^{2}}}$
(95)

In situations where $m\ll p_{\perp}\cosh\eta$, such as those considered in
this dissertation, this reduces to $Y=\eta$, so we can consider them
essentially interchangeable for the processes considered here.

Calculating the momentum fractions proceeds from the law of conservation of
momentum, with the assumption that the plus component of the total momentum is
contributed by the initial quark and the minus and transverse components of
the total momentum are contributed by the gluon:

$\displaystyle x_{p}p_{\textmd{\Pproton}}^{+}$
$\displaystyle=q^{+}=p_{\textmd{\Pphoton}}^{+}+k_{\textmd{\Pquark}}^{+}$ (96a)
$\displaystyle x_{g}p_{\textmd{\HepGenParticle{A}{}{}\xspace}}^{-}$
$\displaystyle=q^{-}=p_{\textmd{\Pphoton}}^{-}+k_{\textmd{\Pquark}}^{-}$ (96b)
$\displaystyle\vec{k}_{\textmd{\Pgluon}\perp}$
$\displaystyle=\vec{q}_{\perp}=\vec{p}_{\textmd{\Pphoton}\perp}+\vec{k}_{\textmd{\Pquark}\perp}$
(96c)

We’ll also need the final state momentum fractions, defined as
$z_{\textmd{\Pphoton}}\coloneqq
p_{\textmd{\Pphoton}}^{+}/p_{\textmd{\Pproton}}^{+}$ and
$z_{\textmd{\Ppi}}\coloneqq p_{\textmd{\Ppi}}^{+}/p_{\textmd{\Pproton}}^{+}$.
Using equations (628) and (93), these can be written as

$\displaystyle z_{\textmd{\Pphoton}}$
$\displaystyle=\sqrt{\frac{p_{\textmd{\Pphoton}\perp}^{2}+M^{2}}{s}}e^{Y_{\Pphoton}}$
$\displaystyle z_{\textmd{\Ppi}}$
$\displaystyle=\frac{p_{\textmd{\Ppi}\perp}}{\sqrt{s}}e^{Y_{\Ppi}}$ (97)

assuming that the pion mass is negligible. Now, we define
$z_{2}\coloneqq\frac{p_{\textmd{\Ppi}}^{+}}{k_{\textmd{\Pquark}}^{+}}$, the
momentum fraction of the pion relative to the quark it fragments from, and
take the pion to be collinear to the quark, so that all components of their
momenta obey this proportionality:
$p_{\textmd{\Ppi}}^{\mu}=z_{2}k_{\textmd{\Pquark}}^{\mu}$. Then in terms of
these variables, the plus component of the conservation equation (96a) leads
to

$x_{p}=z_{\textmd{\Pphoton}}+\frac{z_{\textmd{\Ppi}}}{z_{2}}=\sqrt{\frac{p_{\textmd{\Pphoton}\perp}^{2}+M^{2}}{s}}e^{Y_{\Pphoton}}+\frac{p_{\Ppi\perp}}{z_{2}\sqrt{s}}e^{Y_{\Ppi}}$
(98)

Similarly, we know that $p^{-}=p^{+}e^{-2Y}$ (625), and using that on the
right side of equation (96b), we get

$x_{g}=z_{\textmd{\Pphoton}}e^{-2Y_{\Pphoton}}+\frac{z_{\textmd{\Ppi}}}{z_{2}}e^{-2Y_{\Ppi}}=\sqrt{\frac{p_{\textmd{\Pphoton}\perp}^{2}+M^{2}}{s}}e^{-Y_{\Pphoton}}+\frac{p_{\Ppi\perp}}{z_{2}\sqrt{s}}e^{-Y_{\Ppi}}$
(99)

For the inclusive cross section, we’ll need to integrate over the entire
allowed phase space of the final state particles emerging from fragmentation.
Our code parametrizes that space with the “reduced” photon momentum fraction
$z\coloneqq\frac{p_{\textmd{\Pphoton}}^{+}}{q^{+}}=\frac{z_{\textmd{\Pphoton}}}{x_{p}}$,
which represents the fraction of the total plus momentum carried by the
photon. and the total transverse momentum $\vec{q}_{\perp}$, instead of the
transverse pion momentum $p_{\textmd{\Ppi}\perp}$ and rapidity $Y_{\Ppi}$ and
the momentum fraction $z_{2}$, so we’ll need to remove the latter variables
from the formulas.

To achieve this parametrization, we can use the fact that

$\displaystyle zq^{+}$
$\displaystyle=p_{\textmd{\Pphoton}}^{+}=z_{\textmd{\Pphoton}}p_{\textmd{\Pproton}}^{+}$
$\displaystyle(1-z)q^{+}$
$\displaystyle=k_{\textmd{\Pquark}}^{+}=\frac{p_{\textmd{\Ppi}}^{+}}{z_{2}}=\frac{z_{\textmd{\Ppi}}p_{\textmd{\Pproton}}^{+}}{z_{2}}$
(100)

Combining the equation on the left with (96a) and the definition of
$z_{\textmd{\Pphoton}}$ (97), we find that

$x_{p}=\frac{q^{+}}{p_{\textmd{\Pproton}}^{+}}=\frac{z_{\textmd{\Pphoton}}}{z}=\frac{1}{z}\sqrt{\frac{p_{\textmd{\Pphoton}\perp}^{2}+M^{2}}{s}}e^{Y_{\Pphoton}}$
(101)

For $x_{g}$, we can use the relation
$e^{-Y_{\Ppi}}=\frac{p_{\textmd{\Ppi}\perp}}{p_{\textmd{\Ppi}}^{+}}=\frac{k_{\textmd{\Pquark}\perp}}{k_{\textmd{\Pquark}}^{+}}$,
which follows from equation (628) and the proportionality of
$p_{\textmd{\Ppi}}^{\mu}$ and $k_{\textmd{\Pquark}}^{\mu}$, to express
equation (99) as

$x_{g}=\sqrt{\frac{p_{\textmd{\Pphoton}\perp}^{2}+M^{2}}{s}}e^{-Y_{\Pphoton}}+\frac{k_{\textmd{\Pquark}\perp}^{2}}{k_{\textmd{\Pquark}}^{+}\sqrt{s}}$
(102)

We know from equation (100) on the right that
$k_{\textmd{\Pquark}}^{+}=(1-z)q^{+}$, and from (96a) that
$q^{+}=x_{p}p_{\textmd{\Pproton}}^{+}$. Putting these together with the
assumption that $p_{\textmd{\Pproton}}^{+}=\sqrt{s}$, we find that

$x_{g}=\sqrt{\frac{p_{\textmd{\Pphoton}\perp}^{2}+M^{2}}{s}}e^{-Y_{\Pphoton}}+\frac{k_{\textmd{\Pquark}\perp}^{2}}{(1-z)x_{p}s}$
(103)

Adding in momentum conservation (96c) gets us to the final form of the
expression,

$x_{g}=\sqrt{\frac{p_{\textmd{\Pphoton}\perp}^{2}+M^{2}}{s}}e^{-Y_{\Pphoton}}+\frac{(\vec{q}_{\perp}-\vec{p}_{\textmd{\Pphoton}\perp})^{2}}{(1-z)x_{p}s}$
(104)

### Hybrid Factorization

The final ingredient we will need for the calculation of the Drell-Yan process
is the hybrid factorization formalism. Earlier, in section The Parton Model
and Integrated Parton Distributions, I briefly mentioned that our ability to
split a calculation into the hard scattering process (the low-level
interaction), which is calculable from perturbative QCD, and the parton
distribution, which is not, relies on the existence of factorization theorems.
In fact, factorization can be implemented in several different ways. The
traditional approach [23, 24, 25], the collinear factorization framework,
treats both projectile and target on an equal footing, representing them both
with integrated parton distributions. However, the collinear factorization
approach with integrated PDFs is not well suited to describe dense systems
which are within the saturation regime.

An alternative approach, called hybrid factorization [82], was developed
specifically for asymmetric proton-nucleus collisions. In hybrid
factorization, the target is treated as a color glass condensate and
represented using an unintegrated parton distribution. This framework is
designed for the kinematic region where $Y\gtrsim 2.5$, the larger the better
— this is called the “forward rapidity” region because large values of $Y$,
defined as

$Y=\frac{1}{2}\ln\frac{E+p_{z}}{E-p_{z}}$ (105)

correspond to emission of particles near the direction of the projectile. At
these forward rapidities, the interaction is heavily weighted toward small-$x$
gluons from the target, $x_{g}\sim\num{e-4}$ or even lower, and relatively
large-$x$ partons from the projectile, with $x_{p}\sim\num{e-1}$ or
$\num{e-2}$.

In either approach, factorization necessitates ignoring the dynamics of the
interacting fields above a certain momentum scale, in the same sense as
renormalization group calculations. Physical quantities are, of course,
independent of this factorization scale, designated $\mu$, but when we
calculate a perturbative approximation to fixed order (leading order or next-
to-leading order etc.), we will find that the resulting formula does depend on
the factorization scale. At leading order, the dependence arises entirely
through nonperturbative quantities, the parton distributions and fragmentation
functions, but beyond leading order, the perturbative part of the calculation
can have its own dependence on $\mu$. We’ll see this borne out in the
calculations of chapter NLO Inclusive Pion Production.

There is some loose physical motivation for choosing a particular value for
the factorization scale in any given calculation [83, 84], but in general, the
dependence on the arbitrary scale $\mu$ should be viewed as a theoretical
source of uncertainty in the final result. Incorporating more terms in the
perturbative series tends to reduce the dependence on $\mu$, as we’ll see in
chapter Computational Results and Analysis. This provides a powerful
motivation to extend perturbative quantities beyond leading order.

### Calculation of Drell-Yan Correlation

Let’s now turn to the procedure for CGC calculations described in chapter The
Color Glass Condensate, which we can use to calculate the cross section for
the Drell-Yan process. The calculation begins with the process AA, whose cross
section is calculated in momentum space in reference [85] and in transverse
coordinate space in [38]. However, these calculations only treat the case of a
real photon, $M=0$, whereas we need the result for photon virtuality
$M^{2}>0$. The modifications to the calculation of [38] required to account
for nonzero virtuality are described in reference [86].

Figure 18: Compound diagrams for a quark emitting a photon and interacting
with a nucleus

Following the general procedure described in chapter Hadron Structure From
QCD, we start with the four diagrams contributing to this process at leading
order, as shown in figure 18. They differ only in whether the interaction with
the nucleus occurs before or after the emission of the photon. For the first
diagram, putting together the splitting wavefunction and the Wilson line on
each side of the cut, we get

$\frac{\udc^{6}\sigma}{\udddc\vec{p}_{\Pphoton}\udddc\vec{k}_{\Pquark}}=\iint\frac{\uddc\vec{}x_{\perp}}{(2\pi)^{2}}\iint\frac{\uddc\vec{}x_{\perp}^{\prime}}{(2\pi)^{2}}\iint\frac{\uddc\vec{}b_{\perp}}{(2\pi)^{2}}\iint\frac{\uddc\vec{}b_{\perp}^{\prime}}{(2\pi)^{2}}e^{-i\vec{p}_{\textmd{\Pphoton}\perp}\cdot(\vec{}x_{\perp}-\vec{}x_{\perp}^{\prime})}e^{-i\vec{k}_{\textmd{\Pquark}\perp}\cdot(\vec{}b_{\perp}-\vec{}b_{\perp}^{\prime})}\\\
\times\sum_{\lambda\alpha\beta}\conj{(\psi_{\alpha\beta}^{\lambda})}(\vec{}x_{\perp}^{\prime}-\vec{}b_{\perp}^{\prime})\psi_{\alpha\beta}^{\lambda}(\vec{}x_{\perp}-\vec{}b_{\perp})\braket{\trace
U\bigl{(}z\vec{}x_{\perp}+(1-z)\vec{}b_{\perp}\bigr{)}\herm{U}\bigl{(}z\vec{}x_{\perp}^{\prime}+(1-z)\vec{}b_{\perp}^{\prime}\bigr{)}}_{x_{g}}$
(106)

where $\lambda$ represents the polarization of the photon, and $\alpha$ and
$\beta$ are the helicities of the initial and final quarks, respectively. The
averaging operation $\braket{\cdots}_{x_{g}}$ averages over the unknown color
charge distribution of the nucleus. The exponential factors have their origin
in the LSZ reduction formula.888More detail on this aspect of the derivation
is in references [85, 87], and see also [88]. At the end of this expression,
we see the dipole gluon distribution in the form

$\braket{\trace
U(\vec{}x_{\perp})\herm{U}(\vec{}y_{\perp}))}_{x_{g}}=S_{x_{g}}^{(2)}(\vec{}x_{\perp},\vec{}y_{\perp})$
(107)

Adding in the other diagrams from figure 18 gives the expression

$\frac{\udc^{6}\sigma}{\udddc\vec{p}_{\Pphoton}\udddc\vec{k}_{\Pquark}}=\iint\frac{\uddc\vec{}x_{\perp}}{(2\pi)^{2}}\iint\frac{\uddc\vec{}x_{\perp}^{\prime}}{(2\pi)^{2}}\iint\frac{\uddc\vec{}b_{\perp}}{(2\pi)^{2}}\iint\frac{\uddc\vec{}b_{\perp}^{\prime}}{(2\pi)^{2}}e^{-i\vec{p}_{\textmd{\Pphoton}\perp}\cdot(\vec{}x_{\perp}-\vec{}x_{\perp}^{\prime})}e^{-ik_{\textmd{\Pquark}\perp}\cdot(\vec{}b_{\perp}-\vec{}b_{\perp}^{\prime})}\\\
\times\sum_{\lambda\alpha\beta}\conj{(\psi_{\alpha\beta}^{\lambda})}(\vec{}x_{\perp}^{\prime}-\vec{}b_{\perp}^{\prime})\psi_{\alpha\beta}^{\lambda}(\vec{}x_{\perp}-\vec{}b_{\perp})\bigl{[}S_{x_{g}}^{(2)}\bigl{(}z\vec{}x_{\perp}+(1-z)\vec{}b_{\perp},z\vec{}x_{\perp}^{\prime}+(1-z)\vec{}b_{\perp}^{\prime}\bigr{)}+S_{x_{g}}^{(2)}(\vec{}b_{\perp},\vec{}b_{\perp}^{\prime})\\\
-S_{x_{g}}^{(2)}\bigl{(}\vec{}b_{\perp},z\vec{}x_{\perp}^{\prime}+(1-z)\vec{}b_{\perp}^{\prime}\bigr{)}-S_{x_{g}}^{(2)}\bigl{(}z\vec{}x_{\perp}+(1-z)\vec{}b_{\perp},\vec{}b_{\perp}^{\prime}\bigr{)}\bigr{]}$
(108)

Following reference [38], we make the substitutions
$\vec{u}=\vec{}x_{\perp}-\vec{}b_{\perp}$,
$\vec{u}^{\prime}=\vec{}x_{\perp}^{\prime}-\vec{}b_{\perp}^{\prime}$, and

first term: $\displaystyle\vec{v}$
$\displaystyle=z\vec{}x_{\perp}+(1-z)\vec{}b_{\perp}$
$\displaystyle\vec{v}^{\prime}$
$\displaystyle=z\vec{}x_{\perp}^{\prime}+(1-z)\vec{}b_{\perp}^{\prime}$ second
term: $\displaystyle\vec{v}$ $\displaystyle=\vec{}b_{\perp}$
$\displaystyle\vec{v}^{\prime}$ $\displaystyle=\vec{}b_{\perp}^{\prime}$ third
term: $\displaystyle\vec{v}$ $\displaystyle=\vec{}b_{\perp}$
$\displaystyle\vec{v}^{\prime}$
$\displaystyle=z\vec{}x_{\perp}^{\prime}+(1-z)\vec{}b_{\perp}^{\prime}$ fourth
term: $\displaystyle\vec{v}$
$\displaystyle=z\vec{}x_{\perp}+(1-z)\vec{}b_{\perp}$
$\displaystyle\vec{v}^{\prime}$ $\displaystyle=\vec{}b_{\perp}^{\prime}$

which results in

$\frac{\udc^{6}\sigma}{\udddc\vec{p}_{\Pphoton}\udddc\vec{k}_{\Pquark}}=\iint\frac{\uddc\vec{u}}{(2\pi)^{2}}\iint\frac{\uddc\vec{u}^{\prime}}{(2\pi)^{2}}\iint\frac{\uddc\vec{v}}{(2\pi)^{2}}\iint\frac{\uddc\vec{v}^{\prime}}{(2\pi)^{2}}e^{-i\vec{q}_{\perp}\cdot(\vec{v}-\vec{v}^{\prime})}S_{x_{g}}^{(2)}(\vec{v},\vec{v}^{\prime})\\\
\times\sum_{\lambda\alpha\beta}\conj{(\psi_{\alpha\beta}^{\lambda})}(\vec{u}^{\prime})\psi_{\alpha\beta}^{\lambda}(\vec{u})\bigl{[}e^{-i(\vec{p}_{\textmd{\Pphoton}\perp}-z\vec{q}_{\perp})\cdot\vec{u}}e^{i(\vec{p}_{\textmd{\Pphoton}\perp}-z\vec{q}_{\perp})\cdot\vec{u}^{\prime}}+e^{-i\vec{p}_{\textmd{\Pphoton}\perp}\cdot\vec{u}}e^{i\vec{p}_{\textmd{\Pphoton}\perp}\cdot\vec{u}^{\prime}}\\\
-e^{-i\vec{p}_{\textmd{\Pphoton}\perp}\cdot\vec{u}}e^{i(\vec{p}_{\textmd{\Pphoton}\perp}-z\vec{q}_{\perp})\cdot\vec{u}^{\prime}}-e^{-i(\vec{p}_{\textmd{\Pphoton}\perp}-z\vec{q}_{\perp})\cdot\vec{u}}e^{i\vec{p}_{\textmd{\Pphoton}\perp}\cdot\vec{u}^{\prime}}\bigr{]}$
(109)

In this expression, the integrals over $\vec{v}$ and $\vec{v}^{\prime}$ give
us the momentum space dipole gluon distribution, according to equation (32),
and the integrals over $\vec{u}$ and $\vec{u}^{\prime}$ implement the Fourier
transform of the splitting wavefunction,

$\iint\frac{\uddc\vec{u}}{(2\pi)^{2}}\iint\frac{\uddc\vec{u}^{\prime}}{(2\pi)^{2}}\sum_{\lambda\alpha\beta}\conj{(\psi_{\alpha\beta}^{\lambda})}(\vec{u}^{\prime})\psi_{\alpha\beta}^{\lambda}(\vec{u})e^{-i\vec{}q_{1\perp}\cdot\vec{u}}e^{-i\vec{}q_{2\perp}\cdot\vec{u}^{\prime}}$
(110)

for different values of $\vec{}q_{1\perp}$ and $\vec{}q_{2\perp}$ in each
term.

This is the first point at which the nonzero virtuality becomes relevant. The
explicit form of the splitting wavefunction for a virtual photon is given in
[86] as

$\psi_{\alpha\beta}^{\lambda}(q^{+},p_{\textmd{\Pphoton}}^{+},\vec{u})=2\pi\sqrt{\frac{2}{p_{\textmd{\Pphoton}}^{+}}}\begin{cases}(1-z)MK_{0}(\epsilon_{M}u),&\lambda=0\quad\text{(longitudinal)}\\\
i\epsilon_{M}K_{1}(\epsilon_{M}u)\frac{\vec{u}\cdot\vec{\epsilon}_{\perp}^{(1)}}{u^{2}}\bigl{(}\delta_{\alpha-}\delta_{\beta-}+(1-z)\delta_{\alpha+}\delta_{\beta+}\bigr{)},&\lambda=1\quad\text{(transverse)}\\\
i\epsilon_{M}K_{1}(\epsilon_{M}u)\frac{\vec{u}\cdot\vec{\epsilon}_{\perp}^{(2)}}{u^{2}}\bigl{(}\delta_{\alpha+}\delta_{\beta+}+(1-z)\delta_{\alpha-}\delta_{\beta-}\bigr{)},&\lambda=2\quad\text{(transverse)}\end{cases}$
(111)

where $\epsilon_{M}^{2}=(1-z)M^{2}$. We can plug that into (110) and perform
the sums over $\lambda$, $\alpha$, and $\beta$. Using the identity (658),
equation (110) reduces to

$\iint\frac{\uddc\vec{u}}{(2\pi)^{2}}\iint\frac{\uddc\vec{u}^{\prime}}{(2\pi)^{2}}4\pi^{2}\frac{2}{p_{\textmd{\Pphoton}}^{+}}K_{0}(\epsilon_{M}u)K_{0}(\epsilon_{M}u^{\prime})M^{2}(1-z)^{2}e^{-i\vec{}q_{1\perp}\cdot\vec{u}}e^{i\vec{}q_{2\perp}\cdot\vec{u}^{\prime}}\\\
+\iint\frac{\uddc\vec{u}}{(2\pi)^{2}}\iint\frac{\uddc\vec{u}^{\prime}}{(2\pi)^{2}}4\pi^{2}\frac{2}{p_{\textmd{\Pphoton}}^{+}}K_{1}(\epsilon_{M}u)K_{1}(\epsilon_{M}u^{\prime})\epsilon_{M}^{2}\\\
\times\frac{\vec{u}^{\prime}\cdot\vec{u}-(\vec{u}^{\prime}\cdot\vec{p}_{\textmd{\Pphoton}})(\vec{u}\cdot\vec{p}_{\textmd{\Pphoton}})/p_{\textmd{\Pphoton}}^{2}}{u^{\prime
2}u^{2}}\bigl{[}(1-z)^{2}+1\bigr{]}e^{-i\vec{}q_{1\perp}\cdot\vec{u}}e^{i\vec{}q_{2\perp}\cdot\vec{u}^{\prime}}$
(112)

where the first line comes from the longitudinal component and the second and
third lines from the transverse components.

Performing the Fourier integrals, the cross section for AX becomes [86]

$\frac{\udc^{8}\sigma}{\udc Y_{\Pphoton}\udc
Y_{\Pquark}\uddc\vec{p}_{\textmd{\Pphoton}\perp}\uddc\vec{k}_{\textmd{\Pquark}\perp}\uddc\vec{b}}=\frac{\alpha_{\text{em}}e_{q}^{2}}{2\pi^{2}}(1-z)\mathcal{F}_{x_{g}}(\vec{k}_{\textmd{\Pgluon}\perp})\\\
\times\Biggl{\\{}\bigl{[}1+(1-z)^{2}\bigr{]}\frac{z^{2}k_{\textmd{\Pgluon}\perp}^{2}}{\bigl{[}(\vec{p}_{\textmd{\Pphoton}\perp}-z\vec{k}_{\textmd{\Pgluon}\perp})^{2}+\epsilon_{M}^{2}\bigr{]}\bigl{[}p_{\textmd{\Pphoton}\perp}^{2}+\epsilon_{M}^{2}\bigr{]}}\\\
-z^{2}\epsilon_{M}^{2}\biggl{[}\frac{1}{p_{\textmd{\Pphoton}\perp}^{2}+\epsilon_{M}^{2}}-\frac{1}{(\vec{p}_{\textmd{\Pphoton}\perp}-z\vec{k}_{\textmd{\Pgluon}\perp})^{2}+\epsilon_{M}^{2}}\biggr{]}^{2}\Biggr{\\}}$
(113)

where $e_{q}$ is the charge of the initial-state quark.

To obtain the cross section for
$\HepProcess{\Pproton\HepGenParticle{A}{}{}\xspace\HepTo\Plepton\APlepton\Ppizero\HepGenParticle{X}{}{}\xspace}$,
we need to multiply (113) by the integrated parton distribution, the
fragmentation function, and the factor representing the splitting of the
virtual photon into a lepton pair, which is
$\frac{\alpha_{\text{em}}}{3\pi}\frac{\udc M^{2}}{M^{2}}$. Then we integrate
over all possible values of the unknown momentum fraction $z_{2}$.

$\frac{\udc^{9}\sigma}{\udc Y_{\Pphoton}\udc
Y_{\Ppi}\uddc\vec{p}_{\textmd{\Pphoton}\perp}\uddc\vec{}p_{\Ppi\perp}\uddc\vec{b}\,\udc
M^{2}}=\frac{\alpha_{\text{em}}^{2}e_{q}^{2}}{6\pi^{3}M^{2}}\int_{\frac{z_{\textmd{\Ppi}}}{1-z_{\textmd{\Pphoton}}}}^{1}\frac{\udc
z_{2}}{z_{2}^{2}}\sum_{q}x_{p}q(x_{p},\mu)D_{\Ppi/q}(z_{2},\mu)(1-z)\mathcal{F}_{x_{g}}(\vec{}k_{\Pgluon\perp})\\\
\times\Biggl{\\{}\bigl{[}1+(1-z)^{2}\bigr{]}\frac{z^{2}k_{\Pgluon\perp}^{2}}{\bigl{[}(\vec{p}_{\textmd{\Pphoton}\perp}-z\vec{}k_{\Pgluon\perp})^{2}+\epsilon_{M}^{2}\bigr{]}\bigl{[}p_{\textmd{\Pphoton}\perp}^{2}+\epsilon_{M}^{2}\bigr{]}}\\\
-z^{2}\epsilon_{M}^{2}\biggl{[}\frac{1}{p_{\textmd{\Pphoton}\perp}^{2}+\epsilon_{M}^{2}}-\frac{1}{(\vec{p}_{\textmd{\Pphoton}\perp}-z\vec{}k_{\Pgluon\perp})^{2}+\epsilon_{M}^{2}}\biggr{]}^{2}\Biggr{\\}}$
(114)

The other cross section we need, for the denominator of the correlation, is
the inclusive cross section for AX, which we can get by integrating over all
possible combinations of the pion momentum and everything else (X) that could
come out of the reaction — or, equivalently, we can just ignore the fact that
the final-state quark fragments at all, and just integrate over all possible
momenta of that quark. It’s most straightforward to get this by starting again
from equation (113) and inserting the integrated parton distribution (but not
the fragmentation function) and the lepton splitting factor
$\frac{\alpha_{\text{em}}}{3\pi}\frac{\udc M^{2}}{M^{2}}$.

Then we integrate over the unknown momentum of the quark, represented by
$Y_{\Pquark}$ and $\vec{k}_{\textmd{\Pquark}}$. We’re assuming the quark to be
on-shell, $k_{\textmd{\Pquark}}^{\mu}k_{\textmd{\Pquark}\mu}=m_{\Pquark}^{2}$,
which means there are only three degrees of freedom in this momentum. Now,
since the photon momentum $\vec{p}_{\textmd{\Pphoton}}$ is a constant, we can
change variables from $k_{\textmd{\Pquark}}$ to
$\vec{q}=\vec{k}_{\textmd{\Pquark}}+\vec{p}_{\textmd{\Pphoton}}$, which is a
simple shift. We can then take advantage of the transformations999To derive
the Jacobian determinant for $Y_{\Pquark}$, substitute
$p^{+}=\sqrt{p_{z}^{2}+p_{\perp}^{2}+m^{2}}+p_{z}$ into equation (629) and
differentiate. Only $\pd{Y}{p_{z}}=\frac{1}{E}$ survives into the final
determinant; the other components wind up being multiplied by zero.

$\displaystyle\udc z\uddc\vec{}k_{\Pgluon\perp}$
$\displaystyle=\frac{z}{E_{q}}\udddc\vec{q}$ $\displaystyle\udc
Y_{\Pquark}\uddc\vec{}k_{\Pgluon\perp}$
$\displaystyle=\frac{1}{E_{\Pquark}}\udddc\vec{q}$ (115)

to work out the inclusive cross section to

$\frac{\udc^{6}\sigma}{\udc
Y_{\Pphoton}\uddc\vec{p}_{\textmd{\Pphoton}\perp}\uddc\vec{b}\,\udc
M^{2}}=\frac{\alpha_{\text{em}}^{2}e_{q}^{2}}{6\pi^{3}M^{2}}\sum_{q}\int_{z_{\textmd{\Pphoton}}}^{1}\frac{\udc
z}{z}x_{p}q(x_{p},\mu)\iint\uddc\vec{}k_{\Pgluon\perp}\mathcal{F}_{x_{g}}(\vec{}k_{\Pgluon\perp})\\\
\times\Biggl{\\{}\bigl{[}1+(1-z)^{2}\bigr{]}\frac{z^{2}k_{\Pgluon\perp}^{2}}{\bigl{[}(\vec{p}_{\textmd{\Pphoton}\perp}-z\vec{}k_{\Pgluon\perp})^{2}+\epsilon_{M}^{2}\bigr{]}\bigl{[}p_{\textmd{\Pphoton}\perp}^{2}+\epsilon_{M}^{2}\bigr{]}}\\\
-z^{2}\epsilon_{M}^{2}\biggl{[}\frac{1}{p_{\textmd{\Pphoton}\perp}^{2}+\epsilon_{M}^{2}}-\frac{1}{(\vec{p}_{\textmd{\Pphoton}\perp}-z\vec{}k_{\Pgluon\perp})^{2}+\epsilon_{M}^{2}}\biggr{]}^{2}\Biggr{\\}}$
(116)

which is consistent with other inclusive results derived in references [89,
90, 91, 87].

The correlation is defined as the ratio of the double inclusive to single
inclusive cross sections. For this work, we integrate over transverse momentum
and assume there is no impact parameter dependence, thus considering the
correlation only as a function of rapidity, virtual photon mass, and angle.

$C^{\text{DY}}(Y_{\Pphoton},Y_{\Ppizero},M,\Delta\phi)=\frac{\displaystyle\int_{p_{\textmd{\Pphoton}\perp}>p_{\perp\text{cut}}}\uddc\vec{p}_{\textmd{\Pphoton}\perp}\int_{p_{\Ppi\perp}>p_{\perp\text{cut}}}\uddc\vec{}p_{\Ppi\perp}\frac{\udc^{9}\sigma}{\udc
Y_{\Pphoton}\udc
Y_{\Ppizero}\uddc\vec{p}_{\textmd{\Pphoton}\perp}\uddc\vec{}p_{\Ppi\perp}\uddc\vec{b}\,\udc
M^{2}}}{\displaystyle\int_{p_{\textmd{\Pphoton}\perp}>p_{\perp\text{cut}}}\uddc\vec{p}_{\textmd{\Pphoton}\perp}\frac{\udc^{6}\sigma}{\udc
Y_{\Pphoton}\uddc\vec{p}_{\textmd{\Pphoton}\perp}\uddc\vec{b}\,\udc M^{2}}}$
(117)

### Results

Figures 19 and 20 show the predictions of this model for kinematic conditions
characteristic of RHIC and the LHC, respectively. The RHIC results in figure
19 use $\sqrt{s}=\SI{200}{GeV}$, $A=197$, and
$p_{\perp\text{cut}}=\SI{1.5}{GeV}$, and the LHC results use
$\sqrt{s}=\SI{8800}{GeV}$, $A=208$, and $p_{\perp\text{cut}}=\SI{1.5}{GeV}$.

There are two characteristics that are immediately evident: the double-peak
structure surrounding back-to-back emission ($\Delta\phi\approx\pi$), and the
variation in the parallel-emission behavior around $\Delta\phi\approx 0,2\pi$.
While the double peak is roughly similar for all three models of the gluon
distribution, for most kinematic conditions the BK equation with either fixed
or running coupling gives a large excess for parallel emission relative to the
GBW model. Since the BK solution differs from the GBW model primarily in the
asymptotic behavior of $\mathcal{F}_{x_{g}}$ at large $k_{\Pgluon\perp}$,
power versus exponential respectively, we can see that the parallel emission
comes primarily from large-transverse momentum gluons, whereas the back-to-
back emission comes mostly from low-transverse momentum gluons. This makes
sense intuitively, since if the lepton pair and the hadron are to emerge in
nearly the same direction, they’ll have to have at least
$2p_{\perp\text{cut}}$ of total transverse momentum to be detected, but back-
to-back emission can come from a system with little or no transverse momentum.

Figure 21 shows an explicit breakdown of the momenta
$k_{\textmd{\Pgluon}\perp}$ that contribute to the correlation at each angle,
which confirms this interpretation: the correlation for parallel emission
comes from higher momentum than the central peak. Overall, the correlation at
the LHC receives contributions from higher momenta than at RHIC, due to there
being more total energy $\sqrt{s}$ available in the collision.

Figure 19: The angular correlations between the virtual photons and pions at
RHIC, at medium rapidity, $Y_{\Pphoton}=Y_{\Ppizero}=2.5$. The upper graph
shows the correlation for a virtual photon mass of $M=\SI{0.5}{GeV}$, and the
lower one, for $M=\SI{4}{GeV}$. In each case, the three curves for GBW, fixed
coupling BK, and running coupling BK, exhibit basically the same double-peak
structure around $\Delta\phi=\pi$, but they show differing behavior near
$\Delta\phi=0,2\pi$, the near side correlation. This relates to the
large-$k^{2}$ behavior of the corresponding gluon distributions. Reproduced
from figure 5 of [78].

Figure 20: The virtual photon-pion angular correlations at LHC at rapidity
$Y_{\gamma}=Y_{\pi}=4$. The upper and lower graphs show $M=\SI{4}{GeV}$ and
$M=\SI{8}{GeV}$, respectively. As in figure 19, the key difference is in the
near side correlation, arising due to differences in the high-momentum region
of the gluon distribution. Reproduced from figure 6 of [78].

Figure 21: Breakdown of the contributions to the correlation, split by the
transverse momentum $k_{\textmd{\Pgluon}\perp}$ acquired by the quark as it
interacts with the gluon field of the target. Momentum ranges are in
$\si{GeV}$. The plot on the left corresponds to the top plot of figure 19,
with $M=\SI{0.5}{GeV}$ at RHIC, and the one on the right corresponds to the
top plot of figure 20 with $M=\SI{4}{GeV}$ using LHC parameters. In both
cases, there is a fairly sharp transition between the momenta that contribute
to the central back-to-back emission peak and those that contribute to the
parallel emission peak. For the RHIC plot on the left, that transition is
around $\SI{1}{GeV}$, and on the right, between $\SI{10}{GeV}$ and
$\SI{100}{GeV}$.

#### Peak Characteristics

The double peak of figures 19 and 20 is particularly interesting because it is
unique to the Drell-Yan process; it doesn’t appear in similar calculations for
the correlation between produced hadrons or jets [92]. We can trace the origin
of the central “dip” to the kinematic factor which appears in equations (114)
and (116):

$\bigl{[}1+(1-z)^{2}\bigr{]}\frac{z^{2}k_{\Pgluon\perp}^{2}}{\bigl{[}(\vec{p}_{\textmd{\Pphoton}\perp}-z\vec{}k_{\Pgluon\perp})^{2}+\epsilon_{M}^{2}\bigr{]}\bigl{[}p_{\textmd{\Pphoton}\perp}^{2}+\epsilon_{M}^{2}\bigr{]}}\\\
-z^{2}\epsilon_{M}^{2}\biggl{[}\frac{1}{p_{\textmd{\Pphoton}\perp}^{2}+\epsilon_{M}^{2}}-\frac{1}{(\vec{p}_{\textmd{\Pphoton}\perp}-z\vec{}k_{\Pgluon\perp})^{2}+\epsilon_{M}^{2}}\biggr{]}^{2}$
(118)

This approaches zero as $k_{\Pgluon\perp}\to 0$. If we were analyzing the
parton-level cross section (113) as a function of $p_{\textmd{\Pphoton}\perp}$
and $k_{\Pquark\perp}$, we would expect to see the correlation drop to zero
when $\vec{p}_{\textmd{\Pphoton}\perp}+\vec{}k_{\Pquark\perp}=\vec{0}$. But
the double inclusive result (114) involves an integration over $z_{2}$, so
even when evaluated at specific values of $\vec{p}_{\textmd{\Pphoton}\perp}$
and $\vec{}p_{\Ppi\perp}$, it incorporates a range of values of
$\vec{}k_{\Pgluon\perp}=\frac{\vec{}p_{\Ppi\perp}}{z_{2}}+\vec{p}_{\textmd{\Pphoton}\perp}$.
Thus there won’t be any point
$(\vec{p}_{\textmd{\Pphoton}\perp},\vec{}p_{\Ppi\perp})$ where the correlation
is exactly zero; effectively, the fragmentation blurs the peak structure. We
also integrate over transverse momentum when calculating the correlation
(117), which smooths out the peak even more.

Depending on the kinematic conditions, the central suppression may be small
enough that a true double peak doesn’t appear. This appears in reference [93],
for example, which shows the same analysis for real photons, which is the
$M\to 0$ limit of this calculation, and in figure 22 we can see that the peak
depth decreases as we move toward lower photon masses.

Figure 22: Plots of the half-width of the double peak, on the left, and the
depth of the dip, on the right, as a function of virtual photon mass.

Because the correlation at the back-to-back emission peak comes from low-
momentum regions of the gluon distribution, we can reasonably expect the
characteristics of the peak to reflect the presence and nature of saturation.
In particular, the results in figures 19 and 20 show that the results for the
running coupling BK equation, for which the evolution of the saturation scale
with $\ln\frac{1}{x_{g}}$ is slower, leads to a suppression of the peak height
relative to the fixed coupling. A larger saturation scale $Q_{s}$ seems to
produce a taller peak.

Figure 22 supports the same conclusion by showing that the fixed coupling BK
equation produces a peak which is both shallower and thinner than that of the
GBW model. Looking back at figure 16, we see that the fixed coupling BK
solution has a smaller saturation scale than the GBW model, so again, we see
that a large saturation scale makes the double peak larger.

### Conclusions

There are two main ideas to take away from this chapter. First, the
calculation of the lepton pair-hadron angular correlation serves as an example
of how the CGC and hybrid formalism are used in practice. This lays the
groundwork for the next-to-leading order calculation to be shown in the next
chapter.

Secondly, the correlation is an interesting numerical result in its own right.
It’s clear that it is a sensitive probe of the gluon distribution in both
high- and low-momentum regions. We see an enhancement of the correlation
around parallel emission, which gives us information about the high-momentum
tail of the gluon distribution, and also the unique double peak structure for
near-parallel emission, whose properties can reveal the characteristics of the
gluon distribution in the low-momentum saturation regime.

## NLO Inclusive Pion Production

Leading order calculations, like that in the last chapter, typically give a
decent approximation to measured cross sections. But there are several reasons
to look to the next order of perturbation theory in hopes of improving the
results. Of course, one generally expects that adding more terms to a
perturbation series will give a more accurate approximation of the exact (all-
order) result. Beyond that, though, the cross section depends on the
renormalization and factorization scales $\mu$ through terms which are
subleading in $\alpha_{\text{s}}$. If we are able to calculate some of those
terms explicitly, the remaining dependence on $\mu$, which is unknown and thus
must be considered a theoretical uncertainty, will be reduced.

Calculating the subleading corrections to a cross section is not a trivial
task. Beyond leading order, the number of diagrams that contribute grows
drastically, often roughly exponentially with the order. There are also
additional kinematic degrees of freedom, corresponding to the momenta of
emitted particles, which make the calculation of each diagram more
complicated, and can introduce divergences which need to be regulated.

In this chapter I’ll review the calculations for the inclusive process AX
including all next-to-leading order corrections. The leading order calculation
was performed in reference [94], and [82, 95, 96, 97] added various subleading
corrections. However, the complete LO+NLO result was only derived in [98]
(summarized in [99]), and numerically implemented in [100] with the program
SOLO, “Saturation at One-Loop Order” [101].

### Kinematics

As in section Kinematics of $\mathrm{p}A$ collisions, we take the proton and
nucleon to have momenta in the $\pm z$ directions,

$\displaystyle p_{\textmd{\Pproton}}^{\mu}$
$\displaystyle=(E_{p},0,0,p_{\textmd{\Pproton}})$ $\displaystyle
p_{\textmd{\HepGenParticle{A}{}{}\xspace}}^{\mu}$
$\displaystyle=(E_{A},0,0,-p_{\textmd{\HepGenParticle{A}{}{}\xspace}})$ (119)

leading to
$p_{\textmd{\Pproton}}^{+}=p_{\textmd{\HepGenParticle{A}{}{}\xspace}}^{-}=\sqrt{s}$
(equation (93)). And, again, the fundamental collision involves one parton
from the proton with plus component momentum $x_{p}p_{\textmd{\Pproton}}^{+}$,
and one parton from the nucleus, with a minus component of momentum
$x_{g}p_{\textmd{\HepGenParticle{A}{}{}\xspace}}^{-}$ and transverse momentum
$\vec{k}_{g\perp}$. $p^{\mu}$ represents the momentum of the produced hadron.

#### Leading Order

Applying conservation of momentum at the interaction vertex of figure 23,
assuming that $p_{p}^{-}=p_{\textmd{\HepGenParticle{A}{}{}\xspace}}^{+}=0$
(equations (91b) and (91c)), gives

$\displaystyle x_{p}^{\text{LO}}p_{\textmd{\Pproton}}^{+}$
$\displaystyle=k^{+}$ $\displaystyle
x_{g}^{\text{LO}}p_{\textmd{\HepGenParticle{A}{}{}\xspace}}^{-}$
$\displaystyle=k^{-}$ (120)

From equation (628), assuming the quark mass is negligible
($k^{\mu}k_{\mu}=0$), these become

$\displaystyle x_{p}^{\text{LO}}p_{\textmd{\Pproton}}^{+}$
$\displaystyle=k_{\perp}e^{Y}$ $\displaystyle
x_{g}^{\text{LO}}p_{\textmd{\HepGenParticle{A}{}{}\xspace}}^{-}$
$\displaystyle=k_{\perp}e^{-Y}$ (121)

and then using equation (93) and $k_{\perp}=\frac{p_{\perp}}{z}$ we get

$\displaystyle x_{p}^{\text{LO}}$
$\displaystyle=\frac{p_{\perp}}{z\sqrt{s}}e^{Y}$ $\displaystyle
x_{g}^{\text{LO}}$ $\displaystyle=\frac{p_{\perp}}{z\sqrt{s}}e^{-Y}$ (122)

Evidently, the measured momentum of the detected hadron determines everything
except one unknown momentum fraction. We can choose any of $z$,
$x_{p}^{\text{LO}}$, or $x_{g}^{\text{LO}}$ to parametrize this unknown degree
of freedom, and consequently, the formula for the cross section will involve a
single integral over the chosen variable.

Figure 23: Feynman diagram for AhX at leading order in the quark-quark channel

#### Next-to-Leading Order

At next-to-leading order, the momentum of the emitted gluon is also unknown,
which adds an additional degree of freedom, so the formula will involve an
additional integral. It’s most convenient to express the NLO formulas in terms
of a new variable, $\xi$, representing the fraction of $+$-component momentum
retained by the original parton after it emits a gluon. For example, in figure
25,

$\displaystyle k_{\textmd{\Pquark}}^{+}$ $\displaystyle=\xi k^{\prime+}$
$\displaystyle k_{\textmd{\Pgluon}}^{+}$ $\displaystyle=(1-\xi)k^{\prime+}$
(123)

We can then choose any two of $z$, $x_{p}^{\text{NLO}}$, $x_{g}^{\text{NLO}}$,
and $\xi$ to parametrize the region of integration for this double integral.
The integration region for various choices is shown in figure 24.

Figure 24: Domains of integration over the kinematic variables $z$,
$x_{p}^{\text{NLO}}$, and $\xi$. Any two of the variables are sufficient to
parametrize the integration domain, and the others can be taken as functions
of whichever two are chosen as primary. The different images show the shape of
the integration region for various choices of the two primary variables.

As a momentum fraction, $\xi$ clearly cannot exceed $1$. To determine the
lower bound on this new integral, though, we will need to express $\xi$ in
terms of the other kinematic variables. Starting with momentum conservation,

$\displaystyle x_{p}^{\text{NLO}}p_{\textmd{\Pproton}}^{+}$
$\displaystyle=k^{\prime+}=k_{\textmd{\Pquark}}^{+}+k_{\textmd{\Pgluon}}^{+}$
(124a) $\displaystyle
x_{g}^{\text{NLO}}p_{\textmd{\HepGenParticle{A}{}{}\xspace}}^{-}$
$\displaystyle=k^{\prime-}=k_{\textmd{\Pquark}}^{-}+k_{\textmd{\Pgluon}}^{-}$
(124b) $\displaystyle\vec{q}_{\perp}$
$\displaystyle=\vec{k}_{\perp}^{\prime}=\vec{k}_{\textmd{\Pquark}\perp}+\vec{k}_{\textmd{\Pgluon}\perp}$
(124c)

and plugging the definition (123) into equation (124a) leads to

$x_{p}^{\text{NLO}}p_{\textmd{\Pproton}}^{+}=k^{+}+(1-\xi)x_{p}^{\text{NLO}}p_{\textmd{\Pproton}}^{+}$
(125)

Solving this for $\xi$ and simplifying gives

$\displaystyle\xi$
$\displaystyle=\frac{k^{+}}{x_{p}^{\text{NLO}}p_{\textmd{\Pproton}}^{+}}$
(126)
$\displaystyle=\frac{\sqrt{k_{\perp}^{2}+m_{k}^{2}}}{x_{p}^{\text{NLO}}p_{\textmd{\Pproton}}^{+}}e^{Y}$
equation (628) (127)
$\displaystyle=\frac{k_{\perp}}{x_{p}^{\text{NLO}}p_{\textmd{\Pproton}}^{+}}e^{Y}$
$\displaystyle m_{k}\approx 0$ (128)
$\displaystyle=\frac{p_{\perp}}{zx_{p}^{\text{NLO}}p_{\textmd{\Pproton}}^{+}}e^{Y}$
$\displaystyle k_{\perp}=\frac{p_{\perp}}{z}$ (129)
$\displaystyle=\frac{p_{\perp}}{zx_{p}^{\text{NLO}}\sqrt{s}}e^{Y}$ equation
(93) (130)

Or, inverting the relation, we have

$x_{p}^{\text{NLO}}=\frac{p_{\perp}}{z\xi\sqrt{s}}e^{Y}$ (131)

Now looking at (124b), since the radiated gluon must be on shell
($k_{g}^{2}=0$) we have, for a gluon with nonzero momentum,

$k_{g}^{-}=\frac{k_{\Pgluon\perp}^{2}}{k_{g}^{+}}=\frac{k_{\Pgluon\perp}^{2}}{(1-\xi)x_{p}^{\text{NLO}}p_{\textmd{\Pproton}}^{+}}$
(132)

and assuming that the quark is massless, $m_{k}\approx 0$, in combination with
equation (126) we have

$k^{-}\approx\frac{k_{\perp}^{2}}{k_{+}}=\frac{k_{\perp}^{2}}{\xi
x_{p}^{\text{NLO}}p_{\textmd{\Pproton}}^{+}}$ (133)

Both of these contain the combination
$x_{p}^{\text{NLO}}p_{\textmd{\Pproton}}^{+}$ in the denominator, which can be
rewritten as

$\frac{1}{x_{p}^{\text{NLO}}p_{\textmd{\Pproton}}^{+}}=\frac{z\xi\sqrt{s}}{p_{\perp}}e^{-Y}\frac{1}{\sqrt{s}}=\frac{z\xi}{p_{\perp}}e^{-Y}$
(134)

Then plugging equation (134) into (132) and (133) and those into equation
(124b) gives

$x_{g}^{\text{NLO}}p_{\textmd{\HepGenParticle{A}{}{}\xspace}}^{-}=\frac{zk_{\perp}^{2}}{p_{\perp}}e^{-Y}+\frac{z\xi
k_{\Pgluon\perp}^{2}}{(1-\xi)p_{\perp}}e^{-Y}$ (135)

Since $z=\frac{p_{\perp}}{k_{\perp}}$, this simplifies to

$x_{g}^{\text{NLO}}p_{\textmd{\HepGenParticle{A}{}{}\xspace}}^{-}=k_{\perp}e^{-Y}+\frac{\xi
k_{\Pgluon\perp}^{2}}{(1-\xi)k_{\perp}}e^{-Y}$ (136)