# Cold nuclear matter effects on azimuthal decorrelation in heavy-ion
collisions

Néstor Armesto<EMAIL_ADDRESS>Florian Cougoulic
<EMAIL_ADDRESS>Bin Wu<EMAIL_ADDRESS>Instituto Galego
de Física de Altas Enerxías IGFAE, Universidade de Santiago de Compostela,
E-15782 Galicia-Spain

###### Abstract

The assumption of factorization lies at the core of calculations of medium
effects on observables computable in perturbative Quantum Chromodynamics. In
this work we examine this assumption, for which we propose a setup to study
hard processes and bulk nuclear matter in heavy-ion collisions on the same
footing using the Glauber modelling of heavy nuclei. To exemplify this
approach, we calculate the leading-order corrections to azimuthal
decorrelation in Drell-Yan and boson-jet processes due to cold nuclear matter
effects, not considering radiation. At leading order in both the hard momentum
scale and the nuclear size, the impact-parameter dependent cross section is
found to factorize for both processes. The factorization formula involves a
convolution of the hard cross section with the medium-modified parton
distributions, and, for boson-jet production, the medium-modified jet
function.

##  I Introduction

Factorization in Quantum Chromodynamics Collins:1989gx ; Sterman:1995fz ;
Stewart:2003eft ; Becher:2014oda  plays a central role in hadron collider
physics. While it has been proved in $l^{+}l^{-}$ collisions and Deep
Inelastic Scattering (DIS), in hadron-hadron collisions it has only been
vetted for processes where the final state consists of colorless particles.
Besides, situations are known where factorization fails, see e.g.
Collins:2007nk ; Rogers:2010dm ; Mulders:2011zt  and references therein.

Factorization is also widely assumed for hard processes in heavy-ion
collisions, but it has been analyzed more rarely that in nucleon-nucleon
scattering and most often discussed in the context of DIS and photoproduction,
see, e.g., Qiu:2001hj ; Cao:2020wlm  and refs. therein. In spite of this, its
application lies at the core of studies dedicated to the extraction of
properties of the bulk matter or medium produced in such collisions using hard
probes like jet or large transverse momentum hadron production. Checks of the
(lack of) validity of factorization in collisions involving nuclei require the
description of both hard processes and bulk matter in a theoretically well-
grounded, unified framework. This central issue remains, despite many efforts,
an open question.

Parton saturation, also known as the Color Glass Condensate (CGC), provides a
promising description of bulk matter (soft and semihard partons) in the early
stages of heavy-ion collisions Kovchegov:2012mbw . In this framework, the
production of these partons has been conventionally calculated using classical
color sources, as originally proposed in the McLerran-Venugopalan (MV) model
McLerran:1993ni ; McLerran:1993ka . Alternatively, parton saturation has also
been studied by treating heavy-nucleus states as uncorrelated nucleon states
Mueller:1989st ; Kovchegov:1996ty ; Kovchegov:1998bi . The latter approach
shares some modeling characteristics with the Glauber model Miller:2007ri .
The production of soft gluons at fixed order in the strong coupling
$\alpha_{s}$ has been carried out in refs. Kovchegov:2005ss ; Wu:2017rry ,
revealing qualitative differences from kinetic theory in the coupling
expansion within $\phi^{4}$ theory Kovchegov:2017way . Recently, the
broadening of high-energy partons in the early-stage color fields have been
studied in ref. Li:2021zaw ; Avramescu:2023qvv .

In this work, we investigate whether and how one may formulate hard processes
and bulk matter on the same footing by modelling heavy-nuclei as uncorrelated
nucleons. The focus will be on whether the impact-parameter dependent cross
section, defined in ref. Wu:2021ril , factorizes at fixed order in
$\alpha_{s}$ as well as at leading order in both the hard momentum scale and
the nuclear size. Instead of treating bulk matter and jets differently like in
the previous studies, we carry out a detailed calculation of azimuthal
decorrelation of the Drell-Yan (DY) and boson-jet processes in heavy-ion
collisions. Such an observable has been extensively studied in proton-proton
Mueller:1981fe ; Collins:1984kg ; Becher:2010tm ; Buffing:2018ggv ;
Sun:2018icb ; Chien:2019gyf ; Chien:2020hzh ; Chien:2022wiq  and high-energy
nuclear collisions Mueller:2016gko ; Mueller:2016xoc ; Chen:2016vem ;
Chen:2016cof ; Chen:2018fqu ; Arleo:2020rbm ; Benic:2022ixp ; Wang:2021jgm .
We consider the rescattering of the partons entering and leaving the hard
scattering with partons from different nucleon-nucleon collisions. We do not
include radiation.

The manuscript is organized as follows: In Sec. II we define the collision
parameter-dependent cross section in heavy ion collisions, the associated
observable of interest in this study, namely the azimuthal decorrelation of
the final state pair, which is either a Drell-Yan pair or a photon-jet pair.
In Sec. III we define the Feynman rules used in momentum space and in
coordinate space. In Sec. IV we implement the cold nuclear effects in our
observables and discuss their impact on the initial state distributions. These
are all the modifications necessary to discuss the Drell-Yan pair azimuthal
decorrelation. In Sec. V we implement the cold nuclear effects for the
hadronic final states, which are necessary to complete the discussion of the
photon-jet azimuthal decorrelation. Finally, we discuss the form of the
factorization obtained in this study.

##  II A unified description of hard processes and bulk matter

In this section, we present the main formula to calculate the impact-parameter
dependent cross section for hard processes in heavy-ion collisions. The
constituent nucleons of the colliding nuclei are taken to be uncorrelated,
following the perturbative approach to parton saturation in refs.
Mueller:1989st ; Kovchegov:1998bi  and the Glauber model Miller:2007ri .

###  II.1 The impact-parameter dependent cross section in heavy-ion collisions

In quantum theory, the impact-parameter dependent cross section for any
observable $O$ in the collision of two ultra-relativistic particles can be
generically defined as Wu:2021ril

