Elliptic flow from color-dipole orientation in pp and pA collisions

TL;DR

This paper explores how the orientation of color dipoles in high-energy proton collisions contributes to elliptic flow, revealing non-perturbative effects and matching observed data at the LHC.

Contribution

It introduces a semi-analytic model linking dipole orientation to elliptic flow, emphasizing the role of multiple scattering and non-perturbative gluon exchanges.

Findings

Angular dependence driven by soft gluon exchanges

Multiple scattering alters the sign of v2

v2 peaks at a transverse momentum related to saturation momentum

Abstract

For ultrarelativistic proton-proton and proton-nucleus collisions, we perform an exploratory study of the contribution to the elliptic flow $v_2$ coming from the orientation of the momentum of the produced particles with respect to the reaction plane. Via the CGC factorization valid at high energies, this contribution is related to the orientation of a color dipole with respect to its impact parameter, which in turn probes the transverse inhomogeneity in the target. Using the McLerran-Venugopalan model (with impact-parameter dependence) as an effective description for the soft gluon distribution in the (proton or nuclear) target, we present a semi-analytic calculation of the dipole scattering amplitude, including its angular dependence. We find that the angular dependence is controlled by soft gluon exchanges and hence is genuinely non-perturbative. The effects of multiple scattering turn out to be essential (in particular, they change the sign of $v_2$). We find that sizable values for $v_2$, comparable to those observed in the LHC data and having a similar dependence upon the transverse momenta of the produced particles, can be easily generated via peripheral collisions. In particular, $v_2$ develops a peak at a transverse momentum which scales with the saturation momentum in the target.

Equations

$$\frac{\mathrm{d}{\sigma }^q(qA\to qX)}{\mathrm{d}\eta {\mathrm{d}}^2\mathbf{p}{\mathrm{d}}^2\mathbf{b}}=x_{p}q(x_p)\int \frac{{\mathrm{d}}^2\mathbf{r}}{(2\pi {)}^2}{\mathrm{e}}^{-\mathrm{i}\mathbf{p}\cdot \mathbf{r}}S(\mathbf{b},\mathbf{r},x_g)=\frac{1}{(2\pi {)}^2}{x}_{p}q(x_p)\tilde{S}(\mathbf{b},\mathbf{p},x_g).$$

$$x_p=\frac{p}{\sqrt{s}}{\mathrm{e}}^{\eta },x_g=\frac{p}{\sqrt{s}}{\mathrm{e}}^{-\eta },$$

$$v_n(p,b)\equiv \frac{{\int }_0^{2\pi }\mathrm{d}\varphi \cos (n\varphi )\frac{\mathrm{d}{\sigma }^q(qA\to qX)}{\mathrm{d}\eta {\mathrm{d}}^2\mathbf{p}{\mathrm{d}}^2\mathbf{b}}}{{\int }_0^{2\pi }\mathrm{d}\varphi \frac{\mathrm{d}{\sigma }^q(qA\to qX)}{\mathrm{d}\eta {\mathrm{d}}^2\mathbf{p}{\mathrm{d}}^2\mathbf{b}}}.$$

$$\frac{\mathrm{d}{\sigma }^q(qA\to qX)}{\mathrm{d}\eta {\mathrm{d}}^2\mathbf{p}{\mathrm{d}}^2\mathbf{b}}$$

$${\mathrm{e}}^{\mathrm{i}A\cos \varphi }=\sum _{n=-\infty }^{+\infty }{\mathrm{i}}^n{\mathrm{J}}_n(A){\mathrm{e}}^{\mathrm{i}n\varphi },$$

$$v_2(p,b)$$

$$v_3(p,b)$$

$$\frac{\mathrm{d}{\sigma }^q(qA\to qX)}{\mathrm{d}\eta {\mathrm{d}}^2\mathbf{p}{\mathrm{d}}^2\mathbf{B}}=\frac{1}{(2\pi {)}^2}\int {\mathrm{d}}^2\mathbf{b}{x}_{p}q(x_p,\mathbf{b}-\mathbf{B})\tilde{S}(\mathbf{b},\mathbf{p},x_g),$$

$$v_n(p,B)\equiv \frac{{\int }_0^{2\pi }\mathrm{d}\Phi \cos (n\Phi )\frac{\mathrm{d}{\sigma }^q(qA\to qX)}{\mathrm{d}\eta {\mathrm{d}}^2\mathbf{p}{\mathrm{d}}^2\mathbf{B}}}{{\int }_0^{2\pi }\mathrm{d}\Phi \frac{\mathrm{d}{\sigma }^q(qA\to qX)}{\mathrm{d}\eta {\mathrm{d}}^2\mathbf{p}{\mathrm{d}}^2\mathbf{B}}},$$

$$v_2(p,B)$$

$${\epsilon }_2(p,B)$$

$$\int \frac{{\mathrm{d}}^2\mathbf{p}}{(2\pi {)}^2}\tilde{S}(\mathbf{b},\mathbf{p})=1⟹\int p\mathrm{d}p\int r\mathrm{d}r\mathrm{d}\theta {\mathrm{J}}_0(pr)S(b,r,\theta )=2\pi .$$

$${\epsilon }_2(B)=\frac{\int b\mathrm{d}b\mathrm{d}\alpha b^2\cos (2\alpha )f(|\mathbf{b}-\mathbf{B}|)}{\int b\mathrm{d}b\mathrm{d}\alpha b^{2}f(|\mathbf{b}-\mathbf{B}|)}.$$

$$⟨{\rho }^a(\mathbf{x}){\rho }^b(\mathbf{y})⟩$$

$$\mu (b)$$

$$\tilde{\mu }(\Delta )$$

$$A_a^{-}(\mathbf{x})=\int {\mathrm{d}}^2\mathbf{z}G(\mathbf{x}-\mathbf{z}){\rho }^a(\mathbf{z}),G(\mathbf{b})\equiv \int \frac{{\mathrm{d}}^2\mathbf{q}}{(2\pi {)}^2}\frac{{\mathrm{e}}^{i\mathbf{q}\cdot \mathbf{b}}}{q^2}\simeq \frac{1}{4\pi }\ln \frac{1}{b^2{m}^2},$$

$$⟨A_a^{-}(\mathbf{x})A_b^{-}(\mathbf{y})⟩={\delta }^{ab}\gamma (\mathbf{x},\mathbf{y}),$$

$$\gamma (\mathbf{x},\mathbf{y})\equiv \int {\mathrm{d}}^2\mathbf{z}G(\mathbf{x}-\mathbf{z})G(\mathbf{y}-\mathbf{z})\mu (\mathbf{z})=\int \frac{{\mathrm{d}}^2\mathbf{q}}{(2\pi {)}^2}\frac{{\mathrm{d}}^2{\mathbf{q}}^{\prime }}{(2\pi {)}^2}{\mathrm{e}}^{\mathrm{i}{\mathbf{q}}^{\prime }\cdot \mathbf{x}+\mathrm{i}\mathbf{q}\cdot \mathbf{y}}\frac{\tilde{\mu }({\mathbf{q}}^{\prime }+\mathbf{q})}{q^{\prime 2}{q}^2}.$$

$$N_{2g}(\mathbf{x},\mathbf{y})$$

$$N_{2g}({\bm{b}},{\bm{r}})=\frac{g^{2}C_{F}}{2}\int\frac{{\rm d}^{2}{\bm{\Delta}}}{(2\pi)^{2}}\frac{{\rm d}^{2}{\bm{k}}}{(2\pi)^{2}}\frac{\tilde{\mu}({\bm{\Delta}})}{({\bm{k}}+{\bm{\Delta}}/2)^{2}({\bm{k}}-{\bm{\Delta}}/2)^{2}}\,{\rm e}^{{\rm i}{\bm{\Delta}}\cdot{\bm{b}}}\Big{[}\,{\rm e}^{{\rm i}{\bm{\Delta}}\cdot{\bm{r}}/2}+\,{\rm e}^{-{\rm i}{\bm{\Delta}}\cdot{\bm{r}}/2}-2\,{\rm e}^{{\rm i}{\bm{k}}\cdot{\bm{r}}}\Big{]}.$$

$${\tilde{N}}_{2g}(\mathbf{b},\mathbf{p})$$

$$\frac{1}{(\mathbf{p}+\mathbf{\Delta }/2{)}^2(\mathbf{p}-\mathbf{\Delta }/2{)}^2}=\frac{1}{{\left(p^2+\frac{{\Delta }^2}{4}\right)}^2-(\mathbf{p}\cdot \mathbf{\Delta }{)}^2}\approx \frac{1}{p^4}\left(1-\frac{{\Delta }^2}{2p^2}+\frac{(\mathbf{p}\cdot \mathbf{\Delta }{)}^2}{p^4}\right)+\dots ,$$

$$4\pi R^2\int \frac{{\mathrm{d}}^2\Delta }{(2\pi {)}^2}{\mathrm{e}}^{-{\Delta }^2{R}^2}{\Delta }^i{\Delta }^j{\mathrm{e}}^{i\mathbf{b}\cdot \mathbf{\Delta }}$$

$$-{\tilde{N}}_{2g}(\mathbf{b},\mathbf{p})\simeq \frac{g^2{C}_F}{p^4}\mu (b)\left[1-\frac{b^2}{8p^2{R}^4}\cos (2\varphi )\right].$$

$$v_2(p,b)\simeq \frac{{\int }_{-\pi }^{\pi }\mathrm{d}\varphi \cos (2\varphi )\left[1-\frac{b^2}{8p^2{R}^4}\cos (2\varphi )\right]}{{\int }_{-\pi }^{\pi }\mathrm{d}\varphi \left[1-\frac{b^2}{8p^2{R}^4}\cos (2\varphi )\right]}=-\frac{b^2}{16p^2{R}^4}.$$

$$S(\mathbf{b},\mathbf{r})=\exp \{-N_{2g}(\mathbf{b},\mathbf{r})\},$$

$$N_{2g}(\mathbf{b},\mathbf{r})$$

$$\frac{1}{2}\Big{[}\,{\rm e}^{{\rm i}{\bm{\Delta}}\cdot{\bm{r}}/2}+\,{\rm e}^{-{\rm i}{\bm{\Delta}}\cdot{\bm{r}}/2}-2\,{\rm e}^{{\rm i}{\bm{k}}\cdot{\bm{r}}}\Big{]}\,\simeq\,-{\rm i}{\bm{k}}\cdot{\bm{r}}+\frac{1}{2}({\bm{k}}\cdot{\bm{r}})^{2}-\frac{1}{8}({\bm{\Delta}}\cdot{\bm{r}})^{2}\,,$$

$$N_{2g}({\bm{b}},{\bm{r}})\,\simeq\,\frac{g^{2}C_{F}}{2}\,r^{i}r^{j}\int\frac{{\rm d}^{2}\Delta}{(2\pi)^{2}}\frac{{\rm d}^{2}{\bm{k}}}{(2\pi)^{2}}\frac{\big{(}k^{i}k^{j}-\Delta^{i}\Delta^{j}/4\big{)}\tilde{\mu}(\Delta)}{[({\bm{k}}+{\bm{\Delta}}/2)^{2}+m^{2}][({\bm{k}}-{\bm{\Delta}}/2)^{2}+m^{2}]}\,{\rm e}^{{\rm i}{\bm{\Delta}}\cdot{\bm{b}}}\,.$$

$$G(\mathbf{b})=\int \frac{{\mathrm{d}}^2\mathbf{q}}{(2\pi {)}^2}\frac{{\mathrm{e}}^{i\mathbf{q}\cdot \mathbf{b}}}{q^2+m^2}=\frac{1}{2\pi }{\mathrm{K}}_0(mb),$$

Full text

I Introduction

The unexpectedly large azimuthal asymmetries in hadron production observed in high-multiplicity events in proton-proton ($pp$) and proton(deuteron)-nucleus ($pA$) collisions at the LHC and RHIC Khachatryan et al. (2010); Chatrchyan et al. (2013a); Chatrchyan et al. (2013b); Abelev et al. (2013a, b, 2014); Aad et al. (2013a, b, 2014); Adare et al. (2013, 2015); Adamczyk et al. (2015a, b) have triggered intense debates concerning the physical origin of such phenomena. It is indeed an outstanding problem to understand how a small system like that produced in $pp$ or $pA$ collisions, an order of magnitude smaller than in nucleus-nucleus ($AA$) collisions, can develop a collective behavior which is quite similar to that observed in $AA$ collisions, both in terms of its magnitude and in terms of its dependences upon the transverse momenta, the rapidities, and the masses of the produced hadrons Abelev et al. (2013a, b, 2014); Aad et al. (2013a, b, 2014). Roughly speaking, the associated scientific debate opposes two paradigms. The first of them, which is closer to the generally accepted interpretation of the corresponding phenomena in $AA$ collisions, relates the azimuthal correlations observed in $pp$ and $pA$ collisions to ‘hydrodynamic flow’, i.e. collective effects caused by strong interactions in the final state. Whereas such scenarios may indeed lead to reasonable descriptions of the data (at least for sufficiently small transverse momenta and with suitable choices for the initial conditions) d’Enterria et al. (2010); Bozek (2011); Bozek and Broniowski (2013a, b); Bozek et al. (2013); Qin and M ller (2014); Werner et al. (2014); Bzdak and Ma (2014); Weller and Romatschke (2017), it seems nevertheless difficult to conceive that hydrodynamic flow may develop in such small systems. This motivated the second paradigm, which rather builds upon the ‘initial state’ physics, i.e. the collective phenomena associated with high parton densities in the wavefunctions of (one or both of) the incoming hadrons, prior to their collision Kovchegov et al. (2001); Teaney and Venugopalan (2002); Kovchegov and Tuchin (2002); Dumitru et al. (2008); Gavin et al. (2009); Dumitru et al. (2011); Avsar et al. (2011); Kovner and Lublinsky (2011a, b); Levin and Rezaeian (2011); Iancu and Triantafyllopoulos (2011); Schenke et al. (2012a, b); Dusling and Venugopalan (2012, 2013); Dusling et al. (2016); Kovchegov and Wertepny (2014); Kovner and Lublinsky (2013); Kovner and Rezaeian (2014, 2015); Schenke et al. (2015); Lappi (2015); Lappi et al. (2016); Rezaeian (2016); Schenke and Schlichting (2016); Schenke et al. (2016).

