Multi-dimensional hadron structure through the lens of gluon Wigner distribution

In this review, we present the current status of phenomenological research on constraining the multi-dimensional proton (and nucleus) structure at high energies through studies of the so-called gluon Wigner distributions. We provide a brief pedagogical introduction into the corresponding theoretical definitions and modelling of exclusive and diffractive scattering observables in terms of the Wigner distribution. Also, we present a detailed outlook into the existing and planned experimental measurements that attempt to constrain the Wigner distribution. We briefly overview possible interconnections between various manifestations of the gluon Wigner distribution emerging, for instance, in azimuthal-angle correlations in (semi)exclusive reactions and elliptic flow measurements in inclusive processes. We also summarise the current knowledge on the most important processes that would potentially enable one to constrain the elliptic gluon density in the proton and to separate it from the genuine effect of hydrodynamic evolution in the flow measurements.


I. INTRODUCTION
Over a few decades, a substantial effort of the particle physics community has been devoted to unveiling the proton structure at high energies in both longitudinal and transverse dimensions relative to the collision axis. Since a long ago, it has been well-known that the most detailed information on a given compound quantum system can be unwound through the knowledge of kinematical distributions of its constituents over phase space. In the case of the nucleon target, for instance, this information is contained in the Wigner parton distribution in Quantum Chromodynamics (QCD) [1][2][3][4][5][6][7][8] which represents a comprehensive visualisation of the partonic structure of the nucleon in five dimensions, also referred to as the multi-dimensional parton imaging or tomography. Such an imaging has grown into a paradigm in contemporary studies of hadron structure in high-energy particle collisions [9]. For this reason, the Wigner distribution is also referred to as the "mother of all distributions" [1][2][3] as it is connected to other well known lower-dimensional parton distributions through its integration over one or more dimensions. Its Fourier transform known as the Generalized Transverse Momentum Dependent Distribution (GTMD) [10][11][12][13] is also widely used for modelling of the nucleon structure (for earlier reviews on this topic, see e.g. Refs. [9,14,15]). A comprehensive review of forward small-x physics phenomenology, experimental developments and theoretical models addressing the challenges of the largely yet unknown nucleon structure can be found in Ref. [16].
While the Wigner distribution (and GTMD) encodes all the non-perturbative QCD dynamics of parton constituents inside the nucleon (such as their spatial and momentum distributions as well as the orbital angular momentum [4,17]), it is not calculable in the framework of perturbative QCD and represents the main non-perturbative ingredient of QCD factorization. Consequently, it has been argued in Refs. [18][19][20][21][22][23][24][25][26] that the Wigner distribution can be used to quantify entanglement entropy in the proton wave function. Ongoing developments of phenomenological methods and techniques enabling to constrain the Wigner distribution directly from the experimental data (the co-called "partonometry") is necessary to make further significant advances in nucleon imaging [27][28][29]. A direct measurement of the Wigner distribution appears to be a challenging problem as it generally requires the most detailed knowledge of particle kinematics in the final state in a clean environment (i.e. with the maximal degree of exclusivity).
A large amount of data on various scattering processes are being currently collected by the measurements in hadron and nuclear collisions at the LHC, and it would be desirable to utilise these data to reconstruct the proton structure encoded in the Wigner distribution. Such processes as inclusive DIS or semi-inclusive DIS (SIDIS) are used to probe the integrated (spin-averaged and spin-dependent) collinear Parton Distribution Functions (PDFs) or Transverse Momentum Dependent Distributions (TMDs), respectively, whereas the Generalized Parton Distributions (GPDs) are constrained by the measurements of exclusive reactions, such as the Deeply Virtual Compton Scattering (DVCS) (for a recent review, see Ref. [30]). While these distributions are effectively connected to the Wigner distribution by integrating it over one or more dimensions, none of the above reactions can be utilized as a direct probe for the 5D Wigner distribution itself. Indeed, SIDIS and DVCS are sensitive to the parton longitudinal momentum fraction x and to either the impact parameter b of a given parton inside the nucleon (or nucleus) or to its transverse momentum q.
Diffractive processes such as those in Deep Inelastic Scattering (DIS) of leptons off nucleons or large nuclei play an important role in extracting the crucial information on the structure of the target particle in both coordinate and momentum spaces [31]. All of this information is encoded in the parton Wigner function W (x, q, b) effectively representing the (quasi)probability to find a parton with a given transverse momentum q at the impact parameter b (its transverse separation from the center of the nucleon or nucleus), and with a given momentum fraction x. Among important examples, diffractive vector meson production effectively probes the spatial profile as well as fluctuations of the gluon field in the target [32][33][34].
Throughout this review, we are focused mainly on the gluon Wigner distribution 1 relevant for both hadro-and photoproduction reactions at high energies, i.e. at x ≪ 1. Once the collision axis is fixed and both 2D vectors b and q are in the plane transverse to this axis, we deal with the non-trivial 5D gluon Wigner distribution, W (x, q, b), which is the subject of phenomenological studies of gluon tomography of the nucleon.
As an important example of this research, the pioneering study of Ref. [36] has proposed that the gluon Wigner distribution at a given momentum fraction x ≪ 1 can be probed experimentally in exclusive dijet photoproduction in Deep Inelastic Scattering (DIS), in particular, by measuring the correlation in azimuthal angle between the produced dijet transverse momentum and the recoiled nucleon transverse momentum (see also Ref. [37]). A follow-up work of Refs. [38,39] (see also Ref. [40]) has elaborated on this possibility more quantitatively considering exclusive light-quark dijet and heavy quark-pair photoproduction in ultraperipheral collisions (UPCs). A comprehensive analysis of diffractive dijets production in ep collisions in the Color Glass Condensate (CGC) framework in the context of the Electron-Ion Collider (EIC) has been performed in Ref. [41]. Theoretical uncertainties in the existing models for phase space distributions that actually predict non-trivial angular correlations in impact-parameter space, namely, between the impact parameter b and the dipole separation r, remain large [40]. Hence, it is mandatory to experimentally access such correlations that would enable us to understand the nucleon structure in greater detail.
Given the multi-dimensional character of the Wigner function, rather large-statistics data samples would be needed if one wanted to map all (five) dependencies. Since very low pileup data 2 conditions are preferable, collecting large statistics may require larger integrated luminosities and, hence, longer data-taking periods, depending on the actual average pile-up rate. Furthermore, considering for instance pA UPCs, in the ideal case both, the proton and nucleus on opposite sides should be tagged. If tagging nucleus is problematic or not possible, one could also require a veto in Zero-Degree-Calorimeter (if it exists). If the efficiency of such a veto turns out to be low, one can loosen this requirement by imposing an energy threshold in the forward direction -its value may depend on a given integrated luminosity.
Especially in the case that the double-tagging mentioned above is not possible or inefficient, one could also require a large rapidity gap to accompany the intact nucleon/nucleus on both sides from the interaction point (IP). The rapidity gap is an extremely interesting quantity from both, the theory and experiment points of view. Experimentally, the size of the rapidity gap inevitably depends on the noise level of individual subdetectors and energy thresholds of objects we wish to veto (tracks, calorimeter clusters, jets). A typical lowest p T threshold of jets reconstructed from calorimeter clusters/towers at ATLAS or CMS below which the jet reconstruction is significantly less reliable is 20 GeV, while it can go down to 10 GeV, if tracks are considered. But vetoing jets with such p T thresholds is not considered to be sufficient in reaching exclusivity. One should also be aware that in the context of measurements in UPCs, the incoming quasi-real photon can reveal its structure and hence its remnants then can spoil the gap(s) representing one of the major challenges in precision measurements of exclusive photoproduction processes.
In this review, we summarise the basic properties of the gluon Wigner distribution (and the corresponding GTMD) and the existing methods to access it through detailed measurements of exclusive processes, including the existing phenomenological models and the current experimental situation, as well as their role in our understanding of the nucleon structure. Theoretical modelling of these processes is most conveniently performed in the framework of color dipole approach which has been thoroughly formulated together with the key results from the literature. Besides, we briefly overview the corresponding challenges and published measurements which can be sensitive to the gluon Wigner function, with a potential to probe it either with the existing collected data or with future data.
The review is structured as follows. In Sect. II, we provide key definitions of the phase space parton distributions in QCD and a brief discussion of fundamental properties of the Wigner distributions, with a particular focus on the gluon case, connecting it to the fundamental ingredient of QCD scattering, the dipole S-matrix. Sect. III focuses on observable effects of the dipole orientation with respect to the color background field of the target such as azimuthal angle correlations quantified by the elliptic flow. In Sect. IV, we outline the basics of the phenomenological dipole approach and the dipole orientation effects on phenomenological grounds. Namely, the basic properties of the dipole S-matrix and its widely used phenomenological modelling are discussed in connection to the azimuthal angle correlations. More elaborate QCD evolution based approaches are briefly presented in Sect. V. Sect. VI reviews a few most relevant processes as suitable probes for dipole orientation effects and the associated elliptic gluon Wigner distribution, as well as provides some of the key insights into the corresponding experimental measurements. A brief overview of future measurements that could offer new possibilities for constraining the gluon Wigner distribution is given in Sect. VII. Finally, Sect. VIII concludes the review summarizing its main takeaways.

II. PHASE SPACE DISTRIBUTIONS IN QCD
A remarkable property of QCD is a dramatic difference in complexity of its UV and IR regimes. While in the UV regime one enjoys the simplicity of perturbative description stemming directly from the QCD Lagrangian given in terms of fundamental degrees of freedom, quarks and gluons (partons), in the IR regime one encounters an enormous complexity associated with dynamics of the QCD ground state (quark and gluon condensates), bound states, their interactions and the corresponding nonperturbative QCD phenomena such as confinement and mass gap generation. At short distances, i.e. very far from the "confinement boundary" in the IR regime, an approximation of relatively weakly interacting partons as asymptotic states of the theory works well, at spacetime separations close or larger than this boundary such an approximation breaks down as the asymptotic states can no longer be coloured partons but colorless composite structures of different types (hadrons). A consistent transition between the parton and the hadron levels remains a long-standing theoretical challenge being so far exclusively addressed by lattice QCD simulations in Euclidean spacetime. No single consistent theoretical approach that bridges and unifies both the sub-confined (partonic) and super-confined (hadronic) regimes in a consistent frame-work in Minkowski spacetime has been developed yet (for a recent review on the status of nonperturbative QCD and confinement research, see Ref. [42]). One particularly important aspect of this problem is the hadron structure highly motivated by the practical need for interpretation of collider measurements (such as those at the LHC) and actively explored in ample literature over past few decades, from both theoretical and phenomenological points of view [15].
In scattering of high-energy hadrons, a fast-moving hadron whose internal partonic structure is being probed (or which plays a role of a probe itself) sets the longitudinal, or Light-Cone (LC) direction, such that the motion of its constituents is effectively separated into "longitudinal" and "transverse" parts. At short distances, distributions of asymptotically free partons in their longitudinal momenta in a hadron are described in terms of collinear parton distributions -the basic non-perturbative ingredients of collinear factorisationdefined in a single space (momentum) dimension as functions of parton's dimensionless longitudinal momentum fraction. The collinear factorisation concept works out extremely well in a myriad of different inclusive processes in lepton-hadron and hadron-hadron collisions, and has highlighted the triumph of QCD-improved parton model. However, in certain cases the collinear parton evolution picture does not capture all relevant phenomena, and a more detailed information on parton motion in the transverse plane (to the direction of the parent hadron) is required.

A. Basic definitions
As a suitable choice of the reference frame in studies of a lepton-nucleon or nucleonnucleon scattering process, one typically considers a frame where both colliding particles move fast, e.g. center-of-mass frame, or where one of them is at rest (the target) while the other one moves fast (the projectile), known as the target rest frame. Either way, longitudinal and transverse directions play very different roles enabling the LC decomposition of an incident 4-vector l in terms of its LC components: +/− LC projections l ± = (l 0 ± l 3 )/ √ 2 and 2D transverse momentum l = (l 1 , l 2 ), with its absolute value l ⊥ ≡ |l|. Then, a parton carrying 4-momentum k inside a high-energy nucleon is characterized by its longitudinal momentum fraction x ≡ k + /P + that it takes from the parent nucleon momentum P + , transverse momentum k and the transverse position with respect to the center of the nucleon, b. The phase space distributions separately in k and b are known in the literature as the transverse momentum dependent distribution (TMD) T (x, k) ≡ T (x, k ⊥ ) and the generalized parton distribution (GPD) G(x, b ⊥ ), respectively. The former quantifies transverse momentum of a parton in the considered parent hadron affecting the distribution of produced particles, whereas the latter describes the spatial distribution of partons in the transverse plane which can be found via a Fourier transform over the transverse momentum transferred to the parent hadron (causing, in particular, its deflection at a small angle). While being important ingredients of exclusive cross sections, these distributions encode crucial information on the two-dimensional structure of the nucleon, quantifying it in (x, k ⊥ ) and (x, b ⊥ ) planes of the phase space. In each such case, one probes dynamics of partons at typical hadronic length (or momentum) scales being sensitive to non-perturbative QCD phenomena (such as confinement).
However, in several situations an even more detailed information about distribution of partons in both the k and b variables is required i.e. accounting also for relative (azimuthal) angle between these 2D vectors. One of the most important examples is the orbital angular momentum ∝ b × k particularly relevant for understanding of how the total nucleon spin emerges from the angular momenta of its constituents and their spin-orbital correlations, hence being of significant phenomenological interest.
The most generic object in QCD describing the parton dynamics in the nucleon is the so-called two-parton correlation function. For instance, the manifestly Lorentz-covariant two-quark correlation function is represented as a matrix element of the bilinear quark operator, given in terms of the nucleon state |P ⟩ and its momentum P µ , the momentum transfer in the t-channel ∆ µ = (0, 0, ∆) (|∆| ≡ ∆ ⊥ = √ −t), the quark-antiquark relative separation z µ = (0, z − , z), the staple-shaped Wilson line along z − L introduced to ensure gauge invariance of the corresponding operator, and some Dirac matrix Γ = γ + , γ + γ 5 , . . . defining the spin structure and the relevant twist [44]. The incident momenta satisfying the on-shell nucleon condition, P ∆ = 0 and 4P 2 + ∆ 2 = 4m 2 , with m being the nucleon mass, are depicted in Fig. 1. For a pedagogical discussion of properties of H(k, P, ∆) and its phenomenological significance in different kinematical domains, see e.g. Ref. [43] (see also Ref. [15] and references therein). Considering particular scattering processes, one typically integrates H(k, P, ∆) over one or more components of the parton 4-momentum k giving rise to distribution functions in lower number of dimensions as schematically illustrated in Fig. 2. Starting with integration over k − , one notices that the virtuality of partons is not definite any longer as they are in general off-shell in consistency with confinement. Moreover, the parton evolution in this case is naturally consistent with the parton model in the LC quantisation picture where asymptotically free parton fields are quantised at LC time z + = 0, i.e. just before being probed in scattering. In this case, the k − -integrated parton correlation function gives rise to the six-dimensional GTMD G(x, k, ξ, ∆) describing the emission and absorption of the (anti)quark as well as the emission and/or absorption of a quark-antiquark pair as shown in Fig. 1. Here, the second LC fraction ξ = −∆ + /2P + known as skewness variable determines the total momentum transfer squared as In the forward limit suitable e.g. for exclusive reactions in high-energy lepton-nucleon (lepton-nucleus) or nucleon-nucleon (nucleus-nucleus) collisions, ∆ + → 0 (hence, vanishing ξ), we arrive at the five-dimensional parton GTMD G(x, k, ∆) where the dependence on Relations between different distributions originating from the two-parton correlation function. From Ref. [43].
the skewness parameter ξ is suppressed. Thus, its Fourier-transform in the mixed (impact parameter/momentum space) representation W (x, b, k) describes evolution in both the particle transverse momentum k and its average position in the transverse plane b and is known as the Wigner distribution [1] whose notion has been used in physics in many different contexts. In particular, in the context of exploring the nucleon substructure, the quark Wigner distribution takes the following form [2,3] in terms of the two-quark correlation function defined above. In analogy to quarks, following Ref. [45] (see also Refs. [10,12]), the gluon Wigner distribution in the proton can formally be written as follows in terms of the staple-shaped Wilson line U ± that goes to z + = ±∞ on the LC, and then returns. For further convenience, denoting initial and final proton momenta as p = P − ∆/2 and p ′ = p + ∆, here and below, the brackets ⟨. . . ⟩ denote the off-forward proton matrix element implicitly normalised as ⟨p ′ | . . . |p⟩/⟨p|p⟩.

