Factorization and transverse phase-space parton distributions

We first revisit impact-parameter dependent collisions of ultra-relativistic particles in quantum field theory. Two conditions sufficient for defining an impact-parameter dependent cross section are given, which could be violated in proton-proton collisions. By imposing these conditions, a general formula for the impact-parameter dependent cross section is derived. Then, using soft-collinear effective theory, we derive a factorization formula for the impact-parameter dependent cross section for inclusive hard processes with only colorless final-state products in hadron and nuclear collisions. It entails defining thickness beam functions, which are Fourier transforms of transverse phase-space parton distribution functions. By modelling non-perturbative modes in thickness beam functions of large nuclei in heavy-ion collisions, the factorization formula confirms the cross section in the Glauber model for hard processes. Besides, the factorization formula is verified up to one loop in perturbative QCD for the inclusive Drell-Yan process in quark-antiquark collisions at a finite impact parameter.


I. INTRODUCTION
In hadron and nuclear collisions, cross sections for hard processes are always a blend of short-distance and long-distance physics. QCD factorization allows one to systematically study incoherence between effects at various distance scales [1]. The predictive power of perturbative QCD for hadron collider physics relies on validity of factorization theorems and universality of quantities encoding long-distance physics. Factorization also provides a systematic approach to resum large logarithms of ratios between different scales to all orders in perturbation theory [2].
Long-distance physics in all the processes mentioned above manifests itself either in parton distribution functions (PDFs) [14,15], transverse-momentum-dependent (TMD) PDFs [15] or the beam functions [16] in SCET. They are all defined as matrix elements of gauge invariant operators sandwiched between some momentum eigenstate of colliding hadrons. And they do not possess any information on the spatial distribution of quarks and gluons inside the hadrons, needed for a holistic snapshot of quantum phases-space parton distributions [17].
The spatial distribution of partons in a hadron, say, a proton, could be revealed by studying impact-parameter dependent collisions, as have been extensively studied in nucleusnucleus (AA) collisions. The impact parameter in AA/heavy-ion collisions can be determined via centrality measurements by using the Glauber model [18]. In this model, the cross section for the Drell-Yan production of vector bosons reduces to that in binary nucleon-nucleon collisions, which is consistent with recent measurements at the LHC within experimental uncertainties [19][20][21][22][23]. In order to unambiguously define the impact parameter of protonproton (pp) collisions, one needs first to clarify conceptual difference between small colliding systems like pp collisions and large colliding systems like AA collisions.
The discovery of collectivity in pp collisions [24][25][26][27][28][29], however, blurs the boundary between The transverse momenta of particles A and B are known to be around some design values with an uncertainty ∆p T . Accordingly, their transverse locations x A and x B = x A + b can only be determined with an accuracy ∆x T limited by the uncertainty principle. This accounts for the displacement of their transverse positions across the cut. The product C with momentum p µ C is created at some hard scattering vertex around X, initiated by two partons respectively collinear to A and B. Different collinear partons can communicate with one another via exchange of soft gluons. The misalignment of the hard vertex in the transverse plane across the cut x, that is, its transverse coherence length, is of order 1/p C,T .
large and small colliding systems [30]. Concepts based on (classical) collision geometry and the impact parameter in heavy-ion collisions have been frequently employed to interpret collectivity in pp collisions in many theoretical discussions without scrutiny (see Ref. [31] for a recent review). Nowadays, redefining the boundary of physical concepts respectively applicable to pp, pA and AA collisions is one of the main focuses in high-energy nuclear physics [30,31], which demands a unified theoretical approach in QCD to treat all these collisions on the same footing.
The sole purpose of this paper is to lay the groundwork for a unified description of hard processes in impact-parameter dependent pp, pA and AA collisions based on QCD factorization. As depicted in Fig. 1, we restrict ourselves to a generic inclusive hard process in hadron and nuclear collisions at an impact parameter b: A + B → C + anything else (1) with A and B either hadrons or nuclei, and C some colorless final-state object, such as electroweak gauge bosons and Higgs bosons. Below, we motivate and outline how we proceed to derive the factorized cross section for such a hard process.
In Sec. II, we revisit the fuzzy quantum picture of impact-parameter dependent collisions.
The motivation for this section is two-fold: First, in a typical pp collision at the LHC, the proton size R, the impact parameter b and even the proton transverse spatial dispersion ∆x T 1 could be comparable with one another. In this case, a quantum description of the collision becomes more appropriate. Second, in heavy-ion collisions, the Glauber model [18] has been broadly employed not only in theoretical calculations but also in experimental studies of collision geometry from centrality measurements. Behind the Glauber model underlies a classical picture of the impact parameter, which can even trace back to Rutherford's seminal discovery in Ref. [33] that helped usher in the quantum age. Giving this model a justification in QCD entails a quantum picture instead. However, the formula for the impact-parameter dependent cross section has not been derived in quantum scattering theory [34,35]. This section is aimed to fill this gap.
We consider a collision of two ultra-relativistic particles A and B at the impact parameter b. In collider physics, all we know is that their momenta are determined to be around some design values P µ i with i = A and B. From this information, we do know that the states of the particles can be described by some wave packets in momentum space, with their momentum width constrained by experimental uncertainties. Accordingly, the particles are localized in position space, described by some well-defined spatial wave packets. We aim to define the impact-parameter dependent cross section for the collision, which is intrinsic to the colliding particles and should be independent of the beam particles' wave packets.
We find that such an impact-parameter dependent cross section can be defined, as given in Eq. (27), if the following two conditions are fulfilled. Condition i) is the high-energy limit: |P iz | P iT , ∆p T , ∆p z with ∆p T and ∆p z respectively the particle's transverse and longitudinal momentum dispersions. It is needed to identify the longitudinal momenta of 1 At the LHC, the crossing angle at the interaction point θ c can be measured with an accuracy ∆θ c < 10 µrad [32]. For a E = 7 TeV proton, the uncertainty principle dictates ∆x T ≥ 1 2E∆θc > 1.4 fm.
the wave packets in the amplitude and the conjugate amplitude. Condition ii) is to require that the fuzziness in particles' transverse positions should be smaller than b: b ∆x T . It is needed to identify the transverse momenta of the wave packets in the amplitude and the conjugate amplitude. Therefore, they are both necessary for integrating out the wave packets in order to define the cross section and for unambiguously defining the impact parameter.
Otherwise, the probability for producing any final states in the collision explicitly depends on the wave packets, which are not guaranteed to be the same in different experiments. And, one ends up measuring a quantity which differs across experiments.
With the general formula for the impact-parameter dependent cross section derived in Sec. II, one can calculate it for the hard process in Eq. (1) using perturbative QCD. In general, one needs to deal with the long-distance behavior in perturbative series. We appeal to QCD factorization to factor out long-distance physics in Sec. III.
In Sec. III A, using SCET we first give a detailed derivation of the factorized cross section for the process in Eq. (1). The main features of this factorization formula can be understood based on the following heuristic argument: Imagine that besides the hard scale Q we measure the transverse momentum p C,T Λ QCD of the product C. As illustrated in Fig. 1, C is produced at some hard scattering vertex around X with a transverse coherent length |x| ∼ 1/p C,T 1 fm within a time scale t h ∼ 1/Q. It effectively picks out two partons located around X respectively from A and B with a transverse spatial accuracy of the order of |x|. The distributions of these partons in the colliding particles are described by a new type of PDFs, which are referred to as thickness beam functions in this paper.
The thickness beam functions are expected to be universal as a consequence of the hard-collinear and soft-collinear factorization. Generic hard processes all involve radiation collinear to the two beam directions n µ A and n µ B . The coherent time of a n i -collinear parton with momentum p µ n i is given by where the expansion parameter λ, as a ratio between a soft momentum scale determined by some specific observable to be measured and the hard scale Q, is assumed to be much smaller than unity. That is, the beam collinear radiation takes place too early to interfere with the hard scattering. Different beam collinear partons may communicate, via exchange of soft gluons, with each other or with final-state colored objects since soft gluons have a coherent time t s t h : As long as soft gluons cannot resolve the substructure of collinear splittings, there is factorization between soft and beam collinear partons. And the thickness beam functions should be universal to all the hard processes in which these two types of factorization hold true.
Because their formation time is much longer than t h , the hard scattering occurs too rapidly to radiate soft gluons. That is, soft gluons are detached from the hard scattering vertex. As a consequence, for the process in Eq. (1), soft radiation only couples to the total color charges of collinear partons from A or B and organizes itself into soft Wilson lines.
They act on the vacuum state to define another module of the factorization formula: the soft function. Since it is the vacuum state that is involved in its definition, the soft function is independent of X and does not carry any information on the impact parameter.
Based on the above argument one can expect that the factorization formula, as given in Eq. (44), schematically takes the following form: with the impact parameter b = x B − x A , X denoting the transverse position of the hard scattering and x i being the (average) transverse location of particle i, as shown in Fig.   1. Here, the hard function H j,k→C can be calculated from the partonic process j, k → C + anything else and S is the soft function. They are the same as for conventional cross sections in pp collisions [10,12,13]. The thickness beam functions T j/i are related to the quantum transverse phase-phase PDFs via the Fourier transform. That is, they encode the information on the distribution of parton j carrying a momentum fraction z in both transverse momentum and position spaces inside particle i. A detailed discussion about the properties of T j/i is presented in Sec. III B.
In Sec. III C, we study the connection between the factorization formula in Eq. (4) and the cross section in the Glauber model in heavy-ion collisions. We find that T j/i reduces to a product of the thickness function in the Glauber model and the corresponding beam function, as the Fourier transform of the TMD PDF, when the incoming nucleus is treated as an assembly of uncorrelated nucleons. This gives the success of the Glauber model [36] a QCD justification for such hard processes. On the other hand, the factorization formula allows one to explore refined details about the parton distributions in heavy nuclei. It can be used to optimize the potential of the LHC and HL-LHC, with increased accuracy, in systematically studying cold nuclear effects [37][38][39], which are absent in the aforementioned modelling. Nuclear modifications of PDFs had been first revealed in deep inelastic scattering by the European Muon Collaboration [40] and have been favored by recent measurements in PbPb collisions at the LHC [41]. The factorization formula is important for investigating the LHC/HL-LHC's potential to pin down such effects since the hard and soft functions can be calculated at high accuracy in perturbative QCD.
Using the factorization formula in Eq. At this point, we have two ways to calculate the impact-parameter dependent cross section by either using the general formula in Sec. II or the factorization formula in Sec.
III. This provides a way to verify factorization order by order in perturbation theory: first, calculate the cross section in perturbative QCD using the general formula; then, expand the perturbative QCD results at small λ or, equivalently, large Q; and finally compare the leading-order results in Q to those given by the factorization formula.
In Sec. IV, a verification of the factorization formula is carried out for the inclusive Drell-Yan process qq → γ * . This process involves two scales: the photon virtuality Q and the impact parameter b, which set the expansion parameter λ = 1/(bQ). We only focus on the physically interesting case with λ 1, that is, the impact parameter being much larger than 1/Q. At leading order in λ, the factorization formula is confirmed up to one loop in perturbative QCD. We also calculate the one-loop spatial quark distribution in a fast moving quark, which carries information complementary to the corresponding TMD PDF. The probability to find a quark carrying a momentum fraction z at a transverse distance r from the original transverse location of the incoming quark is found to be inversely proportional to r 2 .
Detailed derivations and calculations backing up the above summary can be found in the ensuing sections. And the interested reader is invited to vet their details.