In practice, the azimuthal asymmetries are most conveniently measured via multi-particle angular correlation. But at a conceptual level, it is often preferable to think in terms of the single-inclusive particle distribution event-by-event and its dependence upon the azimuthal angle $\phi$, as measured w.r.t. the ‘reaction plane’. More precisely, $\phi$ is the angle between the direction of motion of a produced hadron in the transverse plane and its impact factor. Then the azimuthal asymmetries are encoded in the ‘flow coefficients’ $v_{n}(p_{T})$ — the $\cos(n\phi)$ Fourier moments of the single-inclusive distribution in $\phi$ (see e.g. Snellings (2011)). From this perspective, the azimuthal asymmetries reflect a spontaneous breaking of rotational symmetry in the transverse plane, which may have various origins. The best known example is that of non-central $AA$ collisions, where the rotational symmetry is broken by the elliptic shape of the ‘interaction region’ (the overlapping region between the two nuclei) in the transverse plane Snellings (2011). More generally (and including for central $AA$ collisions), azimuthal anisotropies can be generated by fluctuations in the distribution of particles (nucleons, or even gluons inside the participating nucleons) within the incoming nuclei Alver and Roland (2010). Clearly, the typical transverse sizes will be different for nucleon number, respectively, gluon number fluctuations, potentially leading to different laws for the $p_{T}$-dependence of the coefficients $v_{n}(p_{T})$.

In the context of $AA$ collisions, the theoretical ideas concerning the particle (nucleon or parton) number fluctuations are naturally embedded in the initial conditions for the hydrodynamical equations. But such fluctuations can generate momentum-space azimuthal asymmetries already by themselves, that is, even in the absence of interactions in the final state leading to (hydrodynamic) flow. This is particularly interesting for the $pp$ and $pA$ collisions, where the importance of the final-state interactions is far from being established. In this context, most of the calculations associated with the ‘initial-state’ paradigm alluded to above relied on the effective theory for the color glass condensate (CGC) Iancu and Venugopalan (2003); Gelis et al. (2010), which predicts the event-by-event formation of ‘saturation domains’ inside a dense hadronic target. These are regions with a typical transverse size $1/Q_{s}$ ($Q_{s}$ is the target saturation momentum) where gluons have large occupation numbers and are coherent with each other, so they can also be described as condensates of strong chromo-electric fields. These domains introduce a preferred direction in the transverse plane — the orientation of the color fields — thus breaking the rotational symmetry. In turn, this leads to azimuthal correlations in the particle production in the collision between a dilute proton (the ‘projectile’) and the dense target: if 2 partons from the projectile have similar impact parameters (so that they scatter off a same saturation domain) and they are in the same color state, then they will receive similar kicks and hence emerge along nearby angles. The orientations of the saturation domains are of course random, so their effects will be washed out (by the averaging over the events) in the calculation of single-inclusive particle production. But non-trivial correlations survive in the production of two (or more) particles, with very interesting features: the respective spectra are naturally ‘semi-hard’ (the flow coefficients $v_{n}(p_{T})$ are peaked around the saturation scale $Q_{s}$), their strength decreases with increasing $Q_{s}$ (hence, in particular, with increasing energy), and they are suppressed in the limit where the number of colors $N_{c}$ is large (they scale like $1/N_{c}^{2}$). The ‘color field domain’ model proposed in Kovner and Lublinsky (2011b); Dumitru and Giannini (2015); Dumitru and Skokov (2015); Dumitru et al. (2015) can be viewed too a (rather extreme) variant of this scenario: as shown in Lappi et al. (2016), the effects of these ‘color field domains’ can be reproduced by fluctuating color fields in CGC provided non-Gaussian correlations are assumed to be important.

As it should be clear from the above discussion, the CGC-based approaches assume the $pp$ and $pA$ collisions to be of the ‘dilute-dense’ type. This is quite natural when the target is a large nucleus, such as Pb with $A=208$, and it is also justified for $pp$ collisions so long as one considers particle production at very forward rapidities — meaning that the gluon distribution from the target proton has been subjected to the high-energy (or small-$x$) evolution and hence is much denser than that of the projectile proton.

In this paper, we shall remain within the general CGC framework of ‘dilute-dense scattering’, but we shall explore a more elementary mechanism for generating azimuthal correlations: the case where the rotational symmetry in the transverse plane is broken simply by the impact parameter (a 2-dimensional vector ${\bm{b}}$) of an impinging parton from the projectile, i.e. by the very fact that a parton hits the target disk at some point ${\bm{b}}$ which is away from the center ($b\neq 0$). Via its scattering, the parton will acquire some transverse momentum ${\bm{p}}$ and the respective cross-section will generally dependent upon the angle $\phi$ made by the vectors ${\bm{p}}$ and ${\bm{b}}$. (From now on, we suppress the subscript $T$ on transverse momenta or coordinates, to simplify writing.) For this to be the case, the cross-section must depend upon ${\bm{b}}$ in the first place, that is, the target should have some inhomogeneity in the transverse plane. Accordingly, this mechanism naturally generates azimuthal asymmetries which probe the variation of the transverse distribution of matter in the target. These asymmetries do not require non-planar gluon exchanges, hence they admit a non-zero limit when $N_{c}\to\infty$.

The basic idea is not new — it has been originally proposed in Refs. Kopeliovich et al. (2008a, b, c) and more recently revisited in Ref. Levin and Rezaeian (2011); Zhou (2016); Hagiwara et al. (2017). As in these previous works, we shall use the high-energy factorization for particle production in ‘dilute-dense’ collisions in which the cross-section for the production of a parton with transverse momentum ${\bm{p}}$ at ${\bm{b}}$ is related to the Fourier transform (${\bm{r}}\to{\bm{p}}$) of the $S$-matrix $S({\bm{r}},{\bm{b}})$ for the elastic scattering a color dipole with transverse size ${\bm{r}}$ and impact parameter ${\bm{b}}$ (say, a quark-antiquark dipole for the case of quark production, that we shall focus here on, for definiteness). In this framework, the azimuthal asymmetry results from the dependence of the function $S({\bm{r}},{\bm{b}})$ upon the dipole orientation (the angle $\theta$ made by the vectors ${\bm{r}}$ and ${\bm{b}}$).

As compared to the previous literature, we shall use a different theoretical description for the gluon distribution in the target, namely, the McLerran-Venugopalan (MV) model McLerran and Venugopalan (1994), that we here extend to include inhomogeneity in the transverse plane. Our respective extension will be inspired by ‘saturation models’ like IP-Sat Kowalski and Teaney (2003); Kowalski et al. (2008); Rezaeian et al. (2013) and bCGC Iancu et al. (2004); Kowalski et al. (2006); Rezaeian and Schmidt (2013), which include a non-trivial $b$-dependence that has been tested and calibrated via fits to the HERA data for diffractive vector meson production. (See also Refs. Munier et al. (2001); Berger and Stasto (2011, 2013); Armesto and Rezaeian (2014) for other studies of the HERA phenomenology which explore the impact-parameter dependence of the saturation physics.) Note that one cannot directly use those ‘saturation models’ for the present purposes, since they are formulated as parametrizations for the dipole $S$-matrix which do not include any angular dependence. By contrast, the MV model offers a description for the distribution of the ‘valence’ color sources in the target and allows for an explicit calculation of the dipole scattering off the color fields produced by those sources. In the previous literature, this calculation has been performed for a homogeneous target, but here we shall extend it to the case where the ‘valence’ color sources have a Gaussian distribution in impact parameter (including lumpiness effect for a large nucleus). Due to the formal simplicity of the model, we will be able to obtain quasi-analytic results for the dipole $S$-matrix $S({\bm{r}},{\bm{b}})$, including its angular dependence, for both a smooth target (‘a dense proton’) and a lumpy one (‘a large nucleus’). Such a quasi-analytic treatment turns out to be very useful for the physical interpretation of our results.

One of the main conclusions of our analysis is that the azimuthal asymmetry is controlled by soft exchanged momenta, of the same order as the transverse momentum scale introduced by the inhomogeneity of the target. So, strictly speaking, these effects can only be marginally addressed within our present, semi-classical, formalism, which is inspired by perturbative QCD. But this also shows the importance of having a realistic model for the $b$-dependence, that was already tested against the phenomenology. In fact, due to the dominance of soft exchanges, we will also be led to consider the influence of a ‘gluon mass’, namely a mass parameter which controls the exponential decay of the gluon fields at large values of $b$ and mimics confinement. Not surprisingly, this influence appears to be important, on the same footing as that of the other parameters of the model — the transverse scale which characterizes the Gaussian distribution of the color sources and the target saturation scale at $b=0$. These non-perturbative aspects — the transverse inhomogeneity of the target and the gluon confinement — have been differently modeled in the previous related studies Kopeliovich et al. (2008a, b, c); Levin and Rezaeian (2011); Zhou (2016); Hagiwara et al. (2017). This may explain the significant differences that can be observed between our new results and the previous ones in the literature.

A priori, our model produces non-zero values for all the even Fourier coefficients $v_{2n}(p,b)$ with $n\geq 1$, but with a strong hierarchy among them: $|v_{2}|\gg|v_{4}|\gg|v_{6}|$ etc. In order to render the model tractable, we shall perform additional approximations, which will preserve only the information about the largest such a coefficient, the elliptic flow $v_{2}(p,b)$. Using our (quasi)analytic results for the dipole $S$-matrix, it will be quite easy to compute and study this coefficient. We shall thus find that $v_{2}(p,b)$ can be quite large, $v_{2}\gtrsim 0.1$, in peripheral collisions, but it rapidly decreases when moving towards more central collisions. This corresponds to the fact that, in our model, the transverse inhomogeneity is peaked at the edge of the target.

We shall furthermore find that the effects of multiple scattering are truly essential: they can even change the sign of $v_{2}(p,b)$. Namely, $v_{2}(p)$ is found to be negative (but also tiny) for very large values of $p$, where the single-scattering approximation applies, but it turns positive, due to multiple scattering, at lower values — including the most interesting physical regime where $p$ is soft or semi-hard. A positive $v_{2}(p,b)$ means that the preferred direction of motion for a produced particle is along its impact parameter ${\bm{b}}$. In terms of dipole scattering, it means that the scattering is stronger (for a given dipole size $r$) when the dipole is aligned along ${\bm{b}}$ rather than perpendicular on it.

Finally, as a function of $p$, $v_{2}$ shows a maximum at a value proportional to the target saturation momentum $Q_{s}$ (say, as measured at $b=0$). Interestingly, this maximum becomes less pronounced (broader and smaller) when increasing $Q_{s}$, i.e. when the target becomes denser. This $p$-dependence looks acceptable from the viewpoint of the phenomenology and in fact it can even be adjusted to reasonably describe the data in p+Pb collisions at the LHC with reasonable choices for the parameters.

So far, our discussion refers to a fixed value of the impact parameter ${\bm{b}}$: the quantity $v_{2}(p,b)$ characterizes the distribution of the produced particles w.r.t. the reaction plane in a particular event. Since the direction of ${\bm{b}}$ (the ‘reaction plane angle’) is not an observable, it is important to notice that the azimuthal asymmetry under consideration can also be measured via multi-particle correlations. Indeed, as previously mentioned, this asymmetry is sizable only for sufficiently peripheral collisions, in which the interaction region is relatively small. In $pp$ collisions, this region should be much smaller than any of the colliding protons. In $pA$ collisions, it could be as large as the size of projectile proton, but this is still small compared to the relevant impact parameters, of the order of the nuclear radius. The partons from the projectile which participate in such peripheral collisions have similar impact parameters, hence after the scattering they preferentially propagate along nearby directions — namely, along their average impact parameter. This in turn implies the existence of azimuthal asymmetries in the multi-particle correlations; e.g., $v_{2}\{2\}$ — the elliptic azimuthal coefficient extracted from 2-particle correlations — should be non-zero and positive.

The above considerations have consequences not only for the multi-particle correlations, but also for the single inclusive particle spectrum that we shall focus on in this paper. They imply that the strength of the azimuthal asymmetries is also controlled by the geometry of the interaction region. For instance, we shall find that the elliptic flow $v_{2}$ is (roughly) proportional to the eccentricity $\varepsilon_{2}$, which is a measure of the projection of the impact parameters of the participants quarks along the direction of their average impact parameter. Such geometrical aspects are clearly reminiscent of the classical discussion of hydrodynamic flow in $AA$ collisions — in both cases, a flow of particles in the final state is generated via peripheral collisions — but the underlying dynamics is of course different: whereas in peripheral $AA$ collisions the flow is driven by the ‘pressure gradient’ (the final state interactions) associated with the spatial asymmetry of the interaction region, in the new mechanism of interest for us here, the flow is rather a consequence of the angular dependence of the amplitude for dipole scattering.