B. Properties of the gluon Wigner distribution
The parton Wigner distributions contain the most detailed and complete information about the quark and gluon evolution in five-dimensional phase space, and hence encode the parton structure of the nucleon wave function [1][2][3] (see also Refs. [10,12]), being of significant phenomenological interest. In quantum theory, however, such a definition immediately contradicts the uncertainty principle δxδp ≥ ℏ/2 as the particles' position and momentum cannot be determined simultaneously with arbitrary precision. In practice, this fact causes the Wigner distribution to violently oscillate and to become negative in some domains of the phase space. The Wigner distributions (2.3), (2.4) have been evaluated in various models, in particular, for quark one see Refs. [2-5, 28, 46-48]. In simple models without gluons, it turns out to be a positively-definite function. However, once gluons are included in the full QCD description, the Wigner distribution becomes negative for some phase space configurations [5]. For that reason, in QCD the parton Wigner distribution does not have an interpretation of a joint probability distribution in both k and b, but rather of a quasi-probability one, due to not being manifestly positively-definite.
This situation has motivated the authors of Refs. [49,50] to turn to the so-called Husimi distribution as a close analogue to the classical phase-space distribution. In quantum mechanics, it has the meaning of a positively-definite probability density to locate a particle not at precise values of its momentum and coordinate but within an uncertainty band (x ± δx, p ± δp), with δx and δp satisfying δxδp ≥ ℏ/2. In QCD, a transition to Husimi distribution implies a Gaussian smearing of the Wigner distribution in both k and b as follows [49] W(x, k, b) = 1 where the Gaussian widths are related inversely to each other in consistency with the uncertainty relation, hence, ensuring that the resulting distribution is positive semi-definite. Interestingly enough though, by integrating W over one of the transverse vectors yields a proper probability distribution in the remaining one. Indeed, integrating W over k (with an appropriate regularisation procedure) gives rise to the parton distribution f (x, b ⊥ ) in impact parameter space representing a probability density for a parton carrying longitudinal fraction x to be localised at impact parameter b from the nucleon center. A 2D Fourier-transform of f (x, b ⊥ ) represents a GPD in momentum space, G(x, k ⊥ ), for which in general ξ ̸ = 0. On the other hand, integrating W over b results in the corresponding TMD, i.e.
From this point of view, the Wigner distribution is often dubbed as the 'mother distribution' for both TMDs and GPDs. Further possible connections to lower-dimensional distributions are briefly summarised in Fig. 2. As was mentioned above, for a longitudinally polarized nucleon, the Wigner distribution is also related to the canonical orbital angular momentum [3,11] which is an important ingredient in the nucleon spin decomposition (for more details, see also Refs. [51,52]).
C. Gluon Wigner distribution in the dipole picture Traditionally, diffractive particle production processes, for instance, with qq di-jets in the final state, are often regarded as a convenient probe for hadron structure (see e.g. Ref. [53] and references therein). As was demonstrated in Ref. [36], the gluon Wigner distribution at small momentum fractions of the gluons in the target, i.e. x ≪ 1, gets dramatically simplified and can be accessed experimentally in diffractive dijet (electro)production (see also Ref. [37,38]). Indeed, through its connection to the gluon GTMD, which, following Refs. [45,54,55], reads 11) in terms of the coupling constant α s and the number of colors N c = 3 in QCD, one straightforwardly identifies the impact parameter dependent off-forward dipole S-matrix in the eikonal approximation, which determines the dipole survival probability, as a product of two light-like Wilson lines in the fundamental representation U of a quark at b + r 2 and an antiquark at b − r 2 transverse positions, respectively. The Wilson line is an important ingredient of the dipole scattering describing the eikonal propagation of the quark in the colored medium.
In order to see this more explicitly, let us consider the γp scattering in the dipole picture [56,57], where the real photon fluctuates into a qq-dipole that scatters off the target proton p as illustrated in Fig. 3. Here, let us fix the frame in which both initial γ and p are collinear, with the fast proton propagating in the positive z-direction. Then, the process is most conveniently described in the impact parameter representation where the transverse coordinates of the quark and antiquark are set to be respectively. Here, z stands for the longitudinal momentum fraction of the quark taken from the incoming real photon. The lowest Fock state of the photon is described by γ → qq wave function [58], with an emerging qq dipole, whose "center of gravity" is located at impact parameter b ≡ zx 1 + (1 − z)x 2 , which is the same as that of the incoming photon. The transverse separation of the dipole is then found as r ≡ x 1 − x 2 defining the strength of the dipole-target scattering, In the above expressions, the averaging ⟨. . . ⟩ x in the target is defined in the CGC formalism [59,60] (for further details on the gluon Wigner distribution in the CGC framework, see Sect. V B below) yielding the dipole S-matrix as a function of "rapidity" Y ≡ ln(x 0 /x), x ≪ 1, with some adjustable parameter x 0 typically fixed to x 0 = 0.01 unless noted otherwise. Hence, the energy dependence of the dipole S-matrix is effectively encoded through its Y dependence. In the particularly important example of diffractive DIS process, γp → Xp, x is identified with x IP being the LC momentum fraction of the projectile dipole carried by a collective t-channel exchange (often regarded as the Pomeron exchange) between the dipole and the target, and is found as in terms of the transverse mass M ⊥ of the produced diffractive system X. In case of dijet production in diffractive DIS off the proton target, for instance, one has (2. 16) in terms of photon virtuality, Q 2 , the proton mass m p , the c.m. energy squared of the photon-target system, W 2 , and the momentum transfer squared given by the Mandelstam t variable.
Then, it is instructive to switch to momentum space by means of double Fourier trans-form, S Y , as follows [61] where the exchanged gluon transverse momenta, k 1 ≡q + z∆ and k 2 ≡q − (1 − z)∆, flow in opposite directions and are conjugate to the impact parameters, x 1 and x 2 , respectively. The phase factor e −iδ·r recovered in Eq. (2.17) is invariant under interchange of quark and antiquark, i.e. z → 1 − z and r → −r, being the exact symmetry of the dipole amplitude [62] 3 . This phase factor becomes relevant, for instance, in the DVCS with the longitudinally polarized virtual photon [61], while it does not play any role in diffractive dijet production [36]. One then arrives at the gluon GTMD as a Fourier transform of the dipole S-matrix defined in Eq. (2.12). Finally, one finds [36] x Hence, the Wigner distribution is the basic ingredient in the description of the forward cone in various diffractive processes as it naturally incorporates dependence on the transverse momentum transfer between the diffractively produced system and the target, ∆. Coming back to the Wigner distribution, one writes 20) in terms of the scattering matrix T Y found conventionally as S Y ≡ 1 − T Y . Due to the fact that the colorless dipole of a vanishing size ceases to interact with the target -the so-called color transparency property -no scattering off the target gluons occurs leading to T Y → 0, hence, S Y (r, b) → 1 as r ⊥ ≡ |r| → 0. This implies the 'sum rule' while d 2 q S Y (q, ∆) = 0 for ∆ ⊥ ̸ = 0, in consistency with unitarity of the S-matrix.
The above derivation highlights the proof of Ref. [36] showing that the gluon Wigner distribution at small x can be expressed in terms of a Fourier-integral containing the dipole S-matrix defined in the impact parameter space, S Y (r, b). Hence, one may conclude that the framework of Wigner distributions effectively unifies three common descriptions based on the color dipole approach [56,57] as well as on the TMDs and GPDs in the parton distributions approach. It not only gives us access to the transverse momentum distribution of the considered parton at a given rapidity (or x), but also provides an average position of it relative to the center of the target particle.
Such an important connection between the quasi-real Weiszäcker-Williams (WW) type gluon distribution in the nucleon in its rest frame and the dipole scattering matrix probed by a projectile qq dipole propagating through the color background field of the nucleon at small x represents a universality between the corresponding two topologically different gauge-invariant operators [36]. The partial dipole amplitude is an important ingredient of the color dipole picture of QCD scattering where dipoles are considered to be eigenstates of the scattering operator [56] (see below) such that practically any QCD scattering process can potentially be represented as a superposition of different "elementary" dipole scatterings. This concept has been widely explored in a vast literature starting originally from the inclusive DIS in ep collisions, then applying it in inclusive hadron production in pp, pA and AA collisions as well as for a variety of exclusive and diffractive reactions. Typically, the dipole amplitude gets averaged over possible dipole orientations with respect to the background field as most of the measured processes are not sensitive to azimuthal-angle correlations between r and b. From the theory point of view, considerations of Refs. [36,38,39] have pointed out that the non-trivial correlations in angle between the impact parameter and the dipole separation should be sensitive to possible correlations in angle between the gluon's impact parameter and its transverse momentum in the corresponding gluon Wigner distribution. The dependence on the dipole orientation in inclusive processes, for instance, is known to induce flow-like azimuthal angle correlations such as those in production of prompt photons [69], or in two-particle production processes [70]. Below, we will elaborate on the possible origin and significance of such correlations in the small-x regime relevant for high-energy collisions in the framework of the color dipole picture.