$\displaystyle\frac{d\sigma}{d^{2}{\mathbf{b}}dO}=$
$\displaystyle\int\prod\limits_{f}\left[d\Gamma_{p_{f}}\right]\delta(O-O(\\{p_{f}\\})\langle\phi_{1}\phi_{2}|\hat{S}^{\dagger}|\\{p_{f}\\}\rangle\langle\\{p_{f}\\}|\hat{S}|\phi_{1}\phi_{2}\rangle$
$\displaystyle=$
$\displaystyle\int\prod\limits_{f}\left[d\Gamma_{p_{f}}\right]\delta(O-O(\\{p_{f}\\})\text{Tr}\big{[}\hat{S}^{\dagger}|\\{p_{f}\\}\rangle\langle\\{p_{f}\\}|\hat{S}|\phi_{1}\phi_{2}\rangle\langle\phi_{1}\phi_{2}|\big{]},$
(1)

where the wave packages of the colliding particles, denoted as $\phi_{i}$, are
required to be sufficiently localized in space in order to define the impact
parameter $\mathbf{b}$, and the phase-space measure for a particle with mass
$m$ is defined as

$\displaystyle\int
d\Gamma_{p}\equiv\int\frac{d^{4}p}{(2\pi)^{4}}(2\pi)\delta(p^{2}-m^{2})\theta(p^{0}).$
(2)

Note that the unity in $\hat{S}$ should be disregarded, meaning
$\hat{S}=1+i\hat{T}$ is replaced with $i\hat{T}$ to align with the
conventional cross section definition. Here and below, boldface letters denote
two-dimensional transverse vectors.

To calculate the impact-parameter dependent cross section in heavy-ion
collisions using perturbative QCD, it is essential to have detailed
information about the multi-parton distributions within nuclei. As such
information is unavailable, below we express the nuclear multi-parton
distributions in terms of those for their constituent nucleons by employing
expansions based on physical scales.

The first relevant scale is the binding energy per nucleon, denoted as $\Delta
E_{b}$. In the rest frame of heavy nuclei, $\Delta E_{b}$ is known to be of
the order of 10 MeV. This value provides an estimate of the typical reaction
time between nucleons $t_{n}\sim 1/\Delta E_{b}\approx 20$ fm/$c$. In the lab
frame of ultra-relativistic heavy-ion collisions, $t_{n}$ is further dilated.
Since $t_{n}$ is significantly longer than any other timescales, the
interaction and correlation between constituent nucleons in the same nucleus
will be neglected.

The second relevant length scale is the nuclear size $R$, which scales as
$A^{1/3}$ with $A$ the mass number. As the kinetic energy per nucleon is of
the same order as $\Delta E_{b}$, individual nucleons typically possess a
momentum $\Delta p\sim\sqrt{2m_{n}\Delta E_{b}}\sim 140$ MeV, where $m_{n}$
denotes the nucleon mass. Consequently, the nucleons are localized in space
with an uncertainty comparable to the nucleon size $\sim 1/\Lambda_{QCD}$. In
the large $R$ limit, one can conduct an expansion in terms of
$1/(R\Lambda_{QCD})$. Consequently, the size of individual nucleons can be
neglected at leading order in this expansion. Moreover, nucleons can be viewed
as classical particles because $\Delta p$ is also negligible compared to other
relevant momentum scales. The above expansions are, in essence, identical to
the modelling in Mueller:1989st ; Kovchegov:1998bi ; Miller:2007ri .

Below, we assume that nucleus 1 with mass number $A_{1}$ moves along the
positive $z$-axis, while nucleus 2 with mass number $A_{2}$ moves along the
negative $z$-axis. In this case, the density matrix for nucleus 1 can be
written in the form111For brevity, we do not distinguish neutrons from protons
and omit all spin indices.

$\displaystyle|\phi_{1}\rangle\langle\phi_{1}|$
$\displaystyle=\prod\limits_{i=1}^{A_{1}}\int
d\Gamma_{p_{i}}d\Gamma_{p^{\prime}_{i}}|p_{i}\rangle\langle
p_{i}|\phi_{1}\rangle\langle\phi_{1}|p^{\prime}_{i}\rangle\langle
p^{\prime}_{i}|=\prod\limits_{i=1}^{A_{1}}\int\frac{dP_{i}^{+}d^{2}\mathbf{P}_{i}}{(2\pi)^{3}2P_{i}^{+}}$
$\displaystyle\times\frac{1}{2}\int\,db_{i}^{-}d^{2}\mathbf{b}_{i}W_{A_{1}}(P_{i},b_{i})\int\frac{dq_{i}^{+}d^{2}\mathbf{q}_{i}}{(2\pi)^{3}}e^{\frac{i}{2}q_{i}^{+}b_{i}^{-}-i\mathbf{q}_{i}\cdot\mathbf{b}_{i}}|P_{i}+{q_{i}}/{2}\rangle\langle
P_{i}-{q_{i}}/{2}|,$ (3)

where in the final expression on the right-hand side, we have replaced
$1/p_{i}^{+}$ and $1/{p_{i}^{\prime}}^{+}$ in
$d\Gamma_{p_{i}}d\Gamma_{p^{\prime}_{i}}$ by $1/P_{i}^{+}$, the Wigner
distribution function for one nucleon is defined as

$\displaystyle
W_{A_{1}}(P,b)\equiv\int\frac{dq^{+}d^{2}\mathbf{q}}{(2\pi)^{3}2P^{+}}e^{-\frac{i}{2}q^{+}b^{-}+i\mathbf{q}\cdot\mathbf{b}}\langle
P+{q}/{2}|\phi_{1}\rangle\langle\phi_{1}|P-{q}/{2}\rangle,$ (4)

and the light-cone components of a four-vector $v^{\mu}$ read

$\displaystyle v^{\pm}=v^{0}\pm v^{3}.$ (5)

As argued above, the constituent nucleons can be viewed as classical particles
with Kovchegov:2013cva ; Wu:2017rry

$\displaystyle
W_{A_{1}}(p,b)=\hat{\rho}_{A_{1}}(b^{-},\mathbf{b})2(2\pi)^{3}\delta(p^{+}-P_{1}^{+})\delta^{(2)}(\mathbf{p}),$
(6)

where $P^{+}_{1}$ denotes the ”+” momentum of the nucleon, and
$\hat{\rho}_{A_{1}}=\rho_{A_{1}}/A_{1}$ with $\rho_{A_{1}}(b^{-},\mathbf{b})$
the distribution of nucleons (usually with a Woods-Saxon functional form). For
nucleus 2, one only needs to swap the ”+” and ”$-$” variables and replace 1 by
2 in the above formulae.

Figure 1:  Azimuthal decorrelation in heavy-ion collisions. In the frame shown
in this figure, the two colliding nuclei locate respectively at $\mathbf{0}$
and $\mathbf{b}$ in the transverse plane. For the two produced objects, the
one with momentum $p_{4}$ is moving along the negative $x$-axis while the
other one with momentum $p_{3}$ is moving predominantly along the positive
$x$-axis. The azimuthal decorrelation $\delta\varphi$ is determined by the
projection of $p_{3}$ onto the $y$-axis, denoted by $Q_{T}^{y}$.

Expressing the density matrices in terms of the Wigner functions for both
colliding nuclei in eq. (II.1) yields

$\displaystyle\frac{d\sigma}{d^{2}{\mathbf{b}}dO}=$
$\displaystyle\int\prod\limits_{f}\left[d\Gamma_{p_{f}}\right]\delta(O-O(\\{p_{f}\\})$
$\displaystyle\times\prod\limits_{i=1}^{A_{1}}\frac{1}{2P_{1}}\int\,d^{2}\mathbf{b}_{i}\,db_{i}^{-}\hat{\rho}_{A_{1}}(b_{i}^{-},\mathbf{b}_{i})\int\frac{dq_{i}^{+}d^{2}\mathbf{q}_{i}}{(2\pi)^{3}}e^{\frac{i}{2}q_{i}^{+}b_{i}^{-}\,-i\mathbf{q}_{i}\cdot\mathbf{b}_{i}}$
$\displaystyle\times\prod\limits_{j=1}^{A_{2}}\frac{1}{2P_{2}}\int\,d^{2}\mathbf{b}^{\prime}_{j}\,db^{\prime+}_{j}\hat{\rho}_{A_{2}}(b^{\prime+}_{j},\mathbf{b}^{\prime}_{j}-\mathbf{b})\int\frac{dq_{j}^{\prime-}d^{2}\mathbf{q}^{\prime}_{j}}{(2\pi)^{3}}e^{\frac{i}{2}q_{j}^{\prime-}\,b_{j}^{\prime+}\,-i\mathbf{q}^{\prime}_{j}\cdot\mathbf{b}^{\prime}_{j}}$
$\displaystyle\times\langle\\{P_{1}-{q_{i}}/{2}\\},\\{P_{2}-{q^{\prime}_{j}}/{2}\\}|\hat{S}^{\dagger}|\\{p_{f}\\}\rangle\langle\\{p_{f}\\}|\hat{S}|\\{P_{1}+{q_{i}}/{2}\\},\\{P_{2}+{q^{\prime}_{j}}/{2}\\}\rangle,$
(7)

where the momenta of the constituent nucleons within the two nuclei are
respectively given by $P_{1}^{\mu}=P_{1}^{+}n_{1}^{\mu}/2$ and
$P_{2}^{\mu}=P_{2}^{-}n_{2}^{\mu}/2$, with the two beam directions defined as
$n_{1}^{\mu}=(1,0,0,1)$ and $n_{2}^{\mu}=(1,0,0,-1)$. The impact-parameter
dependent cross section in the above form clearly respects translation
invariance. As illustrated in fig. 1, we choose a frame in which nuclei 1 and
2 are located respectively at $\mathbf{0}$ and $\mathbf{b}$ in the transverse
plane.

In the Glauber model Miller:2007ri , the observable $O$ represents
multiplicity or energy deposition in a specific region of the detector. Since
the impact parameter is correlated with this observable, collision geometry
(centrality) is hence determined within the uncertainties dictated by
fluctuations. In this work, we study the impact-dependent cross section for
hard processes with soft particles all integrated over, as discussed in detail
in ref. Wu:2021ril .

###  II.2 The observable: azimuthal decorrelation

To illustrate the utility of eq. (II.1), we carry out a detailed calculation
of the leading-order corrections to the azimuthal decorrelation
$\delta\varphi$ between two final-state objects produced in a hard process due
to the presence of nuclear matter. Such an observable has been proposed to
study the properties of QCD matter in high-energy nuclear collisions
Mueller:2016gko ; Mueller:2016xoc ; Chen:2016vem ; Chen:2016cof ; Chen:2018fqu
; Arleo:2020rbm ; Benic:2022ixp ; Wang:2021jgm . In heavy-ion collisions, the
hard process and bulk matter, conventionally studied separately using
different approaches Cao:2020wlm , are coupled via phenomenological modelling.
Here, we focus on how to formulate a unified description of both hard and soft
sectors of the entire collision.

The definition of azimuthal decorrelation is illustrated in fig. 1. Here, the
momenta for the incoming nuclei and the two final-state objects are
represented by $(A_{1}P_{1},A_{2}P_{2})\to(p_{3},p_{4})$, respectively. The
azimutal decorrelation $\delta\phi$ is determined by the component of
$\mathbf{p}_{3}$ perpendicular to $\mathbf{p}_{4}$, denoted as $Q_{T}^{y}$,
according to

$\displaystyle\delta\varphi=\arcsin(Q_{T}^{y}/|\mathbf{p}_{3}|)$ (8)

in the range $(-\pi,\pi)$. We focus on the kinematic region with
$\delta\varphi\approx Q^{y}_{T}/Q\ll 1$, with the hard scale
$Q=p_{T}\equiv|\mathbf{p}_{4}|\approx|\mathbf{p}_{3}|\gg|Q_{T}^{y}|$ in the
Drell-Yan and boson-jet processes. To justify perturbative calculations,
$Q_{T}^{y}$ is also taken to be semi-hard with $|Q_{T}^{y}|\gg\Lambda_{QCD}$.
In this back-to-back limit at hadron colliders, it is well-known that the
azimuthal distribution is subject to Sudakov suppression as a result of soft
and collinear radiation in the Drell-Yan Mueller:1981fe ; Collins:1985ue ;
Becher:2010tm  and boson-jet Buffing:2018ggv ; Sun:2018icb ; Chien:2019gyf ;
Chien:2020hzh ; Chien:2022wiq  processes. Below, we focus on additional
corrections at tree level due to the presence of nuclear matter.

By exploiting rotation symmetry of the entire collision system in the
transverse plane, the cross section for $\mathbf{b}$ and $\mathbf{p}_{4}$ at
arbitrary orientations can be derived from

$\displaystyle\frac{1}{p_{T}}\frac{d\sigma_{A_{1}A_{2}}}{d^{2}\mathbf{b}d\eta_{3}d\eta_{4}dp_{T}d\delta\varphi}\approx\frac{d\sigma_{A_{1}A_{2}}}{d^{2}\mathbf{b}d\eta_{3}d\eta_{4}dp_{T}dQ_{T}^{y}}$
(9)

with one of their orientations fixed. Here, $\eta_{i}$ is the pseudorapidity
of final-state object $i$. We choose the frame depicted in fig. 1, where
$\mathbf{p}_{4}$ is aligned with the negative $x$-axis, to carry out our
calculations.

At lowest order, there is no imbalance and the cross section for the
underlying partonic process $i(p_{1})+j(p_{2})\to k(p_{3})+l(p_{4})$ is given
by

$\displaystyle\frac{d\hat{\sigma}^{(0)}_{ij\to
kl}}{d\eta_{3}d\eta_{4}dp_{T}^{2}dQ_{T}^{y}}=$
$\displaystyle\frac{1}{16\pi\hat{s}}\frac{\delta(Q_{T}^{y})}{1+\delta_{kl}}|{\overline{M}^{(0)}_{ij\to
kl}(p_{1},p_{2}\to p_{3},p_{4})}|^{2}$
$\displaystyle\times\delta(p_{1}^{+}-p_{3}^{+}-p_{4}^{+})\delta(p_{2}^{-}-p_{3}^{-}-p_{4}^{-}),$
(10)

where the bar over the square of the lowest-order amplitude $M^{(0)}_{ij\to
kl}$ represents the average and sum over the initial- and final-state spins
and colors respectively, $\hat{s}=(p_{1}+p_{2})^{2}=(p_{3}+p_{4})^{2}$, and
$\delta_{kl}$ accounts for the statistical factor for identical particles. The
partonic process involves the following four light-like vectors:

$\displaystyle n_{1}^{\mu}=(1,0,0,1),\qquad n_{2}^{\mu}=(1,0,0,-1)$
$\displaystyle
n_{k}^{\mu}=(1,\sin\theta_{k},0,\cos\theta_{k})=(1,1/\cosh\eta_{k},0,\tanh\eta_{k})$
(11)

where $\theta_{k}$ denotes the polar angle of particle $k$ with $k=3,4$. The
four directions are assumed to be distinct, satisfying $n_{i}\cdot n_{j}\sim
1$ for $i\neq j$ with $i,j=1,\dots,4$. In the following discussion, the lowest
order expression in eq. (II.2) is referred to as the zeroth order, and we
shall evaluate the leading-order (LO) corrections that are enhanced by the
medium length in the Drell-Yan and boson-jet processes.

###  II.3 Perturbative description of hard processes and bulk matter

Figure 2: Unified description of hard processes and bulk matter. The
uncorrelated nucleons are categorized into three groups: the ones initiating
the hard processes (red), the ones participating the semi-hard (soft) process
(orange) and spectators (gray).

For azimuthal decorrelation there are two perturbative scales: the hard scale
$Q$ and the soft scale $Q_{T}^{y}$. As all other final-state particles are
integrated out, the cross section for such an observable is, in principle,
amenable to perturbative calculations. As demonstrated in fig. 2, a generic
diagram can be divided into two coupled parts: the hard process characterized
by $Q$, and the soft part characterized by $Q_{T}$. Even if the correlation
between constituent nucleons is neglected, we do not know a priori whether the
cross section in eq. (II.1) factorizes in the same manner as the corresponding
processes in cleaner types of collisions like those summarized in
Collins:1989gx ; Sterman:1995fz ; Stewart:2003eft ; Becher:2014oda ;
Kovchegov:2012mbw . Therefore, we restrict ourselves to carrying out the
leading-order corrections to the zeroth-order cross section in this work.

Let us now express the cross section in terms of the parton distribution
functions (PDFs) of the constituent nucleons. Following the parton model, we
first define the nucleon PDFs using free bare fields in the interaction
picture, which will be renormalized order by order in $\alpha_{s}$ whenever
necessary. It is convenient to work in coordinate space. As illustrated in
fig. 2, one can take a pair of fields from the vertices located at $x$ and
$x^{\prime}$ in the amplitude and the conjugate amplitude, respectively. After
acting the fields on one of the nucleon states in eq. (II.1), one has

$\displaystyle\int\,d^{2}\mathbf{b}\,db^{-}\hat{\rho}_{A_{1}}(b^{-},\mathbf{b})\int\frac{dq^{+}d^{2}\mathbf{q}}{(2\pi)^{3}}e^{\frac{i}{2}q^{+}b^{-}\,-i\mathbf{q}\cdot\mathbf{b}}$
$\displaystyle\times\sum\limits_{X_{h}}\langle
P_{1}-{q}/{2}|\bar{\psi}^{a}_{\alpha}(x^{\prime})|X_{h}\rangle\langle
X_{h}|\psi^{b}_{\beta}(x)|P_{1}+{q}/{2}\rangle$
$\displaystyle=\int\,d^{2}\mathbf{b}\,db^{-}\hat{\rho}_{A_{1}}(b^{-},\mathbf{b})\int\frac{dq^{+}d^{2}\mathbf{q}}{(2\pi)^{3}}e^{\frac{i}{2}q^{+}(b^{-}-X^{-})\,-i\mathbf{q}\cdot(\mathbf{b}-\mathbf{X})}$
$\displaystyle\times\sum\limits_{X_{h}}\langle
P_{1}-{q}/{2}|\bar{\psi}^{a}_{\alpha}(0)|X_{h}\rangle\langle
X_{h}|\psi^{b}_{\beta}(0)|P_{1}+{q}/{2}\rangle
e^{\frac{i}{2}(P_{1}^{+}-p^{+}_{X_{h}})\Delta x^{-}},$ (12)

where we consider, as an example, quark fields that carry both color (in
superscript) and spinor (in subscript) indices, $X_{h}$ represents the states
of the remnants, $X\equiv(x+x^{\prime})/2$ and $\Delta x\equiv x^{\prime}-x$.
In the last step, we have used the relation $\phi(x)=e^{i\hat{p}\cdot
x}\phi(0)e^{-i\hat{p}\cdot x}$ with $\phi$ denoting a quantum field and
$\hat{p}$ the four-momentum operator, as well as the fact that the nucleon and
the remnants with their total momentum denoted by $p_{X_{h}}$ are moving
predominantly along the positive $z$-axis.

As we expand in $1/(R\Lambda_{QCD})$, we can consistently neglect $q$ in the
hadron state in the above equation. This is justified by the fact that $q$
scales as the inverse of the nucleus size. Consequently, we obtain

$\displaystyle 2\hat{\rho}_{A_{1}}(X^{-},\mathbf{X})\langle
P_{1}|\bar{\psi}^{a}_{\alpha}(\Delta
x^{-}\bar{n}_{1}/2)\psi^{b}_{\beta}(0)|P_{1}\rangle$
$\displaystyle=2\hat{\rho}_{A_{1}}(X^{-},\mathbf{X})\frac{\delta^{ab}}{2N_{c}}P_{1}^{+}\frac{\not{n}_{1\beta\alpha}}{2}\int
d\xi f_{q}(\xi)e^{i\frac{\Delta x^{-}}{2}\xi P^{+}_{1}}$
$\displaystyle=2\hat{\rho}_{A_{1}}(X^{-},\mathbf{X})\int\frac{d\xi}{\xi}f_{q}(\xi)\frac{\delta^{ab}}{2N_{c}}\not{p}_{1}e^{ip_{1}\cdot(x^{\prime}-x)},$
(13)

where $\bar{n}_{1}^{\mu}\equiv(1,-\vec{n}_{1})=n_{2}^{\mu}$ and
$p^{\mu}_{1}=\xi P^{\mu}_{1}$, and $f_{i}$ represents the PDF for parton $i$
within the nucleon. Expanded at the same order, the quark distribution is
consistent with the definition in QCD Collins:1981uw , which is given by

$\displaystyle f_{q}(\xi)=\int\frac{dt}{2\pi}e^{-it\xi P_{1}^{+}}\langle
P|\bar{\psi}(t\bar{n}_{1})\mathcal{P}e^{-ig\int_{0}^{t}dt^{\prime}\bar{n}_{1}\cdot
A(t^{\prime}\bar{n}_{1})}\frac{\not{\bar{n}}_{1}}{2}\psi(0)|P\rangle$ (14)

with $\mathcal{P}$ denoting the path-ordering operator. Applying the above
reasoning to all fields acting on the hadron states within both nuclei, we
finally obtain

$\displaystyle\frac{d\sigma}{d^{2}{\mathbf{b}}dO}=$
$\displaystyle\int\prod\limits_{f}\left[d\Gamma_{p_{f}}\right]\delta(O-O(\\{p_{f}\\})$
$\displaystyle\times\sum\limits_{\\{a_{i},b_{j}\\}}\bigg{(}\prod\limits_{i=1}^{A_{1}}\frac{1}{P_{1}^{+}}\int\frac{d\xi_{i}}{\xi_{i}}f_{a_{i}}(\xi_{i})\bigg{)}\bigg{(}\prod\limits_{j=1}^{A_{2}}\frac{1}{P_{2}^{-}}\int\frac{d\xi^{\prime}_{j}}{\xi^{\prime}_{j}}f_{b_{j}}(\xi^{\prime}_{j})\bigg{)}$
$\displaystyle\times\langle\\{\xi_{i}P_{1}\\},\\{\xi^{\prime}_{j}P_{2}\\}|\hat{S}^{\dagger}|\\{p_{f}\\}\rangle\langle\\{p_{f}\\}|\hat{S}|\\{\xi_{i}P_{1}\\},\\{\xi^{\prime}_{j}P_{2}\\}\rangle$
$\displaystyle\otimes\bigg{(}\prod\limits_{i=1}^{A_{1}}\hat{\rho}_{A_{1}}(X^{-}_{i},\mathbf{X}_{i})\bigg{)}\bigg{(}\prod\limits_{j=1}^{A_{2}}\hat{\rho}_{A_{2}}(Y^{+}_{j},\mathbf{Y}_{j}-\mathbf{b})\bigg{)},$
(15)

where $a_{i}$ and $b_{j}$ iterate over all the parton species, and the
operator $\otimes$ indicates that the incoming partons $i$ and $j$ enter the
diagrams in the amplitude (the conjugate amplitude) at $x_{i}$
($x_{i}^{\prime}$) and $y_{j}$ ($y^{\prime}_{j}$) respectively with
$X_{i}=(x_{i}+x_{i}^{\prime})/2$ and $Y_{j}=(y_{j}+y^{\prime}_{j})/2$. Here,
the initial- and final-state spins and colors are respectively averaged and
summed over in the square of the partonic $S$-matrix element.

Since each of the $A_{i}$ nucleons in nucleus $i$ has an equal chance to
participate the hard collision, at zeroth order, one has

$\displaystyle\frac{d\sigma^{(0)}}{d^{2}{\mathbf{b}}dO}=$
$\displaystyle\int\prod\limits_{f}\left[d\Gamma_{p_{f}}\right]\delta(O-O(\\{p_{f}\\})\frac{1}{s_{NN}}\sum\limits_{ij}\int\frac{d\xi}{\xi}\frac{d\xi^{\prime}}{\xi^{\prime}}f_{i}(\xi)f_{j}(\xi^{\prime})$
$\displaystyle\times\langle\xi
P_{1},\xi^{\prime}P_{2}|\hat{S}^{\dagger}|\\{p_{f}\\}\rangle\langle\\{p_{f}\\}|\hat{S}|\xi
P_{1},\xi^{\prime}P_{2}\rangle\otimes{\rho}_{A_{1}}(X^{-},\mathbf{X}){\rho}_{A_{2}}(Y^{+},\mathbf{Y}-\mathbf{b})$
(16)

with $s_{NN}=P_{1}^{+}P_{2}^{-}$. In the expansion at large $Q$, offshell
propagators shrink to one single spacetime point, and we get

$\displaystyle\frac{d\sigma^{(0)}}{d^{2}{\mathbf{b}}dO}=$ $\displaystyle 2\int
d^{4}X{\rho}_{A_{1}}(X^{-},\mathbf{X}){\rho}_{A_{2}}(X^{+},\mathbf{X}-\mathbf{b})\frac{d{\sigma}^{(0)}_{nn\to
kl}}{dO},$ (17)

where the nucleon-nucleon cross section is defined as

$\displaystyle\frac{d{\sigma}^{(0)}_{nn\to
kl}}{dO}\equiv\sum\limits_{ij}\int{d\xi}{d\xi^{\prime}}f_{i}(\xi)f_{j}(\xi^{\prime})\frac{d\hat{\sigma}^{(0)}_{ij\to
kl}}{dO}.$ (18)

In eq. (II.3), the PDFs are defined in terms of free fields, and infrared
divergences arise at high orders in $\alpha_{s}$. In Deep Inelastic
Scattering, proton-proton and proton-nucleus collisions, it has been
justified, to varying degrees of rigor Collins:1984kg ; Collins:1989gx ;
Sterman:1995fz ; Stewart:2003eft ; Becher:2014oda ; Kovchegov:2012mbw , that
infrared divergences arsing beyond LO in $\alpha_{s}$ can all be absorbed into
universal parton distributions or jet functions defined by corresponding gauge
invariant operators. Accordingly, the cross section in such collisions
factorizes. In the subsequent sections, we shall examine whether the cross
section in eq. (II.3) for the Drell-Yan and boson-jet processes in heavy-ion
collisions factorizes at fixed order.

##  III Expansion of Feynman diagrams at high $Q$

All our calculations in the following sections are to be carried out according
to eq. (II.3). As we are interested in the kinematic region with
$\delta\varphi\ll 1$, we expand all the diagrams in $\delta\varphi$ and keep
only leading-order terms in this expansion. With $Q_{T}^{y}$ held fixed, the
diagrams are equivalently expanded at large $Q$. Furthermore, to be consistent
with the expansion in the nucleus radius as discussed in sec. II.1, only
contributions enhanced by the medium size will be kept among the leading-$Q$
terms. We present the Feynman rules based on this expansion below.

###  III.1 Expansion in momentum space

To organize the expansion of Feynman diagrams, let us first examine the
scaling of momenta in a generic diagram, as illustrated in fig. 2. There are
up to four light-like directions defined by the four light-like vectors
$n_{i}$. At leading order, one only needs to consider the contributions to the
azimuthal distribution due to momentum transfer $\sim\delta\varphi Q$ between
different collinear directions.

The collinear partons engaged in the hard collision, upon interaction with the
soft part of the diagrams, maintain their collinear directions. Consider a
parton $i$ along the light-like direction given by $n_{i}$. The momentum of
this $n_{i}$-collinear parton, denoted as $p$, scales as follows. It satisfies
$n_{j}\cdot p\sim Q$ for $j\neq i$ as all the four light-like directions are
assumed to be distinct. As we look for contributions enhanced by the medium
size, collinear internal momenta should be nearly on-shell in order to travel
for a long distance, so $2n_{i}\cdot p\,n_{j}\cdot p/n_{i}\cdot
n_{j}-|\mathbf{p}_{\perp}|^{2}\simeq 0$. That is, the $n_{i}$-collinear
momentum $p$ scales as

$\displaystyle n_{i}\cdot p\sim\delta\varphi^{2}Q,\qquad n_{j}\cdot p\sim
Q\qquad\text{for $j\neq i$},$ (19)

given the momentum transfer perpendicular to $n_{i}$ and $n_{j}$, denoted
above by $\mathbf{p}_{\perp}$, scales as $\delta\varphi Q$. In the expansion
at high $Q$, one can replace the propagator for a quark or antiquark with

$\displaystyle
S_{F}(p)=\not{p}D_{F}(p)\to\text{sign}(p^{0})\not{p}_{i}D_{F}(p),$ (20)

where the internal momentum $p$ is assumed to be aligned with the direction of
fermion-number flow such that $\text{sign}(p^{0})=1$ for quarks and
$\text{sign}(p^{0})=-1$ for antiquarks, $p_{i}$ is one of the incoming or
outgoing momenta engaging in the hard collision and

$\displaystyle D_{F}(p)\equiv\frac{i}{p^{2}+i\epsilon}.$ (21)

The collinear partons couple to the soft part of the diagrams through the
exchange of gluons, represented by $A^{a\mu}$. The Feynman rules for the
interaction between collinear partons and the exchanged gluon fields have been
widely employed in saturation physics Kovchegov:2012mbw  and jet quenching
physics Casalderrey-Solana:2007knd . Our emphasis here is on extending these
rules to processes that involve multiple distinct collinear directions using
the Feynman gauge. For an incoming $n_{i}$-collinear quark, upon expansion at
high $Q$, one has,

$\displaystyle\begin{array}[]{l}\includegraphics[width=65.04034pt]{image/qqg}\end{array}$
$\displaystyle=-ig\not{p}_{i}D_{F}(p)\not{A}u_{p_{i}}=-2ig(p_{i}\cdot\,A)D_{F}(p)u_{p_{i}}\,,$
(23)

with $A^{\mu}=A^{a\mu}t^{a}$, $t^{a}$ being the generators in the fundamental
representation of $SU(N_{c})$. Similarly, for an incoming antiquark, one only
needs to replace $t^{a}_{bc}$ with $\bar{t}^{a}_{bc}=-t^{a}_{cb}$ in the above
equation, along with a modification of the spinor. Here, the minus sign comes
from the sign of $p^{0}$ in eq. (20). The above rules also apply to outgoing
quarks and antiquarks with the corresponding spinors substitution.

The main complication, in comparison to the rules in physical gauges, arises
from the interaction between collinear gluons and the exchanged gluons. For an
incoming gluon, in the Feynman gauge, one has

$\displaystyle\begin{array}[]{l}\includegraphics[width=86.72267pt]{image/ggg}\end{array}$
$\displaystyle=gf^{abc}A^{a}_{\mu}[g^{\mu\nu}(k+p)^{\rho}+g^{\nu\rho}(-p-p_{i})^{\mu}+g^{\rho\mu}(p_{i}-k)^{\nu}]\epsilon_{\rho}(p_{i})D_{F}(p)$
(25) $\displaystyle\approx-gf^{abc}[2p_{i}\cdot
A^{a}\epsilon^{\nu}(p)-A^{a}\cdot\epsilon(p)~{}p_{i}^{\nu}]D_{F}(p).$ (26)

Given that physical polarization vectors $\epsilon^{\mu}$ are used for all the
incoming or outgoing gluons, the second term, contracted with the rest part of
the diagrams, vanishes (at leading order in $Q$) according to the Ward
identity (see, e.g., Collins:1989gx )

$\displaystyle\langle
M|T\partial_{\mu_{1}}A^{\mu_{1}}(x_{1})\cdots\partial_{\mu_{n}}A^{\mu_{n}}(x_{n})|N\rangle=0,$
(27)

given that $M$ and $N$ are physical states. That is, one can safely take

$\displaystyle\begin{array}[]{l}\includegraphics[width=86.72267pt]{image/ggg}\end{array}$
$\displaystyle=-2gf^{abc}(p_{i}\cdot
A^{a})\epsilon^{\nu}(p_{i})D_{F}(p)=-igT^{a}_{bc}2p_{i}\cdot
A^{a}~{}\epsilon^{\nu}(p_{i})D_{F}(p),$ (29)

in the expansion at large $Q$. This rule also applies to outgoing collinear
gluons.

As illustrated in fig. 2, the soft part of the diagrams can also couple to the
hard part by attaching $A^{a\mu}$ onto the internal propagator in the hard
process. However, these diagrams can be effectively eliminated through an
expansion in $Q$, as such an insertion only splits the off-shell propagator
into two off-shell propagators. Since neither of the resulting propagators are
close to being on-shell, these diagrams are not enhanced by the medium size,
in contrast with the insertions of the exchanged gluon onto the external legs.

As a result, it is sufficient to consider diagrams where the soft part couples
to the external legs of the hard part. As we see above, the spinors and
polarization vectors of all the incoming and outgoing partons pass through the
vertices on the external legs, contracting directly with the off-shell part of
the diagrams. Then, expanding the off-shell part at high $Q$ yields the
zeroth-order amplitudes $iM^{(0)}$. In this way, the medium corrections all
factorize from the amplitude for the hard process on the amplitude level. One
only needs to evaluate the interaction between the external legs and the gluon
fields $A^{a\mu}$ produced by the soft part of the diagrams using the
following rules:

$\displaystyle\begin{array}[]{l}\includegraphics[width=65.04034pt]{image/propagatorp}\end{array}=D_{F}(p),\qquad\begin{array}[]{l}\includegraphics[width=65.04034pt]{image/qqg}\end{array}=-ig\hat{T}^{a}2p_{i}\cdot
A^{a},$ (32)

where all collinear partons are represented by solid lines, and $\hat{T}^{a}$
stands for the color generator: $\hat{T}^{a}_{bc}=t^{a}_{bc}$ for quarks,
$\hat{T}^{a}_{bc}=-t^{a}_{cb}$ for antiquarks, and
$\hat{T}^{a}_{bc}=-if^{abc}$ for gluons.

The rest of this work focuses on interactions between collinear partons via
one gluon exchange. According to the Feynman rules in eq. (32), the
interaction between two collinear partons along the same direction is
negligible at high $Q$. One only needs to consider the interaction between
partons collinear to different directions. Given two different light-cone
directions $n_{i}$ and $n_{j}$, the momentum of the exchanged gluon scales as
$n_{i}\cdot p\sim\delta\varphi^{2}Q\sim n_{j}\cdot p$ according to momentum
conservation. That is, one has $n_{i}\cdot p\,n_{j}\cdot
p\sim\delta\varphi^{4}Q^{2}$, which can always be expanded. Therefore, for a
so-called Glauber gluon Bodwin:1981fv  connecting two collinear partons along
$n_{i}$ and $n_{j}$ respectively, one has

$\displaystyle\begin{array}[]{l}\includegraphics[width=65.04034pt]{image/GGlauber}\end{array}=\frac{ig_{\mu\nu}}{|\mathbf{p}_{\perp}|^{2}}$
(34)

with $|\mathbf{p}_{\perp}|^{2}=-g_{\perp}^{\mu\nu}p_{\mu}p_{\nu}$ and the
transverse metric tensor defined as

$\displaystyle
g^{\mu\nu}_{\perp}=g^{\mu\nu}-\frac{n_{i}^{\mu}n_{j}^{\nu}+n_{i}^{\nu}n_{j}^{\mu}}{n_{ij}}\qquad\text{with
$n_{ij}\equiv\,n_{i}\cdot\,n_{j}$}.$ (35)

Here and below, boldface letters with a subscript $\perp$ denote two-
dimensional vectors perpendicular to two light-like directions $n_{i}$ and
$n_{j}$. To be consistent with previous conventions, the subscript is omitted
if $n_{i}$ and $n_{j}$ are both the two beam directions. That is, the
subscript $\perp$ is reserved as a reminder that at least one of the final-
state collinear partons is involved in the interaction.

For the two beam directions $n_{1}$ and $n_{2}$, the transverse $x$ and $y$
components of a vector are defined as illustrated in fig. 1. When
$\vec{n}_{i}$ is not parallel to $\vec{n}_{j}$, we define the transverse ”1”
and ”2” components using the following two basis vectors:

$\displaystyle
n_{\perp_{1}}^{\mu}=\frac{1}{\sqrt{1-(\vec{n}_{i}\cdot\vec{n}_{j})^{2}}}\bigg{(}1+\vec{n}_{i}\cdot\vec{n}_{j},\vec{n}_{i}+\vec{n}_{j}\bigg{)},\qquad
n_{\perp_{2}}^{\mu}=\bigg{(}0,\frac{\vec{n}_{i}\times\vec{n}_{j}}{|\vec{n}_{i}\times\vec{n}_{j}|}\bigg{)}.$
(36)

As the four three-vectors $\vec{n}_{i}$ are in the same plane,
$n_{\perp_{2}}^{\mu}$ is universal (up to an overall minus sign) for any pair
of $n_{i}$ and $n_{j}$, and it is chosen to be aligned with the $y$-axis.
Accordingly, a four-vector $v^{\mu}$ can be decomposed as

$\displaystyle
v^{\mu}=\frac{v^{n_{i}}n_{j}^{\mu}+v^{n_{j}}n_{i}^{\mu}}{n_{ij}}+v^{1}n_{\perp_{1}}^{\mu}+v^{2}n_{\perp_{2}}^{\mu}\
,$ (37)

with

$\displaystyle
v^{n_{i}}\equiv\,n_{i}\cdot\,v,\qquad\,v^{n_{j}}\equiv\,n_{j}\cdot\,v,$ (38)

and the transverse vector $\mathbf{v}_{\perp}=(v^{1},v^{2})$ with
$|\mathbf{v}_{\perp}|^{2}=(v^{1})^{2}+(v^{2})^{2}$.

###  III.2 Expansion in coordinate space

It turns out to be more convenient to work in coordinate space. As the
amplitude for the hard process only depends on $Q$ and, hence, factorizes in
the high-$Q$ expansion, all the off-shell propagators shrink into one single
spacetime point. Now, one only needs to Fourier transform the building blocks
for the external legs in eqs. (32) and (34) into coordinate space. Or
equivalently, the $S$-matrix elements in eq. (II.3) are calculated directly in
coordinate space: the corresponding parton-level diagrams are evaluated using
the vertices given above in eq. (32) and the free propagators given below in
eqs. (49) and (50), convoluted with

$\displaystyle\begin{array}[]{l}\includegraphics[width=108.405pt]{image/convolution.pdf}\end{array}=\frac{1}{\bar{n}\cdot
P}\int\frac{d\xi}{\xi}f_{i}(\xi)\hat{\rho}_{A}(n\cdot(y+y^{\prime})/2,(\mathbf{y}+\mathbf{y}^{\prime})/2-\mathbf{b})$
(40)

for each incoming or outgoing collinear parton moving along the $n$ direction.
Here, $\bar{n}^{\mu}=(1,-\vec{n})$ and $\mathbf{b}$ stands for the transverse
location of nucleus $A$. All the incoming (outgoing) collinear partons
converge at (originate from) a single spacetime point representing the hard
process, with the coefficient given by $iM^{(0)}$. Note that $iM^{(0)}$ only
depends on the on-shell momenta of the incoming and outgoing particles.

The Feynman rules above are justified as follows. As the momentum transfer is
always between two collinear directions, in the presence of another collinear
direction $n_{j}$, an internal momentum $p$ collinear to the $n_{i}$-direction
scales as $p_{\perp}\sim\delta\varphi Q,$ $n_{i}\cdot
p\sim\delta\varphi^{2}Q$, $n_{j}\cdot p\sim Q$ for $j\neq i$, see eq. (19). In
coordinate space, one has

$\displaystyle x_{\perp}\sim 1/(\delta\varphi Q),\qquad n_{i}\cdot x\sim
1/Q,\qquad n_{j}\cdot x\sim 1/(\delta\varphi^{2}Q).$ (41)

For an $n_{i}$-collinear propagator coupled to another $n_{j}$-collinear
parton, we decompose $x^{\mu}=(x^{n_{i}},x^{n_{j}},x_{\perp})$ according to
eq. (37). From the scalar product $x\cdot
p=\frac{1}{n_{ij}}(x^{n_{i}}p^{n_{j}}+x^{n_{j}}p^{n_{i}})-x^{1}p^{1}-x^{2}p^{2}$,
one can see that the metric tensor in this coordinate system is given by

$\displaystyle g_{\mu\nu}=\left(\begin{array}[]{cccc}0&\frac{1}{n_{ij}}&0&0\\\
\frac{1}{n_{ij}}&0&0&0\\\ 0&0&-1&0\\\ 0&0&0&-1\end{array}\right).$ (46)

With this metric, the Feynman propagator $D_{F}$ in coordinate space can be
expressed as

$\displaystyle D_{F}(x)$
$\displaystyle=\int\frac{d^{4}p}{(2\pi)^{4}}\sqrt{|g|}\frac{ie^{-ip\cdot
x}}{p^{2}+i\epsilon}=\frac{1}{n_{ij}}\int\frac{dp^{n_{i}}dp^{n_{j}}d^{2}\mathbf{p}_{\perp}}{(2\pi)^{4}}\frac{ie^{-ip\cdot
x}}{p^{2}+i\epsilon}$
$\displaystyle=\int\frac{dp^{n_{j}}}{2\pi}\frac{1}{2p^{n_{j}}}e^{-\frac{i}{n_{ij}}p^{n_{j}}x^{n_{i}}}G_{0}(x^{n_{j}},\mathbf{x}_{\perp};p^{n_{j}})[\theta(p^{n_{j}})\theta(x^{n_{j}})-\theta(-p^{n_{j}})\theta(-x^{n_{j}})],$
(47)

where the propagator for a non-relativistic particle in two dimensions is
defined as

$\displaystyle
G_{0}(x^{n_{j}},\mathbf{x}_{\perp};p^{n_{j}})\equiv\int\frac{d^{2}\mathbf{p}_{\perp}}{(2\pi)^{2}}e^{-\frac{i}{n_{ij}}\frac{|\mathbf{p}_{\perp}|^{2}}{p^{n_{j}}}x^{n_{j}}+i\mathbf{p}_{\perp}\cdot\mathbf{x}_{\perp}}=\frac{n_{ij}p^{n_{j}}}{4\pi
ix^{n_{j}}}e^{\frac{i}{2}\frac{n_{ij}}{2}\frac{p^{n_{j}}}{x^{n_{j}}}|\mathbf{x}_{\perp}|^{2}}.$
(48)

As all the coordinate variables are eventually integrated over in our
calculations, $x$ can be so chosen that the corresponding initial- (final-)
state parton flows toward (away from) the hard vertex. In this case, the
$\theta$ functions in eq. (III.2) can be dropped irrespective of parton
species.

In the equation above, $x^{n_{j}}$ plays the role of time and $p^{n_{j}}$
corresponds to the mass (sometimes referred to as Galilean mass in the context
of Light-Cone Perturbation Theory Beuf:2016wdz ; Kogut:1969xa ). That is, the
$n_{j}$-component (resp. $n_{i}$-component) is an analog of the ”+” component
(resp. ”$-$” component) in light-cone coordinates. We have not specified the
relation between $1/p^{n_{i}}\sim 1/(\delta^{2}\varphi Q)$ and the length of
the QCD medium. For simplicity, we assume the path length to be much smaller
than $1/(\delta^{2}\varphi Q)$ in the following discussions, yielding

$\displaystyle
D_{F}(x)=\int\frac{dp^{n_{j}}}{2\pi}\frac{1}{2p^{n_{j}}}e^{-\frac{i}{n_{ij}}p^{n_{j}}x^{n_{i}}}\delta^{(2)}(\mathbf{x}_{\perp})\theta(x^{n_{j}})$
(49)

with $p^{n_{j}}>0$ as per the convention mentioned above.

In the coordinates describe before, the Glauber gluon propagator connecting
the two collinear directions $n_{i}$ and $n_{j}$ in eq. (34) takes the form

$\displaystyle G_{F}^{\mu\nu}(x)$
$\displaystyle=ig^{\mu\nu}n_{ij}\delta(x^{n_{i}})\delta(x^{n_{j}})F_{1}(|\mathbf{x}_{\perp}|),$
(50)

where in dimensional regularization with $d=2-2\epsilon$, $F_{1}$ is defined
as Li:2023vdj

$\displaystyle F_{1}(|\mathbf{x}|)$
$\displaystyle\equiv\mu^{2-d}\int\frac{d^{d}\mathbf{l}}{(2\pi)^{d}}\frac{e^{i\mathbf{l}\cdot\mathbf{x}}}{|\mathbf{l}|^{2}}=\frac{1}{4\pi}\frac{\Gamma(-\epsilon)}{(\pi|\mathbf{x}|^{2}\mu^{2})^{-\epsilon}}=-\frac{1}{4\pi}\bigg{[}\frac{1}{\epsilon}+\gamma_{E}+\ln(\pi|\mathbf{x}|^{2}\mu^{2})+O(\epsilon)\bigg{]}.$
(51)

Here, the infrared divergence results from replacing nucleon states by parton
states within them in eq. (II.3). It reflects the fact that the PDFs in that
equation need to be renormalized. Otherwise, the infrared divergence can be
regularized by modelling nucleons, e.g., with color-singlet dipoles
Mueller:1989st ; Kovchegov:1996ty ; Li:2023vdj .

##  IV Cold nuclear effects on spatial parton distributions

In this section, we study the azimuthal decorrelation in the Drell-Yan process
due to the presence of QCD matter. We first carry out a detailed calculation
of $O(\alpha_{s})$ corrections due to the interaction of the incoming partons
engaged in the hard collision with the nucleons they encounter. Then, we
proceed with the resummation of these terms in an expansion on the nuclear
size.

###  IV.1 Initial-state single scattering

####  IV.1.1 The single-scattering diagrams

Figure 3: Cold nuclear effects before the hard collision associated with
$i(p_{1})+j(p_{2})\to k(p_{3})+l(p_{4})$. The left panel show the real diagram
for single scattering in which one incoming parton exchanges one Glauber gluon
with one encounter nucleon in both the amplitude and the conjugate amplitude.
The right panel shows the corresponding virtual diagrams in which the incoming
parton does not interact with the nucleon either in the amplitude or the
conjugate amplitude.

In ultra-relativistic heavy-ion collisions, the longitudinal sizes of
colliding nuclei experience substantial Lorentz contraction. However, in
contrast to soft gluon production, the hard collision involving high-$x$
partons occurs within a distance even smaller than the contracted sizes of the
nuclei. With the increase in nuclear density compensating for the contracted
path length, the likelihood of an incoming parton interacting with nucleons
from the oppositely moving nucleus remains the same as that in DIS, current-
nucleus and nucleon-nucleus collisions Kovchegov:1998bi .

Let us start with diagrams in which one additional nucleon couples to the two
incoming partons before they engage in the hard collision. The additional
nucleon can be either the constituent of nucleus 1 or 2. The incoming partons
in eq. (II.3) all carry some intrinsic transverse momenta of the order of
$\Lambda_{QCD}$, which are negligible in comparison with $\delta\varphi Q$.
Therefore, one can expect that it is negligible for the hard process to
generate the measured $Q_{T}^{y}$ by simply absorbing one parton from the
addiontal nucleon. In this case, the leading-order correction due to nuclear
effects starts at $O(\alpha_{s}^{2})$ in the fixed-order expansion.

Among all the diagrams at $O(\alpha_{s}^{2})$, it is sufficient to include
those where each incoming hard parton interacts with the nucleon moving in the
opposite direction. Consequently, there are only two types of real diagrams to
be considered. Moreover, only the $t$-channel exchange graphs between the
incoming partons and the additional nucleons need to be kept, as any $s$\- or
$u$-channel exchange would put the ensuing external leg offshell. After
squaring these diagrams, the additional nucleons have to be the same in both
the amplitude and in the conjugate amplitude to ensure color neutrality, as
exemplified by the diagrams in fig. 3. In this figure, the corresponding
virtual diagrams are also included to guarantee unitarity.

The parton from the additional nucleon, after exchanging one gluon with the
hard process, is integrated out. Let us single out this part of diagrams, as
shown in fig. 3 (a). Without loss of generality, we take the parton from the
additional nucleon to be an valence quark in the following discussion.
Assuming the valence quark moving along the negative-$z$ axis, this part of
the diagram reads

$\displaystyle\begin{array}[]{l}\includegraphics[width=104.07117pt]{image/AAreal}\end{array}\equiv\frac{\delta^{ab}}{N_{c}^{2}-1}\langle
A^{c\mu}(x_{1})A^{c\nu}(x^{\prime}_{1})\rangle=\frac{A_{2}-1}{A_{2}}\frac{1}{P_{2}^{-}}\int\frac{d\xi}{\xi}f_{q}(\xi)\int\,d\Gamma_{p}$
(53)
$\displaystyle\times\frac{\text{Tr}(t^{a}t^{b})}{2N_{c}}g^{2}\text{Tr}(\not{p}\gamma_{\alpha}\not{p}_{2}\gamma_{\beta})\int
d^{4}yd^{4}y^{\prime}\rho_{A_{2}}(Y^{+},\mathbf{Y}-\mathbf{b})e^{i\Delta
y\cdot(p-p_{2})}G_{F}^{\mu\alpha}(x_{1}-y)G_{F}^{\nu\beta}(x_{1}^{\prime}-y^{\prime})$

with $p_{2}=\xi P_{2}$, $Y=(y+y^{\prime})/2$ and $\Delta y=y-y^{\prime}$.
Here, the factor $(A_{2}-1)/A_{2}$ arises from the fact that one out of the
$A_{2}$ nucleons participates in the hard process. This factor is taken to be
unity below, as we only consider cases with $A_{1,2}\gg 1$. We then expand
away high order terms in $\delta\varphi$ in the parton-level amplitude
according to the scaling: $p^{-}\sim Q$, $|\mathbf{p}|\sim\delta\varphi$ and
$p^{+}\sim\delta\varphi^{2}Q$ and drop $p^{+}$ everywhere in the above
equation. After using the expression for the Glauber gluon propagator in eq.
(50), one can easily obtain

$\displaystyle\langle A^{c\mu}(x_{1})A^{c\nu}(x^{\prime}_{1})\rangle=$
$\displaystyle
g^{2}C_{F}n_{2}^{\mu}n_{2}^{\nu}\delta(x_{1}^{+}-{x_{1}^{\prime}}^{+})$
$\displaystyle\times\int
d^{2}\mathbf{y}\rho_{A_{2}}(x_{1}^{+},\mathbf{y}-\mathbf{b})\int d\xi
f_{q}(\xi)F_{1}(|\mathbf{x}_{1}-\mathbf{y}|)F_{1}(|\mathbf{x}^{\prime}_{1}-\mathbf{y}|).$
(54)

It is equivalently given by the product of the classical fields in ref.
Kovchegov:1996ty . That is, the valence quark behaves like a classical color
charge moving along the light cone, which is identical to that in the MV model
McLerran:1993ni ; McLerran:1993ka .

We first evaluate the real diagram squared in fig. 3 (a), using eq. (II.3)
along with Feynman rules in coordinate space as presented in sec. III.2. In
this case, the incoming hard parton $i$ is moving along $n_{i}=n_{1}$ while
the additional nucleon is moving along $n_{j}=n_{2}$. Accordingly, we use the
light-cone variables: $v^{n_{i}}=v^{-}$ and $v^{n_{j}}=v^{+}$. In terms of the
above field correlator, one has

$\displaystyle\text{fig.~{}\ref{fig:LO_before}~{}(a)}=$
$\displaystyle\frac{2g^{2}C_{i}}{N_{c}^{2}-1}\int d\xi d\xi^{\prime}\int
d^{4}xd^{4}x^{\prime}d^{4}x_{1}d^{4}x^{\prime}_{1}(p_{1}^{+})^{2}\hat{H}_{ij\to
kl}(X_{1}^{-},\mathbf{X}_{1},\xi;X^{+},\mathbf{X},\xi^{\prime})$
$\displaystyle\times\langle A^{a-}(x_{1})A^{a-}(x_{1}^{\prime})\rangle
D_{F}(x-x_{1})D_{F}^{*}(x^{\prime}-x_{1}^{\prime})e^{-i(x-x^{\prime})\cdot(p_{2}-p_{3}-p_{4})-i(x_{1}-x^{\prime}_{1})\cdot
p_{1}},$ (55)

where $X\equiv(x+x^{\prime})/2$, $X_{1}\equiv(x_{1}+x_{1}^{\prime})/2$,
$C_{i}$ is the quadratic Casimir corresponding to the color representation of
parton $i$, and

$\displaystyle\hat{H}_{ij\to kl}\equiv$
$\displaystyle\frac{f_{i}(\xi)f_{j}(\xi^{\prime})}{2\xi\xi^{\prime}s_{NN}}\rho_{A_{1}}(X_{1}^{-},\mathbf{X}_{1})\rho_{A_{2}}(X^{+},\mathbf{X}-\mathbf{b})|\overline{M}^{(0)}_{ij\to
kl}(p_{1},p_{2};p_{3},p_{4})|^{2}$ (56)

with $s_{NN}=(P_{1}+P_{2})^{2}=P_{1}^{+}P_{2}^{-}$, $p_{1}=\xi P_{1}$ and
$p_{2}=\xi^{\prime}P_{2}$. One can then integrate over $(x^{-},\mathbf{x})$
and $(x^{\prime-},\mathbf{x}^{\prime})$ to obtain

$\displaystyle\text{fig.~{}\ref{fig:LO_before}~{}(a)}=$
$\displaystyle\frac{g^{2}C_{i}}{2(N_{c}^{2}-1)}\int
dx^{+}d{x^{\prime}}^{+}d^{4}x_{1}d^{4}x^{\prime}_{1}\int d\xi
d\xi^{\prime}\hat{H}_{ij\to
kl}(X^{-},\mathbf{X},\xi;X^{+},\mathbf{X},\xi^{\prime})$
$\displaystyle\times\langle
A^{a-}(x_{1})A^{a-}(x_{1}^{\prime})\rangle\theta(x^{+}-x^{+}_{1})\theta({x^{\prime}}^{+}_{1}-{x_{1}^{\prime}}^{+})\bigg{(}\frac{p_{1}^{+}}{p_{3}^{+}+p_{4}^{+}}\bigg{)}^{2}$
$\displaystyle\times
e^{\frac{i}{2}(x^{+}-{x^{\prime}}^{+})(p_{3}^{-}+p_{4}^{-}-p_{2}^{-})-i(\mathbf{x}_{1}-\mathbf{x}^{\prime}_{1})\cdot\mathbf{Q}_{T}+\frac{i}{2}(x_{1}^{-}-{x^{\prime}}_{1}^{-})(p_{3}^{+}+p_{4}^{+}-p_{1}^{+})},$
(57)

where $\mathbf{Q}_{T}\equiv\mathbf{p}_{3}+\mathbf{p}_{4}$, and in light-cone
coordinates
$X^{\mu}\equiv((x^{+}+{x^{\prime}}^{+})/2,(x_{1}^{-}+{x^{\prime}}_{1}^{-})/2,(\mathbf{x}_{1}+\mathbf{x}_{1}^{\prime})/2)$
can be viewed as the spacetime point where the hard collision takes place.

The above result applies to all field correlators. To proceed, one can make
use of the properties of the correlator in eq. (IV.1.1), specifically its
proportionality to $\delta(x_{1}^{+}-{x^{\prime}}_{1}^{+})$ and independence
of $x_{1}^{-}$ and ${x^{\prime}}_{1}^{-}$. Consequently, by integrating over
$x_{1}^{-}-{x^{\prime}}_{1}^{-}$ and $x^{+}-{x^{\prime}}^{+}$, one obtains

$\displaystyle\text{fig.~{}\ref{fig:LO_before}~{}(a)}=$
$\displaystyle\frac{g^{2}C_{i}}{N_{c}^{2}-1}\int
d^{4}X\rho_{A_{1}}(X^{-},\mathbf{X})\rho_{A_{2}}(X^{+},\mathbf{X}-\mathbf{b})\int_{-\infty}^{X^{+}}dX_{1}^{+}d\sigma^{(0)}_{ij\to
kl}(X^{+}-X_{1}^{+})$ $\displaystyle\times\int
dtd^{2}\mathbf{x}e^{-i\mathbf{x}\cdot\mathbf{Q}_{T}}\langle
A^{a-}(X_{1}^{+}+t,\mathbf{X}+{\mathbf{x}}/{2})A^{a-}(X_{1}^{+}-t,\mathbf{X}-{\mathbf{x}}/{2})\rangle,$
(58)

where $t\equiv(x_{1}^{+}-{x^{\prime}}_{1}^{+})/2$, $\mathbf{x}$ is redefined
as $\mathbf{x}_{1}-\mathbf{x}^{\prime}_{1}$ in eq. (56), and the zeroth-order
differential cross section is given by

$\displaystyle d\sigma^{(0)}_{ij\to kl}(X^{+})\equiv$
$\displaystyle\int\frac{d\xi}{\xi}\frac{d\xi^{\prime}}{\xi^{\prime}}\frac{1}{2s_{NN}}f_{i}(\xi)f_{j}(\xi^{\prime})|\overline{M}^{(0)}_{ij\to
kl}(p_{1},p_{2};p_{3},p_{4})|^{2}$
$\displaystyle\times(4\pi)\delta(p_{1}^{+}-p_{3}^{+}-p_{4}^{+})\frac{2\sin[X^{+}(p_{2}^{-}-p_{3}^{-}-p_{4}^{-})]}{p_{2}^{-}-p_{3}^{-}-p_{4}^{-}}.$
(59)

Note that the ”$-$” momentum is not conserved exactly in $d\sigma^{(0)}_{ij\to
kl}(X,X_{1}^{+})$, in contrast to the full result at the zeroth order in eq.
(17). This is a consequence of the uncertainty principle as the nucleon are
confined in the nucleus. However, this should not be taken literally as we
have already expanded away higher-order terms in $1/(RQ)$.

Similarly, for the virtual diagrams in fig. 3 (b), one can first evaluate the
required field correlator to obtain

$\displaystyle\begin{array}[]{l}\leavevmode\hbox to14.44pt{\vbox
to16.87pt{\pgfpicture\makeatletter\hbox{\hskip 7.22191pt\lower-13.53267pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{0.0pt}{0.0pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{\includegraphics[width=216.81pt]{image/AAvirtual}}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-3.8889pt}{-9.36633pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{{+}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\end{array}=$
$\displaystyle-\frac{\delta^{ab}}{N_{c}^{2}-1}\langle
A^{c\mu}(x_{1})A^{c\nu}(x^{\prime}_{1})\rangle,$ (61)

where we have used the relation $\lim\limits_{{x^{\prime}}_{1}^{+}\to
x_{1}^{+}}[\theta(x_{1}^{+}-{x^{\prime}}_{1}^{+})+\theta({x^{\prime}}_{1}^{+}-x_{1}^{+})]=1$.
At the end, including the contribution from virtual diagrams as illustrated in
fig. 3 (b) yields222This result has been obtained by imposing unitarity. We
find that it is necessary to require the Glauber exchange to occur prior to
the ”$+$” time of the hard process in both the amplitude and the conjugate
amplitude for the virtual diagrams.

$\displaystyle\text{fig.~{}\ref{fig:LO_before}}=$
$\displaystyle-\frac{g^{2}C_{i}}{2(N_{c}^{2}-1)}\int
d^{4}X\rho_{A_{1}}(X^{-},\mathbf{X})\rho_{A_{2}}(X^{+},\mathbf{X}-\mathbf{b})\int_{-\infty}^{X^{+}}dX_{1}^{+}d\sigma^{(0)}_{ij\to
kl}(X^{+}-X_{1}^{+})$ $\displaystyle\times\int
dtd^{2}\mathbf{x}e^{-i\mathbf{x}\cdot\mathbf{Q}_{T}}\langle[A^{a-}(X_{1}^{+}+t,\mathbf{X}+{\mathbf{x}}/{2})-A^{a-}(X_{1}^{+}-t,\mathbf{X}-{\mathbf{x}}/{2})]^{2}\rangle.$
(62)

The contributions from the diagrams in which the incoming parton $j$ interacts
with an additional nucleon from nucleus 1 can be obtained by replacing $C_{i}$
with $C_{j}$ and swapping the indices of nuclei 1 and 2 as well as the ”+” and
”-” variables.

####  IV.1.2 The phase-space gluon distributions in heavy nuclei

In hadron-hadron collisions, the only dominate source for a collinear parton
to acquire a transverse momentum of the order of $\delta\varphi
Q\gg\Lambda_{QCD}$ is radiation. Such contributions can be systematically
organized based on Transverse Momentum Dependent (TMD) factorization
Collins:1984kg ; Becher:2010tm ; Buffing:2018ggv ; Sun:2018icb ; Chien:2019gyf
; Chien:2020hzh ; Chien:2022wiq . In the presence of other nucleons in heavy-
ion collisions, the physical picture revealed by the result in eq. (IV.1.1) is
qualitatively different: each incoming hard parton scatters with the phase-
space gluon distributions in the oppositely-moving nucleus before
participating in the hard collision.

In order to substantiate this interpretation, let us evaluate

$\displaystyle G$ $\displaystyle(X_{1}^{+},\mathbf{X},\mathbf{x})\equiv\int
dt\langle[A^{a-}(X_{1}^{+}+t,\mathbf{X}+{\mathbf{x}}/{2})-A^{a-}(X_{1}^{+}-t,\mathbf{X}-{\mathbf{x}}/{2})]^{2}\rangle$
$\displaystyle=2g^{2}C_{F}\int
d^{2}\mathbf{r}\rho_{A_{2}}(X_{1}^{+},\mathbf{X}-\mathbf{r})\int d\xi
f_{q}(\xi)[F_{1}(|\mathbf{r}+\mathbf{x}/2|)-F_{1}(|\mathbf{r}-\mathbf{x}/2|)]^{2}$
$\displaystyle=\frac{\alpha_{s}C_{F}}{2\pi}\int
d^{2}\mathbf{r}\rho_{A_{2}}(X_{1}^{+},\mathbf{X}-\mathbf{r})\int d\xi
f_{q}(\xi)\ln^{2}\bigg{(}\frac{|\mathbf{r}+\mathbf{x}/2|^{2}}{|\mathbf{r}-\mathbf{x}/2|^{2}}\bigg{)},$
(63)

where the vector $\mathbf{r}\equiv\mathbf{X}-\mathbf{y}$ is the transverse
distance vector between the two scattering partons and the impact parameter
$\mathbf{b}$ is omitted. One can easily see that the above integral contains a
large logarithm with an non-perturbative infrared cutoff of $O(1/R)$, as
prescribed by $\rho_{A_{2}}$. Such a large logarithm needs to be absorbed into
the PDFs, and we accordingly introduce a factorization scale $\mu$ to
renormalize the PDFs.333 The factorization scale $\mu$ should be of order
$\Lambda_{QCD}$, which could be set in by modelling the nucleons by color
singlet dipoles Mueller:1989st ; Kovchegov:1996ty ; Li:2023vdj  For the
logarithmic term, one can expand $\rho_{A_{2}}$ around $\mathbf{r}=0$. In this
expansion, one has

$\displaystyle G$
$\displaystyle(X_{1}^{+},\mathbf{X},\mathbf{x})=\pi|\mathbf{x}|^{2}\rho_{A_{2}}(X_{1}^{+},\mathbf{X})\lim\limits_{x\to
0}xG(x,1/|\mathbf{x}|^{2})+\text{nonlogarithmic terms},$ (64)

where the gluon distribution generated by the quark in a nucleon is given by

$\displaystyle
xG(x,1/|\mathbf{x}|^{2})=\frac{\alpha_{s}C_{F}}{\pi}\ln\bigg{(}\frac{1}{\mu^{2}|\mathbf{x}|^{2}}\bigg{)}\int
d\xi f_{q}(\xi,\mu).$ (65)

Although it coincides with the solution of the DGLAP equation at
$\mathcal{O}(\alpha_{s})$ for $x\ll 1$, it represents the Fourier transform of
the transverse-momentum dependent gluon distribution under the same
approximation. Consequently, one can recognize the Fourier transform of the
phase-space gluon distribution, or the thickness beam function, within the
nucleus as defined in ref. Wu:2021ril

$\displaystyle
xf_{g/A_{2}}(X_{1}^{+},\mathbf{X},x,\mathbf{x})=\rho_{A_{2}}(X_{1}^{+},\mathbf{X})xG(x,1/|\mathbf{x}|^{2})$
(66)

at $x\ll 1$.

####  IV.1.3 Azimuthal decorrelation due to single scattering in the Drell-Yan
process

Let us insert the logarithmic term in eq. (64) into the result of the single
scattering in eq. (IV.1.1) and specify the observable $O$ as $\mathbf{Q}_{T}$
at given $p_{T}$, $\eta_{3}$ and $\eta_{4}$ in eq. (II.3). In terms of the
zeroth-order cross section:

$\displaystyle\frac{d\sigma^{(0)}_{ij\to
l^{+}l^{-}}}{d\eta_{3}d\eta_{4}dp_{T}}(X^{+})\equiv\frac{p_{T}}{4(2\pi)^{3}}d\sigma^{(0)}_{ij\to
l^{+}l^{-}}(X^{+})$ (67)

with $d\sigma^{(0)}_{ij\to l^{+}l^{-}}(X^{+})$ given in eq. (IV.1.1) for fig.
3, incorporating additional contributions by swapping nuclei 1 and 2 finally
yields

$\displaystyle\frac{d\sigma^{(1)}_{A_{1}A_{2}\to
l^{+}l^{-}}}{d^{2}\mathbf{b}d\eta_{3}d\eta_{4}dp_{T}d^{2}\mathbf{Q}_{T}}=$
$\displaystyle-\frac{1}{2}\int
d^{4}X\rho_{A_{1}}(X^{-},\mathbf{X})\rho_{A_{2}}(X^{+},\mathbf{X}-\mathbf{b})\int\frac{d^{2}\mathbf{x}}{(2\pi)^{2}}e^{-i\mathbf{x}\cdot\mathbf{Q}_{T}}|\mathbf{x}|^{2}$
$\displaystyle\times$
$\displaystyle\sum\limits_{ij}\bigg{[}\int_{-\infty}^{X^{+}}dX_{1}^{+}\frac{d\sigma^{(0)}_{ij\to
l^{+}l^{-}}}{d\eta_{3}d\eta_{4}dp_{T}}(X^{+}-X_{1}^{+})\hat{q}_{i/A_{2}}(X_{1}^{+},\mathbf{X}-\mathbf{b},|\mathbf{x}|)$
$\displaystyle+\int_{-\infty}^{X^{-}}dX_{1}^{-}\frac{d\sigma^{(0)}_{ij\to
l^{+}l^{-}}}{d\eta_{3}d\eta_{4}dp_{T}}(X^{-}-X_{1}^{-})\hat{q}_{j/A_{1}}(X_{1}^{-},\mathbf{X},|\mathbf{x}|)\bigg{]},$
(68)

where the ”+” and ”-” momenta are interchanged in the zeroth-order cross
section in the last term on the right-hand side, and the jet quenching
parameter is given by

$\displaystyle\hat{q}_{j/A_{i}}(X^{\pm},\mathbf{X},|\mathbf{x}|)=\frac{4\pi^{2}\alpha_{s}C_{j}}{N_{c}^{2}-1}\rho_{A_{i}}(X^{\pm},\mathbf{X})xG(x,1/|\mathbf{x}|^{2}).$
(69)

This parameter is the same as defined in ref. Baier:1996sk , except that the
scale in the gluon distribution is dictated by the observable with
$1/|\mathbf{x}|\sim Q_{T}$.

###  IV.2 Multiple scattering for the Drell-Yan process in large nuclei

Let us now focus on large nuclei and restrict our consideration to leading-
order corrections in the nuclear size $R$. In this case, for typical values of
$X^{+}$ and $X_{1}^{+}$ in eq. (IV.1.1), one can make the following
replacement

$\displaystyle\frac{2\sin[(X^{+}-X_{1}^{+})(p_{2}^{-}-p_{3}^{-}-p_{4}^{-})]}{p_{2}^{-}-p_{3}^{-}-p_{4}^{-}}\to
2\pi\delta(p_{2}^{-}-p_{3}^{-}-p_{4}^{-}).$ (70)

In this approximation, the zeroth-order nucleon-nucleon cross section in eq.
(67) is the same as that in the zeroth-order nucleus-nucleus collisions, as
presented in eq. (17). Combining the zeroth-order and the single-scattering
results results in

$\displaystyle\frac{d\sigma_{A_{1}A_{2}\to
l^{+}l^{-}}}{d^{2}\mathbf{b}d\eta_{3}d\eta_{4}dp_{T}d^{2}\mathbf{Q}_{T}}=$
$\displaystyle 2\int
d^{4}X\rho_{A_{1}}(X^{-},\mathbf{X})\rho_{A_{2}}(X^{+},\mathbf{X}-\mathbf{b})\sum\limits_{ij}\frac{d\sigma^{(0)}_{ij\to
l^{+}l^{-}}}{d\eta_{3}d\eta_{4}dp_{T}}$ $\displaystyle\times$
$\displaystyle\int\frac{d^{2}\mathbf{x}}{(2\pi)^{2}}e^{-i\mathbf{x}\cdot\mathbf{Q}_{T}}\bigg{[}1-\frac{|\mathbf{x}|^{2}}{4}\int_{-\infty}^{X^{+}}dX_{1}^{+}\hat{q}_{i/A_{2}}(X_{1}^{+},\mathbf{X}-\mathbf{b},|\mathbf{x}|)$
$\displaystyle-\frac{|\mathbf{x}|^{2}}{4}\int_{-\infty}^{X^{-}}dX_{1}^{-}\hat{q}_{j/A_{1}}(X_{1}^{-},\mathbf{X},|\mathbf{x}|)\bigg{]}.$
(71)

As each incoming hard parton only interacts with oppositely-moving nucleons,
the cold nuclear effects on the two incoming partons do not interfere.
Therefore, one can straightforwardly iterate our calculations for single
scattering to obtain the results for muliple scatterings. Given that the
average number of scattering is significantly smaller than the mass numbers
$A_{i}$, the summation of multiple scattering results can be substituted by
the exponentiation of single scattering results, namely:

$\displaystyle\frac{d\sigma_{A_{1}A_{2}\to
l^{+}l^{-}}}{d^{2}\mathbf{b}d\eta_{3}d\eta_{4}dp_{T}d^{2}\mathbf{Q}_{T}}=$
$\displaystyle 2\int
d^{4}X\rho_{A_{1}}(X^{-},\mathbf{X})\rho_{A_{2}}(X^{+},\mathbf{X}-\mathbf{b})\sum\limits_{ij}\frac{d\sigma^{(0)}_{ij\to
l^{+}l^{-}}}{d\eta_{3}d\eta_{4}dp_{T}}\int\frac{d^{2}\mathbf{x}}{(2\pi)^{2}}$
$\displaystyle\times$ $\displaystyle
e^{-i\mathbf{x}\cdot\mathbf{Q}_{T}-\frac{|\mathbf{x}|^{2}}{4}\int_{-\infty}^{X^{+}}dX_{1}^{+}\hat{q}_{i/A_{2}}(X_{1}^{+},\mathbf{X}-\mathbf{b},|\mathbf{x}|)-\frac{|\mathbf{x}|^{2}}{4}\int_{-\infty}^{X^{-}}dX_{1}^{-}\hat{q}_{j/A_{1}}(X_{1}^{-},\mathbf{X},|\mathbf{x}|)}.$
(72)

That is, each incoming parton, before participating in the hard collision,
experiences an equivalent amount of momentum broadening as an individual
parton in QCD medium, as studied in refs. Baier:1996sk ; Kovchegov:1998bi .

In Balitsky:2017gis ; Balitsky:2023hmh ; Balitsky:2024ozy , higher twist
corrections to the TMD involved in Drell-Yan were computed in a rapidity
factorization framework. These corrections have the symbolic form

$\langle F^{+\perp}F^{+\perp}\rangle_{A}\langle
F^{-\perp}F^{-\perp}\rangle_{B}+\langle
F^{+\perp}F^{+\perp}F^{+\perp}\rangle_{A}\langle
F^{-\perp}F^{-\perp}F^{-\perp}\rangle_{B}$ (73)

in the case of gluon TMD factorization. It is understood that Wilson lines in
a Drell-Yan staple are involved in the above expression, but are omitted here
for clarity, and the average is over the background fields generated by the
plus-mover hadron $A$ or the minus-mover hadron $B$. The twist two TMDs are
given by the first term, and the second term contains higher twist
corrections. In the present case we focus on the corrections enhanced by
$A^{1/3}$. In the Glauber model, these would take the following symbolic form

$\langle F^{+\perp}F^{+\perp}F^{+\perp}F^{+\perp}\rangle_{A}\langle
F^{-\perp}F^{-\perp}\rangle_{B}+\langle F^{+\perp}F^{+\perp}\rangle_{A}\langle
F^{-\perp}F^{-\perp}F^{-\perp}F^{-\perp}\rangle_{B}\,,$ (74)

which is consistent with Collins:2007nk ; Qiu:2001hj . The higher twist
contributions in eq. (73) do not appear in our case. This simply follows from
the fact that each nucleon is a color singlet and the nucleons are
uncorrelated in the nucleus. As a final remark about the Glauber-gluon
cancellation usually observed at leading twist, the contribution enhanced by
$A^{1/3}$ is a higher twist effect where the cancellation is only partial,
which can be seen in, e.g., eq. (IV.1.1) where the field average would not
vanish between real and virtual contributions.

##  V Cold nuclear effects on jets

In this section, we focus on the azimuthal decorrelation in boson-jet
production due to the presence of QCD matter. The interactions with QCD medium
before the hard collision is identical to the Drell-Yan process as studied in
the previous section. Below, we evaluate the corrections after the hard
collision and the interference between the interactions before and after the
hard collision.

###  V.1 The coordinates for jet production

For a jet collinear to $n_{3}^{\mu}=(1,\sin\theta_{3},0,\cos\theta_{3})$, it
is convenient to choose the four basis vectors according to eq. (37) with
$n_{i}=n_{3}$ and $n_{j}=n_{1}$ or $n_{2}$. When the jet interacts with right-
moving nucleons, one can choose $n_{j}^{\mu}=n_{1}^{\mu}=(1,0,0,1)$, yielding

$\displaystyle
n_{\perp_{1}}^{\mu}=(n_{3}^{+}/\sin\theta_{3},1,0,n_{3}^{+}/\sin\theta_{3}),\qquad
n_{\perp_{2}}^{\mu}=(0,0,1,0).$ (75)

Accordingly, the Feynman propagator in eq. (49) reads explicitly

$\displaystyle D_{F}(x)$
$\displaystyle=\int\frac{dp^{-}}{2\pi}\frac{1}{2p^{-}}e^{-\frac{i}{n_{3}^{-}}p^{-}x^{n_{3}}}\theta(x^{-})\delta^{(2)}(\mathbf{x}_{\perp}),$
(76)

indicating that, in terms of $x^{-}$ and $\mathbf{x}_{\perp}$, the jet parton
follows the trajectory of a classical particle moving in the $n_{3}$
direction. To be more specific, according to the decomposition in eq. (37),
the transverse components of $x^{\mu}$ with respect to the beam direction can
be expressed as

$\displaystyle\mathbf{x}=\frac{x^{-}\mathbf{n}_{3}+x^{n_{3}}\mathbf{n}_{1}}{n_{3}^{-}}+x^{1}\mathbf{n}_{\perp_{1}}+x^{2}\mathbf{n}_{\perp_{2}}=\bigg{(}\frac{x^{-}\sin\theta_{3}}{n_{3}^{-}},0\bigg{)}+\mathbf{x}_{\perp}.$
(77)

This expression is identical to the transverse components of the trajectory of
a classical particle moving along $n_{3}$:

$\displaystyle
x^{\mu}(t)=(0,x^{1},x^{2},0)+n_{3}^{\mu}t=\bigg{(}t,\frac{x^{-}(t)\sin\theta_{3}}{n_{3}^{-}}+x^{1},x^{2},t\cos\theta_{3}\bigg{)}$
(78)

with $x^{-}(t)=t-v_{z}t=n_{3}^{-}t$.

Similarly, for the interaction between the jet and left-moving nucleons, we
choose $n_{j}^{\mu}=n_{2}^{\mu}=(1,0,0,-1)$. Accordingly, the transverse basis
vectors in eq. (36) take the form

$\displaystyle
n_{\perp_{1}}^{\mu}=(n_{3}^{-}/\sin\theta_{3},1,0,-n_{3}^{-}/\sin\theta_{3}),\qquad
n_{\perp_{2}}^{\mu}=(0,0,1,0),$ (79)

and

$\displaystyle D_{F}(x)$
$\displaystyle=\int\frac{dp^{+}}{2\pi}\frac{1}{2p^{+}}e^{-\frac{i}{n_{3}^{+}}p^{+}x^{n_{3}}}\theta(x^{+})\delta^{(2)}(\mathbf{x}_{\perp}).$
(80)

This indicates that, in terms of $x^{+}$ and $\mathbf{x}_{\perp}$, the
transverse location of the jet parton with respect to the beam directions can
be expressed as

$\displaystyle\mathbf{x}=\frac{x^{+}\mathbf{n}_{3}+x^{n_{3}}\mathbf{n}_{2}}{n_{3}^{+}}+x^{1}\mathbf{n}_{\perp_{1}}+x^{2}\mathbf{n}_{\perp_{2}}=\bigg{(}\frac{x^{+}\sin\theta_{3}}{n_{3}^{+}},0\bigg{)}+\mathbf{x}_{\perp},$
(81)

which also coincides with the transverse trajectory of a classical particle
moving along the $n_{3}$ direction.

###  V.2 Single scattering of a jet in QCD medium

####  V.2.1 Single scattering after the hard collision

Figure 4: Cold nuclear effects after the hard collisions. The left panel show
the real diagram for single scattering in which the outgoing parton exchanges
one gluon with one encounter nucleon in both the amplitude and the conjugate
amplitude. The right panel shows the corresponding virtual diagrams in which
the outgoing parton does not interact with the nucleon either in the amplitude
or the conjugate amplitude.

The outgoing jet can interact with nucleons in either of the nuclei. Let us
start with single scattering between the outgoing parton with one additonal
nucleon from nucleus 2, as shown in fig. 4.

The field correlator is expressed in the same way as that before the hard
collision, as given by eq. (53). Here, one only needs to evaluate it in the
coordinate system defined by $n_{2}$, $n_{3}$, and the transverse basis
vectors given in eq. (79). After inserting $D_{F}$ in eq. (49) and
$G_{F}^{\mu\nu}$ in eq. (50) into eq. (53), we find

$\displaystyle\langle
A^{b\mu}(x^{+},\mathbf{x}_{\perp})A^{b\mu^{\prime}}(x^{\prime+},\mathbf{x}^{\prime}_{\perp})\rangle=$
$\displaystyle
g^{2}C_{F}n_{2}^{\mu^{\prime}}n_{2}^{\mu}\delta(x^{+}-{x^{\prime}}^{+})\int
d^{2}\mathbf{y}_{\perp}\rho_{A_{2}}(x^{+},\mathbf{y}-\mathbf{b})$
$\displaystyle\times$ $\displaystyle\int d\xi
f_{q}(\xi)F_{1}(|\mathbf{x}_{\perp}-\mathbf{y}_{\perp}|)F_{1}(|\mathbf{x}^{\prime}_{\perp}-\mathbf{y}_{\perp}|),$
(82)

where the measure $d^{2}\mathbf{y}_{\perp}=dy^{1}dy^{2}$, and $\mathbf{y}$,
the transverse vector with respect to the beams in the nucleon density, is
given by

$\displaystyle\mathbf{y}=\bigg{(}\frac{\sin\theta_{3}}{n_{3}^{+}}x^{+}+y^{1},y^{2}\bigg{)}=\bigg{(}\frac{\sin\theta_{3}}{n_{3}^{+}}x^{+},0\bigg{)}+\mathbf{y}_{\perp}$
(83)

according to eq. (81). Note, for example, in terms of
$(x^{+},\mathbf{x}_{\perp})$, one has
$\mathbf{x}=(\sin\theta_{3}x^{+}/{n_{3}^{+}}+x^{1},x^{2})$ and, consequently,
$|\mathbf{x}_{\perp}-\mathbf{y}_{\perp}|=|\mathbf{x}-\mathbf{y}|$. This shows
that the above basis vectors provide a convenient description of the
interaction between the outgoing parton and nucleons close to its trajectory.

Following the evaluation of fig. 3 (a) with the minus components of
coordinates and momenta replaced with their corresponding $n_{3}$-components,
one has

$\displaystyle\text{fig.~{}\ref{fig:LO_after}~{}(a)}=$
$\displaystyle\frac{g^{2}C_{k}}{N_{c}^{2}-1}\int
d^{4}X\int^{\infty}_{X^{+}}dX_{1}^{+}\rho_{A_{1}}(X^{-},\mathbf{X})\rho_{A_{2}}(X^{+},\mathbf{X}-\mathbf{b})d\sigma^{(0)}_{ij\to
kl}(X_{1}^{+}-X^{+})$ $\displaystyle\times\int
dtd^{2}\mathbf{x}_{\perp}e^{-i\mathbf{x}_{\perp}\cdot\mathbf{Q}_{\perp}}\langle
A^{a-}(X_{1}^{+}+t,\mathbf{X}_{\perp}+{\mathbf{x}_{\perp}}/{2})A^{a-}(X_{1}^{+}-t,\mathbf{X}_{\perp}-{\mathbf{x}_{\perp}}/{2})\rangle,$
(84)

where $\mathbf{X}_{\perp}$ is given by

$\displaystyle\mathbf{X}_{\perp}=(-n_{\perp_{1}}\cdot X,-n_{\perp_{2}}\cdot
X)=\mathbf{X}-(\sin\theta_{3}X^{+}/n_{3}^{+},0),$ (85)

and

$\displaystyle\mathbf{Q}_{\perp}=(-n_{\perp_{1}}\cdot(p_{3}+p_{4}-p_{1}),Q_{T}^{y})=(p_{3}^{1},Q_{T}^{y}).$
(86)

Here, we have used the relation

$\displaystyle
p_{3}^{\mu}=\frac{p_{T}}{\sin\theta_{3}}n_{3}^{\mu}+\frac{\sin\theta_{3}}{2p_{T}}|\mathbf{p}_{3\perp}|^{2}n_{2}^{\mu}+p_{3\perp}^{\mu},\qquad
p_{4}^{\mu}=p_{T}(\cosh\eta_{4},-1,0,\sinh\eta_{4}),$ (87)

and, to obtain $d\sigma^{(0)}_{ij\to kl}$,

$\displaystyle\frac{2\sin[\frac{2}{n_{3}^{+}}(X_{1}^{+}-X^{+})n_{3}\cdot(p_{1}+p_{2}-p_{4})]}{\frac{2}{n_{3}^{+}}n_{3}\cdot(p_{1}+p_{2}-p_{4})}=\frac{2\sin[(X_{1}^{+}-X^{+})(p_{1}^{-}+p_{2}^{-}-p_{4}^{-})]}{p_{1}^{-}+p_{2}^{-}-p_{4}^{-}},$
(88)

where $p_{1}^{+}=p_{3}^{+}+p_{4}^{+}$ and terms
$\propto|\mathbf{p}_{3\perp}|^{2}/p_{T}$ are neglected.

Then, including the virtual diagrams as exemplified by fig. 4 (b) gives

$\displaystyle\text{fig.~{}\ref{fig:LO_after}}=$
$\displaystyle\frac{-g^{2}C_{k}}{2(N_{c}^{2}-1)}\int
d^{4}X\int^{\infty}_{X^{+}}dX_{1}^{+}\rho_{A_{1}}(X^{-},\mathbf{X})\rho_{A_{2}}(X^{+},\mathbf{X}-\mathbf{b})d\sigma^{(0)}_{ij\to
kl}(X_{1}^{+}-X^{+})\int d^{2}\mathbf{x}_{\perp}$ $\displaystyle\times\int
dte^{-i\mathbf{x}_{\perp}\cdot\mathbf{Q}_{\perp}}\langle[A^{a-}(X_{1}^{+}+t,\mathbf{X}_{\perp}+{\mathbf{x}_{\perp}}/{2})-A^{a-}(X_{1}^{+}-t,\mathbf{X}_{\perp}-{\mathbf{x}_{\perp}}/{2})]^{2}\rangle.$
(89)

The expectation value of the gluon field operator

$\displaystyle G$
$\displaystyle(X_{1}^{+},\mathbf{X}_{\perp},\mathbf{x}_{\perp})\equiv\int
dt\langle[A^{a-}(X_{1}^{+}+t,\mathbf{X}_{\perp}+{\mathbf{x}_{\perp}}/{2})-A^{a-}(X_{1}^{+}-t,\mathbf{X}_{\perp}-{\mathbf{x}}_{\perp}/{2})]^{2}\rangle$
(90)

also contains a logarithm in a non-perturbative scale, which needs to be
renormalized. We then change variables from $\mathbf{y}_{\perp}$ in eq.
(V.2.1) to $\mathbf{r}_{\perp}=\mathbf{X}_{\perp}-\mathbf{y}_{\perp}$, and
keep only the logarithmic term in $|\mathbf{x}_{\perp}|$ after integrating
over $\mathbf{r}_{\perp}$ to obtain

$\displaystyle G$
$\displaystyle(X_{1}^{+},\mathbf{X}_{\perp},\mathbf{x}_{\perp})\approx\pi|\mathbf{x}_{\perp}|^{2}\rho_{A_{2}}(X_{1}^{+},\mathbf{X}((X_{1}^{+}-X^{+})/n_{3}^{+}))xG(0,1/|\mathbf{x}_{\perp}|^{2}),$
(91)

where the transverse location of the jet at $t$ is defined as

$\displaystyle\mathbf{X}(t)\equiv\mathbf{X}+{(\sin\theta_{3},0)}t,$ (92)

which is the same as $\mathbf{y}$ in eq. (83) upon expansion at small
$\mathbf{r}_{\perp}$:

$\displaystyle\mathbf{y}=\bigg{(}\frac{\sin\theta_{3}}{n_{3}^{+}}X_{1}^{+},0\bigg{)}+\mathbf{y}_{\perp}\approx\bigg{(}\frac{\sin\theta_{3}}{n_{3}^{+}}X_{1}^{+},0\bigg{)}+\mathbf{X}_{\perp}=\bigg{(}\frac{\sin\theta_{3}}{n_{3}^{+}}(X_{1}^{+}-X^{+}),0\bigg{)}+\mathbf{X}.$
(93)

Accordingly, one has

$\displaystyle\text{fig.~{}\ref{fig:LO_after}}=$
$\displaystyle-\frac{1}{2}\int
d^{4}X\int^{\infty}_{X^{+}}dX_{1}^{+}\rho_{A_{1}}(X^{-},\mathbf{X})\rho_{A_{2}}(X^{+},\mathbf{X}-\mathbf{b})d\sigma^{(0)}_{ij\to
kl}(X_{1}^{+}-X^{+})$ $\displaystyle\times\int
d^{2}\mathbf{x}_{\perp}e^{-i\mathbf{x}_{\perp}\cdot\mathbf{Q}_{\perp}}|\mathbf{x}_{\perp}|^{2}\hat{q}_{k/A_{2}}(X_{1}^{+},\mathbf{X}((X_{1}^{+}-X^{+})/n_{3}^{+})-\mathbf{b},|\mathbf{x}_{\perp}|).$
(94)

####  V.2.2 Incoherence between initial- and final-state effects

Figure 5: Interference between single scattering before and after the hard
collision $ij\to kl$. The left panel show the real diagram for single
scattering in which two different partons engaging in the hard collisions
exchanges one gluon with one encounter nucleon in the amplitude and the
conjugate amplitude. The right panel shows the corresponding virtual diagrams
in which the incoming parton does not interact with the nucleon either in the
amplitude or the conjugate amplitude.

The interference diagrams between single collision before and after the hard
collision, as shown in fig. 5, are power suppressed. For virtual diagrams in
fig. 5 (b), the integration range for $x^{+}$ shrinks to zero. Consequently,
such diagrams can be neglected. For a real diagram such as fig. 5 (a), in
order to reach the same conclusion as the virtual diagrams, we also needs to
require ${x^{\prime}}^{+}>{x}_{1}^{+}$ and ${x^{\prime}}_{1}^{+}>x^{+}$. This
feature is inherent to the time order perturbation theory. The hard process is
defined as a localized event, which, in accordance with the uncertainty
principle, is given by the inverse of the hard scale. Consequently, all other
time scales involved are larger than this time scale, which allows us to
identify both times $x_{1}^{+}$ and ${x^{\prime}}_{1}^{+}$. Further
discussions can be found in Collins:1983ju  and in references therein.

####  V.2.3 Azimuthal decorrelation due to single scattering in the boson-jet
process

Gathering all the above single scattering results, which include single
scattering occurring before the hard collision, as summarized in eqs. (IV.1.3)
for the Drell-Yan process, and single scattering occurring after the hard
collision, as given in eq. (V.2.1), yields

$\displaystyle\frac{d\sigma^{(1)}_{A_{1}A_{2}\to
JV}}{d^{2}\mathbf{b}d\eta_{3}d\eta_{4}dp_{T}dQ_{T}^{y}}=$
$\displaystyle-\frac{1}{2}\int
d^{4}X\rho_{A_{1}}(X^{-},\mathbf{X})\rho_{A_{2}}(X^{+},\mathbf{X}-\mathbf{b})\int\frac{dy}{2\pi}e^{-iyQ^{y}_{T}}y^{2}$
$\displaystyle\times\sum\limits_{ijk}\bigg{[}\int_{-\infty}^{X^{+}}dX_{1}^{+}$
$\displaystyle\frac{d\sigma^{(0)}_{ij\to
kV}}{d\eta_{3}d\eta_{4}dp_{T}}(X^{+}-X_{1}^{+})\hat{q}_{i/A_{2}}(X_{1}^{+},\mathbf{X}-\mathbf{b},|y|)$
$\displaystyle+\int_{-\infty}^{X^{-}}dX_{1}^{-}$
$\displaystyle\frac{d\sigma^{(0)}_{ij\to
kV}}{d\eta_{3}d\eta_{4}dp_{T}}(X^{-}-X_{1}^{-})\hat{q}_{j/A_{1}}(X_{1}^{-},\mathbf{X},|y|)$
$\displaystyle+\int^{\infty}_{X^{+}}dX_{1}^{+}$
$\displaystyle\frac{d\sigma^{(0)}_{ij\to
kV}}{d\eta_{3}d\eta_{4}dp_{T}}(X_{1}^{+}-X^{+})\hat{q}_{k/A_{2}}(X_{1}^{+},\mathbf{X}((X_{1}^{+}-X^{+})/n_{3}^{+})-\mathbf{b},|y|)$
$\displaystyle+\int^{\infty}_{X^{-}}dX_{1}^{-}$
$\displaystyle\frac{d\sigma^{(0)}_{ij\to
kV}}{d\eta_{3}d\eta_{4}dp_{T}}(X_{1}^{-}-X^{-})\hat{q}_{k/A_{1}}(X_{1}^{-},\mathbf{X}((X_{1}^{-}-X^{-})/n_{3}^{-}),|y|)\bigg{]}.$
(95)

Here, we have kept only the $y$-component of $\mathbf{Q}_{T}$ by integrating
out $Q_{T}^{x}$ for single scattering that happens before the hard collision,
and by integrating out $Q_{T}^{1}$ for single scattering that occurs after the
hard collision.

###  V.3 The boson-jet azimuthal decorrelation in large nuclei

In the large $R$ expansion, it is reasonable to make the same replacement as
eq. (70) in $d\sigma^{(0)}_{ij\to kl}$. In this case, one can sum over the
results for multiple scattering and finally obtains

$\displaystyle\frac{d\sigma_{A_{1}A_{2}\to
JV}}{d^{2}\mathbf{b}d\eta_{3}d\eta_{4}dp_{T}dQ_{T}^{y}}=\sum\limits_{ijk}\frac{d\sigma^{(0)}_{ij\to
kV}}{d\eta_{3}d\eta_{4}dp_{T}}\int
dX^{+}dX^{-}d^{2}\mathbf{X}{\rho}_{A_{1}}(X^{-},\mathbf{X})\,{\rho}_{A_{2}}(X^{+},\mathbf{X}-\mathbf{b})$
$\displaystyle~{}~{}~{}~{}\times\int\frac{dy}{2\pi}e^{-iyQ^{y}_{T}-\frac{y^{2}}{4}\int_{-\infty}^{X^{+}}dX_{1}^{+}\hat{q}_{i/A_{2}}(X_{1}^{+},\mathbf{X}-\mathbf{b},|y|)-\frac{y^{2}}{4}\int_{-\infty}^{X^{-}}dX_{1}^{-}\hat{q}_{j/A_{1}}(X_{1}^{-},\mathbf{X},|y|)}$
$\displaystyle~{}~{}~{}~{}\times
e^{-\frac{y^{2}}{4}\int^{\infty}_{X^{+}}dX_{1}^{+}\hat{q}_{k/A_{2}}(X_{1}^{+},\mathbf{X}((X_{1}^{+}-X^{+})/n_{3}^{+})-\mathbf{b},|y|)-\frac{y^{2}}{4}\int^{\infty}_{X^{-}}dX_{1}^{-}\hat{q}_{k/A_{1}}(X_{1}^{-},\mathbf{X}((X_{1}^{-}-X^{-})/n_{3}^{-}),|y|)}.$
(96)

It factorizes into a convolution of a hard cross section and a medium-modified
initial distributions, which is a mixture of the two incoming partons, as well
as a medium-modified jet function. More in detail:

  * •

The initial line of eq. (V.3) represents the probability of selecting a
nucleon within a given nucleus, multiplied by a hard scattering factor that
quantifies the interaction between parton $i$ and $j$ from those nucleons. The
distributions of nucleons in nuclei $A_{1}$ and $A_{2}$ are represented by the
functions $\rho_{A_{1/2}}$.

  * •

The initial-state and final-state distributions are modified by the
introduction of Wilson lines in the last two lines of eq. (V.3). Two of these
Wilson lines are associated with the active partons of the hard process and
are semi-infinite past pointing. These depend on the background field
generated by the other nucleus. In the third line, we observe that the final
state (future-pointing) Wilson lines depend on the background field generated
by both nuclei. We note the relation between the Gaussian expectation value of
correlators of Wilson lines and the form involving $\hat{q}$, as presented in
ref. Baier:1996sk , where $\hat{q}_{k/A_{n}}$ is the jet quenching parameter
of parton $k$ in the medium $A_{n}$ defined in eq. (69) and is related to the
gluon distribution $xG(x,1/|\bf{x}|^{2})$.

###  V.4 Discussion about factorization

The following discussion will present the key points of this analysis. First,
in order to demonstrate the factorization into initial and final state
distributions, we have shown that interference diagrams in which one
background gluon connects to an initial state parton in the amplitude and one
background gluon connects to a final state parton in the conjugate amplitude
are power-suppressed, alike in Collins:1983ju . This relies on two key
features: (a) The hard process in coordinate space is localized to be at the
same coordinate in both the amplitude and the conjugate amplitude. This is
discussed in more detail in the section surrounding eq. (IV.1.1). (b) The
interaction with the background is mediated by a Glauber gluon, which occurs
instantaneously. This feature of the Glauber gluon can be observed in eq.
(50). Consequently, the possibility of interference diagrams is eliminated by
(a) and (b), and the corresponding contributions to the observable are
suppressed. It is anticipated that corrections to this result will be
suppressed by large ratios of the following strongly ordered scales:
$\Lambda_{QCD}\ll Q_{T}\ll Q$. The observable can be described as an
incoherent sum of initial and final state interactions.

Second, a review of the interaction between the active and spectator partons
is now in order. In both the Drell-Yan process and the photon-jet process, the
initial or final distribution involves Wilson lines, which are process
dependent. This process dependence implies the breakdown of collinear
factorization, a concept often discussed in the TMD literature. In this case,
such factorization is labeled as generalized factorization (see, for example,
the references Collins:2007nk ; Rogers:2010dm ; Mulders:2011zt ). For a given
parton from a nucleon, the Wilson line encodes the scattering with a
background field generated by the partons belonging to the other nucleus. This
conclusion is also reached for the Wilson lines in the final state, where the
background is generated by partons from both nuclei. This is in contrast to
the case of Wilson lines that sum collinear gluons belonging to the same
nucleus.

Finally, in order to achieve a comprehensive understanding, we must consider
the influence of spectator-spectator interactions as it is done in the context
of collinear factorization. In collinear factorization, a significant amount
of effort is dedicated to demonstrating how loop momenta involved in
spectator-spectator topologies can be deformed out of the Glauber region. For
a plus moving nucleon, there is a pinch singularity which originates from the
spectator propagator pole and the active parton propagator pole in the minus
momentum complex plane. In order to accommodate the deformation of the
integration contour, it is necessary to remove the spectator’s pole. This can
be accomplished subsequent to the summation of all potential final-state cuts.
Following this removal, the integration contour can be deformed out of the
Glauber region. For large values of transverse momentum flowing from the
spectator to the active parton, and thus passing through the nucleon-active-
spectator vertex, the wave function of the nucleon has the effect of
suppressing the contribution. For further discussion, we direct the reader to
the work Collins:1982wa  and references therein. In the present case, where we
focus on the spectators coming from different nucleons which are those
enhanced by the number of nucleons $A$, the spectator-spectator interactions
do not involve a loop momentum since we assume the nucleons to be
uncorrelated. Nevertheless, the same line of reasoning can be applied. In
particular, in order for a spectator-spectator interaction to contribute to a
significant momentum imbalance, $Q_{T}\gg\Lambda_{QCD}$, this requires a large
momentum to flow through the nucleon-active-spectator vertex. This
contribution is once again suppressed by the nucleon wave function, and thus
is not accounted for in this analysis.

##  VI Summary and Outlook

In this work we have addressed the factorization of hard processes in nuclear
collisions, directly assumed in most studies devoted to this subject, and its
modification due to the interaction with other nucleons than those entering
the hard scattering, called cold nuclear matter effects. We have focused on
two observables: Drell-Yan and $\gamma$-jet. The nuclei have been modelled,
following the Glauber model, as a superposition of independent nucleons,
neglecting any interaction among them.

We have worked at the lowest order in perturbation theory, considering just
one gluon exchange between the active and spectator partons and not including
radiation. Then we have extended the study to multiple gluon exchanges. We
have restricted to the kinematic situation where the final total transverse
momentum of the observed Drell-Yan or $\gamma$-jet pair is much larger than
the imbalance, in turn much larger than any non-perturbative scale. We have
discarded power corrections of these scales, and considered only those
contributions enhanced by the number of nucleons.

We have arrived at a form of factorization for the cross sections where the
parton densities of the projectile and target become dressed by Wilson lines
that depend on the other nucleons in both nuclei, thus being process
dependent. The fields entering the Wilson lines can be described by the
transport coefficient of nuclear matter used in many studies of the field,
leading to broadening. On the other hand, factorization between initial state
and final state effects is preserved at leading power. This kind of
factorization is alike the generalized factorization discussed in the TMD
literature.

The obvious outlook of our work is the inclusion of radiation from active
initial and final partons and the link with the respective evolution equations
and with the phenomenon of radiative energy loss leading to jet quenching.
This further step, of considerable technical difficulty, is left for the
future. As an example, in the present work the field in the Wilson lines that
contain the cold nuclear matter effects, identical to those in DIS and pA
Kovchegov:1998bi ; Kovchegov:2013cva ; Kovchegov:2015zha , results from the
exchange of Glauber gluons with nucleons of both nuclei. It does not include
contributions from two Glauber gluons exchanged from two nucleons from the
same or different nuclei that after merging join a line entering the hard
scattering (see, e.g., Kovchegov:2005ss ; Wu:2017rry  for the resulting
background gluon field). That is, in the background field that couples to the
collinear partons participating in the hard collision (cf. eq. (32)):

$A^{\mu}=A^{\mu}_{1,cl}+A^{\mu}_{2,cl}+A^{\mu}_{11/22/12,CGC}+\cdots,$ (97)

we restrict ourselves to the first two terms on the r.h.s, neglecting
corrections involving the collision of one pair of nucleons, denoted as
$A^{\mu}_{11/22/12,CGC}$.

Figure 6:  Sketch of Wilson lines to be included in an operator definition of
the PDF (in the direction $n_{1}$) containing cold nuclear matter effects.
Gluons in blue are collinear to the $n_{1}$ direction, and gluons in red are
Glauber gluons produced by nucleons moving in the $n_{2}$ direction. The two
bullets indicate the location of the quark field insertions.

One additional point to be clarified is the operation definition of the new
objects, modified TMD PDFs and jet functions, that we find. A tentative idea
would be the inclusion of additional Wilson lines. For example, for DY, on top
of the standard Wilson lines resumming collinear exchanges point to $-\infty$
in the direction of one nucleus, there would be Wilson lines resumming Glauber
exchanges in the direction of the other nucleus, with longitudinal extent
proportional to its length, see fig. 6. The corresponding average should then
be done on the wave functions of both nuclei.

## Acknowledgments

We thank Jianwei Qiu, Carlos Salgado, Matt Sievert and, specially, Yuri
Kovchegov, for most useful comments. This work is part of the project
CEX2023-001318-M financed by MCIN/AEI/10.13039/501100011033. We also
acknowledge financial support from Xunta de Galicia (Centro singular de
investigación de Galicia accreditation 2019-2022, ref. ED421G-2019/05), by
European Union ERDF, and by the Spanish Research State Agency under project
PID2020-119632GBI00. B.W. acknowledges the support of the Ramón y Cajal
program with the Grant No. RYC2021-032271-I and the support of Xunta de
Galicia under the ED431F 2023/10 project. This work has been performed in the
framework of the European Research Council project ERC-2018-ADG-835105
YoctoLHC and the MSCA RISE 823947 ”Heavy ion collisions: collectivity and
precision in saturation physics” (HIEIC), and has received funding from the
European Union’s Horizon 2020 research and innovation programme under grant
agreement No. 824093.

## References

  * (1) J.C. Collins, D.E. Soper and G.F. Sterman, _Factorization of Hard Processes in QCD_ , _Adv. Ser. Direct. High Energy Phys._ 5 (1989) 1 [hep-ph/0409313].
  * (2) G.F. Sterman, _Partons, factorization and resummation, TASI 95_ , in _Theoretical Advanced Study Institute in Elementary Particle Physics (TASI 95): QCD and Beyond_ , pp. 327–408, 6, 1995 [hep-ph/9606312].
  * (3) I.W. Stewart, _Lectures on the Soft-Collinear Effective Theory_ , in _MIT EFT Course 8.851_ , 2003.
  * (4) T. Becher, A. Broggio and A. Ferroglia, _Introduction to Soft-Collinear Effective Theory_ , vol. 896, Springer (2015), 10.1007/978-3-319-14848-9, [1410.1892].
  * (5) J. Collins and J.-W. Qiu, _$k_{T}$ factorization is violated in production of high-transverse-momentum particles in hadron-hadron collisions_, _Phys. Rev. D_ 75 (2007) 114014 [0705.2141].
  * (6) T.C. Rogers and P.J. Mulders, _No Generalized TMD-Factorization in Hadro-Production of High Transverse Momentum Hadrons_ , _Phys. Rev. D_ 81 (2010) 094006 [1001.2977].
  * (7) P.J. Mulders and T.C. Rogers, _Gauge Links, TMD-Factorization, and TMD-Factorization Breaking_ , 1102.4569.
  * (8) J.-w. Qiu and G.F. Sterman, _QCD and rescattering in nuclear targets_ , _Int. J. Mod. Phys. E_ 12 (2003) 149 [hep-ph/0111002].
  * (9) S. Cao and X.-N. Wang, _Jet quenching and medium response in high-energy heavy-ion collisions: a review_ , _Rept. Prog. Phys._ 84 (2021) 024301 [2002.04028].
  * (10) Y.V. Kovchegov and E. Levin, _Quantum chromodynamics at high energy_ , vol. 33, Cambridge University Press (8, 2012), 10.1017/CBO9781139022187.
  * (11) L.D. McLerran and R. Venugopalan, _Computing quark and gluon distribution functions for very large nuclei_ , _Phys. Rev. D_ 49 (1994) 2233 [hep-ph/9309289].
  * (12) L.D. McLerran and R. Venugopalan, _Gluon distribution functions for very large nuclei at small transverse momentum_ , _Phys. Rev. D_ 49 (1994) 3352 [hep-ph/9311205].
  * (13) A.H. Mueller, _Small x Behavior and Parton Saturation: A QCD Model_ , _Nucl. Phys. B_ 335 (1990) 115.
  * (14) Y.V. Kovchegov, _NonAbelian Weizsacker-Williams field and a two-dimensional effective color charge density for a very large nucleus_ , _Phys. Rev. D_ 54 (1996) 5463 [hep-ph/9605446].
  * (15) Y.V. Kovchegov and A.H. Mueller, _Gluon production in current nucleus and nucleon - nucleus collisions in a quasiclassical approximation_ , _Nucl. Phys. B_ 529 (1998) 451 [hep-ph/9802440].
  * (16) M.L. Miller, K. Reygers, S.J. Sanders and P. Steinberg, _Glauber modeling in high energy nuclear collisions_ , _Ann. Rev. Nucl. Part. Sci._ 57 (2007) 205 [nucl-ex/0701025].
  * (17) Y.V. Kovchegov, _Can thermalization in heavy ion collisions be described by QCD diagrams?_ , _Nucl. Phys. A_ 762 (2005) 298 [hep-ph/0503038].
  * (18) B. Wu and Y.V. Kovchegov, _Time-dependent observables in heavy ion collisions. Part I. Setting up the formalism_ , _JHEP_ 03 (2018) 158 [1709.02866].
  * (19) Y.V. Kovchegov and B. Wu, _Time-dependent observables in heavy ion collisions. Part II. In search of pressure isotropization in the $\varphi^{4}$ theory_, _JHEP_ 03 (2018) 157 [1709.02868].
  * (20) M. Li, T. Lappi and X. Zhao, _Scattering and gluon emission in a color field: A light-front Hamiltonian approach_ , _Phys. Rev. D_ 104 (2021) 056014 [2107.02225].
  * (21) D. Avramescu, V. Băran, V. Greco, A. Ipp, D.I. Müller and M. Ruggieri, _Simulating jets and heavy quarks in the glasma using the colored particle-in-cell method_ , _Phys. Rev. D_ 107 (2023) 114021 [2303.05599].
  * (22) B. Wu, _Factorization and transverse phase-space parton distributions_ , _JHEP_ 07 (2021) 002 [2102.12916].
  * (23) A.H. Mueller, _Sudakov Effects and Initial State Interactions in the Drell-Yan Process_ , _Phys. Lett. B_ 108 (1982) 355.
  * (24) J.C. Collins, D.E. Soper and G.F. Sterman, _Transverse Momentum Distribution in Drell-Yan Pair and W and Z Boson Production_ , _Nucl. Phys. B_ 250 (1985) 199.
  * (25) T. Becher and M. Neubert, _Drell-Yan Production at Small $q_{T}$, Transverse Parton Distributions and the Collinear Anomaly_, _Eur. Phys. J. C_ 71 (2011) 1665 [1007.4005].
  * (26) M.G.A. Buffing, Z.-B. Kang, K. Lee and X. Liu, _A transverse momentum dependent framework for back-to-back photon+jet production_ , 1812.07549.
  * (27) P. Sun, B. Yan, C.P. Yuan and F. Yuan, _Resummation of High Order Corrections in $Z$ Boson Plus Jet Production at the LHC_, 1810.03804.
  * (28) Y.-T. Chien, D.Y. Shao and B. Wu, _Resummation of Boson-Jet Correlation at Hadron Colliders_ , _JHEP_ 11 (2019) 025 [1905.01335].
  * (29) Y.-T. Chien, R. Rahn, S. Schrijnder van Velzen, D.Y. Shao, W.J. Waalewijn and B. Wu, _Recoil-free azimuthal angle for precision boson-jet correlation_ , _Phys. Lett. B_ 815 (2021) 136124 [2005.12279].
  * (30) Y.-T. Chien, R. Rahn, D.Y. Shao, W.J. Waalewijn and B. Wu, _Precision boson-jet azimuthal decorrelation at hadron colliders_ , _JHEP_ 02 (2023) 256 [2205.05104].
  * (31) A.H. Mueller, B. Wu, B.-W. Xiao and F. Yuan, _Probing Transverse Momentum Broadening in Heavy Ion Collisions_ , _Phys. Lett. B_ 763 (2016) 208 [1604.04250].
  * (32) A.H. Mueller, B. Wu, B.-W. Xiao and F. Yuan, _Medium Induced Transverse Momentum Broadening in Hard Processes_ , _Phys. Rev. D_ 95 (2017) 034007 [1608.07339].
  * (33) L. Chen, G.-Y. Qin, S.-Y. Wei, B.-W. Xiao and H.-Z. Zhang, _Probing Transverse Momentum Broadening via Dihadron and Hadron-jet Angular Correlations in Relativistic Heavy-ion Collisions_ , _Phys. Lett. B_ 773 (2017) 672 [1607.01932].
  * (34) L. Chen, G.-Y. Qin, S.-Y. Wei, B.-W. Xiao and H.-Z. Zhang, _Dijet Asymmetry in the Resummation Improved Perturbative QCD Approach_ , _Phys. Lett. B_ 782 (2018) 773 [1612.04202].
  * (35) L. Chen, G.-Y. Qin, L. Wang, S.-Y. Wei, B.-W. Xiao, H.-Z. Zhang et al., _Study of Isolated-photon and Jet Momentum Imbalance in $pp$ and $PbPb$ collisions_, _Nucl. Phys. B_ 933 (2018) 306 [1803.10533].
  * (36) F. Arleo and C.-J. Naïm, _Nuclear p ⊥-broadening of Drell-Yan and quarkonium production from SPS to LHC_, _JHEP_ 07 (2020) 220 [2004.07188].
  * (37) S. Benić, O. Garcia-Montero and A. Perkov, _Isolated photon-hadron production in high energy pp and pA collisions at RHIC and LHC_ , _Phys. Rev. D_ 105 (2022) 114052 [2203.01685].
  * (38) S. Wang, J.-W. Kang, W. Dai, B.-W. Zhang and E. Wang, _Probing the in-medium $P_{T}$-broadening by $\gamma$+HF angular de-correlations_, 2107.12000.
  * (39) Y.V. Kovchegov and M.D. Sievert, _Sivers function in the quasiclassical approximation_ , _Phys. Rev. D_ 89 (2014) 054035 [1310.5028].
  * (40) J.C. Collins, D.E. Soper and G.F. Sterman, _Factorization for Short Distance Hadron - Hadron Scattering_ , _Nucl. Phys. B_ 261 (1985) 104.
  * (41) J.C. Collins and D.E. Soper, _Parton Distribution and Decay Functions_ , _Nucl. Phys. B_ 194 (1982) 445.
  * (42) J. Casalderrey-Solana and C.A. Salgado, _Introductory lectures on jet quenching in heavy ion collisions_ , _Acta Phys. Polon. B_ 38 (2007) 3731 [0712.3443].
  * (43) G.T. Bodwin, S.J. Brodsky and G.P. Lepage, _Initial State Interactions and the Drell-Yan Process_ , _Phys. Rev. Lett._ 47 (1981) 1799.
  * (44) G. Beuf, _Dipole factorization for DIS at NLO: Loop correction to the $\gamma^{*}_{T,L}\to q\overline{q}$ light-front wave functions_, _Phys. Rev. D_ 94 (2016) 054016 [1606.00777].
  * (45) J.B. Kogut and D.E. Soper, _Quantum Electrodynamics in the Infinite Momentum Frame_ , _Phys. Rev. D_ 1 (1970) 2901.
  * (46) M. Li, W. Qian, B. Wu and H. Zhang, _Anisotropic flow and the valence quark skeleton of hadrons_ , _JHEP_ 08 (2023) 144 [2304.06557].
  * (47) R. Baier, Y.L. Dokshitzer, A.H. Mueller, S. Peigne and D. Schiff, _Radiative energy loss and p(T) broadening of high-energy partons in nuclei_ , _Nucl. Phys. B_ 484 (1997) 265 [hep-ph/9608322].
  * (48) I. Balitsky and A. Tarasov, _Power corrections to TMD factorization for Z-boson production_ , _JHEP_ 05 (2018) 150 [1712.09389].
  * (49) I. Balitsky, _Rapidity-only TMD factorization at one loop_ , _JHEP_ 03 (2023) 029 [2301.01717].
  * (50) I. Balitsky, _$1/Q^{2}$ power corrections to TMD factorization for Drell-Yan hadronic tensor_, 2404.15116.
  * (51) J.C. Collins, D.E. Soper and G.F. Sterman, _All Order Factorization for Drell-Yan Cross-sections_ , _Phys. Lett. B_ 134 (1984) 263.
  * (52) J.C. Collins, D.E. Soper and G.F. Sterman, _Factorization for One Loop Corrections in the Drell-Yan Process_ , _Nucl. Phys. B_ 223 (1983) 381.
  * (53) Y.V. Kovchegov and M.D. Sievert, _Calculating TMDs of a Large Nucleus: Quasi-Classical Approximation and Quantum Evolution_ , _Nucl. Phys. B_ 903 (2016) 164 [1505.01176].