Although in this paper we shall discuss only the average target geometry, it is quite clear that a similar mechanism should also act when the target inhomogeneity is associated with fluctuations — say, in the gluon distribution produced by the high-energy evolution, or in the distribution of nucleons inside a lumpy nuclear target. In the presence of fluctuations, azimuthal asymmetries can be also generated via more central collisions, but they can be probed only via multi-particle correlations, which are suppressed in the multicolor limit $N_{c}\to\infty$ Dumitru et al. (2008). This discussion suggests that the mechanism to be discussed here is closely connected to that from the ‘glasma’ scenario, where the azimuthal asymmetries are associated with fluctuations leading to ‘saturation domains’ Schenke et al. (2012a, b); Dusling and Venugopalan (2012, 2013); Dusling et al. (2016); Kovchegov and Wertepny (2014); Kovner and Lublinsky (2013); Kovner and Rezaeian (2014, 2015); Schenke et al. (2015); Lappi (2015); Lappi et al. (2016); Rezaeian (2016); Schenke and Schlichting (2016); Schenke et al. (2016). Perhaps a new aspect which is specific to our discussion is the emphasis on peripheral collisions: we show that such collisions can generate sizable azimuthal asymmetries already in the absence of fluctuations. Since related to the (average) target geometry, these asymmetries are expected to factorize in the calculation of multi-particle correlations (e.g. $c_{2}\{2\}\simeq v_{2}^{2}$, where $c_{2}\{2\}$ is the second-order cumulant Snellings (2011)) and also to survive in the large-$N_{c}$ limit.

This paper is organized as follows: In Sect. II, we concisely describe the factorization scheme that we use for quark production in ‘dilute-dense’ collisions and the associated calculation of the azimuthal asymmetry coefficients $v_{n}(p,b)$ in a given event. Sect. III contains our new analytic results. After introducing the (impact-parameter dependent) MV model for the gluon distribution in the dense target in Sect. III.A, we present the calculation of the dipole $S$-matrix with angular dependence, first in the single scattering approximation (in Sect. III.B), next by including the effects of multiple scattering, separately for a proton (in Sect. III.C) and for a large nucleus viewed as a lumpy superposition of independent nucleons (in Sect. III.D). In Sect. IV, we present our numerical results for $v_{2}(p,b)$ and discuss their dependence upon various parameters of the model as well as possible implications for the phenomenology. We summarize our results in Sect. V.

II Color-dipole orientation as the origin of the azimuthal asymmetry

Consider particle production in a dilute-dense collision, say a proton-nucleus ($pA$) collision, for definiteness, but the target could also be another proton provided the produced particle propagates at very forward rapidity. We shall view this process at partonic level to leading order in perturbative QCD at high gluon density (i.e. in the CGC effective theory). For simplicity we shall ignore the fragmentation of the produced parton into hadrons. That is, we shall only compute the cross-section for parton production, with the parton chosen to be a quark. (The discussion of gluon production in this particular set-up would be entirely similar.) To the accuracy of interest, the correct physical picture is as follows: a quark collinear with the projectile proton undergoes multiple scattering off the dense gluon distribution of the target and thus acquires some transverse momentum ${\bm{p}}$. The multiple scattering can be resummed to all orders within the eikonal approximation, which is most conveniently formulated in impact parameter space (since the transverse coordinate of the quark is not modified by the interactions). The cross-section is proportional to the modulus squared of the amplitude and the quark impact parameters in the two amplitudes, direct and conjugate, are different. As a result, one can express the rapidity and ${\bm{p}}$-distribution at fixed impact parameter in terms of an effective $q\bar{q}$ dipole $S$-matrix,

[TABLE]

Here, ${\bm{p}}$ and ${\bm{b}}$ are the transverse momentum and the impact parameter of the produced quark and $\eta$ is its rapidity in the center-of-frame (COM) frame. Furthermore, $x_{p}$ and $x_{g}$ are the longitudinal momentum fractions of the partons participating in the scattering: the ‘collinear’ quark from the proton and a gluon from the wavefunction of the nucleus. Energy-momentum conservation implies

[TABLE]

where $p\equiv|{\bm{p}}|$ and $s$ is the COM energy squared for the scattering between the proton and one nucleon from the nucleus. The quantity $S({\bm{b}},{\bm{r}},x_{g})\equiv S({\bm{x}},{\bm{y}},x_{g})$, with ${\bm{b}}\equiv({\bm{x}}+{\bm{y}})/2$ and ${\bm{r}}\equiv{\bm{x}}-{\bm{y}}$, is the forward $S$-matrix for the scattering between a quark-antiquark dipole (with the quark leg at ${\bm{x}}$ and the antiquark one at ${\bm{y}}$) and the nucleus, for a rapidity separation $Y=\ln(1/x_{g})$. Its Fourier transform $\tilde{S}({\bm{b}},{\bm{p}},x_{g})$ plays the role of a generalized unintegrated gluon distribution (a.k.a. a gluon ‘transverse momentum distribution’, or TMD) in the target. Since the dipole has a finite size and an orientation, its scattering will generally depend upon the angle $\theta$ between ${\bm{b}}$ and ${\bm{r}}$. Via the Fourier transform, this will introduce an anisotropy in the cross-section for quark production, i.e. a dependence upon the angle $\phi$ between ${\bm{b}}$ and ${\bm{p}}$. This anisotropy can be characterized by the ensemble of Fourier components $v_{n}$ (a.k.a. ‘flow coefficients’), defined as

[TABLE]

The Fourier moments involving $\sin(n\phi)$ vanish because the cross-section is symmetric under the parity transformation $\phi\to-\phi$ (the reflection w.r.t. the reaction plane). The above expression can be also evaluated with the dipole amplitude in the coordinate representation (this will be useful e.g. when including the effects of multiple scattering in the eikonal approximation). Rewriting Eq. (1) as

[TABLE]

one can perform the integral over $\phi$ in Eq. (3) with the help of the the following identity,

[TABLE]

where ${\rm J}_{n}(x)$ denotes the Bessel function of the first kind. One thus obtains e.g.

[TABLE]

Notice that the quark distribution function of the proton has canceled in the ratio. But the information about the gluon distribution in the nucleus is still preserved in $v_{n}$, via the dipole $S$-matrix. If one neglects the angular dependence of the latter ($S({\bm{b}},{\bm{r}})\to S(b,r)$), then $v_{n}=0$ for any $n\geq 1$ regardless of the precise shape of the target profile in ${\bm{b}}$. Notice also that a real contribution to $v_{3}$ requires the existence of an imaginary part in the dipole $S$-matrix; hence, a non-zero $v_{3}$ can be related to the odderon contribution to dipole scattering Kovchegov et al. (2004); Hatta et al. (2005).

So far, we have implicitly treated the projectile proton as a pointlike objet (indeed, we assumed that all its valence quarks have the same impact parameter ${\bm{b}}$). As we shall see, this is indeed a good approximation when the target is a large nucleus and for relatively central collisions, where the nuclear matter distributions is quasi-homogeneous. But this is less justified for the case of a proton target, or for peripheral collisions off a nucleus, which are the cases of main interest for what follows. Fortunately, this can be easily remedied (at least at a formal level) by replacing the standard quark distribution $x_{p}q(x_{p})$ with its generalized version (a “generalized parton distribution”, or GPD), which includes impact parameter dependence inside the projectile: $x_{p}q(x_{p})\to x_{p}q(x_{p},{\bm{b}})$, where ${\bm{b}}$ now refers to the position of a quark relative to the center of its parent proton. Then Eqs. (1) and (3) should be then replaced by

[TABLE]

and respectively

[TABLE]

where ${\bm{B}}$ denotes the impact parameter of the proton w.r.t. the center of the target and $\Phi$ is the angle made by the vectors ${\bm{p}}$ and ${\bm{B}}$. A common prescription in the literature, that we here adopt as well, is to assume the factorization of the ${\bm{b}}$-dependence inside the projectile: $x_{p}q(x_{p},{\bm{b}})\simeq x_{p}q(x_{p})f({\bm{b}})$, with $\int{{\rm d}^{2}{\bm{b}}}\,f({\bm{b}})=1$. Under this assumption, the generalization of Eq. (6) to an extended projectile is easily found as

[TABLE]

where $\alpha$ is the angle between ${\bm{b}}$ and ${\bm{B}}$ and we have also assumed that the proton distribution is isotropic: $f({\bm{b}})=f(|{\bm{b}}|)$.

We shall also need the eccentricity $\varepsilon_{2}$ of the interaction region. Writing the 2-dimensional impact parameter of a participating quark as ${\bm{b}}=(x,y)$, where the $x$ and $z$ axes define the reaction plane (that is, the $x$ axis is parallel to the direction of the vector $B$ — the impact parameter of the proton), then $\varepsilon_{2}$ can be estimated as

[TABLE]

where the brackets denote the averaging over the quark impact parameters with weight-function given by the local differential cross-section $\propto\tilde{S}({\bm{b}},{\bm{p}},x_{g})$ and also the angular average over $\Phi$ (i.e. over the direction of propagation of the produced particles). As before, $\alpha$ is the angle between ${\bm{b}}$ and ${\bm{B}}$, hence $x=b\cos\alpha$ and $y=b\sin\alpha$ (see also Fig. 1).

In fact, the quantity that is generally referred to as the ‘eccentricity’ in the literature is the ‘momentum-integrated’ version of Eq. (11), that is, the quantity $\varepsilon_{2}(B)$ which is obtained by separately integrating the numerator and denominator over $p$, with measure $\int p{\rm d}p$. The result of this integration turns out to be very simple and in particular independent of the dipole scattering, because of ‘color transparency’ : a color dipole with zero transverse size cannot scatter, meaning that $S(b,r,\theta)\to 1$ as $r\to 0$. This immediately implies the following ‘sum-rule’ :

[TABLE]

This ‘sum-rule’ is in fact the expression of probability conservation: as manifest from Eq. (8), $\tilde{S}({\bm{b}},{\bm{p}})$ can be interpreted as the probability density for a quark incident at ${\bm{b}}$ to acquire a transverse momentum ${\bm{p}}$. Clearly, the total probability for the quark to emerge with any momentum must be equal to one. Using Eqs. (11) and (12), one finds

[TABLE]

As anticipated, this is a purely geometrical quantity, without any information about the scattering of the dipole: it merely shows how the projectile is ‘seen’ from the center of the target. It is quite obvious that $\varepsilon_{2}(B)$ vanishes for $B=0$ and that it approaches to 1 when $B\to\infty$ (since in that limit the ratio $\langle y^{2}\rangle/\langle x^{2}\rangle$ approaches to zero, meaning that $\alpha\to 0$ as well). This behavior will be confirmed by the explicit calculations to be presented later.

III Elliptic flow from dipole orientation in the MV model

III.1 Dipole scattering in the McLerran-Venugopalan model

So long as the gluon energy fraction $x_{g}$ probed by the scattering is not too low, one can ignore the high-energy evolution of the nuclear gluon distribution and describe the latter within the McLerran-Venugopalan (MV) model McLerran and Venugopalan (1994). In this model, the nucleus is described as a collection of independent color sources (‘the valence quarks’) with a Gaussian color charge distribution in the transverse plane:

[TABLE]

with $\mu({\bm{x}})$ the color charge squared per unit area in the transverse plane. In the original formulation of the MV model, as applying to a (very) large nucleus, there was no explicit impact-parameter dependence: the color charge distribution was assumed to be uniform, $\mu({\bm{b}})=\mu_{0}$, within a large disk with radius $R_{A}\propto A^{1/3}$, with $A$ the atomic number. However, as already mentioned, the inhomogeneity of the target in the transverse plane is essential for the physical effects that we are currently interested in. So, in what follows we shall propose a generalization of the MV model which includes a physically motivated impact-parameter dependence (inspired by the fits to the HERA data). We shall first present this dependence for the case where the target is a single nucleon (say, any of the $A$ nucleons composing a large nucleus), and then extend the model to a nuclear target, in Sect. III.4. But the case of a proton target is also interesting by itself — albeit the applicability of the MV model to this case is of course questionable —, in view of the phenomenology of flow-like effects in high-multiplicity events in $pp$ collisions.

The profile $\mu({\bm{b}})$ of the proton color charge distribution in the transverse plane will be chosen to be a Gaussian, in agreement with saturation fits Kowalski and Teaney (2003); Kowalski et al. (2008); Rezaeian et al. (2013); Rezaeian and Schmidt (2013) to the HERA data on diffraction and vector meson production in DIS (its Fourier transform $\tilde{\mu}(\Delta)$ will be later needed):

[TABLE]

The overall factor $\mu_{0}$, with dimensions of mass squared, is proportional to the total color charge squared of the valence partons; e.g. for a proton target with $N_{c}$ valence quarks, one can write $4\pi R^{2}(N_{c}^{2}-1)\mu_{0}=g^{2}C^{F}N_{c}$, or $\mu_{0}=\alpha_{s}/(2R^{2})$. In practice, this quantity $\mu_{0}$ will be traded for the saturation momentum in the center ($b=0$) of the target, a quantity which is constrained by the fits to the HERA data (see Eq. (35) below). The scale $R$ which fixes the width of the $b$-distribution will be taken too from fits to the HERA data (more precisely from the fits using the IP-Sat model in Refs. Kowalski and Teaney (2003); Rezaeian et al. (2013)). A typical value emerging from these fits is $R\simeq 0.3$ fm. Notice that, with our present conventions, the ‘proton size’ (in the sense of the region in impact parameter space where the valence color charges are distributed) is $2R\simeq 0.6$ fm, and not $R$. Indeed, the exponent of the Gaussian in Eq. (15) becomes equal to one when $b=2R$. Also, for $\Delta=0$, one has $\tilde{\mu}(0)=4\pi R^{2}\mu_{0}$, showing that the natural ‘proton area’ is $4\pi R^{2}$, and not $\pi R^{2}$.

The gauge potential created by the ultrarelativistic color charges is simply the 2-dimensional Coulomb field,

[TABLE]