D. Elliptic gluon Wigner distribution
The angular correlations of the differential cross sections are connected, in particular, to a correlation of the partial dipole amplitude in azimuthal angle between the impact parameter b and the dipole separation r in the course of dipole propagation through the colored medium of the target. One often distinguishes two possible sources of such correlations. The first one associated with the initial-state structure arises from the gluon Wigner distribution (or the corresponding GTMD) and represents a correlation in the angle between the transverse momentum and impact parameter vectors in the transverse plane. The second source is related to dynamics of the produced states in the medium and, to the leading order, is associated with the relative angular dependence of the (anti)quark plane-waves in the qq dipole which probes a gluon in the target nucleon. While both sources of dipole orientation effects may have an impact on multi-partilce correlation observables in a scattering process, a precise prediction of their effect in different kinematical regimes remains a challenging problem under intense studies in the literature. For a comprehensive discussion on the possible origins of azimuthal correlations in small systems, see e.g. Refs. [40,[71][72][73][74][75][76][77][78][79][80] and references therein.
When the initial-state effects encoded in the gluon Wigner distribution are concerned, the dipole orientation effects give rise to a characteristic correlation in azimuthal angle between q and b determined by the so-called "elliptic" gluon Wigner distribution [36,50,81]. Indeed, assuming that the elliptic ∝ cos 2θ dependence of the partial dipole amplitude is dominant, the gluon GTMD decomposition into the angular-independent (isotropic) S Y,0 and elliptic S Y,ϵ components can be written as [36,50,73,81] ignoring small higher harmonics. For hard scattering processes, the elliptic Wigner distribution S Y,ϵ is relatively small (of the order a few percent, see below) compared to the isotropic one, S Y,0 , but exhibits very different dependencies on x and q ⊥ [50]. Indeed, their asymptotic behavior is generically expected as follows [38] Thus, the presence of the elliptic density may lead to distinct experimental signatures such as elliptic flow etc being the subject of intense explorations in the literature (see e.g. Refs. [36,73,81]). The elliptic gluon distribution S Y,ϵ provides crucial information on small-x gluon tomography in the proton and can potentially be constrained through measurements of different processes discussed below, in particular, at the future EIC (for a comprehensive review of the EIC concepts and research, see e.g. Refs. [9,14,82]). The above decomposition (2.22) is typically used in analysis of such observables as, for instance, the differential cross section of exclusive di-jet photoproduction offering interesting opportunities for a direct measurement of the elliptic gluon GTMD. Still, a more generic analysis of feasibility of direct probes for the elliptic gluon Wigner distribution in differential observables accounting for both the initial-and final-state effects on the same footing is still missing, and it remains yet unclear to what extent the elliptic component S Y,ϵ (q, ∆) can actually be accessed in a measurement. For a more detailed discussion on this topic, see Sect. VI A below. Let us first elaborate in more detail on physical origins of such correlations and, further on, how the azimuthal angle correlations can be modelled in the dipole picture of QCD scattering from theoretical and phenomenological points of view.

III. AZIMUTHAL CORRELATIONS AND THE ELLIPTIC FLOW
The experimentally observed fact that the hadron production processes in both protonproton and proton-nucleus collisions at high multiplicities exhibit unexpectedly large azimuthal asymmetries (see e.g. Refs. [83][84][85][86][87][88][89]) has prompted continuous debates in the literature on the origins of such asymmetries. The long-standing theoretical problem is to consistently explain why small systems emerging in pp or pA collisions reveal a very similar collective behavior as order-of-magnitude larger systems formed in AA collisions. Such a universality between small and large systems has been found in various observables differential in rapidities, transverse momenta and masses of final-state hadrons produced in such reactions [87,[90][91][92][93].
Two different approaches are typically put forward to interpret such an universality and confronted to each other. One of them resides on a picture of collectivity due to the finalstate interactions, also known as the hydrodynamic flow -an effect typically considered to be natural in the hydrodynamic description of large systems, at least, for small transverse momenta (see e.g. Refs. [94][95][96][97][98][99][100]). The second approach relies on interpretation of collective phenomena as due to universality in the initial state connected to parton densities in the incoming hadron or nuclear wave functions before they collide (see e.g. Refs. [101][102][103][104][105][106][107][108][109][110][111][112][113][114][115].) In order to measure an asymmetry in the azimuthal angle ϕ reflecting a spontaneous breakdown of the rotational symmetry in the plane transverse to the collision axis, it is of practical relevance to consider an inclusive single-particle distribution as a function of ϕ on event-by-event basis. Experimentally, such an angle is determined in the so-called reaction plane as an angle between the momentum of a produced state and its impact factor, both measured in the transverse plane. Then, the Fourier cos(nϕ)-moments of such a singleparticle distribution characterise, in particular, an ellipticity of the overlapping (interaction) domain of the two non-centrally colliding states. These moments are commonly known as the "flow coefficients" denoted as v n (p ⊥ ) [116] and are traditionally considered as an effective measure of the azimuthal asymmetry. The latter can be relevant even for central collisions and is induced by fluctuations in transverse momentum distributions of incident particles in the compound colliding states such as nucleons in a nucleus or partons in a nucleon [117]. In this sense, v n is sensitive to inhomogeneities in the transverse-plane parton distributions effectively probing variation in the gluon transverse momentum distribution in the target. In the standard hydrodynamic simulations in the case of heavy-ion collisions, the particle number fluctuations are naturally present in the initial conditions, hence, triggering an asymmetry even without final-state interactions. Such an 'initial-state approach' is typically formulated in the effective framework of CGC [60,118].
The high energy limit of QCD (i.e. at very small x) is described in terms of the CGC effective theory which provides a suitable framework for analysis of diffractive processes. According to the CGC picture of the fast moving target, one distinguishes the gluon fields with large occupation numbers ("slow modes"), described by classical Yang-Mills equations of motion, and their sources ("fast modes"). Consequently, the target nucleon in this picture has a sub-structure at very high energies describing the formation of regions of a typical size ∼ 1/Q s , in terms of the semi-hard saturation scale Q s , with large occupation numbers of coherent gluons. Each such region, also known as the "saturation domains", can be interpreted as being in the ground-state (condensate) of strong chromo-electric fields. The orientation of these fields is not symmetric setting up a preferred direction in the transverse plane, thus breaking the rotational invariance and inducing azimuthal correlations in the final-state particles' distribution. The transverse momentum dependence of the flow parameters v n (p ⊥ ) appears to have a maximum at around p ⊥ ∼ Q s as expected from the domain-like substructure of the target. This phenomenon demonstrates the relevance of a realistic saturation model incorporating the correct b ⊥ -dependence in the inhomogeneous transverse-plane distribution of soft gluons in the target that should be tested against the experimental measurements.
Indeed, consider the CGC dipole picture of a very forward particle production (and/or a heavy nucleus in the target) when a dilute projectile probes a dense target at a very low x. In this case, if incoming particles interact with the same saturation domain at a given impact parameter, they would get similar "kicks" giving rise to final states emerging at similar angles. This effect would be washed away in the case of inclusive single-particle production when averaging over scatterings off different domains since their color fields' orientations are assumed to be random. However, very interesting residual effects of such correlations in the CGC approach remain in observables corresponding to production of, at least, two particles in the final state with "semi-hard" spectra decreasing with Q s . Indeed, as Q s grows, v n decrease with energy, while the scaling v n ∼ 1/N 2 c is valid for a large number of colors N c → ∞. Interestingly enough, in the case of non-Guassian correlations the effects of fluctuations in the CGC approach can be reproduced by the color field domain model of Refs. [106,[119][120][121].
For the scattering cross section to have a non-trivial correlation in azimuthal angle between transverse momentum p of the produced hadron, acquired due to scattering, and its impact parameter b, there must be an inhomogeneity in the small-x gluon density in the transverse plane [69,73,81,107,122,123] (see also Ref. [70]). The latter is precisely characterised by the elliptic part of the gluon Wigner distribution in the target in the high-energy (dipole) factorisation valid for 'dilute-dense' collisions. Hence, the flow coefficients v n , more precisely the elliptic flow v 2 , should be sensitive to (or rather determined by) the presence of the elliptic Wigner density of the target which may adopt a finite value even in the large-N c limit due to not relying on non-planar gluon exchanges only. Just like the elliptic density itself, the azimuthal asymmetry is directly connected to the dipole orientation phenomena encoded in the dipole S-matrix.
The simplest illustration of how to probe the azimuthal angle correlations in the target rest frame is via the inclusive single particle production at forward rapidities such as quark production in collinear quark-target scattering, qT → qX, where T = p(A) is the proton (nucleus) target. The scattered quark has impact parameter relative to the center of the projectile proton,b ≡ b − b ′ , and acquires a transverse momentum p due to multiple gluon exchanges with the dense target. The corresponding cross section, differential in both p and the projectile proton impact parameter relative to the center of the target b ′ , can be written in terms of the dipole S-matrixS Y (p, b) ≡S(x g , p, b) in the mixed (transverse momentum/impact parameter) representation as where η is the (pseudo)rapidity of the final-state quark, xq(x, b) is the quark GPD incorporating the quark impact parameter dependence in the dilute projectile proton. In the soft regime, one typically adopts the factorised form of the quark GPD, with an isotropic profile of the projectile, where xq(x) is the collinear quark density in a dilute projectile. The dependence on the transverse size of the projectile proton becomes irrelevant in case of nearly homogeneous large targets in central collisions, the projectile proton can be considered point-like such that f (b) = δ (2) (b) to a good approximation. Above, x q,g are the longitudinal momentum fractions of the projectile quark and exchanged gluon given by respectively, in terms of the total c.m. energy of the (projectile-target) scattering, √ s. Provided that the incoming quark propagates at a fixed transverse position (although different one in the amplitude and its conjugate), multiple gluon exchanges are viewed as effectively resummed to all perturbative QCD (pQCD) orders using the eikonal approximation in the impact parameter space. The latter, after squaring the scattering amplitudes, gives rise to the qq dipole S-matrix S Y at an impact parameter b and the dipole separation r. The transverse-plane positions of the quark and antiquark in the dipole relative to the target can then be written as respectively, while the rapidity separation reads Y = ln(x 0 /x g ) ≫ 1. The unitarity condition (2.21) in this case implies that the total probability for the final-state quark to acquire some momentum p is equal to unity, with the probability densityS Y (p, b) for the quark to emerge at the impact parameter b with the transverse momentum p. Due to a finite size of the dipole, the differential cross section depends on the azimuthal angle ϕ between p and b ′ , thus being sensitive to the dipole orientation with respect to the target color field. In order to measure an anisotropy in ϕ, one introduces a set of flow coefficients, Using the relation (3.1) between the cross section and the dipole S-matrix, the dominant ("elliptic") flow coefficient can be represented as in terms of the dipole S-matrix in the impact parameter representation, S Y (r, b) ≡ S Y (r ⊥ , b ⊥ , θ). Here, the angle α is made by the vectors b and b ′ while θ -by vectors r and b. In the case of a point-like projectile, this expression simplifies to thus, the dependence on the projectile quark distribution is conveniently cancelled in the ratio such that the elliptic flow is determined by the elliptic gluon density in the target.

IV. PHENOMENOLOGICAL MODELLING OF THE DIPOLE S-MATRIX
Being inspired by intricate connections between the Wigner distribution and the partial dipole amplitude discussed above, let us elaborate on the main aspects and theoretical basics of the dipole formulation of QCD scattering processes. Such processes are considered in the so-called dipole (infinite-momentum) frame sometimes regarded simply as the target rest frame (with the target typically being a nucleon or nucleus) corresponding to "dilute-dense" collisions. This is due to the fact that the target gluon density is much higher than that of the projectile hadron as caused by a small-x evolution. The final-state particles in this case originate from a dipole or a system of dipoles and are produced at very forward rapidities. The dipoles consist of fast, (nearly) collinear, partons in a "dilute" projectile that undergo multiple scatterings off gluons in a dense target, thus, acquiring some non-zero transverse momenta. In the eikonal approximation, such gluon exchanges can be effectively resummed to all orders in perturbation theory in the impact parameter representation [124,125]. At the parton level, such a scattering in the leading order may be considered in perturbative QCD as a probe for very high gluon density in the target -subject to continued theoretical and experimental investigations.

A. Universality of the dipole description and DIS
The universality of the dipole picture of QCD scattering resides on the basic argument of Ref. [56] suggesting that the color dipoles are eigenstates of QCD interactions. This promotes the dipole scattering matrix as the universal fundamental ingredient of any QCD scattering process such that any inclusive and diffractive process can be represented as a superposition of "elementary" elastic dipole scatterings. Indeed, color dipoles in this approach are considered to experience elastic scattering only and the corresponding partial amplitudes are characterized only by a transverse dipole separation r at a fixed small value of Bjorken x ≲ 0.01, where the dipole approach is considered to be valid. For illustration, one often considers a projectile hadron state |h⟩ as a subject of excitations due to QCD interactions that can be decomposed into the orthonormal complete set of eigenstates of interactions |α⟩ labeled by a continuous α variable [56,126,127] in terms of the universal elastic dipole scattering operatorf el and its eigenstates f α . Then, elastic h → h and single diffractive h → h ′ scattering amplitudes are represented as, respectively, such that a diffractive excitation is realised as a result of quantum fluctuations in projectile hadron state [127,128]. Finally, the forward single-diffractive and total hadron-proton scattering cross sections read respectively, with the former being expressed as a dispersion of the eigenvalues distribution, σ α is the eigenvalue of the latter known as the universal dipole cross section [56]. Turning to a continuous set of kinematical variables in the dipole (impact parameter) representation, the phase space of dipole scattering is characterised by intrinsic separation between its constituents (quark and antiquark), r, such that in terms of the normalised h → qq transition LC wave function Ψ h (r, z), with z being the incident quark momentum fraction in the projectile hadron h. For small dipoles (compared to the typical hadron size R had ∼ Λ −1 QCD ), σ qq is a function of momentum fraction x of the gluon exchange (i.e. the Bjorken variable of the target gluons). An important example to mention is a heavy system (with total momentum P ) production in pp collisions such as virtual photon in the Drell-Yan process studied in the dipole picture in Refs. [129][130][131] (see also Ref. [132,133] and references therein). Another important example is the inclusive heavy flavour production in pp and pA collisions successfully treated in the dipole picture beyond QCD factorisation [134][135][136][137][138]. The exchanged gluon momentum fraction is found in terms of the phase space variables as in terms of the pp center-of-mass energy s, and M , k T and x F being the invariant mass, transverse momentum and Feynman variable of the produced system, respectively. In the very high energy limit and for soft dipole scattering the dipole-proton collision center-of-mass energy, √ŝ (withŝ = x 1 s), is a more appropriate variable in the dipole cross section than x as advocated e.g. in Refs. [129,139]. Remarkably, it has been demonstrated in Ref. [140] that in the case of the Drell-Yan process the dipole model yields the differential cross sections fully consistent with the results of the next-to-leading order (NLO) collinear factorisation framework. This observation provides a solid argument that the dipole approach effectively incorporates the dominant NLO QCD corrections. Besides, as ample literature suggests, the dipole formalism provides a universal way to quantify such effects as saturation, initial state interactions, nuclear coherence as well as the gluon shadowing (see e.g. Refs. [134][135][136][141][142][143][144] and references therein).
Let us now consider an important example of the inclusive DIS process in γ * p collisions being one of the most precise and rich source of phenomenological information on the dipole scattering in QCD which historically came about due to HERA measurements [145,146]. The optical theorem suggests that the DIS cross section is given in terms of the imaginary part of the amplitude of the elastic virtual photon γ * N → γ * N (with virtuality q 2 = −Q 2 and with longitudinally (L) and transversely (T) polarised photon, λ γ = L, T ) scattering off the nucleon target N in the forward limit, i.e.
where the transverse momentum of the outgoing nucleon is denoted by ∆, and the DIS structure function -by F 2 (x, Q 2 ). Then, a standard illustration of the DIS γ * N → X cross section at small Bjorken x can be seen in Fig. 4 in both Regge (left) and dipole (right) pictures. Indeed, the virtual photon reveals its partonic structure represented as its lowest Fock component to the leading order γ * → qq -the colorless qq dipole characterised by its transverse separation r as introduced above. Such a dipole undergoes scattering off the color field of the target at impact parameter b through a small-x color-singlet compound state exchange, often regarded as a Pomeron exchange in Fig. 4 (left) (see e.g. Refs. [147,148], as well as Ref. [149] and references therein).
In QCD, such a Pomeron is viewed as a colorless multi-gluon exchange in the t-channel such as a di-gluon ladder in Fig. 4 (right), thus, connecting to the unintegrated gluon density (the small-x gluon GTMD, in general kinematics) in the target. The corresponding dipole S-matrix (and hence the gluon Wigner distribution) carries an important information about the gluon saturation effects relevant at small x [60], as well as the dipole orientation effects due to its interactions with the medium. In such a dipole (also known as Light-Front) representation, the elastic γ * N → γ * N amplitude at high energies [63,150], can be viewed as a convolution (over the dipole separation r and a fraction z of the projectile photon LC momentum carried by a quark in the dipole) of three basic ingredients specified as follows. The elastic amplitude of the qq dipole scattering off the proton is denoted as A qq (x, r, ∆), and the amplitudes (also known as the light-front wave functions) of γ * → qq and qq → γ * transitions are labelled as Ψ λγ λλ and (Ψ λγ ) † λλ , for q,q helicities λ,λ, respectively (with the quark flavour index f ). The latter can be found in abundant literature, e.g. in Refs. [63,151] (see also Ref. [58]). Then, combining Eq. (4.7) with Eq. (4.8) one writes the DIS cross section (c.f. Eq. (4.5)), Such a convolution of the kinematical dependence of the LC wave function with the dipole cross section gives rise to the mean dipole size squared determined as an inverse to the quark energy squared, ϵ 2 , depending on the photon virtuality, Q 2 , momentum sharing within the dipole, z, and a constituent quark mass, m q ∼ Λ QCD . The dipole size r ∼ 1/Q is then considered to be "frozen" in the target rest frame in the course of dipole propagation through the color background field of the target. Analogically, one gets the following expression for diffractive DIS γ * N → XN cross section in the forward limit The above well-known results of the dipole model are expressed in terms of the universal dipole cross section which, by means of the optical theorem, is equal to the imaginary part of the amplitude for the elastic qq dipole scattering off the proton, i.e.
The latter is determined by the dipole S-matrix in the impact parameter representation, yielding finally the universal dipole cross section in the form where is known in the literature as the partial (elastic) dipole amplitude off a nucleon (e.g. proton) target N [62]. Then, one deduces the dipole T -matrix in the Glauber approximation for a nucleus with mass number A (A = 1 for a nucleon), For a given interaction potential of a quark with the medium V (b, z), this amplitude adopts an eikonal form, nearly imaginary at high energies [152]. In analogy to Eq. (2.22), the Fourier expansion of the partial dipole amplitude reads where the elliptic part N ϵ quantifies the dipole orientation effect leading to a non-trivial correlation in azimuthal angle ϕ br between the dipole separation r and its impact parameter b.
For small dipoles, i.e. when the dipole separation r ⊥ is small compared to the hadron scale r ⊥ ≪ R had , the universal dipole cross section σ qq realises the color transparency property vanishing quadratically with separation i.e. σ qq (r ⊥ , x) ∝ r 2 ⊥ [56,57,141,153,154]. This is an important consequence of the gauge invariance property of the fundamental QCD theory as well as nonabeliance of strong interactions. Indeed, a colorless object of vanishing size propagating through the target does not interact with its color field, i.e. the color medium is transparent for a point-like color-neutral object. For large dipole sizes, instead, the dipole cross section levels off at the so-called saturation scale Q 2 s which depends either on Bjorken x or on total energy √ŝ in the γ * p c.m. frame. To summarise this part, the universality of the dipole cross section (hence, the dipole S-matrix, S Y , or the related partial dipole amplitude, N Y ) is its most important property making it the main ingredient of a given QCD scattering process, for both nucleon and nuclear targets. Such a colour dipole framework has enjoyed a lot of developments over last few decades, both in theoretical modelling and in phenomenological approaches. Equipped with the universal dipole picture of QCD scattering, one can utilize the same fundamental ingredient, S Y (r, b), as a building block for a vast majority of high-energy scattering phenomena covering the particle production reactions in nucleon-nucleon N N , lepton-nucleon lN (such as the DIS), lepton-nucleus lA, nucleon-nucleus N A, as well as heavy nuclei AA collisions, in case of both inclusive and diffractive topologies of final states. On the other hand, an intricate connection between the gluon Wigner distribution and the dipole scattering enables us to probe factorisation mechanisms at small x which are being actively developed in the literature over past decades.

B. From impact parameter to transverse momentum representation
It is often more convenient and instructive to turn to momentum-space representation for the dipole amplitude in studies of the diffractive reactions as was implemented, for instance, in Refs. [38,39,50,138]. Let us now establish an important connection between the dipole cross section in the impact parameter representation with the intrinsic dipole TMD in the momentum representation, K dip (x, κ 2 ⊥ ). Following Refs. [141,155], one defines such that for large dipoles the cross section saturates to the non-perturbative quantity, The corresponding partial dipole amplitude N Y defined in Eqs. (4.14), (4.15) is connected to an unintegrated (off-forward) gluon density matrix [156] -an off-diagonal generalization of the unintegrated gluon distribution function (UGDF) f (Y, κ 1 , κ 2 ) -as follows: In the two-gluon approximation, each of the phase factors can be straightforwardly related to each of the four possible ways the two gluons can couple to the qq dipole, whose constituents propagate at impact parameters respectively, while the t-channel gluons carry the transverse momenta The genuine gluon GTMD (dipole T -matrix in momentum space) can then be connected to the off-forward gluon density matrix as [40] T In the perturbative regime, k ⊥ ≫ ∆ ⊥ , the δ-functions do not play any role, so up to a factor of two T Y is the same as f . The color transparency property of the dipole amplitude, As was pointed out in Ref. [40], the emergence of δ-functions in Eq. (4.23) means that the Fourier integral does not converge, and hence a regularisation procedure should be advised e.g. via introducing a Gaussian cutoff factor e −εr 2 ⊥ preserving the color transparency [38,39,50]. Then, starting from Eq. (4.17) and utilising the Fourier-Bessel transforms one obtains the regularised isotropic and elliptic components of the gluon GTMD as [36,38,39] Note that the density matrix can also be recast as Then, in the perturbative QCD regime corresponding to a hard gluon transverse momentum in the target, In the forward limit, one recovers In fact, integrating the dipole TMD over κ and considering the perturbative limit of hard Q 2 scales, one recovers in the double logarithmic approximation of the DGLAP equations [155] 1 π thus, establishing an important connection between the dipole formalism and k ⊥ -factorisation approach. Indeed, it has been known for a long time that the forward cross section of the diffractive dissociation of projectile photons (or pions) into di-jets at a perturbatively large jet p T is directly mapped on to the target UGDF (see e.g. Refs. [157,158]). Note, however, that Eqs. (4.29) and (4.31) are approximate and limited to small x and small dipoles (realised e.g. for heavy quarks or high-p T jets production) only. In the soft regime corresponding to large dipole sizes where the saturation and other non-linear QCD effects become relevant, the use of conventional UGDF and the k ⊥ -factorisation formula are not well justified, and the use of phenomenological dipole approach going beyond k ⊥ -factorisation is more sensible due to its universality. Following Ref. [155], one may formally define a two-gluon amplitude in momentum space beyond k ⊥ -factorisation -the dipole TMD -in terms of a known parameterisation of the universal dipole cross section, basically, by inverting Eq. (4.18) yielding in terms of the Bessel function of the first kind, J 0 (x), and σ ∞ qq (x) ≡ lim r ⊥ →∞ σ qq (r ⊥ , x). As was elaborated in detail in Ref. [138], in the framework of dipole approach the dipole TMD is a very convenient object to use for practical calculations directly in momentum representation.

C. Applied dipole parameterisations
A first-principle computation for the dipole cross section (and the corresponding dipole S-matrix) remains a major theoretical challenge and its exact theoretical modelling is still far from complete understanding despite of a three-decade long community effort. A remarkable progress has been made on reconstruction of kinematic and energy dependence of N Y strongly boosted by a vast amount of precision experimental data from the HERA collider. Among important QCD-based advancements, the rapidity (or energy) evolution at high energies can be understood in the CGC formalism [159,160] and is encoded in the Renormalization Group equations represented by an infinite hierarchy of the so-called Balitsky-Jalilian-Marian-Iancu-McLerran-Weigert-Leonidov-Kovner (Balitsky-JIMWLK) equations [161][162][163][164][165][166][167][168][169][170][171][172]. The latter determines how the separation of hard and soft modes evolves towards smaller momentum fractions effectively resuming terms enhanced by large ln(1/x) in leading logarithmic accuracy. In the mean-field approximation, these equations yield the Balitsky-Kovchegov (BK) equation [161,173] whose b ⊥ -dependent solutions for the partial dipole amplitude N Y (r, b) are very difficult to obtain [174]. As the b ⊥ -slope is driven essentially by non-perturbative QCD effects, one often looks for translationally invariant solutions ignoring the b ⊥ -dependence in numerical solutions and then introducing it "by hands" through an extra soft exponential factor to match the observed t-dependence of differential observables for exclusive processes.
The universality of the dipole scattering discussed above enables one to extract it from one process (typically, inclusive DIS) and then to use it in many other hadronic and photoninduced processes in p/l+p, p/l+A and AA collisions. Starting from the fundamental work of Ref. [175], many different phenomenological dipole parameterisations for the partial dipole amplitude have been advised in the literature attempting to account for the saturation phenomena, the QCD-inspired hard-scale and Bjorken x evolution, as well as for impact parameter dependence and the dipole orientation effects. For some of the widely used parameterisations relying on fits to the HERA data [145,146], see e.g. Refs. [63,66,139,150,155,[176][177][178][179][180][181][182][183][184][185].
One of the simplest and well-known phenomenological parametrizations of the saturated dipole cross section satisfying the color transparency property, represents the simplest and phenomenologically consistent way to incorporate the saturation effects and provides a very good description of a number of different observables in ep and pp collisions with both inclusive and exclusive topologies of final states at high energies (with x ≲ 0.01), as well as for processes with nuclear targets. Such a simple parametrization reminds the Glauber model of multiple interactions and is known in the literature as the Golec-Biernat-Wuestoff (GBW) model [139,186]. The GBW parameter values in Eq. (4.33) are found by a fit to the inclusive F 2 data from the HERA collider [145,146] (initially, without accounting for charm quarks) such that the model is found to simultaneously describe the diffractive DIS data (and is denoted in the literature as "GBWold" parameterization). The saturation scale related to the gluon density in the target nucleon is then conveniently given by Incorporating the charm quark contribution (with m c = 1.5 GeV) in the fit yields somewhat different parameter values: σ 0 ≃ 29 mb, x 0 ≃ 4.01 × 10 −5 and λ GBW ≃ 0.277 ("GBWnew" parameterization). It is rather surprising that such a simple phenomenologically inspired parametrization, well justified for small dipoles only, has turned out to be surprisingly successful in giving an acceptable description of HERA data [145,146] in a wide range of Q 2 ∼ [0.1 . . . 30] GeV 2 for both inclusive and diffractive DIS, as well as for a large variety of hard and semi-hard inclusive and diffractive processes in pp, pA and AA collisions (see e.g. Ref. [63] and references therein). For sufficiently hard scales µ 2 = µ 2 (r ⊥ ), it was understood in Refs. [187][188][189] that the saturation scale Q s (x) must depend on the hard scale µ 2 as a manifestation of the QCD evolution and the corresponding scaling violation. Connecting the universal dipole cross section to the leading-order (LO) DGLAP-evolved gluon density in the target, one writes [190] A particular realisation of the phenomenological dipole model, that explicitly incorporates the QCD DGLAP evolution of the collinear gluon density in the target xg(x, µ 2 (r ⊥ )) in the saturation scale, and utilises the same GBW saturated ansatz of Eq. (4.33), is known in the literature as the Bartels-Golec-Biernat-Kowalski (BGBK) model [155]. The target gluon density is found as a solution of the DGLAP evolution equation accounting for the gluon splitting function P gg (z) only, with the following parameterisation for the starting gluon distribution at µ 2 = µ 2 0 : The BGBK model parameters, were extracted through fitting F 2 to the HERA data [145,146] yielding much better fits than the original GBW model, particularly, for large Q 2 (although omitting the charm quarks contribution). A generalisation of the BGBK model incorporating the b ⊥ -dependence of the partial dipole amplitude relevant for a consistent description of the t-dependent observables in exclusive processes at HERA has been proposed in Refs. [63,150,182] and reads where B G is an additional exponential b ⊥ -slope parameter that can be extracted e.g. from the data on t-dependent differential elastic ep scattering. Hence, this so-called Impact-Parameter dependent Saturation (IP-Sat) model -also known as Kowalski-Teaney (KT) model -does not contain any non-trivial correlation in the relative angle between r and b. The IP-Sat parameters found in Ref. [185] by a fit to the precision (both, elastic and inelastic) ep HERA data [145,146] are given by A g = 2.373 , λ g = 0.052 , µ 2 0 = 1.428 GeV 2 , B G = 4.0 GeV 2 , C = 4.0 . (4.41) Turning to momentum-space representation and assuming that Q 2 s (x, µ 2 (r ⊥ )) is either independent of r ⊥ (as in GBW model) or depends on r ⊥ slowly (e.g. logarithmically as in the case of BGBK or KT models) the saturated antsatz of Eq. (4.33) in Eq. (4.32) leads to the dipole TMD in a nearly-Gaussian form and hence to the diagonal "dipole UGDF" F(x, k 2 ⊥ ) via Eq. (4.29). Then, accounting for the diffractive slope the off-forward gluon density matrix defined in Eq. (4.28) can be found as explicitly incorporating the exponential dependence on ∆ 2 ⊥ ≡ −t (with the slope B G = 4.0 GeV 2 ) in a factorised form. The latter has been used in the analysis of diffractive vector meson production in Ref. [191] (see also Refs. [69,192]).
Imposing translational invariance and, hence, ignoring the b ⊥ -dependence of the partial dipole amplitude N Y (r, b) → N (r, x), one often looks for the corresponding solution of the BK equation describing its evolution in rapidity Y ≡ ln(x 0 /x) as formulated in Refs. [193][194][195].
where r 2 = r − r 1 . Including the Renormalization Group running of the strong coupling, the BK kernel reads [195] K(r, r 1 , in terms of the number of active quark flavours N f and an adjustable parameter C. In this "rcBK" formulation, the initial r ⊥ -shape of the dipole amplitude at the starting value x = 0.01 is chosen to take a McLerran-Venugopalan (MV) [159,196] form, with Λ QCD = 0.241 MeV, Q 2 s0 = 0.165 GeV 2 , γ = 1.135 and C = 2.52 extracted from the fit yielding σ 0 = 32.895 mb and assuming that α s (r 2 ⊥ ) freezes at r ⊥ > r 0⊥ , determined by α s (r 2 0⊥ ) = 0.7 [195]. A further variation of the dipole model can be obtained utilising the collinearly improved kernel [197] given by, (4.47) Here, the parameter A 1 is fixed to 11/12, while the positive sign corresponds to r ⊥ < min(r 1⊥ , r 2⊥ ), and the definitionᾱ s ≡ α s (min(r 2 ⊥ , r 2 1⊥ , r 2 2⊥ )) Nc π is introduced. Then, variable number of flavours scheme is applied [195] such that Λ QCD becomes flavour-dependent and is computed using Similarly to the rcBK model, one employs the MV initial condition which in the collinear formulation of Ref. [197] takes the form with the definitionᾱ sat ≡ Nc π α sat ,ᾱ s (r 2 ⊥ ) ≡ Nc π α s (r 2 ⊥ ), where α sat can be fixed to unity. Fitting the normalisation of the corresponding dipole cross section σ qq (r ⊥ , x) = σ 0 N (r ⊥ , x) to the HERA data leads to the following parametrization, also known as "colBK": σ 0 = 31.4055 mb, Q 2 s0 = 0.4 GeV 2 , C = 2.586 and p = 0.807 [197]. Iancu, Itakura and Munier in Ref. [176] suggested yet another parametrization (denoted as "IIM", in what follows), in terms of an effective anomalous dimension γ eff (r ⊥ , x), the saturation scale Q 2 s (x) in the standard parametrization and the normalisation of the dipole cross section σ 0 = 2πR 2 p . Above, a continuity of the dipole cross section at r ⊥ Q s (x) = 2 is ensured by a specific choice of the parameters, while the remaining parameters were fitted to the HERA data in Ref. [183] yielding the values: γ = 0.7376, λ = 0.2197, x 0 = 0.1632 × 10 −4 , κ = 9.9, N 0 = 0.7, and σ 0 = 27.33 mb. A modification of the IIM model to incorporate the b ⊥ -dependence explicitly has been achieved in Refs. [63,182] through a modified saturation scale (4.53) The corresponding parametrization adjusted to fit the HERA data in Ref. [66] For light hadron scattering at high energies off a given target driven by the soft-scale QCD dynamics (i.e. when no hard scale is involved) the total gluon-target collision energy squared in the c.m. frameŝ ≡ x 1 s is a more appropriate variable than Bjorken x, such that σ qq (r ⊥ , x) → σ qq (r ⊥ ,ŝ) is implied as was suggested by Ref. [129]. This is the case for instance of photoproduction such that Q 2 → Λ 2 QCD where even for moderate or low energies one gets to very small values of x. The latter indeed becomes inappropriate for soft reactions at high energies, and the saturation scale becomes a function of ratherŝ than x. Other examples of such processes are the pion-proton scattering, the diffractive Drell-Yan and the gluon radiation processes. Retaining the saturated ansatz of Eq. (4.33) and following the definition of Ref. [69,129], the parameterization for the dipole cross section known in the literature as the Kopeliovich-Schafer-Tarasov (KST) model is given by in terms of the pion-proton total cross section σ πp tot (ŝ) = 23.6(ŝ/s 0 ) 0.08 mb [198], with s 0 = 1000 GeV 2 and the mean pion charge radius squared ⟨r 2 ch ⟩ π = 0.44 fm 2 [199]. It is worth noting here that despite of its soft origin this parametrization provides an acceptable description of the pion-proton scattering cross section for scales reaching as high as Q 2 ∼ 20 GeV 2 .   In order to illustrate typical uncertainties associated with a spread between different phenomenological parametrizations for the dipole cross section outlined above, in Fig. 5 adopted from Ref. [200] we show σ qq (r ⊥ , x) for each model as a function of the dipole separation for two distinct Bjorken-x values: x = 10 −2 (left) and x = 10 −4 (right). One notices big variations in both the magnitudes and shapes of the dipole cross section, particularly, for large dipoles in the saturation domain, as well as for very small dipole sizes in the perturbative (color transparency) domain r ⊥ ≲ 0.05 ÷ 0.06 fm. The corresponding theoretical uncertainties in the soft domain grow with decreasing Bjorken x. So the uncertainties in the treatment of saturation effects are rather large even before incorporating angular dependence in the corresponding gluon GTMD. Other elaborate dipole parameterisations will be briefly reviewed below in connection to the dipole orientation effects.

D. Elastic dipole amplitude and dipole orientation
As will be discussed in more detail below, one of the most important examples of phenomenological probes for the elliptic gluon density (i.e. the processes that are sensitive to it) is the fully differential cross section of the exclusive di-jet photoproduction [36]. Indeed, the angular dependence of it originates from a non-trivial dependence of the partial dipole amplitude on the dipole orientation. As was mentioned above, the qq dipole (with net color charge zero) interacts with the target only due to a non-zero relative difference (separation) r between the quark and anti-quark impact parameters relative to the scattering "center of gravity". As the "center of gravity" of the dipole is located at the impact parameter b, the interaction strength encoded in the partial dipole amplitude is expected to disappear for r ⊥ b, while it gets maximized for r ∥ b. This is one of the generic model-independent properties of the dipole S-matrix. Below, we wish to summarise some of the basic theoretical and phenomenological models for dipole orientation effects employed in the literature and discuss their basic features and phenomenological consequences.
As was pointed out recently in Ref. [40], in the impact-parameter representation a nontrivial correlation of the dipole amplitude in angle between r and b emerges even if the corresponding GTMD is fully isotropic (i.e. independent of the angle between k and ∆). This can be illustrated starting, for instance, from the isotropic gluon density matrix in the form of Eq. (4.43) and substituting it into Eq. (4.20). This way, one recovers the dipole amplitude as follows [40]: where we notice an emergence of the dipole orientation effect ∝ cos ϕ br controlled by σ 0 (x) being the dipole cross section in the non-perturbative limit of large dipoles given in Eq. (4.19), and Such a correlation has a straightforward geometric interpretation as coming due to contributions where only q orq interact with the target at impact parameters b q or bq given in Eq. (4.21), respectively. The corresponding elliptic contribution to N naturally vanishes in the perturbative limit of small r ⊥ ≪ Λ QCD as the gluon density in the target can be considered constant (homogeneous) over distances ∼ r ⊥ probed by such a very small dipole. Any genuine correlation between k and ∆ in the gluon GTMD representing the dynamical elliptic density in the target may, in principle, survive in the hard regime and comes on top of the above effect of the final-state (anti)quark wave-functions. Opportunities for disentangling both contributions unambiguously and hence for accessing the genuinely dynamic elliptic Wigner (or GTMD) gluon density experimentally constitutes yet unsolved problem being a subject of intense studies in the literature.
At Born level, considering for simplicity a symmetric dipole where quark and antiquark carry the same longitudinal momenta (i.e. having the same distances from the dipole "center of gravity"), the partial elastic amplitude of qq dipole interacting with a target quark reads [69] (see also Ref. [62]) with the effective gluon mass m g accounting for non-perturbative (e.g. confinement) effects, and the modified Bessel function K 0 (x). Indeed, as the soft gluon exchanges are dominant in the process, the gluon mass parameter encodes the exponential gluon fields' decay at large impact parameters effectively mimicking the confinement phenomenon. In Ref. [69], m g has been fixed to the pion mass enabling to reproduce typical hadronic cross sections. As anticipated above, the terms in the last line of Eq. (4.57) cancel each other for b · r = 0, thus, yielding a non-trivial angular correlation between r and b, already in the simplistic Born approximation. One should also note here that, similarly to the example above, only the first source of dipole orientation effects -due to a combination of the phase factors representing the asymptotic plane-wave states of the dipole constituents -emerges at Bornlevel, and no elliptic GTMD contribution is generated at the leading order. One therefore concludes that an angular correlation between the k and b vectors of a probed target gluon should be a higher-order (or non-perturbative) effect. Following the footsteps of Ref. [56], one can perform an analogical calculation and derive the partial dipole amplitude off the nucleon target N in the Born two-gluon exchange approximation and for a symmetric dipole resulting in |Ψ N ⟩ are the nucleon and two-quark nucleon form factors, respectively, represented in terms of the three-quark nucleon wave function Ψ N (ρ 1 , ρ 2 , ρ 3 ).
As was extensively discussed in the previous subsection, precision measurements of the proton structure function at the HERA collider at high energies (small x) have enabled to reconstruct the dipole cross section σ qq (r ⊥ , x) in wide kinematic ranges of r ⊥ and x. The Born approximation illustrated above is, of course, too naive as it does not reproduce the experimentally well established Bjorken x (or energy) dependence of σ qq (r ⊥ , x), particularly manifested at high energies. The experimentally observed rise of σ qq (r ⊥ , x) at x ≪ 1 is due to a steep evolution of the unintegrated gluon density F(x, κ 2 ⊥ ) in the target nucleon as follows from their relation given by Eq. (4.18) in the leading-logarithmic approximation, or equivalently, In order to recover a correct x-dependence of the dipole partial amplitude in the twogluon exchange model, in Ref. [69] a generalised off-diagonal unintegrated gluon density F(x, q, q ′ ) has been introduced, such that in terms of the longitudinal light-cone momentum fractions of anti-quark 1 − β and quark β, such that the dipole "center of gravity" gets effectively shifted towards the fastest q orq, andᾱ s ≡ α s (q 2 ⊥ )α s (q ′ 2 ⊥ ). The Born approximation of Eq. (4.58) is then recovered for Thus, Eq. (4.60) advises us that the dipole S-matrix should be sensitive to dipole separation not only in the transverse plane, but also in longitudinal direction due to, in general, unequal sharing of energy between q orq in the dipole [69,129]. More specifically, we notice from the structure of Eq. (4.60) that the dipole orientation effect in the partial dipole amplitude appears to be sensitive to both the angle between the transverse momenta of the exchanged gluons q and q ′ and the longitudinal momentum sharing between q orq quantified by the β variable.
Due to a remarkable simplicity of the GBW parameterisation (4.33), it enables to analytically derive both the form of the corresponding off-diagonal unintegrated gluon density [69], 62) and the phenomenological ansatz for the related partial (elastic) dipole amplitude satisfying Eq. (4.59) in the two-gluon exchange approximation (4.60), The latter contains a non-trivial dependence on azimuthal angle θ between r and b in the transverse plane. Thus, it manifestly accounts for not only the relative dipole orientation effect, but also an explicit dependence on the q andq light-cone fractional momenta. Here, are the elastic slope and its part R 2 N originated from the Pomeron-proton form factor F p IP (k 2 ⊥ ) = exp(−k 2 ⊥ R 2 N /2) and determined in terms of the proton mean charge radius squared, ⟨r 2 ch ⟩ (the current value is ⟨r 2 ch ⟩ ≃ 0.875 fm [201]), respectively. Note that the angular correlation in the phenomenological ansatz for the unintegrated gluon density (4.62) takes the following simple form, which can, in principle, be used to extract the elliptic gluon density by mapping the above expression to the ansatz of Eq. (2.22). Starting from the dipole elastic amplitude off the nucleon target f N qq , we can derive the elastic cross section of dipole-N scattering -the main building block of the dipole picture of the QCD scattering. Ignoring the real part which is typically small at high energies, one writes [69] dσ (qq)N el The latter relation can be used to determine the slope of the differential elastic cross section which has been measured at HERA collider in diffractive ρ-meson electroproduction process at large Q 2 to be [202] The latter is in rough agreement with its relation to the proton-proton elastic slope, B (qq)N el (r ⊥ ) ≈ B pp el /2. One may straightforwardly generalise the phenomenological formalism for the dipole scattering off a nucleon outlined above to the case of nuclear target, T ≡ A, with atomic mass number A. Correspondingly, the partial dipole amplitude for a qq dipole scattering off A at impact parameter b ⊥ is found in the Glauber picture as in terms of the effective nuclear thickness [203] T Here, the standard nuclear thickness function T A is found as an integral of the conventional nuclear density ρ A (b ⊥ , z) taken along the trajectory of the particle. Provided that the mean value ⟨s 2 ⊥ ⟩ ≈ 0.4 fm 2 found from Eq. (4.67) is much smaller than the nuclear radius squared (for a heavy nucleus), i.e. ⟨s 2 ⊥ ⟩ ≪ R 2 A , one may expand and apply the latter in Eq. (4.68) giving rise to the approximate formula for the partial dipole amplitude off a nucleus target, where We notice that for a large nucleus, being the density at the center of the nucleus) are small for b ⊥ ≪ R and are peaked at the nuclear periphery. Besides, in the eikonal approximation applicable for heavy nuclei one may use T A (r, b) ≈ T A (b ⊥ ) such that the angular correlations are naturally small.
The formula (4.63) (in both GBW and KST parameterisations as well as its generalisation to nuclear targets in Eq. (4.71)) represents the first known phenomenologically driven model for the 5D partial dipole amplitude that has been fit to the experimental data and that can be connected to the dipole S-matrix and hence to the gluon Wigner distribution using Eqs. (4.15) and (2.20), respectively. The θ-averaged part of the Wigner distribution is directly related to the isotropic unintegrated gluon density in the target F, the basic nonperturbative ingredient of k ⊥ -factorization. The difference between the Wigner distribution and its θ-averaged part is proportional to the elliptic gluon distribution that can be extracted from the exponential parametrization of the dipole amplitude in Eq. (4.63), or in more elaborate models discussed below.
In the following, we wish to illustrate how the azimuthal angle correlations of the gluon GTMD can be understood at small x in the framework of QCD evolution based models, also leading to emergent effects of saturation and elliptic flow.

V. QCD EVOLUTION BASED MODELS FOR DIPOLE SCATTERING
Both angular-dependent and angular-independent components of the Wigner distribution are highly sensitive to the soft and saturation physics, and thus are a subject of intense research in the literature [50,70]. Here, we would like to overview the widely used models for the Wigner distribution (dipole S-matrix) based upon QCD evolution.

A. Azimuthal angle correlation in the BFKL model
In the BFKL approximation valid at x ≪ 1 (or high energies), the T -matrix of a dipole scattering off another dipole at transverse position x in the impact parameter space reads [204][205][206] (see also Ref. [36]) valid for a symmetric energy sharing between the quark and antiquark in the dipole. Here, the holomorphic eigenfunction of the BFKL operator and the corresponding characteristic function read respectively, with h = 1−n 2 + iν andh = 1+n 2 + iν. At x ≪ 1, the dominant contribution comes from the n = 0 term. In order to disentangle a non-trivial correlation in azimuthal angle ϕ br between r and b, one can apply the saddle point approximation in Eq. (5.1) in a vicinity of ν = 0 and consider the limit of x ⊥ ≪ r ⊥ , b ⊥ , yielding the result for the dipole T -matrix [36] T in terms of Thus, one observes that, in consistency with the discussion in the previous section and with Ref. [69], the partial dipole amplitude is enhanced for a dipole configuration with parallel b and r (i.e. ϕ br = 0) compared to the case of b ⊥ r. Analogically, beyond the saddle point approximation, i.e. close to the nonlinear saturation domain where the partial dipole amplitude reaches a constant value T (r, b, Y )| r ⊥ →1/Qs → const at the saturation momentum scale Q s , one finds Resolving the latter with respect to Q s at large impact parameters (small momentum transfers) corresponding to b ⊥ ≫ r ⊥ ≃ 1/Q s , one gets [36] in overall agreement with a numerical analysis of the nonlinear small-x dynamics of Refs. [174,207]. Note that a non-trivial angular correlation is observed also for configurations with r ⊥ ∼ b ⊥ and hence is a generic phenomenon that should be inherent for realistic saturation models. Switching to momentum-space amplitude by means of a Fourier transform T Y F = ⇒ T Y over r and b and averaging the result over all possible target dipole x orientations, one obtains [36] T .

(5.8)
Thus, the angular correlation is manifested also in momentum space, particularly, for even harmonics, starting from a sizable elliptic ∼ cos 2ϕ part, while the contributions of the odd harmonics vanish. Furthermore, in the linear (saddle-point) approximation and for finite momentum transfers ∆, one may write the relation between the partial dipole amplitude and the Wigner function in momentum space following Ref. [36] as

B. Elliptic gluon distribution in the Color Glass Condensate
Let us turn now to a discussion of the gluon Wigner distribution in a high-energy nucleon or nucleus in the saturation regime at small-x accounting for the gluon saturation phenomenon. The latter can be incorporated in the dipole S-matrix formalism through small-x QCD evolution described by the Balitsky-Kovchegov (BK) equation with impact parameter dependence [161,173], providing a rapidity evolution for S Y in the leading-log accuracy and large-N c limit. Here, x ≡ b + r 2 and y ≡ b − r 2 are the transverse positions of a quark and an antiquark in the dipole, respectively.
A computation of both the isotropic and elliptic components of the gluon Wigner distribution (and the corresponding GTMD) through numerical solution of the BK equation (5.10) has been performed by Hagiwara, Hatta, and Ueda in Ref. [50] (we denote it as HHU model, in what follows). Here, a simplifying assumption that the considered BK solutions are invariant under an SO(3) subgroup of the maximal symmetry of the BK equation, the conformal Möbius group in the transverse plane, has been applied [208]. Namely, the BK solutions were sought in the following SO(3)-symmetric form, and with symmetric initial condition, in terms of an arbitrary length-scale parameter R. The above ansatz for the considered solution relies on a conjecture that the BK equation may dynamically restore SO(3) symmetry at the level of solutions despite its breaking by realistic initial conditions as has been observed numerically at small or large r ⊥ in Ref. [174] (see also Ref. [209]). Note, the ansatz (5.11) is symmetric with respect to the change of variables, r ⊥ → r 2 m /r ⊥ , with r m = 2 b 2 ⊥ + R 2 ≃ 2b ⊥ for b ⊥ ≫ R, thus, recovering the observation of Ref. [174] that the partial dipole scattering amplitude T Y is maximised for the dipole separation of r ⊥ = 2b ⊥ .
With these considerations, the BK equation (5.10) can be recast into the following SO(3)symmetric form [50], with d 2 given in Eq. (5.11). The numerical results of Ref. [50] for S Y (d 2 ) and T Y (r ⊥ , b ⊥ = 1, cos(ϕ b − ϕ r ) = 0) are shown in Fig. 6 in left and right panels, respectively, fixing the units by setting R = 1, as well as taking N c α s /π = 0.2. As expected above, the peak in T Y appears at r ⊥ = r m at any value of Y .
As was further pointed out in Ref. [50], there is an artefact of the BK equation whose kernel maintains the perturbative Coulomb interaction at large separations which is a nonphysical behavior due to confinement. In real QCD, however, the effect of confinement is expected to lead to the black disc limit corresponding to T Y (r ⊥ ≫ R) → 1 such that the contribution to the Fourier integral in the Wigner distribution (2.20) from the perturbative tail at large r ⊥ ≫ r m domain should vanish and cannot affect the observables. A natural way to eliminate non-physical contributions is to use the Husimi distribution (2.5) where the integral over r would be cut off at large distances by means of a Gaussian factor. A similar dumping factor has been introduced in Ref. [50] also in the Wigner distribution (2.20) such that where a natural choice of ϵ = 1/4 has been made, thus, effectively cutting off the nonphysical contribution from large r ⊥ ≳ 2. As follows from Eq. (5.11), the Fourier transform of W ′ g contains only even harmonics such that the first ∝ cos 2ϕ (elliptic) term provides the leading angular dependence [36]. Isolating the isotropic and elliptic contributions in Eq. (5.15), one finds and shown as functions of k ⊥ (at fixed b ⊥ = 1) for different rapidities Y = [0, . . . , 10] in Fig. 7, left and right panels, respectively, and ϕ br is as usual an angle between b and r vectors. We notice that the position of the peak in xW ′ g,0 found at the saturation momentum k peak = Q s (Y, b ⊥ ) grows with rapidity Y and falls with impact parameter b ⊥ . The latter behavior is in consistency with the geometric scaling, T Y (r, b) ∝ (r ⊥ Q s ) 2γ for small r ⊥ . Here, γ is weakly dependent on r ⊥ and is found to be between γ = 1 at r ⊥ → 0 and γ ≃ 0.63 at r ⊥ ≲ 1/Q s [210,211]. The peak value in the elliptic part xW ′ g,1 is only a few percent of that of xW ′ g,0 , and evolves much slower with Y since there is no geometric scaling in the elliptic contribution [50].

C. Azimuthal correlations in the McLerran-Venugopalan model
As was mentioned earlier, the color dipole orientation with respect to its impact parameter b probes the transverse inhomogeneity of the gluon distribution in the target. A robust semianalytical computation of the partial dipole amplitude incorporating the dipole orientation effects has been performed in Ref. [70]. For this purpose, in the latter work an effective   McLerran-Venugopalan (MV) model [159] has been employed and extended to incorporate inhomogeneities in the transverse-plane distribution of soft gluons in the target. Such an extension procedure is in overall consistency with the phenomenological dipole models that feature a non-trivial impact-parameter dependence and saturation such as the bCGC [63,66,176] and IP-Sat [150,185,212] parameterisations being successfully matched to the HERA data as was discussed above (see also Refs. [65,207,213,214]). However, they do not contain a non-trivial azimuthal angle dependence essentially assuming the rotational invariance of the dipole S-matrix.
Being equipped with the MV model embedding an inhomogeneous Gaussian-like valence quark distribution of the color-field sources in the impact parameter b ⊥ , in Ref. [70] Iancu and Rezaeian have computed explicitly the angular dependence of the partial amplitude for dipole scattering off such sources in the target (this model will then be called MV-IR in what follows). This has been worked out for both smooth and dense targets such as nucleon and for a lumpy target like a nucleus. An important observation here is that the elliptic gluon density (and hence the azimuthal asymmetry) has been shown to originate as a genuinely non-perturbative phenomenon significantly impacted by soft or semi-hard gluon multiple scattering effects at characteristic momentum scales of order of the scale of inhomogeneities in the target (i.e. the saturation scale). By focusing on valence quark distributions as a starting point, the MV model relies on a semi-classical pQCD-motivated approach that may not reflect the full picture of the soft QCD phenomena. It is believed, however, that it, at least, captures their most pronounced characteristics signifying the relevance of the realistic b ⊥ -dependence in the dipole scattering process. Important physical parameters of the MV-IR model [70] are the effective gluon mass setting up the gluon confinement scale, the width of the Gaussian source distribution in the transverse plane as well as the saturation scale characterising the transverse inhomogeneity of the target. A significant model dependence and variations among predictions (e.g. for v 2 ) existing in the literature typically come from different modelling of these intrinsically non-perturbative ingredients.
It has been demonstrated in Ref. [70] that the inhomogeneity in the transverse plane is largest at the edge of the target, hence, causing a large v 2 (p ⊥ ) (peaked at p ⊥ ∼ Q s ) for peripheral collisions, while multiple scatterings strongly impact v 2 causing a change in its sign. Let us elaborate on this phenomenon in more detail as it would help us to shed light on the common origin of elliptic flow in nucleon and nuclear collisions.
In the original MV model [159] formulated for the case of very large (hence, homogeneous) target nucleus with atomic mass number A, the nuclear gluon density is assumed to be saturated in the high-energy limit, such that its small-x QCD evolution is ignored. In this regime (which can, in principle, be considered as an initial condition for the high-energy QCD evolution), one views the target as a collection of 'valence' color-field sources ρ a (b) whose two-point correlation function reads in terms of the uniform color-charge squared of the valence quarks per unit area µ(b) → µ 0 = const within a disk of radius R ∝ A 1/3 . As was advocated in Ref. [70], inhomogeneities in the transverse color-charge distribution are responsible for azimuthal asymmetries in particle production. To this end, a phenomenologically inspired impact-parameter dependence of the isotropic color-charge distribution µ(b) ≡ µ(b ⊥ ) has been introduced in the MV-IR model in a simple Gaussian form suggested by successful fits to HERA data on diffractive production processes as well as on vector meson production in DIS of Refs. [66,150,185,212]. Such a generalised MV-IR model has been formulated for both the proton and nuclear targets as an attempt of a unified description of flow-type effects in events with high-multiplicity final states in pp and heavy-ion collisions. While µ 0 scale characterises the saturation momentum squared at b ⊥ = 0, in Ref. [70] the Gaussian width R ≃ 0.3 fm (setting up the soft scale 1/R ∼ Λ QCD ) and the semi-hard saturation scale squared, (here, C F = (N 2 − 1)/2N for SU (N )) have been taken from the phenomenological fits to the HERA data in the framework of IP-Sat model [150,185]. Employing the eikonal approximation, in the MV model it is assumed that the color dipole scatters off the sources of the color field independently. For a single such scattering, the forward amplitude corresponds to an exchange of two gluons in the t-channel and reads, in the impact parameter representation, where γ(x, y) = γ(y, x) represents a partial elastic amplitude, with one of the two gluons coupled to the quark at transverse position x and another one -to the antiquark at position y. The latter can be found in terms of the two-point correlation function of the two-dimensional Coulomb-like color fields A − a (x) of the ultrarelativistic color charges, 22) in terms of the two-dimensional effective Coulomb propagator, where the modified Bessel function of the second kind K 0 exponentially falls at large separations m g b ⊥ ≫ 1. One should notice here that in the case of m g → 0 the small-x gluon field of the color source would extend far beyond the typical size of the source itself due to a power-like decay of the Coulomb propagator, in immediate contradiction to confinement. Hence, the presence of the effective gluons mass m g ∼ Λ QCD manifestly introduces confinement into the picture which is an important physical property of the gluon propagator even though the amplitude N 2g (x, y) is well-behaved in the infrared limit in perturbative m g → 0 regime of the theory. Using (5.18) and (5.22), one obtains , (5.24) in terms of the Fourier transform of the color-charge distribution,μ(q) ≡μ(q ⊥ ), withμ(0) = 4πR 2 µ 0 , such that the area of valence charge distribution of the proton is 4πR 2 and hence the transverse proton size is 2R. Introducing the average transverse momentum transfer between the quark and the target, k = (q ′ − q)/2, and a difference between the momentum transfers between the amplitude and its conjugate characterising the inhomogeneity of the target, ∆ = q ′ + q, one finally arrives at [70] N 2g (r, b) = g 2 C F 2 with the standard dipole variables -the dipole impact parameter with respect to the target centre, b = (x + y)/2, and the dipole separation, r = x − y. As the integral above is controlled by predominantly soft (or semi-hard) gluon exchanges, ∆ ⊥ ≲ 1/R, one could apply the eikonal approximation where multiple gluon exchanges exponentiate yielding the Glauber approximation for the dipole S-matrix, S = exp(−N 2g ). For small dipoles, r ⊥ ≪ b ⊥ , corresponding to the perturbative limit with the hard scale Q 2 ∼ 1/r 2 ⊥ ≫ Q s (b ⊥ ) concerned, the leading-order single-scattering approximation S ≃ 1 − N 2g can be used to a good accuracy. Being an even function of both r and b, N 2g and hence S are functions of the azimuthal angle θ only through powers of (b·r) 2 . Hence, the flow parameters vanish (v n = 0) for any odd n. This can be readily seen in the single-scattering approximation, where the Fourier-transformed amplitude reads for an exchanged single hard gluon with momentum (having the same momentum as the final-state quark) p ⊥ ≡ k ⊥ , hence resolving a localized valence substructure in the target.
In this regime, we are not sensitive to a gluon mass. The angular correlation emerges at a sub-leading 1/p 6 ⊥ order such that the cross section reaches a maximum at π/2. As a result, the produced quark is more likely to move perpendicular to its impact parameter vector leading to a negative leading (elliptic) flow coefficient, 27) while the higher harmonics appear in the higher-order terms in 1/p ⊥ -expansion and hence are suppressed relative to v 2 . As expected, the correlation for very small dipoles r ⊥ ≪ R and for central collisions b ⊥ → 0 should vanish as the dipole orientation makes no impact on the production cross section in those cases. Besides, the cos(2ϕ) term in the partial dipole amplitude (5.26), and hence the v 2 coefficient, are proportional to derivative of the valence color-charge distribution µ(b ⊥ ) over the impact parameter b ⊥ . This illustration provides a clear picture of azimuthal correlations that can be relevant only for soft processes and in the case of non-central impacts, and are generated essentially due to the presence of inhomogeneities in the transverse gluon distribution in the target [70]. For softer gluon exchanges, i.e. for p ⊥ ≲ Q s (b ⊥ ), particularly relevant for understanding the collective phenomena in soft particle production in pp and pA collisions, the singlescattering approximation is no longer valid, and the Glauber-like exponentiated dipole Smatrix S(b, r) = exp (−N 2g (b, r)) →S(b, p) should be considered. As the multiple gluon exchange is concerned, the momentum p acquired by final-state quark in the course of its propagation through the target can be much larger than the momentum k exchanged in a single scattering. Moreover, as has been emphasised in Ref. [70], the angular correlations are controlled by soft momentum exchanges, i.e. with 1/r ⊥ ≫ k ⊥ ∼ ∆ ⊥ , so one can no longer expand the amplitude in power series over ∆ ⊥ /k ⊥ . Instead, expanding (5.25) in the limit of small r ⊥ and keeping the leading order terms only, where spatial indices j, l = 1, 2 are summed, the dependence on the gluon mass m g becomes relevant restricting the phase space allowed for soft gluons, k ⊥ ≲ Λ QCD ∼ m g , and typically |N ϵ | ≪ |N 0 |. As a result, the elliptic flow coefficient (3.7) takes the form, . (5.29) where I 0 and I 1 are the modified Bessel functions of the first kind. Due to the inhomogeneities in the target, the quantity ∆ ⊥ effectively plays a role of an infrared cutoff in the k ⊥ -integration, and no infrared divergencies are developed. In the UV limit, a logarithmic divergence is effectively cut off by the dipole separation, i.e. at k ⊥ ∼ 1/r ⊥ , such that the leading term at large k ⊥ can be subtracted from the integrand in Eq. (5.28), and then the corresponding contribution to the integral can be replaced by a regularised term, (5.30) with the saturation scale defined in Eq. (5.20). Here, the logarithmic term emerges due to k ⊥ -integration over k ⊥ ⊂ [m g , 1/r ⊥ ], while the constant term fixes the renormalization scheme. As a result of the MV-IR model, the isotropic N 0 (r ⊥ , b ⊥ ) and elliptic N ϵ (r ⊥ , b ⊥ ) components of the amplitude in Eq. (5.28) for the proton target can be represented in the following integral form [70]: where the k ⊥ -integrals are restricted to k ⊥ ⊂ [m g , ∆ ⊥ ]. The first term in Eq. (5.31) corresponds to one in the original MV model [159]. Taking ∆ ⊥ -integral analytically one arrives at [39] N 0 (r ⊥ , b ⊥ ) = Q 2 0,s r 2 which are rather straightforward to use in practical calculations. The functions N 0 (r ⊥ , b ⊥ ) and N ϵ (r ⊥ , b ⊥ ) are illustrated in Fig. 8 (left and right panels, respectively). For completeness, we also highlight the MV-IR result of Ref. [70] for the nucleus target, 36) given in terms of the nuclear thickness function, is the normalised nuclear density distribution in the Woods-Saxon parametrization [215], with the nuclear radius R A = (1.12 fm)A 1/3 and δ = 0.54 fm. Notably, in the limit of a homogeneous target,μ(∆) ∝ δ (2) (∆), only the first term in N 0 (r ⊥ , b ⊥ ) remains while the elliptic component N ϵ (r ⊥ , b ⊥ ) vanishes. Hence, in the considered MV-IR model formulation the azimuthal angle correlations are indeed due to transverse inhomogeneities in the source distribution of the target. Besides, the logarithmic divergence for m g → 0 in N 0 (r ⊥ , b ⊥ ) is cancelled by the one coming from the second term yielding infrared stability of the result. Yet, m g ∼ 0.3 GeV is there on physical grounds effectively cutting off too soft gluon momenta, k ⊥ ≲ ∆ ⊥ , not consistent with confinement. As was shown in Ref. [70] starting from Eq. (5.29), multiple scattering effects turn v 2 positive and large for peripheral collisions, i.e. b ⊥ ≳ R, whereas it is negative in the single scattering approximation, see Eq. (5.27). The MV-IR model enables one to study quantitatively the dipole orientation phenomena encoded in the dipole S-matrix. Using Eqs. (5.31) and (5.32) one obtains the behavior of S(r ⊥ , b ⊥ , θ) and N 2g (r ⊥ , b ⊥ , θ) shown in Fig. 9 (left and right, respectively). One immediately notices that, as expected, the dipole survival probability is largest for b ⊥ r and smallest for b||r, while the target inhomogeneities become irrelevant and the dipole orientation effect disappears for small b ⊥ ≲ R = 0.3 fm.
For illustration, Ref. [40] shows the isotropic part of the gluon GTMD in the forward kinematics, T 0 (k ⊥ , ∆ ⊥ = 0.01 GeV), for a fixed x ≡ x IP , obtained using Eq. (4.26) for different selected models for the partial dipole amplitude N briefly reviewed above. Such a comparison reveals a significant theoretical uncertainty in modelling of the gluon GTMD in the soft regime of k ⊥ ≲ m c GeV, even in its isotropic part. Ref. [40] further illustrates how such uncertainties in the gluon GTMD modelling are translated into the uncertainties in determination of the differential cross sections of the exclusive diffractive cc photoproduction. Theoretical uncertainties in modelling the elliptic component of the Wigner distribution (or GTMD) are found to be even larger as was discussed extensively in Refs. [38][39][40] posing a natural question of constraining it directly from the experimental data. Fig. 10 illustrates the forward dipole T -matrix components -the isotropic T 0 (upper panels) and elliptic T ϵ (lower panels) ones -computed in the MV-IR model using Eqs. (4.26) and (4.27), respectively, as functions of k ⊥ for different values of ∆ ⊥ and for different targets -the proton target (left panels) and lead nucleus target (right panels) [39]. The elliptic effect in nuclear targets is much more pronounced than in the proton in general. Nevertheless, it get quite significant at k ⊥ ∼ ∆ ⊥ even for the proton case, while the ballpark of sensitivity to variations in small ∆ ⊥ occurs for k ⊥ ≲ 2 GeV. To conclude this part, one typically implies that the effects of hydrodynamic evolution become increasingly important and even dominate for large collision systems (especially, at large particle multiplicities), while at low particle multiplicities, often probed in collisions of small systems, initial-state effects are expected to dominate. This is, for example, illustrated in Fig. 11 summarizing measurements of the particle multiplicity dependence of v 2 using 2, 4, 6 and 8-particle cumulants, and v 3 and v 4 using 2-particle cumulants for pp, pPb, PbPb and XeXe collisions. Comparisons of the data with prediction of the hydrodynamic model MUSIC using IP-Glasma for initial conditions and UrQMD for hadronic rescattering [76,216] may indicate that the overall description could be improved once measurements of the gluon Wigner function become available.