II. THE IMPACT-PARAMETER DEPENDENT CROSS SECTION
In this section we revisit the concepts of the impact parameter and the impact-parameter dependent cross section in quantum field theory (QFT). These two concepts are well-defined in classical physics, as exemplified by the original derivation of the Rutherford scattering formula [33]: the deflection angle of a charged particle scattering off a Coulomb potential is given by a unique function of the impact parameter. In contrast, in the textbook derivation of the conventional cross section (see, e.g., Chapter 4 of [35]), there is no unique relation between the deflection angle and the impact parameter. Yet, Rutherford's formula can be easily reproduced. Below, we redo the derivation of the formula for the cross section in QFT in order to restore its impact-parameter dependence.
Consider a collision of two ultra-relativistic particles A and B respectively from two counter-moving beams. Before the collision, the momentum of particle i with i = A or B is accelerated to be around some design value P i with an uncertainty ∆ p i . So, all we know is that its wave packet peaks about P i with a momentum width equal to or smaller than ∆ p i , which can be generically written in the form 2 where φ i is normalized, that is, the phase factor in front of φ i ( p) accounts for the spatial translation in the transverse plane and x i is the transverse position vector of particle i, to be determined later. Here and below, we denote two-dimensional vectors in the transverse plane by bold letters and threedimensional vectors by letters with an arrow overhead.
The impact parameter of the collision is defined as once x A and x B are given. However, the particles are, at best, known to locate somewhere with an uncertainty limited by the uncertainty principle. In order to quantitatively define x i , one needs the particle's spatial wave packet, which is given bỹ with the position eigenstate [42] As long as x is not measured with a resolution better than the particle's de Broglie wavelength, |φ i | 2 admits interpretation as the probability density to find particle i at x. Givenφ i , x i can be chosen to be the average transverse position vector of particle i. Obviously, there are alternative choices for x i , such as the center of mass or the transverse peak location of Below, we show that such an ambiguity in defining b ≡ |b| is negligible when one is allowed to define an impact-parameter dependent cross section.
In principle, one can predict the probability for producing any final state |{p f } in the collision at the impact parameter b according to Since we are only interested in the cross section, we can replace the S-matrix element by Then, plugging the wave packets in Eq. (5) into Eq. (10) and using one of the delta functions from the above replacement to integrate outp Az ,p Bz andp B yields where the measure n j=1 dp j 2π is denoted by for brevity, the off-diagonal cross section is defined as The impact-parameter dependent probability P b generally depends on the wave packets.
However, the beam particles' wave packets are not measured in collider physics. Moreover, they are not guaranteed to be the same in different experiments. The cross section, on the other hand, is intrinsic to the colliding particles and, therefore, should be independent of the wave packets in order to allow comparison across experiments. Below, we show that the following two conditions are sufficient for defining the impact-parameter dependent cross section from P b . Here, ∆x T , ∆p T and ∆p z are respectively the transverse spatial, transverse and longitudinal momentum dispersions of the colliding particles.
Condition i) is the high-energy limit in which both particles are moving predominantly along the beam (±z) directions. In this limit, the solutions to Eq. (14) are given bȳ for P Az > 0 and P Bz < 0. Since these terms are small 3 , we will drop them and takē Following [43], we define transverse Wigner functions In terms of W i , P b can be expressed as with large P iz . For collider phenomenology, one only needs to keep leading-order terms in ∆/P iz with ∆ = P i , ∆p T or ∆p z , which is equivalent to replacing P i by 0 and p iz by their design values P iz in the off-diagonal cross section 4 . |q| ∼ 1/b could also be much smaller than |P iz |.
However, as it will become clear in the following sections, we do not expand the amplitude squared in the off-diagonal cross section about |q|/P iz = 0 here. We only do so when it is needed in detailed calculations, as exemplified in Sec. IV. As a result, we have where the Mandelstam variable the incoming momenta are given by with masses being neglected, and for any four-vector V µ we define Note that one has V 2 T = −|V| 2 . Here, the light-like vectors are chosen to be and the transverse metric is defined as As a result, the integrand on the right-hand side of Eq. (19) depends on P i only through the transverse Wigner functions. In order to integrate out the wave packets and get unity, one needs to drop X i from its phase factor. Condition ii) is sufficient to justify such an 4 Another justification of this approximation is that modern detectors typically have limited resolution which can not resolve such a small variation of P i [35].
approximation. If one goes back to Eq. (12), this, equivalently, means that the difference between p i andp i is negligible compared to their average. Obviously, this is also the condition needed to eliminate the ambiguity in defining the impact parameter in Eq. (7) due to alternative choices of x i .
At the end, one can integrate out the transverse Wigner functions in Eq. (19) and obtain the impact-parameter dependent differential cross section for producing any observable O dσ where O({p f }) defines the observable O as a function of the final-state momenta {p f }, the incoming momenta are given in Eq. (23) and the phase-space measure in d-dimensional spacetime for a particle of mass m with y the particle's rapidity.