where $m$ is an infrared cutoff, physically associated with confinement: $m\sim\Lambda_{{}_{\rm QCD}}$. (So, the second estimate for $G({\bm{b}})$ given above applies only for sufficiently small distances $b\lesssim 1/m$.) This implies that the distribution of the color fields is Gaussian as well, with 2-point correlation

[TABLE]

where

[TABLE]

In the MV model and in the eikonal approximation, the projectile dipole independently scatters off the color charges in the nucleus. Accordingly, its multiple scattering exponentiates, as in the Glauber approximation: $S=\exp\{-N_{2g}\}$, where $N_{2g}({\bm{x}},{\bm{y}})$ is the amplitude for a single scattering via the exchange of two gluons:

[TABLE]

Note that, unlike the valence color charges, which are effectively confined in the transverse plane within a disk with radius $R$, cf. Eq. (15), the color fields created by these charges (the small-$x$ gluons) can be delocalized over much larger distances, due to the slow decay of the 2-dimensional Coulomb propagator at large distances. In particular, it is easy to check that for very large impact factors $x\sim y\gg R$ (with $x\equiv|{\bm{x}}|$ and $y\equiv|{\bm{y}}|$), the dipole amplitude predicted by this model shows a power tail: $N_{2g}({\bm{x}},{\bm{y}})\propto({\bm{b}}\cdot{\bm{r}})^{2}/b^{4}$, with ${\bm{r}}={\bm{x}}-{\bm{y}}$ and ${\bm{b}}=({\bm{x}}+{\bm{y}})/2$.

At this stage, it is convenient to change variables, from ${\bm{x}},{\bm{y}}$ to ${\bm{r}},{\bm{b}}$, and from ${\bm{q}}^{\prime},{\bm{q}}$ to ${\bm{\Delta}},{\bm{k}}$, with ${\bm{k}}=({{\bm{q}}^{\prime}-{\bm{q}}})/{2}$ and ${\bm{\Delta}}={\bm{q}}^{\prime}+{\bm{q}}$. For the physical discussion to follow, it is useful to keep in mind the physical meaning of the momenta ${\bm{k}}$ and ${\bm{\Delta}}$ from the viewpoint of our original problem, that of quark production: (i) ${\bm{k}}$ is the average transverse momentum transmitted by the nucleus to the quark via a single collision; (ii) ${\bm{\Delta}}$ is the difference between the transverse momenta acquired by the quark in the direct amplitude and the complex conjugate amplitude, respectively; as such, it is a measure of the additional momentum transfer associated with the inhomogeneity of the target. Clearly, ${\bm{\Delta}}$ is a soft momentum, $\Delta\sim 1/R\sim\Lambda_{{}_{\rm QCD}}$, whereas is generally semi-hard, that is, it is either comparable to the final momentum $p$ of the produced quark, or to the saturation momentum $Q_{s}(b)$ of the target. Yet, soft values for ${\bm{k}}$ will be important too, when discussing the flow coefficients in the presence of multiple scattering.

We thus obtain,

[TABLE]

The first two terms within the square brackets, which are independent of ${\bm{k}}$, represent ‘tadpole’ contributions where the two gluons exchanged with the target are attached to a same fermion leg (the quark or the antiquark). The final term, which is negative, refers to ‘exchange’ contributions, where one gluon is attached to the quark leg and the other one, to the antiquark.

Since $\tilde{\mu}({\bm{\Delta}})$ is truly a function of $\Delta\equiv|{\bm{\Delta}}|$, it is quite obvious that $N_{2g}({\bm{b}},{\bm{r}})$ is an even function of ${\bm{b}}$ and also an even function of ${\bm{r}}$; hence, it depends upon $\theta$ (the angle between ${\bm{b}}$ and ${\bm{r}}$) only via the squared dot product $({\bm{b}}\cdot{\bm{r}})^{2}$. This in turn implies that all the odd ‘flow coefficients’, like the ‘radial flow’ $v_{1}$ and the ‘triangular’ one $v_{3}$, must vanish. In what follows, we shall compute the elliptic flow $v_{2}$. For pedagogy, we shall first present the respective calculation in the single scattering approximation.

III.2 The single scattering approximation

The single scattering approximation $S({\bm{b}},{\bm{r}})\!\simeq\!1-N_{2g}({\bm{b}},{\bm{r}})$ applies so long the dipole is small enough for its transverse resolution $Q^{2}\equiv 1/r^{2}$ to be much larger than the (local) saturation momentum $Q_{s}({\bm{b}})$ at its impact parameter. Equivalently (since $pr\sim 1$ by virtue of the Fourier transform in Eq. (1)), the produced quark is relatively hard, with a transverse momentum $p\gg Q_{s}({\bm{b}})$. The saturation scale $Q_{s}({\bm{b}})$ in the MV model will be more precisely defined in the next subsection, where we discuss multiple scattering. Here, we anticipate that this is a semi-hard scale, comparable to, or larger than, the momentum scale $1/R$ introduced by the impact-parameter distribution $\mu(b)$.

To compute $v_{2}$ in the single scattering approximation (SSA), it is convenient to first perform the Fourier transform of the dipole amplitude, $N_{2g}({\bm{b}},{\bm{r}})\to\tilde{N}_{2g}({\bm{b}},{\bm{p}})$, and then use Eq. (3). It is quite clear that the ‘tadpole’ pieces in Eq. (20) do not significantly contribute in the kinematics of interest: via the Fourier transform, the respective exponentials ${\rm e}^{\pm{\rm i}{\bm{\Delta}}\cdot{\bm{r}}/2}$ select $\Delta=2p$, but the function $\tilde{\mu}(2p)$ is exponentially suppressed for $2p\gg 1/R$. As for the Fourier transform of the ‘exchange’ piece in Eq. (20), this is simply obtained by replacing ${\bm{k}}\to{\bm{p}}$. We deduce

[TABLE]

Physically, the fact that ${\bm{k}}={\bm{p}}$ means that the momentum carried by the final quark must be acquired via its only collision with the target.

Eq. (21) can be further simplified by using the fact that $pR\gg 1$, whereas the integral is controlled by softer values $\Delta\lesssim 1/R$. Accordingly, one can expand the integrand in powers of $\Delta/p$ and keep only the leading order piece,

[TABLE]

where the dots stand for terms of order $({\Delta}/{p})^{4}$. After also using Eq. (15), we are led to a Gaussian integral

[TABLE]

Putting everything together and using the trigonometric identity $2\cos^{2}\phi-1=\cos(2\phi)$, we finally deduce

[TABLE]

This holds up to terms suppressed by higher powers of $1/(pR)^{2}$. In this approximation, the dipole amplitude is proportional to $\mu(b)$, hence it is as localized in $b$ as the ‘valence’ color charges from the target. This is so because the scattering involves the exchange of a hard gluon, with momentum $p$, and this exchange is quasi-local.

The leading-order contribution at large $p$, proportional to $1/p^{4}$, is independent of $\phi$. This is recognized as the standard result for the particle spectrum produced via a single, hard, scattering. The angular dependence enters via the subleading term $\propto 1/p^{6}$, whose sign is quite remarkable: this is such that the cross-section for quark production (which in the present approximation is proportional to $\tilde{S}({\bm{b}},{\bm{p}})\simeq-\tilde{N}_{2g}({\bm{b}},{\bm{p}})$) is largest when $\theta=\pi/2$. Physically, this means that a quark produced via a single scattering has more chances to propagate along a direction which is perpendicular on its impact parameter (${\bm{p}}\perp{\bm{b}}$), rather than parallel to it (${\bm{p}}\parallel{\bm{b}}$). In turn, this implies that the elliptic flow coefficient $v_{2}$ is negative in this regime. Namely, by inserting Eq. (24) into Eq. (3), one finds

[TABLE]

Except possibly for its sign, which is somewhat unexpected, the above result for $v_{2}$ shows the expected trends: it vanishes when $b\to 0$, since for such central collisions the orientation of the incoming dipole plays no role, and it decreases quite fast when increasing the momentum $p$ of the produced quarks, as this corresponds to exploring dipoles with very small sizes $r\ll R$.

The above calculation also illustrates another generic feature of the $v_{2}$ (more generally, of the azimuthal anisotropy) generated by the current mechanism: this is directly related to the target inhomogeneity in the transverse plane, i.e. it is proportional to the derivatives of the $b$-distribution $\mu({\bm{b}})$. It should be furthermore clear that the higher azimuthal harmonics $\cos(2n\phi)$ with $n\geq 2$ would be generated via the higher-order terms in the large-$p$ expansion, hence the corresponding Fourier coefficients are parametrically suppressed — by powers of $b^{2}/(p^{2}R^{4})\sim 1/(pR)^{2}$ when $b\sim R$ — compared to the elliptic flow $v_{2}$.

Notice that, in this single-scattering approximation, the overall normalization $\mu_{0}$ of the charge-charge correlator, cf. Eq. (15), and also the coupling constant $g^{2}C_{F}$, drop out from the calculation of $v_{2}$. This last feature will be of course modified by the inclusion of multiple scattering, which becomes compulsory for softer momenta $p\lesssim Q_{s}(b)$ and will be discussed in the next section.

III.3 Adding multiple scattering

The multiple scattering between the quark projectile and the target becomes important when the transverse momentum of the produced particle $p$ is comparable to, or smaller than, the nuclear saturation momentum $Q_{s}(b)$. This is actually the most interesting situation for the phenomenology of flow in $pp$ and $pA$ collisions at the LHC, as we shall see. In that case, we must return to the general expression for the dipole $S$-matrix (within the framework of the MV model, of course), namely

[TABLE]

with $N_{2g}$ as given in Eq. (20). Due to the exponentiation, the Fourier transform ${S}({\bm{b}},{\bm{r}})\to\tilde{S}({\bm{b}},{\bm{p}})$ is more complicated. Physically this reflects the fact that the momentum ${\bm{p}}$ of the produced quark gets accumulated via several scatterings and hence needs not be identified with the momentum ${\bm{k}}$ transferred by a single collision. The typical situation, to be referred to as soft multiple scattering, is such that the number of quasi-independent scatterings is quite large, so that the typical value of $k$ is much smaller than the final momentum $p$.

In order to isolate the angular dependence of the $S$-matrix, one may be tempted to perform the small-$\Delta$ expansion as in Eq. (22) already before performing the Fourier transform. However, this manipulation, which corresponds to an expansion in powers of $\Delta^{2}/k^{2}$, would generate infrared divergences, leading to a result which is meaningless except for the leading order term, which has no angular dependence. For instance, to first order in $\Delta^{2}/k^{2}$, one finds

[TABLE]

Here we have assumed that $k\gg\Delta$, yet if one attempts to compute the above integral over ${\bm{k}}$ (for fixed ${\bm{\Delta}}$), one faces strong infrared ($k\to 0$) divergences, showing that this expansion is not really justified. To better see these divergences, notice that for sufficiently soft $k$ and $\Delta$, the ${\bm{r}}$-dependence within Eq. (27) can be expanded out:

[TABLE]

where the linear term in the r.h.s. vanishes (by parity) after the ${\bm{k}}$-integration.

Using the above, one sees that the dominant term $\propto 1/k^{4}$ in the large-$k$ expansion in Eq. (27) gives rise to a logarithmic integration for momenta ${\bm{k}}$ within the range $\Delta\ll k\ll 1/r$. This is a well-known result Iancu and Venugopalan (2003): the (angle-averaged) scattering amplitude $N(b,r)$ for a small dipole in the MV model is logarithmically sensitive to all transferred momenta within the interval $m\ll k\ll 1/r$, where $m\sim\Lambda_{{}_{\rm QCD}}$ is the infrared cutoff introduced in Eq. (16). In the present context, where the target is inhomogeneous, there is no genuine infrared divergence in the calculation — the associated momentum $\Delta$ effectively acts as an infrared cutoff on $k$ —, but we recover the logarithmic enhancement of the amplitude averaged over dipole orientations.

However, the second-order terms in the expansion in Eq. (27), which in particular carry an angular dependence, appear to develop a quadratic infrared divergence as $k\to 0$. This shows that this particular effect — the angular dependence of the dipole amplitude — is in fact controlled by soft exchanged momenta, $k\sim\Delta$, whose contribution cannot be computed via the expansion in powers of $\Delta^{2}/k^{2}$. Importantly, this also means that, for semi-hard momenta $p\lesssim Q_{s}(b)$, one cannot perform a reliable calculation of $v_{2}$ from first principles — not even within the limits of the MV model. Indeed, the soft momenta $k\sim\Delta\lesssim 1/R$ lie within the realm of the non-perturbative, confinement, physics, so their description within the MV model is not really justified. This being said, this model offers a convenient set-up for at least approaching the physics of the dipole orientation, while at the same time being consistent with the pQCD description of the angular-averaged amplitude. In that sense, we believe that the results of our subsequent analysis are still useful for qualitative and even semi-quantitative studies of the phenomenology.

We thus conclude that, for the present purposes, one cannot expand the double integral in Eq. (20) in powers of $\Delta^{2}/k^{2}$. Yet, the above discussion points towards another simplification, which is quite useful in practice: within the interesting regime of soft multiple scattering, all the relevant contributions come from relatively small transferred momenta $k\ll 1/r$, for which one can expand the $r$-dependence as shown in Eq. (28). This yields

[TABLE]

We have introduced here the infrared cutoff $m$ as a ‘gluon mass’ in the 2-dimensional Coulomb propagator,

[TABLE]