VI. PHENOMENOLOGICAL PROBES FOR DIPOLE ORIENTATION EFFECTS
Here, we would like to elaborate on several phenomenologically important examples of the processes particularly relevant for exploration of the nucleon or nucleus structure through measurement of the azimuthal angle correlations in the particle production processes.  Ref. [218].

A. Exclusive light-quark dijet photoproduction
A multitude of new experimental data on various exclusive diffractive processes in hadronic and nuclear collisions are now becoming available from the LHC and making use of these data for detailed studies of the Wigner distribution gains its importance. One of the important classes of diffractive processes that can be utilized for probing non-trivial correlations in phase space parton distributions is the meson-pair photoproduction. For instance, in Ref. [219] it has been shown that the azimuthal cos 4ϕ-asymmetry in exclusive di-pion production close to the ρ 0 -resonance in UPCs recently observed by STAR Collaboration can receive sizeable contributions from the elliptic gluon Wigner distribution. Another important class of processes for probing the gluon tomography concerns diffractive electroor photoproduction of light-and heavy-quark jets which will be elaborated upon in detail below.
Typically, an extraction of GTMD directly from the data is a rather challenging task due to the fact that a typical cross section, such as that of a semi-inclusive DIS process or an exclusive vector meson production, is computed in terms of an integral convolution of the GTMD with other kinematical functions and an inversion of such an integral is, in general, a very difficult problem. As a way out, one could consider other possible reactions where such an inversion can be realised in a simpler way. A recent proposal to reconstruct both the isotropic and elliptic components of the gluon GTMDs (and hence the gluon Wigner distribution) at low-x through a measurement of exclusive diffractive dijet photoproduction process γT → (qq)T illustrated in Fig. 3, where the target T can be either proton p or a nucleus A, has been made in Ref. [36,38] (see also a recent work of Ref. [220] for more details). In the same way as in the DIS process discussed above, the emerging colorless qq state, the color dipole, is a lowest-Fock fluctuation of the projectile photon γ → qq that plays an important role of the probe of the color background field (and hence substructure) of the target T . For comprehensive analysis of this process in the framework of color dipole approach, see e.g. Refs. [37,41,221] and references therein.
Such a process can be realised e.g. in ep and eA collisions as well in ultra peripheral pA and AA collisions (UPCs) where an electron, a proton or a heavy projectile nucleus of mass number A is considered as a source of quasi-real photons while another proton or nucleus plays a role of the target whose internal structure is being probed. As will be discussed below, this type of processes can be efficiently studied in UPCs at the LHC, as well as at RHIC and a future EIC providing plausible opportunities for probing the elliptic component theoretically expected to be at the level of only a few percent effect in the perturbative regime. For instance, the EIC [9,14,82] is expected to enable proton tomography with unprecedented precision through reconstructing partons' positions in the impact parameter space, their transverse momenta and spins in the proton and nucleus. This may be achieved via precision measurements of diffractive dijet production observables as functions of target's recoil momentum and dijet's total and relative transverse momenta.
The colliding states in UPCs such as relativistic nucleons and/or nuclei scatter at large impact parameters interacting only electromagnetically through an exchange of quasi-real (and hence mostly transversely polarised) WW photons [222,223]. Let us consider this process in the nucleus-target center-of-mass frame where the WW photon is collinear to the beam axis, with energy ω. In the light-cone decomposition, its momentum reads: q = ( √ 2ω, 0, 0). One then introduces the quark and antiquark x 1,2 momentum fractions of the light-front plus momentum of the projectile nucleus, respectively, written in terms of their transverse momenta k 1,2⊥ and center-of-mass rapidities y 1,2 as while the quark and antiquark light-cone momentum components read where z is the quark longitudinal momentum fraction of the parent photon in the γ → qq transition. Neglecting the quark mass, it reads z = k 1⊥ e y 1 k 1⊥ e y 1 + k 2⊥ e y 2 . (6.2) In the quark massless limit, momentum conservation then provides the standard expressions for the invariant mass squared of the produced dijet and the photon energy, respectively, in terms of the azimuthal angle ϕ 12 between k 1 and k 2 . Note, M 2 jj plays a role of the hard scale of the process enabling the use of QCD perturbation theory for sufficiently large values of quark transverse momenta. The kinematics of the multi-gluon exchange in the t-channel is determined in terms of its longitudinal momentum fraction x taken by such an exchange from the target as defined in Eq. (2.15), or equivalently, x = M 2 jj /(x γ s). At high energies and in the forward scattering domain, one requires a large rapidity gap Y ≡ log(x 0 /x) ≫ 1 (with x 0 = 0.01 as often adopted in the literature) stretched between the recoiled target and the produced dijet.
The WW photon carries the fraction x γ ≡ x 1 + x 2 whose flux can be written as in terms of the atomic number (total charge) of the projectile particle Z, the fine structure constant α em , the radius of the target (projectile) R T (R A ), respectively, ξ is defined in such a way to effectively eliminate an overlap between T and A in the transverse plane, and m p is the proton mass. The proton and nuclear radii are typically fixed as R p = 0.8 fm and R A = (1.12 fm)A 1/3 , respectively. Consequently, in the c.m. frame, the differential cross sections of qq production processes in AT UPCs and γT collisions are related as, One of the important advantages of using a large nucleus as a source of WW photons is that a suppression due to the electromagnetic coupling α em can be largely eliminated by means of a large prefactor Z 2 . Moreover, in the target rest frame, the projectile hadron or nucleus appears as a source of WW photons whose spectrum is rather broad, and the peak photon energy scales linearly with the Lorentz γ-factor in the target rest frame. This process has been studied earlier at the next-to-leading order in QCD in Ref. [224] where it has been shown that the corresponding dijet differential cross section shows a large extension in photon-proton invariant mass and is very sensitive to small momentum fractions of the exchanged gluons.
The shear advantage of pA and AA UPCs is not only a strong reduction of QCD backgrounds and a cleaner hadronic environment due to photon-induced reactions and a large relative separation between the scattered systems, but also an enhancement of the WW photon flux from the projectile nucleus A due to its scaling as the square of the nuclei electric charge, Z 2 . The measurements in UPCs do not feature complications related to a significant pileup and enable a good subtraction of non-exclusive backgrounds as previously elaborated in e.g. Refs. [225,226]. The LHC forward proton detectors (FPDs) such as Roman pots in TOTEM [227], ALFA [228] and AFP [229,230] at ATLAS, and CT-PPS [231] in the CMS provide the necessary means for precision analysis of exclusive processes in pA and AA UPCs enabling full kinematic reconstruction by detecting the forward final-state protons and possibly nuclei, in addition to reconstruction of the diffractive dijet system.
A proper extraction of the "elliptic" gluon Wigner distribution (or its GTMD counterpart S(q ⊥ , ∆ ⊥ )) requires a detailed reconstruction of kinematics of each final-state jet. Even this may not be sufficient due to a potentially large momentum leakage into the rapidity gap between the dijet and the final-state nucleon (or nucleus). So an extra vetoing of any hadronic activity between the dijet and the forward nucleon at the detector level would enable to suppress the QCD backgrounds related to the Pomeron and/or the nucleon breakup. In order to ensure the full exclusivity of the process and the maximal reduction of the dominant backgrounds, it necessary to also reconstruct the full momentum of the forward nucleon by simultaneously measuring its energy loss and its full transverse momentum. This is in principle available by using forward proton detectors as will be discussed below. In this case, momentum conservation helps veto any additional hadronic activity in the rapidity gap and, hence, reliably probe the Pomeron-exchange contribution relevant for the Wigner distribution measurement. It may be relevant to explore a large class of inclusive and exclusive processes sensitive to the Wigner distribution (or GTMD) to make its extraction more constrained and reliable, and we will further discuss several promising examples below.
As was mentioned above, in order to reconstruct the gluon GTMD (and hence the corresponding dipole S-matrix), one has to identify a full set of observables sensitive to both q and ∆. In the considered example of exclusive dijet photoproduction in the forward region, one aims at precision measurement of the jet transverse momenta, k 1 and k 2 , enabling to reconstruct the total transverse momentum transfer with the target, ∆ = −(k 1 + k 2 ). This momentum is transferred across the large rapidity gap in elastic qq-dipole scattering off the target. The same ∆ is carried by the target in the final state which remains intact and can be, in principle, independently measured in the forward detector. The latter would substantially reduce a possible unobservable momentum leakage, both into the rapidity gap (due to a Pomeron break-up) and off the jets (parton energy loss due to a momentum leakage out of the jet cone), thus ensuring the exclusivity of the process. The full exclusivity then also requires that the projectile nucleus (source of the WW photons) also remains intact and whose kinematics is measured to precisely reconstruct the WW spectrum. The second independent kinematic variable is the relative transverse momentum between jets, P = 1 2 (k 1 −k 2 ), whose measurement provides an information about the q-dependence of the GTMD.
Focusing for simplicity on the transversely polarised projectile WW photon and con-sidering the direct photon (not resolved) process only, the differential cross section of the exclusive di-jet production in γT -scattering reads [36], Here, the quark energy squared ϵ 2 q = z(1 − z)Q 2 + m 2 q ≈ 0 in the case of photoproduction of light quarks. Turning to the dipole representation [56,57], it is instructive to write down the γT → (qq)T cross section in terms of the impact parameter space ingredients, namely, the imaginary part of the partial dipole amplitude, N Y (r, b) defined in Eq. 4.15, and the light-front wave function of the γ → qq fluctuation, Ψ λγ λλ (z, r) [58]. The cross section is explicitly averaged over the photon polarisation λ γ and summed over the quark λ/2 and antiquarkλ/2 helicities. Using Eq. (4.20), one gets [40] A λγ λλ qq where the impact factor describes the coupling of two t-channel gluons to the γ → qq amplitude and vanishes when κ = ±∆/2. In Ref. [36] an extraction of the gluon GTMD has been proposed through a measurement of light quark jets neglecting the quark mass m q and considering Q 2 → 0 (photoproduction limit). In this case, the momentum integral in Eq. (6.6) gets dominant contributions for q ∼ P and can be analytically inverted [38]. Indeed, taking into account that the dominant angular correlation in the Wigner distribution is encoded in its elliptic part as in Eq. (2.22) (see also Refs. [36,50]) and performing the angular integral in the amplitude analytically, in the massless quark limit one arrives at [38] dσ AT →A(qq)T where the azimuthal angles of P and ∆ vectors are ϕ P and ϕ ∆ , respectively, summation over the light quark flavours q is made explicit, and the following two structure functions 10) (6.11) are introduced that contain the basic information about the transverse-momentum dependent isotropic and elliptic gluon densities in the target, respectively. Note that in Eq. (6.9) only the linear terms for S Y,ϵ ≪ S Y,0 are kept. A direct reconstruction of A and B structure functions through a precision measurement of the fully differential cross section (6.9) would then enable one to extract the isotropic and elliptic components of the gluon GTMD as follows [38] S Y,0 (P ⊥ , It is worth noticing here that, when considering non-zero virtualities of the projectile quasireal photon, one should also add a contribution from the longitudinally polarized virtual photons [232,233]. The latter may somewhat impact the extraction of the Wigner distribution discussed above. With the two commonly used models for the gluon Wigner distribution -BK and MV-IR models introduced above -the predictions for the structure functions A and B have been obtained in Ref. [38]. In particular, Fig. 12 illustrates the results for rapidity dependence of A (left) and B (right) found as functions of P ⊥ using a numerical solution of the impactparameter dependent BK equation [161,173] for the proton target following Ref. [50] (here, ∆ ⊥ = 1 and the scale parameter R = 0.4 fm). The transverse momenta P ⊥ and ∆ ⊥ are defined in units of 1/R. The function A(P ⊥ ) exhibits a peak at the saturation scale Q s (Y ) which grows with Y following the geometric scaling, while the peak in B(P ⊥ ) moves much slower with Y . Considering a possibility to probe the Wigner distribution in the nucleus in AA → A(jj)A UPCs, Ref. [38] also shows the predictions for A and B for a large nucleus in the target (such as lead with atomic mass number A = 208). For this purpose, the MV-IR model for the nuclear gluon distribution has been utilized taking into account its inhomogeneities in the transverse plane as formulated in Ref. [70] -see Figs. 13 and 14 for smaller ∆ ⊥ ≤ 0.5 and larger ∆ ⊥ ≥ 0.7, respectively. The corresponding analytical expressions for the isotropic S Y,0 (r ⊥ , b ⊥ ) and elliptic S Y,ϵ (r ⊥ , b ⊥ ) densities in the nuclear targets in the small-x limit (not incorporating an explicit x-dependence) are found in terms of the nuclear thickness function T A (b) and its derivatives in Ref. [70] (see also Ref. [81]).
We notice that the peaks of A and B in the case large nuclear targets for small ∆ ⊥ < 0.3 exceed those for the proton target evaluated in the BK model by one-two orders of magnitude. While the nuclear A drops very quickly with ∆ ⊥ and eventually falls below that of the proton, the MV-IR model predicts the nuclear B being of an order of magnitude larger than that of the proton at ∆ ⊥ ≥ 0.7 where both A and B structure functions feature a node in their P ⊥ dependence. This observation of Ref. [38] signifies that a measurement of the exclusive dijet production in AA UPCs has a better chance to constrain the elliptic Wigner distribution.

HERA measurements
Dijets produced diffractively in photon-induced processes have been measured rather thoroughly at HERA. From comparing data with theory predictions based on the QCD factorization and hence on the universality of diffractive proton PDFs extracted from previous HERA measurements, one usually judges the level of the QCD factorization breaking (the farther from unity, the bigger the breaking is). This factorization breaking is usually explained as a consequence of additional interactions between incoming photon and proton. The H1 Collaboration published results of three analyses, making use of data with large rapidity gap [234,235] as well as with a very forward proton spectrometer [236], the latter allowing one to control better the non-diffractive background and contributions from proton dissociation. The ZEUS Collaboration published one analysis relying on large rapidity gaps [237]. While all three H1 analyses report the above ratio to be around 0.6, the ZEUS experiment measures it to be close to unity and hence a different level of the factorization breaking, if at all, is claimed by the two collaborations. We conclude that this is rather a small tension given the fact that all three H1 analyses confirm each other after carefully evaluating many sources of uncertainties and using several experimental approaches. As far as the potential information about the Wigner function is concerned, we have to conclude that no cos(2ϕ) or even jet ϕ measurement was done. And since the data taking at HERA finished in 2007, chances for any supplemental future measurements of this kind are small.

LHC measurements
One of the main conclusions of Ref. [38] is that the Wigner function can be observed at rather low jet transverse momenta p T by reconstructing the full kinematic dependence of the structure functions A(P ⊥ , ∆ ⊥ ) and B(P ⊥ , ∆ ⊥ ) and then inverting them using Eqs. (6.13). As explained above, for LHC experiments the usual p T threshold for jets based on calorimeter information is 20 GeV. As was evident from Figs. 12-14 above, theoretical predictions show a non-trivial dynamics (e.g. the appearance of the node in P ⊥ dependence of the structure functions, and a growing relative importance of the elliptic term) in the low ∆ ⊥ and P ⊥ < 4 GeV domains while they are monotonously falling at larger transverse momenta such that their features are very difficult to probe at the LHC. To get to lower hard scales via jet production, it was suggested to use track-based jets at the LHC or to use data from RHIC. Low-p T jets can also be probed at EIC as will be explained in the section VII.
A very useful measurement of the cos(2ϕ) distribution (where ϕ is the azimuthal angle between the vector sum and the vector difference of momenta of two highest E T jets) in PbPb collisions has been performed by CMS [238]. It is based on central detector only, so no nucleus tagging is used and exclusivity requirements are reduced to requiring empty regions in the tracker and/or calorimeter, depending on the noise of individual subdetectors. The transverse momenta of the two jets are required to be above 30 GeV for the first jet and 20 GeV for the second jet and at the same time it is required that the absolute value of the difference is larger than the absolute value of the sum. As Fig. 2 (right) of Ref. [238] documents the slow increasing trend for the mean value of cos(2ϕ) as a function of ∆ ⊥ observed in the data is described decently by calculations [239] based on perturbative soft QCD gluon radiation in the final state up to ∆ ⊥ < 15 GeV, although Wigner function effects are mainly of non-perturbative origin. This is probably not unexpected given the relatively large hard scale given by the jet p T and the absence of forward detectors, more specifically Roman Pots (RPs).
Indeed, the advantage of RPs is manifold: they can be used to i) tag nuclei, ii) measure their fractional energy loss, ξ, iii) time-of-flight (ToF) and iv) transverse momenta. The first three would allow one to use more stringent requirements on exclusivity than those based on rapidity gaps used in Ref. [238], by requiring that the value of an observable measured in the central detector is identical to the value measured by FPDs within resolutions of subdetectors used in the above measurements. In this way, quantities based on ξ (dijet mass or rapidity) obtained in the central detector are required to match those measured in Silicon trackers of FPDs, while z vtx of a primary vertex obtained in the central tracker is required to match that obtained by the ToF detector in FPDs. Finally, the measurement of the nucleus transverse momentum would then be significantly more precise than the measurement of the sum of the two jet momenta.