COLLISIONS
Based on the heuristic argument outlined in the introduction, one can expect that a new type of PDFs, which describe the parton distributions in transverse phase space, can be universally defined for inclusive hard processes in impact-parameter dependent hadron and nuclear collisions. In this section, we justify this argument using SCET [7][8][9][10][11].

A. Factorization for inclusive hard processes with colorless final states
In this subsection, we derive a factorized form of Eq. (27) for the process in Eq. (1). In hadron and nuclear collisions, there are additional length scales that one needs to consider: the sizes R i of the colliding particles. Accordingly, the impact-parameter dependent cross section for hadron and nuclear collisions is defined in the range: modes with p 2 f Λ 2 QCD that contribute to the observable O({p f }) to be measured. Since these non-perturbative modes are collinear to the colliding hadrons or nuclei, we take them as submodes of the corresponding beam collinear modes.

Basics of SCET
We first briefly review the elements of SCET that are relevant to our discussion. With the modes of q T ∼ 1/b taken as submodes of the corresponding beam collinear modes, all the relevant infrared degrees of freedom for the process under study are the same as those for the conventional cross section in Refs. [10,12,13]: where λ 1 is an expansion parameter and we have defined for any V µ in terms of a pair of light-like vectors n i andn i with n i ·n i = 2.
The S-matrix element at leading order in λ can be expressed generically in the form with {p X } standing for momenta of unmeasured infrared partons and the amplitude operator M being a convolution of relevant SCET operators and corresponding Wilson coefficients.
Let us combine everything but the collinear and soft fields in the SCET operators with the Wilson coefficients and denote their sum by the coefficient C. In this way, irrespective of the species of the product C,M can be always written in the form where the coefficient C is to be determined by the matching procedure after the species of the product C being specified, a i and α i are respectively the color and Lorentz/spinor indices, S n i is the soft Wilson line along the collinear direction n i and the collinear building blocks φ n i stand for [44] respectively for n i -collinear quarks, antiquarks or gluons with D µ n i T ≡ ∂ µ T − ig s A µ n i T and W n i the n i -collinear Wilson line.
Both the hard-collinear factorization and the soft-collinear factorization (at leading order in λ) are implemented throughM , which are independent of the initial states of the colliding particles. Soft radiation decouples from the hard scattering encoded in the coefficient C, which has been proved for certain processes based on infrared power counting in perturbative QCD [1,2]. Soft gluons can couple to collinear partons only through their unphysical polarizations n i · A s , which gives rise to the soft Wilson lines inM . Soft and collinear fields decouple in the SCET Lagrangian for both SCET I (after decoupling transformation [9]) and SCET II . The soft Wilson lines in coordinate space take the form 6 where φ n i ± respectively stand for the positive and negative energy parts of φ n i and A s ≡ A a s T a with T a being SU (N c ) generators in the corresponding color representation. The collinear Wilson lines take the form The inclusion of collinear Wilson lines in the collinear building blocks is mandated by the collinear gauge invariance [9].