where ${\rm K}_{0}(x)$ is the modified Bessel function of the second kind. After this modification, the propagator $G({\bm{b}})$ shows an exponential decay at large transverse separations $mb\gg 1$, which mimics confinement. As already stressed, the insertion of this ‘mass’ $m$ is not required by infrared divergences: the integral over ${\bm{k}}$ in Eq. (29) is well-defined in the ‘infrared’ ($k\to 0$) even when $m=0$; and indeed, we shall later study the limit $m\to 0$ of our results. Rather, the ‘gluon mass’ $m$ is needed in order to restrict the phase-space allowed to very soft momenta $k\sim\Lambda_{{}_{\rm QCD}}$, which control the physics of the dipole orientation. (In real QCD, this phase-space would be of course restricted by confinement.) On the other hand, the integral over $k$ in Eq. (29) develops a logarithmic ‘ultraviolet’ ($k\to\infty$) divergence; it is understood that this divergence is cut off at the scale $k\sim 1/r$ (see Eq. (III.3) below for details).

It is also interesting to notice that the expansion (28) in powers of ${\bm{k}}\cdot{\bm{r}}$ ‘does not commute’ with the single scattering approximation studied in the previous section: in the latter, the exchanged momentum ${\bm{k}}$ is identified (via the Fourier transform) with the final momentum ${\bm{p}}$, hence ${\bm{k}}\cdot{\bm{r}}={\bm{p}}\cdot{\bm{r}}\sim\mathcal{O}{(1)}$ and a finite-order expansion in powers of ${\bm{k}}\cdot{\bm{r}}$ would be incorrect. Accordingly, a calculation using $S=\exp\{-N_{2g}\}$ together with Eq. (29) for $N_{2g}$ cannot reproduce the value of $v_{2}$ at very large momenta $p\gg Q_{s}(b)$ previously obtained in Eq. (25). More precisely, such a calculation would correctly reproduce the leading-order contribution $\propto 1/p^{4}$ to $\tilde{N}_{2g}({\bm{b}},{\bm{p}})$ in Eq. (24), which is independent of $\phi$, but not also its subleading piece $\propto 1/p^{6}$, which carries the interesting $\phi$-dependence.

The double integral in the r.h.s. of Eq. (29) has a relatively simple tensorial structure, which immediately implies that its result must be written as a linear combination of the following two rank-2 tensors: $\delta^{ij}$ and $b^{i}b^{j}/b^{2}$. Equivalently, the ensuing approximation for $N_{2g}({\bm{b}},{\bm{r}})\equiv N_{2g}(b,r,\theta)$ has the following generic structure

[TABLE]

without higher Fourier components. This is easily verified via direct calculations of the angular integrals in Eq. (29), which can be analytically completed. This is detailed in Appendix A, from which we here quote the final results:

[TABLE]

and respectively

[TABLE]

The above expression for $\mathcal{N}_{0}(b,r)$ has been obtained from Eq. (61) in Appendix A via the following manipulations: we have first subtracted the dominant behavior of the integrand at high $k$ and then replaced the subtracted piece via its following, regularized, form:

[TABLE]

with the impact-parameter dependent ‘saturation momentum’ $Q_{s}(b)$ defined as

[TABLE]

$Q^{2}_{0s}=\alpha_{s}C_{F}\mu_{0}$ is the central value of the saturation momentum at $b=0$. The coefficient of the logarithm $\ln(1/r^{2}m^{2})$ in the r.h.s. of Eq. (34) unambiguously follows from the logarithmic integration over the range $m\ll k\ll 1/r$, whereas the constant term under the log specifies our ‘renormalization’ scheme. Notice that all the results throughout this paper depend upon the QCD coupling $\alpha_{s}$, the fundamental Casimir $C_{F}$ and the 2-dimensional density $\mu_{0}$ of color charge squared only via this quantity $Q^{2}_{0s}$, to be treated as a free parameter of our model. In spite of our notations, $Q^{2}_{s}(b)$ is not exactly the saturation scale in the present model, but it is comparable to it, as we shall shortly argue.

The first piece in the r.h.s. of Eq. (III.3), proportional to $Q_{s}^{2}(b)$, would be the only one to survive in the case of a homogeneous target, i.e. when $\tilde{\mu}({\bm{\Delta}})\propto\delta^{(2)}({\bm{\Delta}})$. This piece has an apparent logarithmic divergence in the limit $m\to 0$. However, in the present context, where the target is inhomogeneous, this divergence is compensated by a corresponding divergence generated by the second, integral, term in Eq. (III.3). This is demonstrated in Appendix A, where we will also show that, when $m\to 0$, the mass parameter $m$ gets replaced by $1/R$ within the argument of the logarithm. This being said, the insertion of a non-zero ‘gluon mass’ $m$ is still necessary, on physical grounds.

The saturation momentum is more precisely defined by the condition that the scattering becomes strong: $N_{2g}(b,r,\theta)\sim 1$. This condition is controlled by the orientation-averaged piece $\mathcal{N}_{0}(b,r)$, which is numerically (much) larger than the piece $\mathcal{N}_{\theta}(b,r)$ encoding the angular dependence. This is manifest for sufficiently small dipoles $r\ll 1/m$, when the first piece in $\mathcal{N}_{0}(b,r)$ is enhanced by the large logarithm $\ln(1/r^{2}m^{2})$, but it is generally true for all the values of $r$ and $m$ of relevance for this work (see e.g. Fig. 2). The actual saturation momentum in the present set-up, to be denoted as $\mathcal{Q}_{s}(b)$, is conveniently defined by the condition $\mathcal{N}_{0}(b,r)=1$ when $r=2/\mathcal{Q}_{s}(b)$. This could be numerically extracted (as a function of $R$ and $m$), if needed, but for the present purposes it will be sufficient to use the following, qualitative, estimate, which strictly holds to leading-logarithmic accuracy:

[TABLE]

We have previously argued that the angular dependence of the dipole amplitude comes from relatively soft transferred momenta $k\sim\Delta$. It is interesting to check that at the level of Eq. (33). To this aim, let us take the limit $m\to 0$ in that equation. (The corresponding limit for $\mathcal{N}_{0}$ will be discussed in Appendix A.) Using

[TABLE]

it is easy to see that the expression within the square brackets inside the integrand becomes

[TABLE]

so that the whole contribution to $\mathcal{N}_{\theta}(b,r)$ indeed comes from soft momenta $k\leq\Delta/2$. As a matter of facts, the ensuing integral over $k$ is dominated by its upper limit $\Delta/2$ and the final result for $m=0$ takes a rather simple form:

[TABLE]

But albeit formally well defined, the limit $m\to 0$ of $\mathcal{N}_{\theta}(b,r)$ is physically meaningless, since very soft momenta $k\lesssim m$ are not allowed by QCD confinement. In that sense, the ‘massive’ case in Eq. (33) is more useful in practice, albeit our current treatment of confinement is merely heuristic. To illustrate the uncertainty introduced by this treatment, we represent in the left panel of Fig. 2 the result of the double integral in Eq. (33) as a function of $b$ for several values of $m$, including $m=0$. As one can see there, the $m$-dependence become stronger with increasing $b$, a feature which is easy to understand: the integral over $k$ is effectively restricted to values $m\lesssim k\lesssim\Delta$ and larger values for the impact parameter $b$ correspond to smaller values for $\Delta$.

In Fig. 2, one also sees that the function $\mathcal{N}_{\theta}(b,r)$ develops a maximum at a value of $b$ which is proportional to $R$ and roughly independent of $m$. For $m=0$, Eq. (39) shows that $\mathcal{N}_{\theta}\propto b^{2}$ at small $b\ll R$ and we expect this to remain true for any value of $m$. Another interesting aspect of the dipole amplitude in Eq. (39) is the fact that it exhibits a power tail $\propto 1/b^{2}$ at sufficiently large distances $b\gtrsim R$. This is in agreement with the discussion after Eq. (19): it reflects the fact that the angular dependence of the dipole amplitude is controlled by soft gluon exchanges, for which there is no confinement in the limit $m\to 0$. For a non-zero ‘gluon mass’ $m$, this power law tail will of course be replaced by the decaying exponential ${\rm e}^{-mb}$, which mimics confinement. This behavior too is visible in Fig. 2.

For more clarity, we also plot the angular-averaged amplitude $\mathcal{N}_{0}(b,r)$, under the same assumptions as for $\mathcal{N}_{\theta}(b,r)$ (see the right panel of Fig. 2). The fact that for small values of the ‘gluon mass’ $m\leq 2$ GeV, the maximum of $\mathcal{N}_{0}(b,r)$ as a function of $b$ appears to be displaced at non-zero values for $b$ is probably an artifact of the model. But this also shows that the second, integral, term in the r.h.s. of Eq. (III.3) is indeed important (by itself, this contribution is negative for sufficiently small values of $b$, but it becomes positive at $b\gtrsim R$). This feature will have no incidence on our subsequent numerical studies of $pp$ collisions, where we shall restrict ourselves to larger values $m\geq 0.3$ GeV. For such values, the maximum of $\mathcal{N}_{0}(b,r)$ is located at $b=0$, as expected on physical grounds.

Using the above results for $N_{2g}$ together with $S=\exp\{-N_{2g}\}$ and the representation (6) for the elliptic flow coefficient $v_{2}$, we finally deduce the following estimate for the latter:

[TABLE]

where ${\rm I}_{0}$ and ${\rm I}_{1}$ (the modified Bessel functions of the first kind) have been generated via

[TABLE]

${\rm I}_{1}(z)$ is an odd function which has the same sign as its argument. In fact, the quantity $\mathcal{N}_{\theta}$ is numerically small in the physical regime of interest (see Fig. 2), hence one can use the approximation ${\rm I}_{1}\left(\mathcal{N}_{\theta}(b,r)\right)\simeq\mathcal{N}_{\theta}(b,r)/2$. This shows that $v_{2}(p,b)$ is significantly large only for peripheral collisions, i.e. for impact parameters $b\gtrsim R$, where lies the peak of the function $\mathcal{N}_{\theta}(b,r)$. It furthermore shows that the elliptic flow generated via multiple scattering is positive 111So long as $p$ is not too large, $p\lesssim Q_{s}(b)$, the above integrals over $r$ are effectively cut off by the Gaussian ${\rm e}^{-\mathcal{N}_{0}(b,r)}$ at a value $r\sim 1/Q_{s}(b)$. Hence, within the range for $r$ which is relevant for the integration, the Bessel functions ${\rm J}_{2}(pr)$ and ${\rm J}_{0}(pr)$ remain positive, meaning that the sign of $v_{2}$ coincides with that of the function $\mathcal{N}_{\theta}(b,r)$ defined in Eq. (33). — that is, it has the opposite sign as compared to the case of a single hard scattering discussed in Sect. III.2.

Via Eq. (31), the sign of $v_{2}$ can be related to properties of dipole scattering. Namely, the fact that $\mathcal{N}_{\theta}(b,r)$ is positive implies that the scattering is stronger when the dipole orientation is (anti)parallel to its impact parameter ($\theta=0$ or $\theta=\pi$) than for a dipole perpendicular on ${\bm{b}}$ ($\theta=\pi/2$). Equivalently, the $S$-matrix $S({\bm{b}},{\bm{r}})$, which measures the dipole survival probability, is larger when ${\bm{r}}\perp{\bm{b}}$ than for ${\bm{r}}\parallel{\bm{b}}$ (see Fig. 3). This property is studied in more detail in Fig. 4: in the left panel, we show the dipole $S$-matrix as a function of $\theta$ (for two different dipole sizes and a fixed impact parameter); on the right, we show the scattering amplitude $N_{2g}(b,r,\theta)$ as a function of $b$ for a fixed value of $r$ and two extreme orientations: $\theta=0$ and $\theta=\pi/2$. As one can see, the difference between ‘parallel’ and ‘perpendicular’ scattering increases with the dipole size $r$ and also with the impact parameter $b$. These features are intuitively understandable since a point-like dipole should not be sensitive to its orientation. Besides for very small impact factors $b\lesssim R=0.3$ fm, the target looks quasi-homogeneous and then the dipole orientation is irrelevant. We therefore expect the associated $v_{2}$ to follow a similar trend. This will be confirmed by the numerical results to be presented in Sect. IV.

Returning to the case of the single scattering approximation, as applying at high $p\gg Q_{s}(b)$, it might be tempting to interpret the negativity of $v_{2}$ in that case as an opposite trend for the dipole scattering, namely $S({\bm{r}}\parallel{\bm{b}})>S({\bm{r}}\perp{\bm{b}})$. However, we believe that such an interpretation is truly misleading: in that case, the sign of $v_{2}$ follows from an analysis that was performed fully in momentum space. Such an analysis gives one information about the unintegrated gluon distribution in the target (proportional to $-\tilde{N}_{2g}({\bm{b}},{\bm{p}})$, cf. Eq. (24)), but not about the dipole scattering as a function of ${\bm{r}}$. To compute the latter, i.e. the function $N_{2g}({\bm{b}},{\bm{r}})$, one needs its Fourier transform $\tilde{N}_{2g}({\bm{b}},{\bm{p}})$ for all values of ${\bm{p}}$, and not just for the relatively hard values for which Eq. (24) applies. In fact, even for small values of $r$, the angular dependence of $N_{2g}({\bm{b}},{\bm{r}})$ is controlled by relatively soft values of $p$ within the inverse Fourier transform $\tilde{N}_{2g}({\bm{b}},{\bm{p}})\to N_{2g}({\bm{b}},{\bm{r}})$ (cf. the discussion following Eq. (27)).

Finally, let us generalize the previous results to the case where the proton projectile itself has a Gaussian distribution in the transverse plane, $x_{p}q(x_{p})\to x_{p}q(x_{p})f({\bm{b}})$, with $f({\bm{b}})={\rm e}^{-b^{2}/4R^{2}}/(4\pi R^{2})$. Using this Ansatz for $f({\bm{b}})$ together with the expression (31) for the dipole amplitude, one can easily perform the angular integrations in Eq. (10) for $v_{2}$ and thus obtain (the identity $\int_{0}^{2\pi}{\rm d}\alpha\,{\rm e}^{-z\cos\alpha}\,\cos(2\alpha)=2\pi{\rm I}_{2}(z)$ is also useful)