B. Diffractive photoproduction of heavy quarks
As was elaborated above, the elliptic component of the the gluon Wigner distribution becomes relatively more pronounced compared to its isotropic component when the transverse momenta of final-state quark and antiquark approach the saturation scale of the gluon density in the target. As there is no other hard scale in the γ + (gg) → qq matrix element for light quarks, with gg in the color singlet state, except for the (anti)quark transverse momentum (being translated into the jet-p T ), pushing into the soft domain poses natural questions about the validity of perturbation theory. Besides, the existing experimental capabilities at ATLAS are limited by the detector acceptance imposing with a lower cut-off on jet at 20 GeV or so. These constraints, together with substantial backgrounds, severely restrict one's capability to probe the elliptic density in the kinematics domain where it is typically enhanced.
As a way out of this problem, over the past years it has been suggested in Refs. [39,40] (see also Ref. [41]) that exclusive photoproduction of heavy (Q = c, b) quark pairs in highenergy diffractive eT (with the target T = p, A) collisions as well as in pT and AA UPCs offers a powerful tool for probing the gluon Wigner distribution in the target T in both perturbative and soft kinematics domains. This is due to the fact that the quark mass m Q provides a naturally hard scale (instead of a large P ⊥ in the case of light-quark jets). The latter enables us to go effectively to much lower transverse momenta for final-scale heavy-quark jets (i.e. a jet that contains an open heavy flavour meson reconstructed in a detector) and, hence, to probe the gluon GTMD at transverse momenta P ⊥ ≳ ∆ ⊥ down to a GeV scale, or even lower. The analysis of Refs. [39,40] was inspired by an earlier work of Refs. [36,38] on light-quark jets discussed above in the previous subsection, and below we outline some of its basic details.
Besides providing a natural hard scale for the perturbation theory to more reliably work at low quark transverse momenta, this process is also advantageous from the experimental point of view as the reconstruction of open heavy-flavor states (such as B and D mesons) can be achieved at lower transverse momenta than those for light-quark jets. Finally, despite of a smaller cross section, this process is expected to benefit from a better control of QCD backgrounds. However, there is a price to pay for these gross advantages: the integral in the diffractive amplitude is not analytically invertible for a large m Q with respect to the components of the Wigner distribution as was shown above for massless quarks. Thus, so far in the literature only predictions for the differential observables were given, for a set of distinct theoretical models, and no efficient reconstruction procedure for the elliptic density directly from the data has been advised yet in this case.
The formalism of exclusive heavy quarks' photoproduction in A + T (T = p, A) UPCs is very similar to that of light quarks, upon an appropriate introduction of the finite quark mass terms [39,40]. The differential cross section for the diffractive γT → (QQ)T process can be represented as where the matrix elements can be represented in terms of the gluon GTMD in the target, T (Y, k, ∆), absorbing all the ∆-dependence as follows [40] In this representation, the last k-independent term in both integrands effectively drops out upon the use of the "sum rule" (4.25). Then, one ends up with the "formfactor representation" for the matrix elements [39] leading to a generalization of Eq. (6.9): where the incident WW photon flux is defined in Eq. (6.4). Note that two additional formfactors C and D emerge in the heavy quark case whose expressions are given explicitly in Ref. [39]. The impact of the additional C and D formfactors due to the finite-mass corrections on the angular correlation of the differential cross section is illustrated in the case of exclusive cc-pair photoproduction in Fig. 15. Here, the ratio of the angular cosine-weighted average cross section defined as to the angle-integrated one is shown as a function of P ⊥ at fixed ∆ ⊥ = 0.2 GeV (upper panels) and of ∆ ⊥ at fixed P ⊥ = 10 GeV (lower panels). Here, MV-IR model discussed above has been used. Such an observable appears to be a good probe for the angular correlation since, indeed, it tends to more positive the more collinear (parallel or antiparallel) P and ∆ become (while negative case corresponds to the vectors turning perpendicular). Naturally, the quark mass effect diminishes at larger transverse momenta but becomes quite important for softer P ⊥ and in the forward limit ∆ ⊥ ∼ Λ QCD . For more details on the quark mass effects, see Ref. [39], while a comparison of exclusive photoproduction cross section of charm quarks computed for different models for the dipole T -matrix can be found in Ref. [40]. In the literature, there are also other relevant exclusive diffractive processes involving heavy quarks that are sensitive to the elliptic component of the gluon Wigner distribution. An earlier complementary analysis of exclusive charm production in ep collisions at the EIC (at large Q 2 ) has been performed in Ref. [41] in the framework of CGC. Here, the gluon Wigner distribution has been modelled by solving the JIMWLK equations for Wilson lines using the initial conditions constrained through the HERA data. Besides, the quarkonia pair production process and its connection to the gluon Wigner distribution has been studied in Ref. [240].