Derivation of the factorization formula
With the S-matrix element given in Eq. (31), we are ready to derive the factorized cross section for the process in Eq. (1). We impose the two conditions in Eq. (15) and identify the impact-parameter dependent probability P b with the cross section at the outset: Here, for definiteness we measure the rapidity y C and the transverse momentum p C of the object C.
The following derivation is pretty much the same as that in the previous section, except that we keep X unintegrated. It tells us the transverse location of the hard scattering vertex as shown in Fig. 1. One can first integrate out X 0 and X z in Eq. (36) by using the momentum operator and then the longitudinal momenta associated with the wave packets in the conjugate amplitude. After that, one identifies the longitudinal momenta in the amplitude and the conjugate amplitude, as in Eq. (17). And, in terms of the transverse Wigner functions in Eq. (18), the cross section can be expressed as where in terms ofM , the off-diagonal cross section in Eq. (21) takes the form and the incoming momenta are given by with q µ i = p µ i −p µ i orthogonal to n i andn i . Finally, the two conditions in Eq. (15) allow one to integrate out the transverse Wigner functions and one has Since the collinear and soft modes decouple, the right-hand side of the above equation can be further written in a factorized form. The coefficient C inM is independent of the initial states of the colliding particles, which is combined into the hard function together with its complex conjugate. The soft Wilson lines only act on the vacuum state to define the soft function. As a result, the soft function is independent of X, as shown below. That is, the soft and hard functions are the same as for the conventional cross section [10,12,13].
The collinear fields act on the incoming states to give a new type of PDFs, which are referred to as thickness beam functions in this paper. The n i -collinear sector in Eq. (40) takes the form with r i ≡ X − x i and the spin of particle i being implicitly averaged over. Given some projector P α α , we define the corresponding thickness beam function T j/i as with d c i the dimension of the color representation of φ n i and j the parton species corresponding to φ n i . The most common projectors for unpolarized collisions are given by the spin/polarization average, which take the form for gluons (43) with d the spacetime dimension.