[TABLE]

We recall that the ‘dummy’ variable $b$ is the impact parameter of a participating quark, whereas the external variable $B$ refers to the center of the projectile. The integral over $b$ in the numerator of Eq. (42) is restricted by the support of the function $\mathcal{N}_{\theta}(b,r)$, cf. Fig. 2, hence it receives most of its contribution from relatively large values $b\gtrsim R$. For the nearly central proton-proton collisions with $B\ll R$, the overall elliptic flow is negligible, by rotational symmetry: the individual contributions from various (peripheral) values of $b$ can have any orientation, so they compensate each other. Indeed, using ${\rm I}_{2}(x)\simeq x^{2}/8$ for $x\ll 1$, it is easy to see that $v_{2}(p,B)$ vanishes as $B^{2}$ when $B\to 0$. But for larger impact parameters $B\gtrsim R$, the rotational symmetry of the interaction region is badly broken (recall Fig. 1) and one expects a non-trivial net result. Geometrical considerations suggest that $v_{2}(p,B)$ should be proportional to the eccentricity of the overlapping region, as defined in Eqs. (11) or (13), which can be more explicitly written as

[TABLE]

and respectively

[TABLE]

Indeed, one can understand this eccentricity as the expectation value $\varepsilon_{2}=\langle\cos(2\alpha)\rangle$, where we recall that $\alpha$ is the angle made by the impact parameter ${\bm{b}}$ of an individual quark with respect to that, ${\bm{B}}$, of the center of the projectile (see Fig. 1). Hence larger values for $\varepsilon_{2}$ imply that all participating quarks have similar impact parameters, hence the produced particles are preferentially produced along a common direction — that of ${\bm{B}}$ —, thus generating a sizable value for the elliptic flow. And indeed, by inspection of the equations above, it is clear that both $v_{2}(p,B)$ and $\varepsilon_{2}(p,B)$ [or $\varepsilon_{2}(B)$] are proportional to $B^{2}$ so long as $B\lesssim R$, hence they are proportional to each other. This relation between the elliptic flow and the eccentricity will be further investigated in Sect. IV, where all these quantities will be numerically computed.

III.4 Dipole-nucleus scattering: the case of a lumpy target

The most straightforward generalization of the previous set-up to the case where the target is a large nucleus with atomic number $A\gg 1$ would be obtained by assuming that the valence color charges (and hence the associated gluon distribution) are uniformly distributed throughout the nuclear volume — the so-called “smooth nucleus”. Experience with nuclear physics at lower energies suggests that a reasonable approximation for the 3-dimensional distribution of the nuclear matter within a large nucleus is provided by the Woods-Saxon distribution $\rho_{A}(\vec{r})$. By boosting this distribution and assuming that it also applies to the valence color charges, we conclude that the case of a “smooth nuclear target” can be obtained by replacing the 2-dimensional density $\mu(b)$ in Eq. (15) as follows

[TABLE]

where $T_{A}(b)$ is the nuclear thickness function,

[TABLE]

where $R_{A}=(1.12\,{\rm fm})A^{1/3}$ is the nuclear radius and $\delta=0.54$ fm is the width of the “nuclear edge” (the radial distance across which the nuclear density is rapidly dropping). The quantity $\mu_{0}$ has the same meaning as before — the color charged squared for the valence quarks of the nucleon per unit transverse area —, hence it is independent of $A$. The overall factor of $A$ visible in Eq. (45) reflects the fact that the density $\rho_{A}$ is normalized to unity: $\int{\rm d}^{3}\vec{r}\rho_{A}(\vec{r})=1$. This in turn implies that the normalization factor $N_{A}$ scales like $1/A$, hence $T_{A}\propto 1/A^{2/3}$ and the color charge density therefore has the canonical scaling with the number $A$ of nucleons: $\mu_{A}(b)\propto A^{1/3}$.

Under the above assumption, the formal calculation of the dipole $S$-matrix would proceed in the same way as before, leading to expressions similar to those already presented in Eqs. (26), (31), (III.3) and (33). The ensuing numerical evaluation however would likely lead to considerably smaller values for $v_{2}$, due to combined effect of the larger value for the nuclear radius $R_{A}$ and to the fact that the Woods-Saxon profile is less rapidly varying with $b$ than the Gaussian.

This being said, it is quite clear that a real nucleus is not homogeneous; rather, it is a lumpy superposition of $A$ distinct nucleons and this lumpiness is known to have important consequences for the phenomenology. In particular, it can generate a privileged direction of motion for the produced particles (for a given impact parameter), via the following mechanism: the effective dipole, with a given orientation and size ${\bm{r}}$, will scatter off the nucleon which happens to be located at the dipole impact parameter ${\bm{B}}$. (From now on, we shall use ${\bm{B}}$ to denote the impact parameter of the dipole w.r.t. the center of the nucleus, and keep ${\bm{b}}$ for its impact parameter w.r.t. the struck nucleon.) As a result, the produced quark will preferentially move along the direction of the local impact factor ${\bm{b}}={\bm{B}}-{\bm{b}}_{i}$, where ${\bm{b}}_{i}$ is the position of the struck nucleon w.r.t. the center of the nucleus. If nucleons are randomly distributed around the given ${\bm{B}}$, then the information about the orientation of the produced particle will be washed out after averaging over the nucleon distribution. For large values of $A$, this will likely be the case at impact parameters ${\bm{B}}$ deeply inside the nucleus, where the nuclear distribution is quasi-homogeneous. But even in that case, this cannot happen for impact parameters close to the periphery ($B\sim R_{A}$), which will therefore generate nonzero contributions to $v_{2}$. These qualitative considerations will be confirmed via an explicit calculation to which we now turn.

For a given configuration of the $A$ nucleons inside the nucleus and assuming the dipole to independently scatter off any of them, the dipole $S$-matrix should read (see Kowalski and Teaney (2003) for a more complete discussion)

[TABLE]

For simplicity, we shall further assume that the various nucleons are distributed independently from each other; for each of them, its central position ${\bm{b}}_{i}$ is distributed according to the Woods-Saxon thickness function $T_{A}({\bm{b}}_{i})$. The physical observable is then obtained by averaging over all possible configurations of the nucleons, as follows

[TABLE]

The most interesting regime, including for the phenomenology of $pA$ collisions at the LHC, is such that the scattering between the dipole and a single nucleon is weak, $N_{2g}({\bm{b}},{\bm{r}})\ll 1$, yet the overall scattering can be strong (meaning that the $S$-matrix can be small compared to unity: $S_{A}({\bm{B}},{\bm{r}})\ll 1$). Under these assumptions, one can expand the exponential ${\rm e}^{-N_{2g}}$ to lowest non-trivial order, perform the integral over ${\bm{b}}$ and then re-exponentiate the result, to finally obtain (recall the normalization condition $\int{\rm d}^{2}{\bm{b}}\,T_{A}({\bm{b}})=1$)

[TABLE]

with the following definition for the dipole-nucleus scattering amplitude in the two-gluon exchange approximation (divided by the number $A$ of nucleons)

[TABLE]

The above integral over ${\bm{b}}$ is effectively restricted (by the support of the dipole-proton amplitude $N_{2g}({\bm{b}},{\bm{r}})$) to the area of the proton disk, which is small compared to that of the nucleus. In other terms, one has $b\ll B$ for the most interesting values $B\sim R_{A}$. In view of this, one may be tempted to approximate $T_{A}({\bm{B}}-{\bm{b}})\simeq T_{A}({\bm{B}})$, as often done in the literature Kowalski and Teaney (2003). However, this approximation would wash out the information about the dipole orientation, which is important for us here. To keep trace of this information, one needs to go one step further in the small $b/B$ expansion, namely up to quadratic order (the linear term does not contribute to the integral over ${\bm{b}}$, by parity). We thus write

[TABLE]

Plugging the above expansion and the generic form of $N_{2g}(b,r,\theta)$ given in Eq. (31) into Eq. (50), one can easily perform the integral over the angle $\theta$ between ${\bm{b}}$ and ${\bm{r}}$ and thus obtain (from now on, we use $\theta$ to denote the angle made by the vectors ${\bm{r}}$ and ${\bm{B}}$)

[TABLE]

where

[TABLE]

The $\theta$-dependent piece $\mathcal{N}^{A}_{\theta}(B,r)$ is proportional to the (first and second) derivatives of the thickness function $T_{A}(B)$, hence its support is limited to values of $B$ near the edge of the nucleus, within a distance $\Delta B\sim\delta$ around $B=R_{A}$ (see Fig. 5). This is in agreement with our previous physical discussion and confirms that the mechanism under consideration can generate a sizable $v_{2}$ only in peripheral $pA$ collisions. As also illustrated in Fig. 5, the special combination $T^{\prime\prime}_{A}(B)-T^{\prime}_{A}(B)/B$ which enters $\mathcal{N}^{A}_{\theta}(B,r)$ is positive for most values of $B$ within its support. (It can become slightly negative at intermediate values of $B$, but the corresponding value $v_{2}(p,B)$ is anyway very small, as we shall see.) Together with the positivity of the respective proton amplitude $\mathcal{N}_{\theta}(b,r)$, as numerically observed in Sect. III.3, this implies $\mathcal{N}^{A}_{\theta}(B,r)>0$. That is, as in the case of a proton target, the scattering is stronger when the dipole orientation ${\bm{r}}$ is (anti)parallel to its nuclear impact factor ${\bm{B}}$, rather than perpendicular on it.

At this stage, one could use the integral representations for the functions $\mathcal{N}_{0}(b,r)$ and $\mathcal{N}_{\theta}(b,r)$, as given in Eqs. (III.3) and respectively (33), to numerically perform the integrals in Eq. (53). This would amount to computing a sequence of three radial integrations, with integrands involving the oscillatory Bessel functions. This is indeed possible in practice, but rather tedious and very time-consuming. It turns out that this whole calculation can be efficiently reorganized, in such a way to provide fully analytic results for the nuclear amplitudes $\mathcal{N}^{A}_{0}(B,r)$ and $\mathcal{N}^{A}_{\theta}(B,r)$. This is explained in Appendix B, from which we show the final results (whose general structure is indeed consistent with Eq. (53)):

[TABLE]

The first line in the r.h.s. of Eq. (54) for $\mathcal{N}^{A}_{0}$, which is proportional to the large logarithm $\ln({1}/{r^{2}m^{2}})$, represents the dominant contribution to the dipole amplitude. Its present calculation within the MV model is indeed under control (at least for sufficiently small dipole sizes $r\ll 1/m$), since this contribution is dominated by relatively large exchanged momenta, $m\ll k\ll 1/r$. Within that contribution, the dominant piece is the one proportional to $T_{A}(B)$. This argument shows that the nuclear saturation momentum $Q_{sA}$ (the inverse dipole size for which the dipole amplitude $AN_{2g}^{A}({\bm{B}},{\bm{r}})$ becomes of order one) can be estimated as

[TABLE]

This scales like $A^{1/3}$ and has the same impact-parameter dependence as the nuclear thickness function $T_{A}(B)$.

The $\theta$-dependent piece $\mathcal{N}^{A}_{\theta}$ of the amplitude, which is the most interesting one for the present purposes, is proportional to $1/m^{2}$, which demonstrates its non-perturbative origin: it has been generated by integrating over soft momenta $k\sim m$. In that sense, the above calculation is merely heuristic and in particular strongly dependent upon our recipe for implementing the infrared cutoff $m$. At this point, one may wonder about the difference between the small-$m$ behaviors observed in $pA$ and respectively $pp$ collisions: when $m\to 0$, the $\theta$-dependent piece of the dipole amplitude remains finite for $pp$ collisions, cf. Eq. (39), whereas it is quadratically divergent in the case of $pA$ collisions, cf. Eq. (55). This difference can be traced back to the integral over $b$ which needs to be performed when passing from $pp$ to $pA$ collisions, cf. Eq. (53). When $m\to 0$ and for large $b\gtrsim R$, the respective amplitude for $pp$ collisions has a power tail $\mathcal{N}_{\theta}(b,r)\propto 1/b^{2}$, as visible in Eq. (39). Therefore, the integral $\int{\rm d}b\,b^{3}\,\mathcal{N}_{\theta}(b,r)$ which enters Eq. (55) for $\mathcal{N}^{A}_{\theta}$ would be quadratically divergent in the absence of confinement. After adding the latter in the form of a ‘gluon mass’ $m$, this integral is cutoff at $b\sim 1/m$, thus yielding $\mathcal{N}^{A}_{\theta}\propto 1/m^{2}$.

We are finally in a position to compute the elliptic flow coefficient $v_{2}$ for $pA$ collisions: by inserting the dipole $S$-matrix obtained according to Eqs. (49) and (52) into our master formula (6), we obtain, similarly to Eq. (40),

[TABLE]

This can be numerically computed using Eqs. (54) and (55), with results to be discussed in the next section. The generalization of Eq. (57) to an extended projectile is straightforward and will be considered too in Sect. IV.

IV Numerical results for $v_{2}$ and physical discussion

In this section we present the numerical results for $v_{2}$ in $pp$ and $pA$ collisions (with $A=208$) as emerging from our model. For more clarity, in the following (and in all plots) we shall denote the transverse momentum with $p_{T}$. We shall limit ourselves to the scenario which includes the effects of multiple scattering, as discussed in Sect. III.3 for $pp$ collisions and in Sect. III.4 for the $pA$ collisions. For both cases, $pp$ and $pA$ collisions, we shall exhibit $v_{2}$ as a function of the transverse momentum $p_{T}$ of the produced quark, for various choices of the impact-parameter, the central saturation momentum in the proton $Q_{0s}$, and the infrared cutoff $m$. The only other parameter of our model, i.e. the width $R$ of the proton color charge distribution in the transverse space, is fixed to the average value emerging from a fit to $t$-distribution of diffractive vector meson production at HERA, that is $R^{2}=2\,\text{GeV}^{-2}$. Strictly speaking, such a fit is based on a different ‘saturation model’, namely IP-Sat Rezaeian et al. (2013) but this difference is not essential for the subsequent discussion, which will be mostly qualitative. Note also that the value of $R^{2}$ extracted using the bCGC model in a fit to the same data Rezaeian and Schmidt (2013) is slightly larger.