C. Photon bremsstrahlung
A particularly clean and important phenomenological source of new information on the nucleon and nucleus structure is direct (or prompt) production of photons. Such a reaction is typical, for instance, in DIS processes well explored at the HERA collider. As photons do not take part in strong interactions, they transmit unaffected information about the hard scattering in pp or pA collisions, and can be measured in the final state, either real or virtual, through a detection of a Drell-Yan lepton pair.
In the parton model, the direct photon emission can be thought of as a Compton scat-tering gq → γq. On the other hand, in the target rest frame the same process should be considered as photon bremsstrahlung off the projectile quark scattered off the target T (nucleon, T = N , or nucleus, T = A). Its amplitude is represented by an interference of two simple diagrams at tree-level as is shown in Fig. 16. The corresponding color-dipole description in the case of direct production of real photons has been thoroughly discussed e.g. in Refs. [69,143,[241][242][243], while the process of virtual (massive) photon production with dilepton (Drell-Yan) final states have been analysed both in pp and pA collisions e.g. in Refs. [131,144,[244][245][246][247][248][249][250]. As was discussed above, the elliptic gluon Wigner distribution in 16. Photon radiation by a projectile quark in the target rest frame. From Ref. [132].
the target encodes a non-trivial dependence on the dipole orientation with respect to the color background field of the target. The latter can be probed by measuring the azimuthal angular correlation of the final-state photon and a jet or a leading hadron originating from the recoil quark after hadronisation. Here, we briefly elaborate on the main results and implications of this process for phenomenology of the elliptic gluon distribution.
In order to highlight the relevance of this process for gluon tomography of the nucleon, we follow the analysis of Ref. [69] based on phenomenological dipole framework. Following the discussion of Ref. [143], the factorized dipole formula for the transverse momentum distribution of the direct (unpolarised) photons radiated by the projectile quark at a certain impact distance b from the centre of the target T reads α em 2π 2 m 2 q α 2 α 2 K 0 (αm q r 1⊥ )K 0 (αm q r 2⊥ ) (6.18) in terms of the constituent quark mass m q ≈ 0.2 GeV that can be viewed as an adjustable IR cutoff parameter [129,251], the final-state photon transverse momentum p, its LC momentum fraction α ≡ p + γ /p + q , the Bjorken x variable of soft gluons in the target, the LC wave function of the lowest-Fock q → qγ transition, ϕ γq (α, r), at a transverse separation r, and the effective dipole amplitude off a target T , F T (x, αr 1 , αr 2 , b). At high energies, the latter can be represented as a superposition of partial dipole amplitudes accounting for dipole-target scatterings of colorless dipoles with different separations at impact parameter b ⊥ , namely, As was elaborated above, in order to quantify the azimuthal asymmetry in the direct photon production one could track the correlation in azimuthal angle ϕ p of the transverse momentum p with respect to a fixed impact parameter vector b represented by the elliptic flow coefficient, where the averaging can be performed analytically. For this purpose, it is convenient to switch from the integration of the differential cross section (6.18) over direction of vector p at fixed b to the integration over direction of vector b at fixed p, such that ϕ p → ϕ b . This is a helpful step due to a particularly simple (Gaussian) form of the phenomenological b ⊥ -dependence in the effective amplitude F T stemming from the dipole parameterisations in Eqs. (4.60) and (4.71) for the nucleon and nucleus targets, respectively. Now, we have all the ingredients to write down the azimuthal asymmetry of direct photons emitted in quark-target collisions as [69] v where In particular, for a heavy nuclear target T ≡ A, using Eq. (4.71) one obtains [69] Imf A qq (x, r ⊥ , b ⊥ ) = This result explicitly demonstrates that the azimuthal asymmetry is non-zero only if a nontrivial correlation between the transverse vectors b and r is present in the partial dipole amplitude f N qq (r, b) (otherwise, g(r ⊥ ) = h(r ⊥ ) = 0). Besides, the asymmetry is strongly enhanced at the periphery of the nucleus, b ⊥ ∼ R A , as expected. Fig. 17 illustrates the sensitivity of the azimuthal angular correlations of the prompt photon production in quark-nucleon qN (left panel) and quark-nucleus qA (right panel) collisions represented by the corresponding elliptic flow parameter, v qN 2 (p ⊥ , b ⊥ , α), to the impact parameter b ⊥ and the photon transverse momentum p ⊥ for a fixed α = 1 [69].
Here, a few cases with different center-of-mass collision energies (at RHIC and LHC) and impact parameters b ⊥ are shown. One notices an increase of the azimuthal anisotropy with the rise of b ⊥ while it vanishes with p ⊥ . As can be anticipated on general grounds, the azimuthal anisotropy reaches maximal values at the nucleon and nucleus peripheries where the sensitivity to the dipole orientation effects is largest, while the latter diminishes for smaller dipoles. It was demonstrated in Ref. [69] that the azimuthal asymmetry of direct photons is strongly sensitive to the behavior of the nuclear thickness function at the nuclear edge. When the quark propagates through the nucleus at small values of b ⊥ , it effectively interacts with nucleons placed at different azimuthal angles with respect to the quark path such that their contributions largely cancel each other in v qA 2 making it smaller compared to v qN 2 . This cancellation is not exact, hence the azimuthal symmetry is broken, due to a non-trivial b ⊥ -dependence of T A .

Exclusive and semi-exclusive mesons and photons at HERA
The quasi-elastic production of vector mesons, production of photons in the Deeply Virtual Compton Scattering (DVCS) or prompt photon production belong to the simplest processes at HERA. The vector meson production was thoroughly measured in DIS as well as in inclusive or diffractive photoproduction. Thanks to the wealth of data and the level of details, the measurements form a legacy. Properties of vector mesons of the photoproduction origin are mapped in Refs. [252][253][254][255][256][257][258][259][260][261][262][263][264][265]. The vector meson species including ρ, ω, ϕ, ρ ′ , J/Ψ, Ψ ′ and Υ had been detected via mainly two-prong decays by reconstructing charged decay products and requiring no further activity beyond the noise levels in the relevant subdetectors. Due to this requirement of exclusivity, we have |t − t min | = p 2 T , where the vector meson transverse momentum p T is obtained from its decay products and t min is the kinematically allowed minimum which is usually negligible. Thanks to the high-precision trackers, the t variable is measured precisely without the need to tag the intact proton. Other variables whose distributions were usually measured are Q 2 and W . The Φ angles have not been measured, so the elliptic part can not be constrained but all the HERA results can be used in global analyses to extract the radial part (i.e. TMD).
The two processes with the photon in the final state are attractive due to the fact that hadronization effects do not need to be considered and thanks to a high precision of the energy measurement of an isolated object in the calorimeter. This, however, is balanced by a necessity to demand a good level of isolation whose aim is to separate the signal photon from background photons (e.g. from decays of π 0 and η 0 ).
Photoproduction of the prompt photon in association with a jet was a subject of numerous HERA analyses, see Refs. [266][267][268][269][270][271][272]. Besides the distributions of E γ T , η γ and ξ γ of the prompt photon shown in all the above papers, in Refs. [266,269] there are distributions of the difference between the azimuthal angle of the jet and prompt photon.

Inclusive γ+jet production
Globally speaking the inclusive production of prompt photons in association with jets in proton-proton collisions is measured routinely and with a high precision. Precise measurements are useful for many reasons. They provide a testing ground for perturbative QCD (for example, the t-channel quark exchange) and serve to tune MC models or validate generators used to estimate backgrounds in searches for BSM physics signals which involve photons. The results include also measurements of angular correlations between the photon and the jet which, on one hand, can be used to probe the dynamics of the hard-scattering process but on the other, have a potential to give crucial information about the elliptic part of the WF, as we stress throughout this review. To separate the prompt-photon in a final state with accompanying hadrons one has to require a tight isolation of the photon to avoid the large contribution from neutral-hadron decays into photons. The production of isolated photons in association with jets in pp collisions at 7, 8 and 13 TeV was studied by the ATLAS [273][274][275][276] and CMS [277][278][279][280] Collaborations. This process is also extremely useful in the procedure of determining jet energy scales at hadron colliders (see detailed descriptions by ATLAS [281] and CMS [282]). Once the reconstructed jets are corrected to the particle level, the so called in-situ calibration needs to be done whose purpose is to account for differences between the jet response in data and simulations. The jet response is evaluated by balancing the p T of a jet against that of a well-calibrated reference object, the Z 0 -boson or a photon.

D. Deeply virtual Compton scattering
In the previous subsection we discussed the prompt photon production as a bremsstrahlung process off the projectile quark in the dipole picture where the angular correlation between the produced photon and the leading hadron or jet (emerging from the recoil final-state quark) has been shown to probe the elliptic Wigner contribution. In lepton-nucleon (and lepton-nucleus) colliders, the real photon can be produced either in the elecromagnetic bremsstrahlung (Bethe-Heitler) process, which dominates in fixed-target experiments (such as COMPASS) for not-so-hard photon virtualities, or through the deeply-virtual DVCS process [283,284].
Let us discuss how the elliptic gluon Wigner distribution impacts the DVCS process providing an important complementary way to analyse the inner structure of the proton, both theoretically and experimentally, being of particular relevance for lepton-nucleon colliders (see e.g. Refs. [9,14,[285][286][287][288]). Among particularly important implications of DVCS, it enables to access the information about orbital angular momentum of partons inside a nucleon [289][290][291][292]. A thorough evaluation of the DVCS cross section accounting for the cos 2ϕ term in the dipole picture at small x realised in the framework of CGC formalism [59] has been performed in detail in Ref. [61], and here we briefly outline the basic results. The differential DVCS cross section in ep-scattering takes the following standard form: in terms of the standard leptonic and hadronic tensors, L µν and M µν , respectively. The DVCS kinematics in the dipole approach is analogical to that of exclusive dijet production discussed above (and depicted in Fig. 18). In particular, denoting the lepton momenta in the initial and final states as l and l ′ , while keeping the same notation for initial and final proton momenta, p and p ′ = p + ∆, respectively, one finds the momentum of the collinear initial virtual photon, q = l − l ′ , q 2 = −Q 2 . The standard DIS variables then read, x Bj = Q 2 /(2q · p), y = q · p/(l · p), t = −∆ 2 ⊥ and W 2 = (q + p) 2 ≈ Q 2 /x Bj , such that the lepton transverse momentum, l ⊥ = l ′ ⊥ , satisfies l 2 ⊥ = 1−y y 2 Q 2 . In the dipole factorization approach valid at small x, the DVCS amplitude, illustrated in Fig. 18 in impact-parameter (left) and momentum (right) representations, can be found as [62,63,293,294] A DVCS ∼ d 2 be ib·∆ dzd 2 rΨ γ * (z, r)Ψ * γ (z, r)T Y (r, b) , (6.24) in terms of the dipole T -matrix T Y ≡ 1 − S Y discussed above valid at x ≪ 1. In this process, x is related to the Bjorken variable as x ≃ x Bj /2 being the same as the skewness ξ = (p + − p ′+ )/(p + + p ′+ ). Indeed, at high energies, or small x, the dipole amplitude is dominated by its imaginary part which to the leading order is found in terms of GPDs at ξ = x. Besides, Ψ γ * (z, r) and Ψ * γ (z, r) are the light-cone wave functions for the virtual (T, L-polarised) photon in the initial state and for real (T -polarised) photon in the final state, respectively (found e.g. in Refs. [61,63,151]).
In the considered light-front kinematics, the gluon GPDs, in particular, are defined as such that the conventional unpolarised gluon PDF is recovered in the forward limit as xH g (x, ∆ ⊥ → 0) ≡ xG(x), while E T g stands for the gluon transversity (or helicity-flip) GPD normalised as in Ref. [295]. Then, using the relation of the above matrix element to the dipole S-matrix found in Ref. [36] in the form (c.f. Eq. (2.19)) and utilising the gluon GTMD S Y decomposition in Eq. (2.22), one arrives at the following important relations connecting the gluon GPDs to the isotropic and elliptic components of the small-x gluon Wigner distribution as follows [61] xH g (x, ∆ ⊥ ) = 2N c α s d 2 q q 2 ⊥ S 0,Y (q ⊥ , ∆ ⊥ ) , (6.26) such that the elliptic gluon Wigner distribution is responsible for the helicity-flip contribution to the DVCS amplitude. Hence, the azimuthal cos 2ϕ-correlation in DVCS induced by the helicity-flip gluon GPDs as initially predicted in Refs. [295][296][297][298][299] can be considered as an important probe for the CGC dynamics and the gluon tomography at small x. The full result for the differential DVCS cross section in the dipole formulation reads [61] dσ(ep → e ′ γp ′ ) dx Bj dQ 2 d 2 ∆ ⊥ = α 3 em πx Bj Q 2 1 − y + y 2 2 (A 2 0 + A 2 2 ) + 2(1 − y)A 0 A 2 cos(2ϕ ∆l ) where the correlation in the azimuthal angle ϕ ∆l of the final state photon with respect to the lepton plane is concerned. Here are the UV-finite leading (in the collinear Q 2 ≫ q 2 ⊥ limit) contributions to the helicity-conserving and helicity-flip amplitudes for T -polarised photons, as well as the finite amplitude for L-polarised photons where ϵ 2 q = z(1−z)Q 2 in the massless quark limit. In fact, as was pointed out in Ref. [61] the latter amplitude describing γ * L → γ transition proportional to cos ϕ is often neglected in the dipole model calculations and vanishes unless one incorporates the correct phase factor e −iδ·r [62]. The results of the CGC dipole-model analysis have demonstrated the full consistency with those of the collinear factorisation framework at Q 2 ≫ q 2 ⊥ , in particular, for the cos ϕ and cos 2ϕ correlations [286,287]. Due to a sensitivity to the transverse-momentum gluon distribution in the target at small x, the dipole approach becomes more appropriate for lower Q 2 .
Hence, a possible reconstruction of A 2 (∆ ⊥ ) and A L (∆ ⊥ ) from the measurement of cos ϕ ∆l and cos 2ϕ ∆l angular correlations in the fully-differential DVCS cross section would potentially enable us to extract the elliptic gluon Wigner distribution S ϵ,Y . Since the latter can also source the cos 2ϕ correlation in inclusive hadron production, hence, giving rise to the elliptic flow in high energy pp and pA collisions [71,73], a combination of both measurements would certainly help understand the physical origin of the elliptic flow and the role of hydrodynamic evolution in inclusive hadron production in pA collisions as well as better constrain the nucleon and nucleus structure (gluon tomography) at high energies [81]. The DVCS analysis can straightforwardly be generalised to the diffractive vector meson (such as J/ψ, ρ and ϕ) production in DIS [14,63,150,300,301]. For more recent works along these lines, see e.g. Refs. [32-34, 200, 302-309] and references therein.
As discussed in Ref. [310], HERA data are sensitive to GPDs at low x. With respect to the vector meson production, the advantage of DVCS is the absence of the wave function but the disadvantage is a lower production cross section (due to the coupling of the final-state photon). Published results from HERA about DVCS [311][312][313] show a good precision of the Q 2 , W and t measurements. In Ref. [312] the azimuthal Φ angle between the plane of incoming and outgoing lepton and the plane of virtual and real photon was measured (the incoming lepton is the beam electron/positron and the virtual/real photon is the incoming/outgoing photon).

VII. FUTURE MEASUREMENTS
In this section, we briefly discuss proposals for future measurements, with a potential to probe the gluon Wigner function.
Very large integrated luminosities assumed to be achievable at future hadron colliders promise to give high-statistics data samples that would allow us to map as many as all five kinematical dependencies in the Wigner distribution with a relatively fine binning. We stress the importance of measuring the ϕ-angle of the objects in the pair required to be present in the final state, respectively their difference, if we want to access the elliptic part of the Wigner distribution. Thanks to precise trackers present in large experiments, the precision of the ϕ-angle measurement is usually rather high. The t-variable can also be potentially measured using the central tracker if the final-state objects are individual tracks (and not jets) and if the tracker stays more-or-less empty. This is impossible to reach at future hadron colliders that are designed to run at high instantaneous luminosities and hence with large amounts of pile-up events that lead to overwhelming combinatorial backgrounds.
There are two approaches to cope with the pile-up problem: i) to take data using the so-called special runs, and ii) to utilize forward proton detectors equipped with the so-called time-of-flight (ToF) detectors and to try to impose exclusivity criteria. The advantage of the first approach is low pile-up rates, while the disadvantage is low statistics of the collected data. The asset of the second approach (besides the measurement of the t-variable directly from forward protons, see Section VI A 2) is an access to very-high-statistics samples but we pay for that by a limited acceptance of forward detectors and a limited efficiency of the exclusivity criteria including a finite resolution of the ToF detector. The requirement of exclusivity practically means requiring a matching between measurements of the same variable by the forward proton detectors and by the central detectors, within resolutions of the considered detectors.
At the HL-LHC with an assumed pile-up rate of 200 events per bunch crossing and assumed integrated luminosity of 4 ab −1 (or 140 pile-up events and 3 ab −1 ), a 10 ps ToF resolution with a fine granularity (of the order of 1 mm) seems to be imperative. The latest Run 2 results indicate about 25-30 ps resolution achieved with a 3-5 mm granularity of quartz bars, obtained at relatively low pile-up rates, but it is believed that the design characteristics above are conceivable if proper R&D studies are performed, especially with respect to performance and life-time of photo-multipliers expected to be exposed to extremely high radiation doses during the nominal HL-LHC running.
In the case of Central Diffraction process of interest (all (semi)exclusive processes considered in this review), the dominant part of the combinatorial background comes from a sum of three independent events: one hard-scale event and two soft SD events, each giving a proton on one side of the forward proton detector. The virtue of the ToF detector lies in the capability to suppress this combinatorial background. If our signal is selected using the double-tag signature (tagging nucleon/nucleus on each side of the IP), combining the time information from both ToF detectors enables us to constrain the longitudinal position of the IP, z vtx . From a direct comparison of such measurements with a very precise z vtx measurement from the central tracker (the latter of the order of 100 µm), we can suppress the combinatorial background coming from vertices more distant from the primary vertex than the resolution of the ToF method. It is thus evident that the better resolution of the ToF detector, the better suppression of the combinatorial backgrounds can be expected (for more performance studies, see e.g. Ref. [314]. The role of ToF detector granularity and µ-dependence is studied in Ref. [315]). The impact of the ToF resolution on significance to extract various CD processes from the backgrounds can be found e.g. in Refs. [316,317].
At very high levels of pile-up rate, for example at the HL-LHC, more handles will most likely be necessary. One possibility is to make use of timing detectors in the central part of the main detector. Such detectors (HGTD [318] in ATLAS, MTD [319] in CMS) have been approved as upgrades for the HL-LHC period and are currently being built. With a design resolution for both of about 30 ps, one can get an additional information about timeof-flight of tracks in the central tracker and hence suppress tracks and vertices belonging to pile-up. Another handle is the so-called track veto (or z-vertex veto). By requiring no (or very few) tracks and vertices in a very narrow region around the primary vertex, one can significantly reduce the effect of pile-up. This requirement is an alternative to requiring large rapidity gaps in an environment where the gaps are obscured by tracks from pile-up. One should, however, be aware that the efficiency of all suppressors of the combinatorial background (exclusivity (or matching including ToF of forward protons), ToF of central tracks and track veto) decrease with the increasing level of pile-up. It will therefore be advisable to use some combination of them, with the exact composition depending on the signal and background processes under study. In this context it is useful to note that not necessarily the double-tag signature should be required to get the optimum statistical significance. Requiring both forward detector sides to be tagged leads to a better control over the combinatorial background but unavoidably leads to smaller signal yields and hence sometimes to smaller significances. A single-tag requirement naturally gives more signal but leads to a more dangerous background contamination. This should not discourage us from considering this option, if we are able to model the combinatorial background properly, for example by data-driven methods such as event-mixing or machine learning. The event mixing technique has successfully been used in modeling the combinatorial background in searches for Axion-like particles with AFP based on the single-tag signature [320], or in searches for New Physics with CT-PPS based on the double-tag signature [321].
The experimental environment in detectors to be used at the future lepton colliders, such as FCC-ee [322] and EIC [14], is believed to be more friendly compared to future hadron colliders due to a much lower pile-up and underlying event contaminations. This promises to provide very "clean" jets with access to low p T and masses which, together with a sensitivity to fragmentation reduced compared to final-state hadrons, makes them an excellent probe of the gluon Wigner function properties [323]. On the other hand, the jets are going to contain significantly fewer particles than at hadron colliders, and also the fragmentation contributions are more visible at low jet masses, which may entail jetalgorithm and fragmentation-model dependencies. The ability to polarize beams at the EIC allows one to access spin-dependent observables in which jets can play a leading role. Diffractive dijets (as elaborated above) and jet-lepton correlations in DIS [324,325] probe the 3D unpolarized and polarized TMD, and serve as important roads to the nucleon tomography at the EIC.

VIII. SUMMARY AND CONCLUSIONS
In summary, in this review we have collected some of the basic theoretical and phenomenological developments and experimental considerations regarding the fundamental properties and principal possibilities of eventual measuring the 5D gluon Wigner distribution in QCD through scattering processes at high energies. The characteristic cos 2ϕ-correlations in the azimuthal angle of the produced particles, in both inclusive and (semi)exclusive reactions in high-energy pp, pA and AA collisions provide important means for probing the nucleon and nucleus structure in the transverse plane and at x ≪ 1, i.e. the gluon tomography, encoded in the elliptic Wigner distribution. This is especially pronounced in the dipole scattering picture in the target rest frame where multi-particle correlations are influenced directly by the dipole orientation relative to the color field of the target nucleon or nucleus, and can be theoretically estimated e.g. in the Color Glass Condensate framework.
While the fully differential cross sections for (semi)exclusive processes are most sensitive to the kinematical dependencies of the Wigner distribution, their measurement is very challenging as the statistics are very low, and no reliable measurement directly accessing the elliptic density is yet available. Multi-particle correlations in inclusive processes, that manifest themselves through, for instance, the elliptic flow effect, have been measured in a number of different reactions for both small (pp and pA) and large (AA) systems. However, information about the initial-state phenomena encoded in the elliptic Wigner distribution is convoluted with hydrodynamics effects impacting the hadronic final states, especially for high-multiplicity events and large collision systems. A comprehensive global analysis combining the yet-to-be-measured azimuthal-angle correlations in (semi)exclusive reactions with the elliptic flow measurements in inclusive processes in several different channels (such as DVCS, diffractive dijet and prompt photon production) would be mandatory for the reliable reconstruction of the elliptic Wigner distribution and for constraining the genuine effect of hydrodynamic evolution in the flow measurements. We believe that our review will be useful for QCD theory and phenomenology research in these interconnected areas.