Inserting the expression ofM in Eq. (32) into Eq. (40) and making the replacement in
Eq. (42) eventually yields the factorization formula: where we have taken x A = 0 and x B = b, P µ i =n i ·P i 2 n µ i , the soft function is defined by (45) with x ± ≡ X ± x/2, and the hard function, given by the partonic process j(z A P A ) + k(z B P B ) → C(p C ) + anything else({p f }), takes the form withC given bỹ For the spin-averaged projectors in Eq. (43), the hard function with |M | 2 the square of the amplitude averaged over initial-state colors, spins or polarizations [45].

B. Thickness beam functions and transverse phase-space parton distributions
The thickness beam functions are universal. Their emergence only relies on the hardcollinear and soft-collinear factorization. Therefore, they should universally show up in all inclusive hard processes in hadron and nuclear collisions as long as the processes admit of these two types of factorization. And hard processes with colorless final states like that in Eq. (1) provide the cleanest way to measure the thickness beam functions. The discussion in the previous subsection can be easily generalized to deep inelastic scattering in electronproton and electron-ion collisions. Therefore, they could also be measured using the future Electron-Ion Collider [46] once the impact-parameter dependence of the collisions can be determined experimentally.
The thickness beam functions are related to transverse phase-space PDFs (TPS PDFs) 7 via the Fourier transform with respect to x: as for the beam functions [16] and TMD PDFs [15].
They are two-dimensional analogues of quantum phase-space distributions in the rest frame of a proton defined in Ref. [17]. This can be easily seen from the definition of T j/i in Eq. (42), which generally takes the form with Γ α α to be determined after P α α are chosen. By using the momentum operator, one can write In the first line, the right-hand side of this equation can be interpreted as the "probability" 8 of measuring a collinear field of type j at r inside particle i and producing a final state |m with total momentum p µ m . Accordingly, the corresponding TPS PDF is which, after being coarse-grained, is the probability of finding a collinear parton at r carrying a momentum fraction z and transverse momentum p.
In the aforementioned coarse-grained sense, thickness beam functions and TPS PDFs admit interpretation respectively as the beam functions [16] and TMD PDFs [15] at r since one typically has |r| |x| for hard processes. For inclusive hard processes with |x| ∼ 1/Q, x can be dropped from the thickness beam functions as a result of multipole expansion.
And T j/i (r, z, 0) can be viewed as the corresponding conventional parton distribution functions [14,15] at r, that is, the transverse spatial PDFs. In this case, the off-diagonal matrix element needed for defining T j/i (r, z, 0) is the same as generalized parton distributions (GPDs) [47] with the incoming momenta differing only in the transverse plane.
At the end, we explicitly write down the expressions for the thickness beam functions corresponding to the spin-averaged projectors in Eq. (43): for quarks, for antiquarks, and for gluons For gluons, there are also other possible projectors as combinations of g µν T , x µ T and r µ T . We will not exhaust all the possible forms of the projectors in this paper, whose relevance will depend on specific processes under consideration.