Concerning $Q_{0s}$ — the proton saturation momentum at $b=0$ —, we shall consider a rather wide range of values, from $Q_{0s}^{2}=0.2-0.4\ \text{GeV}^{2}$ up to $Q_{0s}^{2}=2.4\ \text{GeV}^{2}$. The lowest values emerge from phenomenological analyses based on the Balitsky-Kovchegov equation with running coupling (rcBK) to either the HERA data Albacete et al. (2011); Iancu et al. (2015a); Albacete (2017), or to the $pp$ data at RHIC and the LHC Jalilian-Marian and Rezaeian (2012); Albacete et al. (2013); Rezaeian (2013). The highest value could in principle be reached in high-multiplicity events characterized by large fluctuations Hatta et al. (2006); Avsar et al. (2011). (Notice that the fits to HERA data in Iancu et al. (2015a); Albacete (2017) use a more complete version of the BK equation which besides a running coupling, also includes collinear improvement Beuf (2014); Iancu et al. (2015b, a).)

We are now prepared to present our numerical results, starting with $pp$ collisions. As stressed in the Introduction, we have in mind an asymmetric situation, where one of the protons (‘the target’) looks dense and can be described by the MV model, while the other one (‘the projectile’) is dilute. This might be the case for particle production at very forward rapidities and also for rare, high-multiplicity, events in which the target proton develops ‘hot spots’ via fluctuations in the high-energy evolution Hatta et al. (2006); Avsar et al. (2011). We first show our results for the idealized case of a point-like projectile, cf. Eq. (40), and then for the more realistic case of a projectile which has a non-trivial extent in the transverse plane, cf. Eq. (42).

In Fig. 6 we show the azimuthal asymmetry $v_{2}$ computed according to Eqs. (III.3), (33) and (40) for different choices of the central saturation scale $Q_{0s}$. These plots illustrate the scaling of the peak position $p_{\rm max}$ with $Q_{0s}$ : when we plot $v_{2}$ as a function of $p_{T}/Q_{0s}$, the peak position $p_{\rm max}/Q_{0s}$ is quasi independent of $Q_{0s}$ and rather close to 1/2. This scaling property indicates the importance of the saturation physics. A larger saturation scale shifts the unintegrated gluon distribution (the integrand of $v_{2}$) to higher transverse momenta. In Fig. 7 we show the dependence of the azimuthal asymmetry $v_{2}$ upon the impact parameter $b$ (left panel) and upon the infrared cutoff $m$ (right panel). As expected, the strength of $v_{2}$ is increasing with $b$. Remarkably though, one see that quite large peak values $v_{2}(p_{\rm max},b)\gtrsim 0.15$ are obtained already for not so large impact parameters, $b\lesssim 0.3$ fm, that is, for collisions which are peripheral, but not ultra-peripheral. (Recall that the typical transverse size of the color charge distribution in the target is $2R\sim 0.6$ fm.) It is also interesting to notice that, albeit the height of the peak is rapidly increasing with $b$, its position $p_{\rm max}$ changes only slightly when going from rather central ($b=0.1$ fm) to more peripheral ($b=0.3$ fm) values. This observation should be correlated with the fact that, as manifest on Eq. (39), the piece $\mathcal{N}_{\theta}(b,r)$ of the amplitude which is responsible for the angular dependence is proportional to the central value $Q_{0s}$ of the saturation scale, and not to its local value at the actual impact parameter. As anticipated, the $m$-dependence is quite strong: a slight increase in $m$, from 0.3 GeV to 0.4 GeV, reduces the peak value of $v_{2}$ by a factor 3.

Turning now to an extended projectile with a Gaussian distribution in impact parameter space, the corresponding $v_{2}$ is shown in Fig. 8, for various values for $Q_{0s}^{2}$ and $B$. (We have checked that the $m$-dependence of the results is similar to that observed for a point-like projectile, cf. Fig. 7.) One may expect the strength of the azimuthal asymmetry to be reduced, perhaps even significantly, after averaging over the surface of the projectile, but this is actually not the case: as visible in Fig. 8, the peak value of $v_{2}$ remains as large as for a point-like projectile. To see such a sizable $v_{2}$, however, one needs to go to larger values for the impact parameter $B$, which now refers to the center of the projectile (recall Fig. 1). This is in agreement with the discussion at the end of Sect. III.3, which also suggests that the value of $v_{2}$ should be correlated to the eccentricity $\varepsilon_{2}$ of the interaction region.

To check this conjecture, we have numerically computed $\varepsilon_{2}(p_{T},B)$ and $\varepsilon_{2}(B)$ according to Eqs. (43)–(44), with the results shown in Fig. 9 (left panel). These results should be compared to the $B$-dependence of $v_{2}$, as exhibited in the right panel of the same figure. These plots confirm that $v_{2}$ and $\varepsilon_{2}$ show a similar trend with $B$: they monotonously increase with $B$ — actually, they are both proportional to $B^{2}$ so long as $B$ is small enough, $B\lesssim R$. On the other hand, they show rather different behaviors with $p_{T}$. The plots for $v_{2}(p_{T},B)$ in the right panel of Fig. 9 are in agreement with those in the left panel of Fig. 8: $v_{2}$ vanishes as $p_{T}\to 0$ and has a pronounced peak at $p=p_{\rm max}$ with $0.5<p_{\rm max}<1$ GeV. On the other hand, $\varepsilon_{2}$ has a rather weak dependence upon $p_{T}$: the curves corresponding to different values for the momentum are rather close to each other, and also to the curve representing the integrated eccentricity $\varepsilon_{2}(B)$. This reflects the fact that the quantity $\varepsilon_{2}(p_{T},B)$ is only weakly sensitive to the dipole scattering, since mostly controlled by the geometry.

We now turn to the case of $pA$ collisions, for which the present formalism is somewhat better justified. The respective $v_{2}$ is computed by numerically integrating Eq. (57) with the dipole amplitude given by the analytic results in Eqs. (54) and (55). The two plots in Fig. 10, which exhibit $v_{2}$ as a function of $p_{T}$ (left panel) and $p_{T}/Q_{0s}$ (right panel), for different values of the proton saturation scale $Q_{0s}$, are quite similar to the corresponding plots for $pp$ collisions, cf. Fig. 6. In particular, the peak position appears to respect the expected scaling with the nuclear saturation momentum $Q_{sA}=A^{1/6}Q_{0s}$ : indeed, the maximum occurs at, roughly, $p_{\rm max}/Q_{0s}\simeq 1.2$, which is larger by a factor $A^{1/6}\simeq 2.4$ (for $A=208$) than the respective value observed for $pp$ collisions. However, in order to reach values for $v_{2}$ which are comparable to those in $pp$ collisions, one now needs to go up to much larger values of the impact parameter $B\gtrsim R_{A}$, there the inhomogeneity in the nuclear distribution is located (cf. the discussion in Sect. III.4). The $B$-dependence of the function $v_{2}(p)$ is illustrated in the left panel of Fig. 11. This is controlled by the combination $T^{\prime\prime}_{A}(B)-T^{\prime}_{A}(B)/B$, cf. Eq. (53), and the results in Fig. 11 are indeed in agreement with the previous discussion of Fig. 5. Namely, $v_{2}$ is seen to be sizable and positive for all values $B\geq R_{A}\approx 6.5$ fm.

Notice that in the present approximations, the dipole amplitude (hence our estimate for $v_{2}$) for the case of a nuclear target depends upon the two scales $R^{2}$ and $Q_{0s}^{2}$ mostly via their product $R^{2}Q_{0s}^{2}$. (This becomes obvious by inspection of Eqs. (54) and (55).) Accordingly, the effect of increasing $Q_{0s}^{2}$ at fixed $R^{2}$, as visible in Fig. 10, can alternatively be associated with increasing $R^{2}$ for a fixed value $Q_{0s}^{2}$. In the right panel of Fig. 11, we show the dependence of $v_{2}$ in $pA$ collisions upon the infrared cutoff $m$. Similarly to the case of $pp$ collisions, one finds that this dependence is rather strong: by decreasing $m$ from the ‘confinement’ value $m=0.3$ GeV to the pion mass $m=0.14$ GeV, one increases the peak value of $v_{2}$ by a factor of 3.

In Fig. 12 we illustrate the effect of using an extended proton projectile, with Gaussian distribution in impact parameter. The corresponding formula for $v_{2}(p_{T},B)$ is the straightforward generalization of Eq. (42), obtained by replacing $\mathcal{N}_{0}(b,r)\to A\mathcal{N}^{A}_{0}(b,r)$ and $\mathcal{N}_{\theta}(b,r)\to A\mathcal{N}^{A}_{\theta}(b,r)$. As visible in Fig. 12, the effect is quite small — at most a change of 20% in the value of $v_{2}$ at its peak.

The systematics of the above results for $v_{2}$ can be physically understood as follows. First of all, we found that $v_{2}(p_{T})$ is small for both very small and very large values of $p_{T}$, but has a maximum at some intermediate value $p_{\rm max}$. In particular, $v_{2}=0$ for $p_{T}=0$, as already obvious by inspection of Eq. (6). These features are easy to understand: the angular orientation cannot play any role when either the momentum $p_{T}$, or the dipole size $r$, are too small. Since typically $r\sim 1/p_{T}$, the second argument explains the rapid decrease of $v_{2}$ that we observe at high $p_{T}$. But the detailed shape of the function $v_{2}(p_{T})$ — in particular the position, the width, and the height of its maximum — are strongly dependent upon the impact parameter and also upon the values of the 3 parameters $Q_{0s}$, $m$, and $R$.

Specifically, as visible in the left panels in both Fig. 7 and Fig. 11, $v_{2}$ is negligible for relatively small impact parameters (in particular, it vanishes as $b,B\to 0$), but it becomes large — in the sense that it reaches a peak value $v_{2}(p_{\rm max})\gtrsim 0.1$ — when the impact parameter is comparable to the typical size for inhomogeneity in the target, that is, $b\sim R\gtrsim 0.2$ fm for a proton and respectively $B\sim R_{A}\gtrsim 6.5$ fm for a large nucleus. This is understandable, given that the angular orientation would play no role for a target which is homogeneous in impact-parameter space. We recall that, for the mechanism under consideration, the elliptic flow is driven by the sensitivity of the color-dipole orientation to the variation in the gluonic or nuclear distribution in the transverse plane.

We furthermore found that the peak in $v_{2}(p_{T})$ moves towards larger values of $p_{T}$ and becomes broader when increasing $Q_{0s}^{2}$, see Fig. 6 and Fig. 10. This is as expected: the larger saturation momentum in the target, the larger is the typical momentum of the produced parton and the wider is its distribution in $p_{T}$. Interestingly, for both $pp$ and $pA$ collisions we found that the position $p_{\rm max}$ of the peak in $v_{2}(p_{T})$ is proportional to $Q_{0s}$. A similar observation was recently made in Ref. Hagiwara et al. (2017). When $v_{2}(p_{T})$ is plotted as a function of $p_{T}/Q_{0s}$, the peak position $p_{\rm max}/Q_{0s}$ is quasi-independent of $Q_{0s}$, albeit its height and shape are still strongly dependent (see the right panels in Fig. 6 and Fig. 10). Specifically, the maximal value at the peak $v_{2}(p_{\rm max})$ appears to increase when decreasing $Q_{0s}$, i.e. when the target becomes more dilute. This may seem counter-intuitive since, as already stressed, the multiple scattering represents an essential ingredient of the mechanism under consideration (it even changes the sign of $v_{2}$ as compared to the single-scattering approximation). However, the importance of the dipole orientation depends in a crucial way upon the balance between the dipole size and the size of its impact parameter. The dipole size is fixed by the transverse momentum of the produced quark, $r\sim 1/p_{T}$, which in turn is determined by the target saturation momentum: $p_{T}\sim Q_{s}(b)\propto Q_{0s}$. Hence, if one keeps increasing $Q_{0s}$, the dipole size eventually becomes much smaller than $b$ and the dipole orientation plays no role anymore. A similar effect is seen when the saturation momentum increases as a consequence of the high-energy evolution Kovner and Lublinsky (2013); Lappi (2015); Lappi et al. (2016).

Finally, given the importance of soft, non-perturbative, exchanges for the angular-dependence of the dipole amplitude, is should be no surprise that our results for $v_{2}$ are rather strongly dependent to the ‘confinement’ scale $m$: the anisotropy is enhanced when decreasing $m$, since the phase-space for soft exchanges is rapidly increasing, see the right panels in Fig. 7 and Fig. 11. For the angular-dependent piece of the dipole amplitude and for a proton target, this dependence has been already exhibited in Fig. 2, whereas for a nuclear target, it is directly visible by inspection of Eq. (55) for $\mathcal{N}^{A}_{\theta}$.

From the previous considerations in this paper, it should be clear that our current analytic description for the mechanism under consideration is too crude to allow for quantitative predictions, or realistic applications to the phenomenology. This being said, we would like to show via an example that this scenario is not excluded by the current data. Namely, we will show that, by appropriately choosing the values of the impact parameter $B$ and of the free parameters of the models, one can give a reasonable description of the $p_{T}$-dependence of the elliptic flow $v_{2}(p_{T})$ extracted from multi-particle azimuthal correlations in p+Pb collisions at the LHC, in a given multiplicity class. This should not be confounded with a genuine fit to the data — it is merely an exploratory comparison. Given the uncertainties inherent in our model, we shall adopt a rather crude strategy for relying the predictions of this model to the phenomenology.