C. The Glauber model for hard processes in heavy-ion collisions
In heavy-ion collisions, one has That is, r i -dependence of thickness beam functions lies deep in the non-perturbative regime.
In principle, one could also quantitatively study such non-perturbative degrees of freedom using the effective field theory approach. We, instead, content ourselves with connecting our factorization formula in Eq. (44) to the cross section in the Glauber model [18] by modelling large nuclei.
impact parameter in Eq. (7) can be determined better than 1 fm as long as one is not aimed at high accuracy in the determination of the beam particles' transverse momenta and keeps ∆p T 1 GeV 9 . Since the quantum fuzziness is much smaller than the impact parameter, the classical concept of collision geometry used in the Glauber model is indeed a reasonable approximation to the underlying quantum picture.
Let us approximate the thickness beam functions by appealing to a commonly used model for heavy nuclei as in the Glauber model [18]. We also consider the difference between protons and neutrons [48]. Large nuclei are known to be loosely bound with the binding energy per nucleon ∆E ≈ 8 MeV in their rest frame. This means that the typical time scale for internal nucleon interactions is of order ∆t =n ·P i 2m i 1 ∆E ≈ 25γ fm in the lab frame. Such a time scale is much longer than any other time scales in the problem. Therefore, nucleons in the nuclei are to be treated as free particles with a normalized distributionρ i proportional to that of electric charges [49]. Accordingly, the probability to find a nucleon per unit transverse area around r i is given bŷ And, the thickness beam functions of nucleus i, which is made of Z i protons and N i neutrons, can be replaced by where the nuclear thickness function is defined as [18] T  9 In this paper, we are not concerned about how to determine the impact parameter of a collision experimentally but only about what is speakable and unspeakable about the impact parameter in the quantum picture.
with the binary nucleon-nucleon cross section σ nn given by the factorization formula for colliding two nucleon beams with neutron-to-proton ratios respectively equal to N i /Z i . That is, the nuclear modification factor 10 with D. Impact-parameter dependent pp collisions Collective phenomena among produced soft particles have been studied in pp collisions [24][25][26][27][28][29], which are presumptively related to collision geometry according to some models [31]. Inclusive hard processes in impact-parameter dependent pp collisions are worth being explored experimentally as well. Based on the factorization formula in Eq. (44), such hard processes can be potentially used to measure transverse phase-space parton distributions inside protons.
Caution is, however, needed when one studies impact-parameter dependent pp collisions.
For hard processes, one has Therefore, like heavy-ion collisions, T j/i (r, z, x) can be viewed as the distribution of parton j with a transverse size ∼ |x|, i.e., the beam function B j/i (z, x), located at r inside proton i. On the other hand, in contrast to heavy-ion collisions, the two conditions in Eq. (15) are not always fulfilled in pp collisions. As discussed in Sec. II, in order to measure universal quantities across experiments, one needs to maintain ∆p T , however, can not be too large on modern colliders. For example, the crossing angle at the interaction point θ C ∼ 100 µrad at the LHC [32], which is one of the crucial parameters 10 In experiments, one measures R AA with both its numerator and denominator averaged over a range of impact parameter ∆b corresponding to some centrality bin. Since ∆b 1/|p C |, one can identify the average impact-parameter dependent cross section with the average number of hard processes per collision.
for achieving high luminosity. θ C can be determined with an accuracy ∆θ C = 10 µrad, which gives us for E = 7 TeV proton beams. If the measurements like those in Refs. [24][25][26][27][28][29] are subject to a similar constraint (with lower beam energies), the required theoretical calculations always involve the wave packets of colliding protons, as shown in Eqs. (19) and (37), which, however, are not measured 11 .

COLLISIONS
The modes associated with q T ∼ 1/b are non-perturbative in hadron and nuclear collisions. In order to make a model-independent verification of the factorization formula in Eq. (44), in this section we study impact-parameter dependent qq collisions with b 1/Λ QCD .
Since the factorization formula is valid at all orders in α s and at leading order in λ, its validity in such collisions can be verified order by order in perturbation theory by comparing to the results of the general formula in Eq. (27). For this task, the factorization formula is required to reproduce the perturbative QCD results, expanded to leading order in λ.

A. The impact-parameter dependent Drell-Yan cross section
We calculate the impact-parameter dependent cross section for the Drell-Yan process Let us first derive the impact-parameter dependent total cross section from the general formula in Eq. (27). One first integrates over the observable O and then singles out one of the final-state particles as γ * to obtain the general formula for the total cross section with p C = p A + p B − p f and the incoming momenta given by Eq. (23). For unpolarized collisions, it can be written in the following compact form For the factorization formula, one needs to integrate over y C and p C in Eq. (44). This can be easily done by using Eq. (28) to introduce δ(p 2 C − Q 2 ) and then integrate out the four-momentum p µ C instead. We only consider the production threshold: Q 2 ∼ s. In this case, one can replace δ(p 2 C − Q 2 ) by δ(n A · p CnB · p C − Q 2 ) and the p C -integral gives δ (2) (x). Then, integrating out x gives withn i · p C = z ini · P i − n i · p f . Here, we have used the fact that the soft function becomes unity at x = 0 and the hard function H is equal to |M | 2 for the partonic process j, k → C + anything else, as given in Eq. (48). We define which are referred to as transverse spatial PDFs. In heavy-ion collisions, the above equation gives the factorization formula for the inclusive Drell-Yan production of γ * , W ± and Z 0 with a single hard scale Q, valid for all orders in α s . The interested reader is referred to Refs. [37,48] for fixed-order calculations and phenomenological studies. Below, we only focus on the verification of the factorization formula by a detailed calculation for qq collisions at leading order in λ.