First, we have not attempted to compute the consequences of our mechanism for multi-particle azimuthal correlations; rather, we shall make the simplifying assumption that the final particles are correlated with each other only through the flow correlations with the reaction plane (that is, we neglect possible ‘non-flow correlations’). Under this assumption, we can write $c_{2}\{2\}\simeq v_{2}^{2}$, and similarly $c_{2}\{4\}\simeq-v_{2}^{4}$, where $c_{2}\{2\}$ and $c_{2}\{4\}$ are the second-order two-particle and respectively four-particle cumulants, as defined e.g. in Borghini et al. (2001); Snellings (2011), and $v_{2}$ is the usual elliptic flow coefficient, as discussed throughout this paper. More precisely, our present estimates for $v_{2}$ refer to a fixed impact parameter $B$, that is, we have a prediction for the function $v_{2}(p_{T},B)$, whereas the data for p+Pb (and also p+p) collisions are rather classified according to the particle multiplicity (the number of reconstructed tracks) in the final state. Due to the large multiplicity fluctuations, the correlation between the multiplicity classes and the cuts in impact parameter is rather loose and not well under control (see e.g. the discussion in Adam et al. (2015)). To cope with that and by lack of any better alternative, we shall simply select a value of $B$ for which the predictions of our model appear to reasonably agree with the data in a given high-multiplicity class. We shall similarly proceed with the free parameters of our model, which in the case of $pA$ collisions are the dimensionless product $R^{2}Q_{0s}^{2}$ and the infrared cutoff $m$. It is understood that the values for $B$, $m$ and $R^{2}Q_{0s}^{2}$ which will emerge from this procedure must also reflect the influence of other uncertainties or simplifications inherent in our present approach, like the omission of gluons. (Gluon production within the present set-up should give rise to a similar $p_{T}$-dependence in $v_{2}$, but only modify its overall magnitude.)

Finally, to have a better comparison with the $p_{T}$-dependence of the data, which refer to hadrons, one must take into the account the effect of the quark (or gluon) fragmentation into hadrons. This too will be implemented in a rather heuristic way, by assuming that $p_{\text{hadron}}=\langle z\rangle p_{\text{quark}}$ where $\langle z\rangle$ is the average value of the splitting fraction and will be treated as a free parameter.

In Fig. 13, we compare our results for $v_{2}$ in peripheral $pA$ collisions (with $A=208$) with the experimental data from the CMS experiment for the second-order elliptic harmonic $v_{2}\{4\}$ in p-Pb collisions, as extracted from a four-particle cumulant analysis in events where the number of reconstructed tracks lies in the range $120<N_{\rm trk}^{\rm off}<150$ Chatrchyan et al. (2013b). The two theoretical curves are obtained by using an extended projectile (so the impact parameter $B$ refers to the center of this projectile, as in Fig. 12), together with the following values for the free parameters: $B=7$ fm, $R^{2}=2\,\,\text{GeV}^{-2}$, $Q_{0s}^{2}=3$ GeV2, $m=0.1$ GeV, and two values for $\langle z\rangle$: 0.6 and respectively 0.7. Notice that this value $Q_{0s}^{2}=3$ GeV2 is in the ballpark of the theoretical expectations for the proton saturation momentum at the LHC energies. It is quite remarkable that, in spite of the many simplifications and the crude assumptions, our mechanism appears to be able to produce the correct $p_{T}$-dependence of $v_{2}$ and also its correct size, with reasonable values for the free parameters.

V Conclusions and perspectives

In this paper, we have analyzed a less explored, albeit not totally new, mechanism for generating azimuthal asymmetries in particle production in ‘dilute-dense’ collisions: the dependence of the cross-section for single-inclusive particle production upon the azimuthal orientation of the momentum of the produced particle w.r.t. its impact parameter. As compared to previous related studies in the literature Kopeliovich et al. (2008a, b, c); Levin and Rezaeian (2011); Zhou (2016); Hagiwara et al. (2017), we have considered a different model for the gluon distribution in the dense target — an extension of the McLerran-Venugopalan (MV) model — which combines the proper pQCD tail at high transverse momenta with a Gaussian profile (inspired by fits to the HERA data) for the distribution of the color charges in the impact-parameter space and a ‘gluon mass’ which mimics confinement. A realistic, or at least physically motivated, treatment of the non-perturbative aspects related to the transverse inhomogeneity and to confinement is indeed essential, since, as demonstrated by our analysis, these aspects do actually control the azimuthal asymmetries produced by this mechanism.

Within this set-up, we have given a systematic, semi-analytic, study of the angular dependence of the dipole scattering amplitude and we have used the result to compute the elliptic flow coefficient $v_{2}$, separately for proton-proton and proton-nucleus collisions (the main difference being the inclusion of lumpiness effects in the case of a nuclear target). We thus found that, as a function of the impact parameter, $v_{2}$ is rather strongly peaked near the edge of the target, where the transverse inhomogeneity is more pronounced. Hence the present mechanism will significantly contribute to azimuthal asymmetries only via peripheral collisions. Furthermore, as a function of the transverse momentum of the produced particle, $v_{2}$ shows a maximum at some intermediate, semi-hard, value which is proportional to the target saturation momentum at $b=0$. This maximum broadens and lowers with increasing the saturation momentum. The overall $p_{T}$-dependence is quite similar to that observed in the phenomenology of p+Pb collisions at the LHC. In fact, a semi-quantitative agreement with the data can be obtained with reasonable choices for the model parameters, but this agreement should be taken with a grain of salt, because our model is quite crude and additional approximations are performed when comparing with the data.

To convincingly demonstrate this mechanism, further studies are necessary. First, one would like to understand its evolution with increasing energy. To that aim, one should solve the Balitsky-Kovchegov (BK) equation Balitsky (1996); Kovchegov (1999) with initial conditions at low energy provided by our present model. In the respective solutions, one must keep trace of the impact-parameter dependence (including the angular dependence) of the dipole amplitude and one must enforce confinement within the soft gluon emission kernel, preferably by using the same infrared regulator (‘gluon mass’) $m$ as in the initial condition. Similar solutions have been already considered in Motyka and Stasto (2009); Berger and Stasto (2011, 2013).

Furthermore, since the orientation of the reaction plane cannot be experimentally measured, one must compute the imprint of the azimuthal asymmetries generated by the present mechanism on multi-particle correlations. For instance, one can study the simultaneous production of 2 quarks in $pA$ collisions. The calculation can be simplified by assuming that the two quarks (originally collinear with the projectile) scatter independently of each other and by taking the multicolor limit $N_{c}\to\infty$. (We recall that the multi-particle correlations generated by the present mechanism survive in the large-$N_{c}$ limit.) As anticipated in the Introduction, we expect the multi-particle correlations to be important only for sufficiently peripheral collisions and to lead to flow — a collective motion of particles which are produced independently from each other, but which are all correlated with the reaction plane defined by their average impact parameter. Note that even for independent partons in the projectile at large $N_{c}$, the double-patron-scattering (DPS) Hanbury Brown and Twiss (HBT) correlations may lead to anisotropy Kovner and Rezaeian (2017). This effect is not accounted for in the present mechanism.

Given the prominence of the peripheral physics for the physical problem at hand, we furthermore expect that the fluctuations in the shape of the colliding hadrons should play an important role: they should amplify the inhomogeneity in impact parameter space and thus enhance the azimuthal asymmetries. The importance of such fluctuations is supported by a recent analysis of the HERA data for incoherent exclusive diffractive vector meson production in deep inelastic scattering M ntysaari and Schenke (2016). It would be therefore interesting to redo our present analysis of the angular dependence of the dipole scattering for the case where the (projectile and/or target) proton has strong shape fluctuations, say as described by the models used in M ntysaari and Schenke (2016) and which are constrained by the HERA data.

Acknowledgements.

We would like to thank Y. Hagiwara, Y. Hatta, B.-W. Xiao, and F. Yuan for sharing with us a preliminary version of their recent work Hagiwara et al. (2017) and for related discussions. We are grateful to T. Lappi, J.-Y. Ollitrault and L. McLerran for comments on the manuscript. This research was supported in part by ECOS-Conicyt C14E01. The work of E.I. was supported in part by the Agence Nationale de la Recherche project ANR-16-CE31-0019-01. The work of A.R. was supported in part by Fondecyt grant 1150135, Anillo ACT1406 and Conicyt PIA/Basal FB0821.

Appendix A Computing the dipole amplitude in the two-gluon exchange approximation

In this Appendix, we present details on the calculation leading to Eqs. (III.3) and (33) for the dipole scattering amplitude in the two-gluon exchange approximation. The starting point is Eq. (29) for $N_{2g}({\bm{b}},{\bm{r}})$, which, we recall, has been obtained from the general MV-model formula (20) by expanding to second order in ${\bm{k}}\cdot{\bm{r}}$.

Specifically, let $\beta$ and $\theta_{k}$ denote the angles between ${\bm{\Delta}}$ and ${\bm{k}}$ and, respectively, between ${\bm{b}}$ and ${\bm{k}}$. Then we can rewrite Eq. (29) as

[TABLE]

We first perform the integral over the angle $\theta_{k}$ using the identity in Eq. (5), to obtain

[TABLE]

The ensuing integrals over $\beta$ can be now performed by using the following identities

[TABLE]

with $A=k^{2}+\Delta^{2}/4+m^{2}$ and $B=k\,\Delta$. We thus confirm the general structure in Eq. (31) and at the same time obtain more explicit expressions for the two functions $\mathcal{N}_{0}(b,r)$ and $\mathcal{N}_{\theta}(b,r)$, as double radial integrals to be numerically computed. Specifically, we find

[TABLE]

together with Eq. (33) for $\mathcal{N}_{\theta}(b,r)$. It is easy to check that the double integral yielding $\mathcal{N}_{\theta}(b,r)$ is well defined as is stands, so its numerical evaluation poses no special problem. On the other hand, the integral over $k$ giving $\mathcal{N}_{0}(b,r)$ has a logarithmic ultraviolet divergence and needs to be cut off at $k\sim 1/r$, as already explained. The precise implementation of this UV cutoff introduces some ambiguity, that we shall fix by replacing Eq. (61) with Eq. (III.3), as discussed in the main text.

As also explained in Sect. III.2, the limit $m\to 0$ of $N_{2g}$ is well defined, because of the explicit transverse inhomogeneity in the target. We have already studied this limit for the piece $\mathcal{N}_{\theta}(b,r)$ which controls the angular dependence, cf. Eq. (39). In what follows we shall perform the corresponding study for the other piece $\mathcal{N}_{0}$. To that aim, it is more convenient to use Eq. (61), in which we set $m=0$ and simply cut the integral over $k$ in the ultraviolet at $k\sim 1/r$. Using Eq. (37), we find

[TABLE]

This is strictly true for $1/r\gg\Delta$, meaning $r\ll R$ (recall that the integral over $\Delta$ is restricted to $\Delta\lesssim 1/R$ by the transverse inhomogeneity in the target). At this point, it is convenient to separate the logarithm as $\ln({1}/{r^{2}\Delta^{2}})=\ln(R^{2}/r^{2})+\ln(1/R^{2}\Delta^{2})$, where the first piece in the r.h.s. is the ‘large logarithm’ (since typically $r\ll R$), whereas the argument of the second logarithm is of order one. One then easily finds

[TABLE]

where in the second line we kept only the power-law tail at large $b$, as developed by the integral over $\Delta$ from the first line. (This power law tail can be easily obtained by letting $\tilde{\mu}(\Delta)\to\tilde{\mu}(0)$ inside the integrand, as appropriate for large $b\gg R$.) Notice that for very large impact parameters $b\gg R$, the first term proportional to the local saturation momentum $Q_{s}^{2}(b)$ is exponentially suppressed and the scattering amplitude (which is small anyway) is controlled by the power-law tail. However, this is not the situation that we have considered throughout this work; indeed, we have looked at impact parameters $b\lesssim R$, where the scattering amplitude is dominated by the first term $\propto Q_{s}^{2}(b)$.

Appendix B More details on the case of a lumpy nuclear target

In this Appendix, we shall present more details on the derivation of Eqs. (54) and (55) for the dipole scattering off a lumpy nucleus in the two-gluon exchange approximation. Namely, we would like to compute the integral over ${\bm{b}}$ in Eq. (50) with the nuclear thickness function given by the second-order expansion in Eq. (51). To that aim, it is useful to introduce the following Fourier transform

[TABLE]

To the order of interest, $N_{2g}({\bm{b}},{\bm{r}})$ is given by Eq. (29), which immediately implies

[TABLE]

The interesting integrals over ${\bm{b}}$ can be related to the behavior of $\tilde{N}_{2g}({\bm{\Delta}},{\bm{r}})$ near $\Delta=0$. We first have

[TABLE]

We have also used here $k^{l}k^{m}\to(k^{2}/2)\delta^{lm}$ together with the renormalization prescription in Eq. (34) and $\tilde{\mu}(0)=4\pi R^{2}\mu_{0}$. Not surprisingly, the integral over ${\bm{b}}$ has generated a result proportional to the proton area $\sim R^{2}$.

For the terms quadratic in $b$, one similarly finds

[TABLE]

where in going from the second to the third line we have used the fact that, under the integral over ${\bm{k}}$, one can replace

[TABLE]

By combining the above equations (66) and (67) with the expansion in Eq. (51), one immediately finds the results exhibited in Eqs. (54) and (55).