B. Factorization at the Born level
At zeroth order in α s , the off-diagonal amplitude squared is given by with e q the electric charge of the quark and antiquark. Expand it around λ = 0, and at leading order one has where s =n A · P AnB · P B and we have used the following identities Inserting Eq. (72) into Eq. (68) and expanding the δ function in Q as well 12 gives That is, the hard scattering is initiated by the quark and antiquark only when they pass very close to one another at a distance ∼ 1/Q, which defines "point-like" in our calculation.
and the hard function is given by Plugging them into Eq. (69) gives the same result as Eq. (75), hence confirming the validity of factorization at the Born level. 12 As an exercise to illustrate the validity of such an expansion, one can work out by using test functions in both b and x.
It also contains collinear and soft divergences.
We use the method of regions (see Ref. [45] for an introduction) to expand the real (Eq. (78)) and virtual (Eq. (80)) integrals in each relevant region of the gluon momentum k µ in order to verify the factorization formula at one loop. The relevant regions include: 1. The hard region: k µ ∼ Q In this region, upon expanding the off-diagonal amplitude squared, at leading order in λ it equals the conventional amplitude squared with p i andp i both replaced by their design values P i . As a result, the expansion of the virtual and real integrals in this region gives with σ (1) the conventional one-loop cross section for qq → γ * , which can be found in, e.g., Ref. [51]. The double poles in virtual and real contributions cancel out but the collinear divergences do not cancel, which produces the 1/ pole in the final result of σ (1) . If one inserts into the factorization formula the zeroth-order transverse spatial PDFs in Eq. (76) and the spin-averaged virtual and real amplitude squared with incoming momenta equal to P i as the one-loop hard function, one evidently reproduces the same result as Eq. (81).
2. Two collinear regions: k µ ∼ Q(λ 2 , 1, λ) n ini Expanded in these regions, the virtual diagrams become scaleless integrals and, hence, vanish. As a result, one only needs to consider the real diagrams as shown in Eq. (79).
Let us take for example the n A -collinear region in which With the above scaling in mind, expanding the δ function in Eq. (80) gives with (1 − z A ) ≡n A · k/n A · P A . Then, expand the real amplitude M (1) r in λ as well. After some algebra, we have One only needs to keep terms of O(λ −2 ) in M. After expanding it according to the scaling in Eq. (82), one can easily obtain The one-loop transverse spatial quark distribution in a fast-moving quark can be calculated analytically. One combines the denominators on the right-hand side of Eq. (93) by introducing a Feynman parameter x: Then, by changing variables tok, one can easily integrate out k in Eq. (87) and obtain where we have used the integral in Eq. (A6) and replaced the bared coupling g 2 s by with µ the MS renormalization scale and γ E the Euler constant. At the end, by using we have The singularity in arises from the fact that the virtual-gluon radiation can only contribute to the prefactor in front of δ(r) and there is no real-virtual cancellation at a nonzero transverse distance r.
Similarly, one can also identify the contribution from the transverse spatial antiquark distribution function of the incomingq by expanding the integrand of Eq. (80) in the n B -collinear region. One hence confirms the contributions from one-loop transverse spatial PDFs in the factorization formula.
3. The soft region: k µ ∼ λQ As a consistency check, we show that the contribution from the soft region vanishes. We expand both the diagram for the one-loop spatial quark distribution function in Eq. (88) and the virtual and real integrals for the general formula given respectively by Eqs. (78) and (80) in this region. The former corresponds to the zero-bin contribution to the spatial PDFs while the latter corresponds to the one-loop correction to the soft function. In both cases, one ends up with scaleless integrals, which vanish in dimensional regularization. Therefore, the soft region is indeed irrelevant.
In summary, we have verified the validity of the factorization formula at one loop: the correction from the one-loop hard function is given by the expansion of the virtual and real integrals for the general formula in the hard region, as given in Eq. (81); the correction from one-loop transverse spatial PDFs is equivalent to expanding these virtual and real integrals respectively in the two collinear regions, as given in Eq. (86) for the quark collinear region; and the correction from the soft function vanishes. In this appendix, we evaluate the integral .
One can choose r to align with the positive x-axis and split the integral into two pieces: Here, with the contour C running from −∞ to +∞. This integral is well-defined only for the range 1 2 < Re < 1, where one can deform the contour from C to C 1 +C 2 as shown in Fig.  2. On C 1 the arguments of x − i0 + and x − i0 − are −i 3π 2 and i π 2 respectively while on C 2 they are both equal to i π 2 . Using this fact, one has I y is defined as With I x and I y evaluated above, one finally has