## Abstract

We study the effects on multigluon production at mid-rapidity in the Color Glass Condensate of the non-eikonal corrections that stem from relaxing the shockwave approximation and giving the target a finite size. We extend previous works performed in the dilute-dilute approximation suitable for proton–proton collisions, to the dilute-dense one applicable to proton–nucleus. We employ the McLerran–Venugopalan model for the projectile averages. For the target averages, we use the Golec-Biernat–Wüsthoff model and restrict to the leading contributions in overlap area that allow a factorization of ensembles of Wilson lines into products of dipoles. We make the connection with the jet quenching formalism and compare with previous results in the literature, providing a parametrization of the so-called decorated dipoles. We show that the non-eikonal effects on single inclusive particle production, contrary to what happens in jet quenching, are only sizable, of a few percent, for modest energies \(\sqrt{s_{NN}}\le 100\) GeV and central rapidities. On the other hand, we find that the effects on double inclusive gluon production are larger for the same kinematics. We show that, as found previously in the dilute-dilute situation, non-eikonal corrections break the accidental symmetry in the CGC, allowing for the existence of non-vanishing odd azimuthal harmonics.

### Similar content being viewed by others

Avoid common mistakes on your manuscript.

## 1 Introduction

The understanding of multiparticle production is a central issue in high energy physics. Multiparticle production is a direct consequence of the interaction and also acts as background for new phenomena that may eventually appear. It is of particular interest in collisions at high energies or involving extended objects like nuclei. There, the number of produced particles is so large that a microscopic treatment is defying and investigating the possibilities of simplified treatments becomes compulsory.

In heavy ion collisions evidence for the creation of a new state of matter, the Quark Gluon Plasma (QGP), is found. After an initial short time, smaller than 1 fm/c, such state is amenable to a macroscopic description in terms of relativistic hydrodynamics (see e.g. [1, 2]). One of the key findings at the Large Hadron Collider (LHC) at CERN is the observation that many of the behaviors of observables that in heavy ion collisions are considered as evidence of QGP formation and as tools to analyse its properties, have also been observed in collisions of smaller systems, proton–proton (*pp*) and proton–nucleus (*pA*) [3,4,5,6,7]. In these small collision systems the success of the hydrodynamic explanation based on final state interactions seems difficult to justify and alternatives based on the dynamics in the initial state have also been essayed.

In *pA* collisions at high energies, particle production is commonly studied under the approximations that consider the projectile as a dilute probe which scatters off a dense target described by an intense color field. This process can be computed within the Color Glass Condensate (CGC) effective theory [8, 9] where the dilute projectile probes the small-*x* glue inside the target. This framework relies on the eikonal approximation where the collision is described by infinitely boosted partons in the projectile that traverse a medium in the target with infinitesimal width – a situation referred to as the shock-wave approximation. In this approximation, fluctuations in the projectile wave function are long-lived and can be considered as frozen throughout the interaction time with the target. Indeed, the eikonal approximation amounts to considering the target field of a left-moving nucleus \(A^\mu (x^+,x^-,{\textbf{x}})\propto \delta ^{\mu -} \delta (x^+) A^-({\textbf{x}})\), where light-cone coordinates \(x^\pm =(x^0\pm x^3)/\sqrt{2}\) are used and transverse coordinates are denoted by boldface letters \({\textbf{x}}=(x^1,x^2)\). The infinite boost is then present in three different ways in the expression for the target field: the \(\delta \)-function in \(x^+\) which sets the target width to zero, the independence of \(x^-\) which is equivalent to having recoilless interactions, and keeping only the minus component which is enhanced by the boost. Each of these features can be relaxed separately, thus leading to sub-eikonal corrections.

The CGC offers a framework where particle production at the very early phases of the collision [8,9,10] and initial state correlations (see the review [11]) can be analysed. The origin of such initial state correlations lies in the quantum effects (Bose enhancement in the case of gluons) in the wave function of the colliding hadrons [12,13,14,15,16,17,18], but problems still exist to include subleading density effects (see recent efforts in [19,20,21]), and to understand the behaviour of correlations with centrality or multiplicity in the event [22, 23]. Furthermore, odd coefficients of the azimuthal Fourier decomposition of the distribution of produced particles are absent in standard CGC calculations, and can only be generated by including density corrections [24, 25], or a modification of the usual isotropic field ensembles in the target [26, 27], or non-eikonal corrections.

Sub-eikonal effects scale with the inverse of the beam energy. Therefore, the eikonal approximation is well justified at the kinematics of the Large Hadron Collider (LHC) at CERN where *pA* collisions are performed at centre-of-mass energies per nucleon, \(\sqrt{s_{\textrm{NN}}}\), of order \(10^3\) GeV. It is only approximate, with the corrections being around a few percent [28, 29], at the top energies, \(\sqrt{s_{\textrm{NN}}}=200\) GeV, of the Relativistic Heavy Ion Collider (RHIC). Therefore, with the upcoming Electron-Ion Collider (EIC) [30, 31] where electron-nucleus collisions will be performed at \(\sqrt{s_\textrm{NN}}\sim 20\div 100\) GeV, these sub-eikonal contributions will become relevant, as they already are for RHIC energies. For this reason studies aimed at including non-eikonal corrections, such as finite width effects, have experienced an intense growth in recent years.

Further, spin physics requires going beyond the eikonal approximation where spin-flip terms and sub-leading exchanges in the *t*-channel are absent. Sub-eikonal corrections are introduced in the quark and gluon propagators by first assuming a target with finite width and then relaxing the other mentioned approximations [32,33,34,35,36,37,38]. In order to compute quark and anti-quark helicity transverse momentum distributions and parton distribution functions, modifications to the JIMWLK evolution equations to include helicity-dependent effects are introduced through polarized Wilson lines [39,40,41,42,43,44,45], with several numerical analyses of the impact of small-*x* evolution on the comparison with experimental data done recently [46, 47]. A different method in order to study the effects of the non-eikonal corrections was pursued in [48,49,50] by including the longitudinal momentum exchange between the projectile and the target during the interaction. Moreover, studies including the effects of an \(x^-\)-dependent, i.e., dynamical, target field in the quark propagator are being carried out [38, 51, 52].

As mentioned before, it has been shown by the inclusion of finite width effects in double gluon production that the so-called accidental symmetry, i.e., an artificial azimuthal symmetry that appears in the leading order multigluon spectrum within the CGC and results in vanishing odd azimuthal harmonics [11], is broken [28, 29]. Thus, finite width effects open a new window for explaining, within the initial state framework provided by the CGC, the away- and near-side asymmetry seen in small collision systems at RHIC where non-eikonal corrections are sizeable.

Finite width effects appear naturally in jet quenching where, in order to study the modification of jets traversing the QGP, one assumes that a high-energy parton traverses a finite colored medium and loses its energy through the induced emission of gluons [53,54,55]. One is usually interested in hard partons created in the initial collision which have not had time to radiate before interacting with the medium and therefore would not have any modification if one keeps the strict shock-wave approximation. In order to have non-trivial contributions from medium effects it is necessary to consider a finite longitudinal extent of the target, thus including sub-eikonal corrections which are resummed through in-medium propagators of gluons and quarks in a background field. Thus, inspired by jet quenching calculations, the approach of this manuscript will be to include the non-eikonal effects to the multigluon spectrum coming from considering a target with finite width through a modification of the eikonal gluon propagator. It is worth noting that, even though the techniques are the same and the calculations resemble each other to some degree, we are indeed calculating very different things in the jet quenching case from particle production in *pA*. What appears as the main contribution for in-medium emission (transverse motion along a finite longitudinal extent) is only a subleading correction when the incoming particles are allowed to develop long-lived fluctuations before scattering with a background medium.

The main goal of this work is to generalize the results of [28, 29], where some of us have relaxed the shock-wave approximation and analysed the effects of including such non-eikonal corrections in multigluon production in *pp* collisions, to the case in which one of the participants is dense. We study the *finite width effects* in multigluon production in the dilute-dense limit of the Color Glass Condensate framework which is suitable for *pA* collisions, while still assuming that all momentum exchanges are only transverse and the only non-zero component of the background field is the minus component. We note that we work at a semiclassical level, and no rapidity evolution is considered. This manuscript is organised as follows. In Sect. 2 we present a short introduction of the framework, similar to that employed for jet quenching, which will be used along this manuscript. In Sect. 3 we generalize the results of [56], where we have computed the target averages of Wilson lines within a model [57,58,59] – named the area enhancement model – that keeps those contributions that are leading in the product overlap area of the collision times the squared saturation momentum of the target, to the non-eikonal case in which the Wilson line has to be substituted by the scalar gluon propagator. In Sect. 4 we compute the single gluon spectrum beyond the eikonal accuracy – in the sense of finite width effects. We analyze the effects of the non-eikonal corrections in Sect. 4.2. In Sect. 4.3 we compare our results with the next-to-next-to-eikonal expansion performed in [32, 33] obtaining a parameterisation, within the Gaussian approximation for target averages, for the decorated dipole functions defined in those works. In Sect. 5 we introduce a general framework for computing the multigluon spectrum for any number of gluons. In Sect. 5.1 we apply this framework to double inclusive gluon production and we analyse the dependence of the double gluon spectrum on the non-eikonal effects. Finally, in Sect. 6 we conclude with a summary and outlook. Technical details are provided in Appendices A and B.

## 2 Brief review of the theoretical framework

In this section, we briefly summarize the theoretical framework that we adopt to study multigluon production in *pA* collisions beyond eikonal accuracy. In this framework, the dilute projectile, which is a highly boosted right-moving proton, is composed of large-*x* partons that act as a color source with color charge \(\rho _p^a(\textbf{x})\), with index *a* denoting color and \(\textbf{x}\) being transverse position. The left-moving dense target is defined through the target field \(A^-(x^+,\textbf{x})\) which has a finite support from 0 to \(L^+\) along the \(+\) direction, and thus goes beyond the eikonal description of the target where the target fields are localized around \(x^+=0\) – the shockwave limit. In this setup, by restricting ourselves to the light-cone gauge \({A^+=0}\), the production amplitude of a gluon with transverse momenta \(\textbf{k}\), longitudinal (plus) momentum \(k^+\), color *a* and polarization \(\lambda \) can be obtained by using the LSZ reduction formula and reads^{Footnote 1}

where \(\rho ^b_p(\textbf{q})\) is the Fourier transform of the projectile color charge density and the reduced matrix amplitude is given by^{Footnote 2}

The three terms in the reduced matrix amplitude, Eq. (2.2), describe emission of the gluon from the projectile color charges after, before or inside the medium respectively and are illustrated in Fig. 1.

The propagation of the gluon inside the medium is described by the scalar retarded propagator

This path integral, from transverse position \({\textbf{y}}\) at \(y^+\) to \({\textbf{x}}\) at \(x^+\) along trajectory \({\textbf{z}}(z^+)\), encodes both the color rotation and the motion of the gluon in the transverse plane while traversing the medium. The propagation of the gluon outside the medium is given by the free scalar propagator

(which can be obtained by setting \(U^{ab}\) to \(\delta ^{ab}\) and computing the remaining integral). Finally, \(U^{ab}(x^+,y^+; \textbf{x})\) is the Wilson line which accounts for multiple gluon exchanges between the projectile and the target and it is given by the path-ordered exponential

In this approach the single inclusive gluon spectrum is given by

where \(\langle \cdots \rangle _{p,T}\) accounts for the average over the color configurations of the projectile and the target respectively. The explicit expression of this function is derived later in this manuscript and can be found in (4.1) and (4.3).

Within the same setup, single inclusive gluon production was studied in Refs. [32, 33] by developing a systematic method to compute higher order corrections to the eikonal approximation (or equivalently to the shockwave approximation) that are due to the non-zero longitudinal width of the target. Specifically, the method was developed to expand the scalar retarded propagator given in Eq. (2.3). The expansion was performed at next-to-next-to-eikonal accuracy and the result was applied to the single inclusive gluon production process in *pA* collisions. In the present work, we follow a different approach to account for the non-eikonal effects that stem from the finite longitudinal thickness of the target. Contrary to the studies performed in Refs. [32, 33], in the present work we discuss how to perform target averaging of non-eikonal objects that appear at the level of the squared amplitude without adopting the aforementioned eikonal expansion.

As mentioned briefly in Sect. 1, inclusive multigluon spectra is an essential observable for the study of particle correlations from the initial state point of view for small size systems such as *pp* or *pA* collisions. Within the CGC framework, these studies are usually performed in the “glasma graph approximation” [12, 13, 61,62,63] which amounts to allowing the emission of a gluon from a single color source in the projectile and neglecting the contributions where more than one gluon is emitted from the same color source in the projectile. One important aspect of the glasma graph calculations is that it is only valid for dilute-dilute scattering processes, i.e., *pp* collisions, because it only considers two gluon exchanges with the target. However, its extension that accounts for the multiple interactions with the target which corresponds to dilute-dense scatterings, or equivalently *pA* collisions, was performed in [59, 64]. In a recent study [56], the extension of the glasma graph calculations to compute *n*-gluon production in *pA* collisions was developed and the four-gluon production spectrum was explicitly computed in the strict eikonal limit.

Here, we extend the eikonal calculation of the *n*-gluon spectrum to the non-eikonal case where the effects of the finite longitudinal width of the target are accounted for. The non-eikonal generalization of the eikonal *n*-gluon spectrum computed in [56] can be written as

where a shorthand notation \({\underline{k}} \equiv (k^+,\textbf{k})\) is introduced. Due to the finite longitudinal extend of the target, the reduced matrix amplitude is a more complicated object compared to its eikonal limit. It includes not only the standard Wilson lines as in the case of strict eikonal limit but also the retarded scalar propagator \({\mathcal {G}}_{k^+}(x^+,\textbf{x} ; y^+,\textbf{y}) \) (see Eq. (2.2)). Therefore, at the squared amplitude level one gets not only eikonal multipoles (dipoles, quadrupoles, etc.) but their non-eikonal generalizations which include scalar propagators. In the next section, we discuss how to evaluate the target averaging of such new objects that appear beyond eikonal accuracy.

## 3 Non-eikonal target averaging

In order to obtain an analytical result for the inclusive *n*-gluon spectra given in Eq. (2.7), the target average of the 2*n* reduced matrix amplitudes has to be computed. It is well known that in the strict eikonal limit of the single inclusive gluon spectrum, one ends up with the eikonal dipole function which is usually evaluated by using some model like McLerran–Venugopalan (MV) [65, 66] or Golec-Biernat–Wüsthoff (GBW) [67, 68]. However, for two or more gluons multipole functions such as quadrupoles, sextuples, etc., appear. In order to easily compute the target averaging of these multipole functions, the Area Enhancement (AE) model^{Footnote 3} is introduced in Refs. [57,58,59]. This model is based on the following ansatz: Any multipole can be written in terms of dipole functions through a Wick expansion for those configurations of multipole functions that maximize the phase space integration, up to terms that are suppressed by the collision area.^{Footnote 4} In other words, in the AE model, after the integration over the phase space, the Wilson lines can approximately be described by a Gaussian distribution up to the corrections that are of the order of the inverse of the phase space area. In this section, we will generalize the target averaging of eikonal two point functions (eikonal dipole functions) to compute the non-eikonal two point functions (i.e., non-eikonal dipole functions) that appear in the single inclusive gluon production when one includes the target finite longitudinal width effects. Then, these results will be used together with the AE argument to study the non-eikonal multigluon spectra.

Non-eikonal single inclusive gluon spectrum is given by setting \(n=1\) in Eq. (2.7), which requires the evaluation of the following three objects:

where we have written \(U_\textbf{x}(x^+,y^+) \equiv U(x^+,y^+;\textbf{x})\) for simplicity. It is straightforward to note that \(d^{(0)}(x^+,y^+|\textbf{y},\bar{\textbf{y}})\) given in Eq. (3.1) can be identified with the eikonal dipole function evaluated over a finite longitudinal extent \(z^+\in [y^+,x^+]\) of the target. The functions defined in Eqs. (3.2) and (3.3) are non-eikonal generalizations of the dipole function.

For simplicity, we assume that the eikonal dipole function is described by the GBW model which is valid as long as the dipole size is much smaller than \(\Lambda _{\textrm{QCD}}^{-1}\). In this case, the eikonal dipole function reads

where \(q_s(x^+,y^+)\) can be identified as the effective saturation momentum in a longitudinal slice \([y^+,x^+]\) of the target. Within the MV model, saturation momentum can be defined as^{Footnote 5}

with \(\mu ^2\) the color density of the target and \(C_R\) the Casimir of the projectile parton interacting with the target. With this definition of the saturation momentum, the effective saturation momentum \(q_s(x^+,y^+)\) can be defined as

Thus, assuming that the color density \({\tilde{\mu }}^2(z^+)\) is constant along the target longitudinal extent, i.e., \({\tilde{\mu }}^2(z^+) = \mu ^2/L^+\), the effective saturation momentum takes the following form:

The object defined in Eq. (3.2) is one of the generalizations of the eikonal dipole function which stems from the finite longitudinal extent of the target and we refer to it as the 1st order NE dipole function. It requires averaging over the scalar propagator \({\mathcal {G}}_{k^+}(x^+,\textbf{x} ; y^+,\textbf{y})\) defined in Eq. (2.3) and therefore implies solving a path integral. In order to do so, we discretize the longitudinal space into *N* slices where each discretized point is labeled as \(z_i^+\). In this case, the scalar propagator can be written as [32, 33, 60]

where we have defined \({{\underline{x}}} \equiv (x^+,\textbf{x})\). This equation describes the discrete random walk of the emitted gluon through the transverse points \(\textbf{z}_n(z_n^+)\) inside the target. The emitted gluon is propagated from \(\textbf{z}_0(z_0^+)=\textbf{y}\) with \(z_0^+=y^+\) to \(\textbf{z}_N(z_N^+)=\textbf{x}\) with \(z_N^+=x^+\). The only dependence on the target configuration in Eq. (3.8) appears through the discretized Wilson lines, \(U^{ab}_{\textbf{z}_n}(z_{n-1}^+,z_n^+)\), that account for the eikonal propagation of the gluon in each longitudinal slice. Thus, the 1st order NE dipole function defined in Eq. (3.2) can be written as

This expression can be further simplified by realizing that, due to the properties of path-ordered exponentials, the Wilson line evaluated over some longitudinal extent \([y^+,x^+]\) can be factorized into a product of independent contributions at each slice of the discretized axis:

where \(x_0^+=y^+\) and \(x_n^+=x^+\). Moreover, noting that the MV model is local in the longitudinal direction which implies that the average of Wilson lines evaluated at different points on the longitudinal axis factorizes into independent averages, we simplify the target average in Eq. (3.9) to

Finally, by using the local GBW model given in Eq. (3.4) and noting that \(z_i^{+}-z_{i-1}^{+}=L^+/N\), the 1st order NE dipole function takes the following form:

In the continuum limit, by defining \(\textbf{r}=\textbf{z}-\bar{\varvec{y}}\), Eq. (3.12) can be written as

Two comments are in order here. First, the result given in Eq. (3.13) can be obtained in a straightforward manner by using the expression for the scalar propagator \({\mathcal {G}}_{k^+}(x^+,\textbf{x} ; y^+,\textbf{y})\) given in Eq. (2.3) in the definition of the \(1^\textrm{st}\) order NE dipole function and adopting the GBW model directly without introducing the discretization along the longitudinal axis. However, the way that Eq. (3.13) is obtained above justifies the use of the locality argument which plays a crucial role in our discussion. Second, Eq. (3.13) is the well known path integral of the harmonic oscillator with “mass” \(k^+/2\) and imaginary “frequency” \(\sqrt{-iQ_s^2/2L^+k^+}\) that is frequently used in jet quenching calculations (see for example [54]). For completeness, we compute this integral in Appendix A with the result given in (A.22). Thus, after performing the path integral, the 1st order NE dipole function reads

where we have defined \(\Delta ^+=x^+-y^+\), \(\textbf{a}_0=\textbf{y}-\bar{\varvec{y}}\), \(\textbf{a}_N=\textbf{x}-{\bar{\varvec{y}}}\) and

Note that \(\epsilon _i^2\) is a dimensionless parameter that vanishes in the eikonal limit (\(k_i^+ \rightarrow \infty \) and \(L^+ \rightarrow 0\)).^{Footnote 6} Its role is to weight the strength of the finite longitudinal width effects and we will refer to it as the *non-eikonal parameter*.

The next two-point function that we consider is defined in Eq. (3.3). It is referred to as the \(2^\textrm{nd}\) order NE dipole function which can be computed in the same manner as the 1st order one. By using the GBW model and the locality of the MV model we can write it as

where \(\textbf{r}(y^+)=\textbf{y}\), \(\textbf{r}(x^+)=\textbf{x}\), \({\bar{\textbf{r}}}(y^+)={\bar{\varvec{y}}}\) and \({\bar{\varvec{r}}}(x^+)={\bar{\varvec{x}}}\). This equation can be identified with the path integral of two coupled harmonic oscillators which can be decoupled. The solution is also derived in Appendix A and the final result is given in Eq. (A.44). Thus, the 2nd order NE dipole function reads

where we have defined for simplicity \(\epsilon _{\pm }^2=\epsilon _1^2\pm \epsilon _2^2\), \(\textbf{r}_0=\textbf{y}-\bar{\varvec{y}}\), \(\textbf{r}_N=\textbf{x}-\bar{\varvec{x}}\), \(\textbf{b}_0=(k_1^+ \textbf{y}+k_2^+ \bar{\varvec{y}})/(k_1^++k_2^+)\) and \(\textbf{b}_N=(k_1^+ \textbf{x}+k_2^+ \bar{\varvec{x}})/(k_1^++k_2^+)\).

We would like to emphasize that this solution is novel and only in the limit \(k_2^+=k_1^+\) it simplifies and leads to the known results [53, 69,70,71].

which is required for single inclusive gluon production or medium-induced gluon radiation. But the more general result for \(k_2^+\ne k_1^+\) is required for double inclusive gluon production, that we analyse below (see also Appendix B).

With Eqs. (3.4), (3.14) and (3.17) we have determined the non-eikonal dipole functions evaluated within a longitudinal extent \(\Delta ^+\) of the target medium. We note that these functions are extensively used by the jet quenching community where the effects of the parton propagation within a dense medium is computed. In that case, the GBW model is usually referred to as the *harmonic oscillator approximation* and the effective saturation scale defined in Eq. (3.4) is written

where \({\hat{q}}(z^+)\) is the medium transport coefficient. For a static homogeneous medium in which \({\hat{q}}(z^+)\) is constant, we get \({\hat{q}}=Q_s^2/L^+\).

Before we conclude this section, we would like to comment on the computation of the higher order functions, i.e., multipoles, where one has to evaluate the non-eikonal target average of multiple Wilson lines and scalar propagators. As stated shortly at the beginning of this section, in order to compute the NE multipoles we use the AE model. This model is based on the chromo-electric domain structure in the target transverse plane and leads to a similar result when compared with the MV model. The advantage of this approach is that, upon integration over the transverse phase space, the Wilson lines follow approximately a Gaussian distribution – with the non-Gaussian corrections suppressed by the collision area – and therefore allows one to adopt Wick’s theorem. However, in the case of a target with a finite longitudinal extent, one should be more careful since the chromo-electric domains may decohere in the longitudinal direction.

To justify the application of the AE model in target ensembles with non-zero width, \(L^+\), we examine the space-time kinematics of the interaction in the center-of-mass (CoM) frame. In this frame the Lorentz contracted target width is \(L^+ \sim A^{1/3}/\sqrt{s_\textrm{NN}}\). On the other hand, the chromo-electric fields are defined by the low-*x* gluons which have the following spread in the longitudinal direction: \(\Delta x^+ \sim 1/q^- = A/(xQ_T^-)\), where \(Q_T^- = A \sqrt{s_{\textrm{NN}}}/\sqrt{2}\) and \(q^-\) are the longitudinal momentum (in the CoM frame) of the target and the gluon respectively. Thus, as long as \(\Delta x^+ \gg L^+\) the domain structure of the target will not decohere within its extent and the AE model is justified. Therefore, for \(L^+/\Delta x^+ \sim x A^{1/3} \ll 1\) the Area Enhancement model should be a good approximation in the non-eikonal case.^{Footnote 7}

## 4 Non-eikonal single inclusive gluon production

Before computing the non-eikonal multi gluon production, we first consider the single inclusive case as a warm up. In this section, we first derive the analytical solutions of the non-eikonal single inclusive spectrum, then perform the numerical analysis of these results and finally compare our results with the ones that exist in the literature (specifically with the results of Ref. [32, 33]).

### 4.1 Analytical results

The non-eikonal single inclusive gluon spectrum is given by Eq. (2.7) for \(n=1\), which reads

We employ the MV model to compute the averaging over the projectile color sources:

where the factor \(1/(N_c^2-1)\) is introduced by convenience but will not be relevant in the analysis performed in this manuscript.^{Footnote 8} Since this contribution is proportional to \(\delta ^{b_1 b_2}\), it leads to a color trace of the reduced amplitudes in Eq. (4.1). By using the expression for the reduced amplitude given in Eq. (2.2), the target average of the reduced amplitudes for single inclusive gluon spectrum can be written as^{Footnote 9}

where we have used the convolution property of the scalar propagator

in order to break the target average of the scalar propagator into its local pieces where the number of objects is fixed.

It is worth mentioning that since we are using the GBW model (i.e., a Gaussian form) for the non-eikonal dipole functions, all the integrals over the transverse coordinates or transverse momenta appearing in Eq. (4.3), or in general in this manuscript, are of the form

and, therefore, although very tedious in some cases, it is straightforward to compute them. Moreover, the argument of the integrals over \(\textbf{y}\) and \(\bar{\varvec{y}}\) is translational invariant under these coordinates and thus the integral is proportional to \(\delta ^{(2)}(\textbf{q}_1-\textbf{q}_2)\). Then, the single inclusive gluon spectrum given in Eq. (4.1) can explicitly be written as

where we have used the fact that \(\mu ^2(\textbf{q},-\textbf{q})=\pi B_p\), with \(B_p\) the gluonic transverse area of the projectile. We have also changed our notation from the longitudinal momentum to rapidity: \(dk^+/k^+=d\eta \). The terms that appear in the argument of the integral in Eq. (4.6) result from squaring the contributions in Eq. (2.2) where the gluon is emitted before, after or inside the target medium (see Fig. 2).

By using the expressions for the eikonal dipole function and the two non-eikonal dipole functions, explicit expressions of each contribution in Eq. (4.6) can be computed. Let us start with the \({\mathcal {I}}_{\text {aft-aft}}\) term. This term is the contribution where the gluon is emitted after the interaction of the source with the target on both sides of the cut. It reads

Note that we have dropped the dependence on the impact parameter \(\textbf{b}\) because the eikonal dipole function is translational invariant.

The next contribution in Eq. (4.6) is the \({\mathcal {I}}_{\text {bef-bef}}\) term which corresponds to the emission of the gluon before the interaction of the source with the target on both sides of the cut. It reads

We would like to mention that we keep the impact parameter \(\textbf{b}\) dependence in Eq. (4.8) for consistency of the notation but, after performing the Fourier transforms, the final result is translationally invariant, i.e., independent of \(\textbf{b}\) (see also the comment before Eq. (4.6)). Moreover, the Fourier transform of the 2nd order NE dipole when the two longitudinal momenta are equal (as in the case of Eq. (3.18)) reads

which shows that in this case the 2nd order NE dipole function does not provide any non-eikonal correction. As we shall see later, this is not the case when the two longitudinal momenta are different. Thus, the result for the \({\mathcal {I}}_{\text {bef-bef}}\) contribution reads

The \({\mathcal {I}}_{\text {aft-bef}}\) term corresponds to the case where the gluon is emitted before the interaction of the source with target on one side of the cut and after on the other side. This term can be written as

Noting that the Fourier transform of the 1st order NE dipole given in Eq. (3.14) reads

\({\mathcal {I}}_{\text {aft-bef}}\) takes the following form:

It is straightforward to realize that from this equation one can recover the results coming from the GBW model by considering the eikonal limit \(|\epsilon ^2| \rightarrow 0\) (see the dependence on \({\textbf {k}}\) and **q** in Eq. (4.7)).

\({\mathcal {I}}_{\text {aft-in}}\) corresponds to the case where the gluon is emitted after the interaction of the source with the target on one side of the cut and inside the target on the other side. This contribution can be written

where \(\textbf{y}=\textbf{b}+\frac{\textbf{r}}{2}\), \(\bar{\varvec{y}}=\textbf{b}-\frac{\textbf{r}}{2}\) and \(y^+\) is the longitudinal coordinate where the gluon is emitted inside the target. Using Eq. (4.12) and defining \(\alpha =\frac{L^+-y^+}{L^+}\) as the fraction of the target that the gluon traverses, we can write this contribution as

Here, we also introduced the notation \({\tilde{\alpha }}=1-\alpha \) that corresponds to the longitudinal extent of the target traversed by the source before emitting the gluon. It is worth mentioning that when the trigonometric functions are expanded to the leading order in \(\epsilon \), \({\mathcal {I}}_{\text {aft-in}}\) is \({\mathcal {O}}(\epsilon ^2)\) while the first three contributions (\({\mathcal {I}}_{\text {aft-aft}}\), \({\mathcal {I}}_{\text {bef-bef}}\), \({\mathcal {I}}_{\text {aft-bef}}\)) are \({\mathcal {O}}(\epsilon ^0)\). This implies that \({\mathcal {I}}_{\text {aft-in}}\) term is the first genuine non-eikonal contribution which is absent in the shockwave approximation.

\({\mathcal {I}}_{\text {in-bef}}\) term is another interference contribution that corresponds to emission of the gluon inside the target on one side of the cut and before on the other side and it can be written

This integral can be performed in the same manner as before, using Eq. (4.9) and the same definition of \(\alpha \). The result reads

Finally, \({\mathcal {I}}_{\text {in-in}}\) is the last contribution to the single inclusive gluon spectrum given in Eq. (4.6) and it accounts for the case in which the gluon is emitted inside the target on both sides of the cut. This contribution can be written

Here the gluon on the left side of the cut is emitted at longitudinal position \(y^+\) while the one on the right side of the cut is emitted at position \(\bar{y}^+>y^+\). By defining \(\xi =(\bar{y}^+-y^+)/L^+\) as the fractional longitudinal length traveled by the first gluon before the second one is emitted, \({\tilde{\alpha }}=y^+/L^+\) and using Eqs. (4.9), (3.14) and (3.4), we can write

In Eq. (4.19) we introduced a new variable \(\gamma =1-\xi -{\tilde{\alpha }}\) that corresponds to the longitudinal fraction of the target which the two gluons traverse at the same time. We would also like to mention that by expanding the trigonometric functions at leading order in powers of \(\epsilon \), \({\mathcal {I}}_{\text {in-in}}\) term is \({{{\mathcal {O}}}}(\epsilon ^4)\) and therefore this non-eikonal correction is sub-leading with respect to the other contributions.

Summing up, the non-eikonal single inclusive gluon spectrum is given by Eq. (4.6), with the explicit expressions for the contributions inside the integral provided in Eqs. (4.7), (4.10), (4.13) (4.15), (4.17) and (4.19) respectively. In the two following subsections, we investigate these results further.

### 4.2 Numerical results

This subsection is devoted to a numerical analysis of the results derived in Sect. 4.1 to illustrate the effects of the non-eikonal corrections on single inclusive gluon production in *pA* collisions. These non-eikonal effects, that stem from the finite longitudinal width of the target, are encoded in the so-called *non-eikonal parameter* \(\epsilon \) defined in Eq. (3.15). This non-eikonal parameter depends on the saturation momentum \(Q_s\), the gluon longitudinal momentum \(k^+\) and the longitudinal width of the target \(L^+\). For convenience of the analysis, the longitudinal momentum of the gluon \(k^+\) and the width of the target \(L^+\) can be written as (see Ref. [29] for a detailed discussion)

Thus, non-eikonal corrections depend on the values of the mass number *A* of the target, pseudorapidity \(\eta \) and transverse momentum \(p_\perp \equiv |\textbf{k}|\) of the produced gluon, and center of mass energy \(\sqrt{s_{\textrm{NN}}}\). We analyze the dependence of the single inclusive gluon spectrum on these variables.

In our analysis, in order to regulate the infrared divergences that appear, we introduce a mass in the denominators that lead to these divergences,^{Footnote 10} i.e., \(1/\textbf{q}^2 \rightarrow 1/(\textbf{q}^2 + m^2)\) in Eqs. (4.10),(4.13) and (4.17). We fix the value of the regulator to \(m = 0.6\) GeV and check the dependence of our results on this value. Our analysis shows a mild dependence on the value of the regulator: non-eikonal corrections become slightly smaller for smaller values of the infrared regulator. Moreover, we fix the value of the saturation scale to \(Q_s=\sqrt{2}\) GeV and take \(A=197\) (gold nucleus as target). Note that the non-eikonal parameter depends mildly, \(\propto A^{1/3}\), on *A*.

In Fig. 3, the ratio of the non-eikonal yield given in Eq. (4.6) to the strict eikonal result (where \(\epsilon =0\)) is plotted as a function of transverse momentum, \(p_\perp \), for three values of the center-of-mass energy per nucleon and for \(\eta =0\). The results show that while non-eikonal corrections decrease the yield at low \(p_\perp \), they cause an enhancement at \(p_\perp \gtrsim 1\) GeV. Moreover, the effect of non-eikonal contributions are about \(\sim 10\) % when \(\sqrt{s_{\textrm{NN}}} = 50\) GeV and \(p_\perp \sim 3\) GeV, but they are almost negligible for the highest center-of-mass energy \(\sqrt{s_{\textrm{NN}}} = 200\) GeV at RHIC. The results suggest that the impact of the finite width effects on phenomenology studies in single inclusive gluon production at relatively high energies is almost negligible.

In Fig. 4, we plot the same ratio as a function of the center-of-mass energy per nucleon for three values of pseudorapidity and for \(p_\perp =2\) GeV. The results indicate that for \(\eta =0\) the finite width effects are changing roughly between 10 % to 3 % up to between \(\sqrt{s_{\textrm{NN}}}\sim 40\) GeV to \(\sqrt{s_{\textrm{NN}}}\sim 100\) GeV, for \(\eta =0.5\) they lead to up to 4 % effect between \(\sqrt{s_{\textrm{NN}}}\sim 40\) GeV and \(\sqrt{s_{\textrm{NN}}}\sim 50\) GeV and for \(\eta =1\) they are almost negligible. Note that even though the results are plotted starting from \(\sqrt{s_{\textrm{NN}}} = 20\) GeV, at this low energy Bjorken-*x* of the target parton probed by the projectile is \(x \sim \frac{1 \ \textrm{GeV}}{\sqrt{s_{\textrm{NN}}}} e^{-\eta }\), i.e., not small, and therefore our approach cannot not be considered reliable. Therefore, we consider our results starting from \(\sqrt{s_{\textrm{NN}}}\sim 40\) GeV where our approach is still valid.

Finally, in Fig. 5, we plot the ratio as a function of the pseudorapidity of the produced gluon for three values of \(\sqrt{s_{\textrm{NN}}}\) and for \(p_\perp =2\) GeV. The conclusions are analogous to the ones extracted from previous figures: the non-eikonal corrections yield to up to \(\sim 6\) % effect for \(\eta <1\) and when \(\sqrt{s_{\textrm{NN}}} = 50\) GeV but it is almost negligible at higher energies at all values of the pseudorapidity.

### 4.3 Connection with earlier results in the literature

To finalize our study on the non-eikonal single inclusive gluon production, we make the explicit connection between our results and the ones obtained in Refs. [32, 33]. We would like to emphasize that the results derived in this manuscript do not adopt any kind of expansion in the non-eikonal parameter \(\epsilon ^2\), and are therefore valid to all orders. On the other hand, in Refs. [32, 33] a systematic expansion of the scalar propagator is proposed and the corrections are obtained up to next-to-next-to-eikonal accuracy. The expansion proposed in Refs. [32, 33] for the scalar propagator is performed in two steps. First, an expansion of the scalar propagator around its classical trajectory, i.e., the one that minimizes the path integral and leads to a linear path instead of the Brownian motion. Then, a second geometric expansion by assuming that the angle between the classical trajectory and the longitudinal axis with constant transverse coordinate is small since it is proportional to \(L^+/k^+\), and that the target field is weak far away from the classical solution. The corresponding result is an expansion in powers of \(\frac{L^+}{k^+}=\frac{2 i \epsilon ^2}{Q_s^2}\). It is effectively performed for the function \({\mathcal {R}}_{{\underline{k}}}(L^+,0;{\textbf {y}})\) which encodes the finite width target effects and is related to the scalar propagator via

The result of the expansion at next-to-next-to-eikonal accuracy reads [33]

where the color objects \(U^{i \cdots j}_{[\alpha ,\beta ]}(L^+,0; {\textbf {y}})\)^{Footnote 11} are referred to as *decorated Wilson lines*. They can be identified with the non-eikonal partners of the standard Wilson lines.

In order to make the connection between the two results, one can immediately realize that the Fourier transform of the non-eikonal dipole functions that we have introduced in Eqs. (3.3) and (3.2) can be written in terms of \({\mathcal {R}}_{{\underline{k}}}(L^+,0;{\textbf {y}})\) as

Therefore, by using the non-eikonal expansion given in Eq. (4.23) we can write the Fourier transform of the \(1^\textrm{st}\) order NE dipole function as

where the tensors \({\mathcal {O}}^{i \cdots j;l\cdots m}_{[\alpha ,\beta ];[\gamma ,\delta ]}({\textbf {r}})\) are referred to as *decorated dipoles* and are the leading eikonal corrections to the eikonal dipole function \({\mathcal {O}}({\textbf {r}})=d^{(0)}(\textbf{r})\). Furthermore, we note that in order to write Eq. (4.26), we assume translational invariance of the decorated dipoles and write \(\textbf{r}=\textbf{y}-\bar{\varvec{y}}\). The definitions of these dipole functions are:

On the other hand, the Fourier transform of the 2nd order NE dipole function can be written in the same fashion, reading

We can now make a one-to-one comparison between our result Eqs. (4.12) and (4.26) (and also between Eqs. (4.9) and (4.30)) by expanding our results to the appropriate order in order to obtain a parametrization of the decorated dipoles within the GBW model, i.e., assuming a Gaussian form. Expanding the Fourier transform of the 1st order NE dipole function given in Eq. (4.12) at next-to-next-to-eikonal accuracy, i.e., in terms of \(\epsilon ^2\) up to order \(\epsilon ^4\), we get

A term by term comparison of Eqs. (4.26) and (4.31) leads to the following parametrization of the decorated dipoles in our Gaussian ansatz:^{Footnote 12}

On the other hand, by using Eq. (4.9), the Fourier transform of the 2nd order NE dipole function can be written as

As mentioned earlier, the 2nd order NE dipole function does not encode any finite width effect in the case in where both longitudinal momenta are equal and therefore it is of \({\mathcal {O}}(\epsilon ^0)\). Thus, comparing Eqs. (4.30) and (4.38) and keeping in mind that the non-eikonal corrections in Eq. (4.30) vanish in the Gaussian approximation, we obtain the following constraints for the decorated dipoles:

We conclude that by comparing the model independent non-eikonal expansion given in Refs. [32, 33] with our approach, which relies on the Gaussian form of the dipole functions, we are able to obtain a parameterization of the decorated dipoles in the Gaussian ansatz that is valid for small size dipoles. This parameterization may be useful for including finite target width effects in future phenomenological studies.

## 5 Non-eikonal multigluon production

In this section we study the general case that corresponds to the non-eikonal *n*-gluon production whose generic expression is given in Eq. (2.7). Throughout this study, we adopt the Area Enhancement approximation in order to write the multipoles as non-eikonal two point functions (or dipoles). The validity of the AE argument when the finite width target effects are included in the computation is discussed at the end of Sect. 3. We would like to also mention that in Ref. [56] we investigated multigluon production in the strict eikonal limit where we introduced a diagramatic approach to study the 2*n*-point correlators within the AE argument. Our aim, in this section, is to generalize that approach to the non-eikonal case where the finite longitudinal width of the target is accounted for.

Let us start by writing the *n*-gluon spectrum given in Eq. (2.7) in terms of the simplified notations introduced in Ref. [56]. Within the AE model, the non-eikonal *n*-gluon spectra can be written as

where \(\chi =\{1,2,\dots ,2 n\}\) represents the total number of color charges in the amplitude and in the complex conjugate amplitude (i.e., on both sides of the cut) and \(\Pi (\chi )\) represents the set of partitions of \(\chi \) with disjoint pairs. The projectile is represented by function \(\Omega _i\) which is given by

with the property that for odd *i* the color source sits on the left side of the cut and for even *i* it sits on the right side. The two point correlator of the function \(\Omega _i\) is defined as^{Footnote 13}

On the other hand, the object \(\Lambda _\alpha \) represents the target and its interaction with the projectile. It is defined as^{Footnote 14}

with odd values of \(\alpha \) corresponding to the left side of the cut and even values of \(\alpha \) to the right side of the cut. Therefore, for even values of \(\alpha \) one takes the complex conjugate of the reduced amplitude.

Similar to the strict eikonal case, we impose the following constraints on the produced gluons:

namely the color, polarization and both transverse and longitudinal momenta of the produced gluons have to be the same on both sides of the cut. We would like to remind that the arguments introduced at the end of Sect. 3 to justify the use of AE model in the non-eikonal case where the finite longitudinal width effects are accounted for, are based on the fact that the fluctuation time of the target fields is much larger than its longitudinal extent. This implies that that chromo-electric domains do not decohere during the interaction of the projectile and the target. Moreover, it also allows one to argue that each domain has to be color neutral and the target average has to be a singlet. Therefore, we can write the two point correlators of functions \(\Lambda _\alpha \) as

Using the definition of the reduced amplitude given in Eq. (2.2), and the definitions of the eikonal and non-eikonal dipole functions introduced in Eqs. (3.3), (3.2) and (3.1), the 2-point correlator becomes

where we have introduced the notation \(\pm k_{a} \equiv (-1)^{a+1} k_a\). It corresponds to stating that when the respective momentum (either longitudinal or transverse) is on the left side of the cut, i.e., with *a* odd, it is multiplied by \((+1)\) and when it is on the right side it is multiplied by \((-1)\), changing its sign due to the Fourier transform. The explicit expressions for the functions *I* are given in Appendix B.

Thus, in order to compute the *n*-gluon spectra, one needs to evaluate the Wick expansions given in Eq. (5.1). On the projectile side, since we are adopting the MV model for the projectile configurations, projectile averages factorize into 2-point correlators defined in Eq. (5.3). On the target side, we are using the AE model to perform the averages over the target fields which also factorize into 2-point correlators of the reduced matrix amplitudes, given in Eq. (5.9). Expanding the \(\left( \frac{(2n)!}{2^n n!} \right) ^2\) terms in Eq. (5.1) and evaluating the 2-point correlators (on both the projectile and the target sides), we are able to compute the *n*-gluon spectra for any value of *n* neglecting, of course, contributions that are subleading in powers of \(\rho _p\). Since multigluon spectra contain a very large number of terms, it is convenient to use the diagrammatic approach introduced in Ref. [56].

Finally, let us note that the single inclusive gluon spectrum can be computed in the approach presented in this section by just writing

It is straightforward to see that if we use Eqs. (5.3) and (5.9), this equation reduces to Eq. (4.6). In the next subsection, we study the double inclusive gluon production by using the setup introduced here.

### 5.1 Non-eikonal double inclusive gluon production

Now we focus on the non-eikonal double inclusive gluon production – the \(n=2\) case in the multigluon spectrum in Eq. (5.1). In our set up, the double inclusive spectrum reads

where the 2-point functions for the projectile and the target sides are given by Eqs. (5.3) and (5.9) respectively.

Double inclusive gluon production has been studied extensively in the eikonal approximation (see Refs. [12,13,14,15,16,17,18]). In the eikonal limit, double inclusive gluon production spectrum contains not only a fully uncorrelated piece, but also Bose Enhancement and Hanbury-Brown-Twiss-like (HBT) correlations encoded in the terms that are suppressed by powers of \((N_c^2-1)\) with respect to the fully uncorrelated one. Moreover, it has been shown that the eikonal results give non-zero even harmonics but vanishing odd harmonics. This is due to an accidental symmetry of the CGC which is encoded in the double inclusive spectrum through its \(\textbf{k}_3\rightarrow -\textbf{k}_3\) symmetry, with \(\textbf{k}_3\) being the transverse momenta of the second gluon. Over the last decade, there have been many suggested ways to break this symmetry (see Ref. [11] for a review). One of the ways to break this symmetry and generate non-zero odd harmonics through CGC calculations is to include non-eikonal corrections, as has been shown in Refs. [28, 29] for proton–proton collisions.

The non-eikonal double inclusive spectrum given in Eq. (5.11) resembles its eikonal limit in many features. First, the term proportional to \(\big \langle \Omega _1 \Omega _2 \big \rangle _{p} \big \langle \Omega _3 \Omega _4 \big \rangle _{p} \big \langle \Lambda _1 \Lambda _2 \big \rangle _{T} \big \langle \Lambda _3 \Lambda _4 \big \rangle _{T}\) is the fully uncorrelated piece and it is the leading term in the \(N_c\) counting. Bose Enhancement and HBT-like correlations are encoded in the other terms in Eq. (5.11) and they are subleading in powers of \((N_c^2-1)\) and collision area (actually in \(B_p Q_s^2\)) just as in its eikonal limit. The main difference between the non-eikonal double inclusive spectrum given in Eq. (5.11) and its eikonal limit is that the non-eikonal corrections break the accidental symmetry of the CGC. Note that Eq. (5.11) is symmetric under the exchange of \(({{\underline{k}}}_3 \rightarrow - {\underline{k}}_3)\) where \({{\underline{k}}}_3\equiv (k_3^+, \textbf{k}_3)\),^{Footnote 15} but not under \(\textbf{k}_3 \rightarrow -\textbf{k}_3\). Therefore, non-eikonal corrections introduce an asymmetry in the spectrum with respect to the azimuthal angle.

In order to see this asymmetry explicitly, we compute the non-eikonal double inclusive gluon spectrum given in Eq. (5.11) numerically. For that purpose, we need to expand the terms inside the parenthesis in Eq. (5.11) and use Eqs. (5.3) and (5.9) as the definitions of the 2-point correlators for the projectile and the target averages respectively. Similar to the treatment in the single inclusive case, the infrared divergences that appear in the limit \(\textbf{q}_i \rightarrow 0\) in Eq. (5.9) are regulated by introducing a mass term \(m^2\) in the denominators, i.e., we perform the change \(1/\textbf{q}_i^2 \rightarrow 1/(\textbf{q}_i^2+m^2)\), where the value of the regulator is fixed to \(m = 0.4\) GeV. We also fix in the numerical evaluations \(N_c=3\), \(A^{1/3}=6\), \(B_p=6\) GeV\(^{-2}\) and \(Q_s^2 = 2\) GeV\(^{2}\).

The results are presented in Fig. 6, where we plot the solution of Eq. (5.11) and its respective eikonal limit as a function of the azimuthal angle \(\Delta \phi = \arccos \frac{\textbf{k}_1 \cdot \textbf{k}_3}{|\textbf{k}_1||\textbf{k}_3|}\). In this plot, we used \(\sqrt{s_{\textrm{NN}}}= 100\) GeV, \(\eta _1=\eta _3=0.2\) and \(|\textbf{k}_1|=|\textbf{k}_3|=1\) GeV. As discussed above, a clear asymmetry between the near- and away-side peaks appears in the non-eikonal double gluon spectrum. This observation confirms the breaking of the accidental symmetry of the CGC by introducing the non-eikonal corrections that here we check in the dilute-dense (*pA*) situation. Moreover, Fig. 6 shows that the non-eikonal spectrum differs from the eikonal one by \(4 \%\) in the near-side and \(8 \%\) in the away-side peak. In this kinematics, this is a slightly larger correction than the one in the single inclusive spectrum discussed in Sect. 4.2.

These results clearly indicate that non-eikonal corrections provide non-vanishing odd harmonics since the accidental symmetry of the CGC is broken by their inclusion. We leave a detailed analysis of the azimuthal harmonics for a future work.

## 6 Conclusions and outlook

In this work we generalized to proton–nucleus collisions the framework proposed in Refs. [28, 29] for proton–proton collisions. We computed multigluon spectra including the non-eikonal corrections that stem from the finite longitudinal width of the target. We use the Area Enhancement model that allows for an expansion of multipoles in terms of dipoles by neglecting contributions suppressed by powers of the overlap area of projectile and target. The formalism adopted in this work is inspired by the one used in jet quenching studies [53,54,55] where a non-zero longitudinal extent of the dense medium is essential to study parton propagation in the Quark Gluon Plasma. We adopt and extend such formalism to compute particle production and correlations in the Color Glass Condensate for the first time, where the effect of considering a finite medium width is small as compared to the jet quenching case where this is the main contribution to emission processes.

We have shown that the non-eikonal corrections in the target 2-point correlators can be computed analytically using the GBW approximation that assumes that the eikonal dipole function has a Gaussian form. This approximation is justified as long as the dipole size is \(\ll \Lambda _{\textrm{QCD}}^{-1}\). The results for the non-eikonal dipole functions are given in Eqs. (3.4), (3.14) and (3.17), where the finite width target effects are accounted for, encoded in the dimensionless non-eikonal parameter \(\epsilon ^2=\frac{Q_s^2 L^+}{2 i k^+}\). Using these results, we computed the single inclusive gluon spectrum beyond eikonal accuracy. The results of the single inclusive gluon production are summarized in Figs. 3, 4 and 5. They show that including the non-eikonal corrections to the single inclusive gluon spectrum results in between 2 % to 10 % effect for center-of-mass energies per nucleon smaller than 100 GeV and for pseudorapidities of the produced particle smaller than 1. Thus, we conclude that the non-eikonal corrections are not sizable for phenomenological studies in single inclusive gluon production at top RHIC or LHC energies. In Sect. 4.3, we compare the next-to-next-to-eikonal approximation of our results, i.e., an expansion up to order \(\epsilon ^4\), with the one presented in [32, 33]. We obtained a parameterization of the so-called decorated dipoles within the Gaussian ansatz that we have used here.

In Sect. 5 we generalize the framework introduced in [56] to compute multigluon production, to the non-eikonal case where the corrections to the eikonal approximation stem from the finite longitudinal with of the target. Like in its eikonal limit, our framework is based on the Area Enhancement model that we argue should be valid for a target with finite longitudinal size since the coherence length of the chromo-electric domains is much larger than the longitudinal extent of the target. We then focused on double inclusive gluon production and performed its numerical analysis (see Fig. 6). The results suggest two important outcomes. First, for not very high energies and rapidities, the non-eikonal corrections to the double inclusive gluon spectrum are slightly larger than those to the single inclusive case. Second, the non-eikonal corrections induce an azimuthal asymmetry in the double inclusive spectrum which demonstrates that the inclusion of such corrections breaks the accidental symmetry of the CGC, as in *pp* [28, 29]. Therefore, we expect to obtain non-vanishing odd harmonics from the non-eikonal double inclusive gluon spectrum in *pA* collisions. A detailed study of the effects of non-eikonal corrections to both even and odd azimuthal harmonics demands a dedicated numerical effort that we leave for the future.

## Data Availability Statement

This manuscript has no associated data or the data will not be deposited. [Authors’ comment: The results shown in the manuscript have been obtained by numerical implementation of the formulae that are provided. The corresponding code can be obtained from the authors by request.]

## Notes

Hereafter we employ the notation \(\int _{\textbf{x}}\equiv \int d^2 \textbf{x}\), \(\int _{\textbf{q}}\equiv \int d^2 \textbf{q}\).

This model has been called dipole approximation in [21], where its validity is also discussed.

See Appendix A of Ref. [56] for a detailed discussion of the validity of the AE model and its numerical comparison with the MV model for double dipole and quadrupole operators in the fundamental representation. Specifically, the difference between the Fourier transform for the double dipole in both models is just a few % and clearly decreasing with increasing collision area. But for the Fourier transform of the quadrupole such difference is considerably larger, up to 30 %, and the convergence with increasing collision area much slower than for the double dipole. Further numerical checks and discussions can be found in [21].

This definition of the saturation momentum is arbitrary and depends on the representation of the target parton that interacts with the projectile. However, it is convenient since it makes the dipole function independent of the value of the Casimir of the representation.

The meaning of the non-eikonal parameter in (3.15) is the following: \(Q_s\) is the typical transverse momentum of the fluctuations of the projectile that are resolved by the target, and \(k^+\) their plus momentum. Then \(\epsilon \) becomes the ratio of \(L^+\) to the formation time of such fluctuations, the latter being infinite (very large) in the strict eikonal limit. The deviations from eikonality are then given by \(\epsilon \) becoming larger and larger, which happens when the formation time becomes of the order or even smaller that the length of the target.

Note that its application in jet quenching calculations, where the length of the medium is not Lorentz contracted in this frame but also the longitudinal spread of the target gluons may be larger, seems more delicate.

Compared to other definitions in the literature, this color factor goes into the definition of \(\mu ^2\).

As explained before, this equation is very similar to the medium-induced radiation spectrum used in jet quenching calculations [72,73,74]. The main difference between this equation and those used in jet quenching calculations is that the reduced amplitude given in Eq. (2.2) and therefore the trace given in Eq. (4.3) include interference terms with the contribution where the gluon is emitted before the interaction of the source with the target, which are absent in the jet quenching framework [54, 75].

A comment is in order here. In the eikonal limit, the reduced matrix amplitude defined in Eq. (2.2) can be further factorized into the Lipatov vertex and the standard Wilson lines that describe the interaction between the projectile and the target (see for example [59]). In the eikonal limit, the mass term introduced here regulates the infrared divergence appearing in the Lipatov vertex. In the case of a non-eikonal treatment of the target one can not factorize the Lipatov vertex from the reduced matrix amplitude. However, the mass term introduced here regulates the infrared divergence appearing the kinematic factor that is the analog of the Lipatov vertex of the eikonal limit.

For simplicity of the expressions we have dropped the arguments of these functions. Moreover, since the definition of these color objects is lengthy we refer to Ref. [33] for their explicit expressions.

These parametrizations can also be obtained in a straightforward way from the definition of the decorated Wilson lines, assuming a local version of the GBW model with Eq. (3.7). We thank Carlota Andrés and Alexis Moscoso for this observation.

We would like to note that both the definition of the function \(\Omega _i\) (Eq. (5.2)) and the definition of its 2-point correlator (Eq. (5.3)) differ from the definitions introduced in Ref. [56] for the following reason. In that reference, the transverse momenta of each projectile color source was denoted as \(\textbf{k}-\textbf{q}\). In this manuscript, transverse momenta of each color source is denoted as \(\textbf{q}\). Therefore, compared to Ref. [56], we perform a change of variable \(\textbf{k}-\textbf{q}\rightarrow \textbf{q}\). The factor \((-1)^{i+1}\) in the definition of the 2-point correlator is introduced in order to account for the change of the sign in the transverse momentum when the source is on the right side of the cut, i.e., for even

*i*.The definition of \(\Lambda _{\alpha }\) introduced in this manuscript differs from its definition introduced in Ref. [56]. Here, since we account for the finite longitudinal width effects of the target, the reduced amplitude \(\overline{{{{\mathcal {M}}}}}^{a_\alpha b_\alpha }\) and therefore the object \(\Lambda _{\alpha }\) is more involved compared to its eikonal limit which was used in Ref. [56].

Two comments are in order here. First, in the notation that is adopted to write Eq. (5.11) the momenta of the second gluon is represented by \(k_3\). Second, the exchange of \(({{\underline{k}}}_3 \rightarrow - {{\underline{k}}}_3)\) in Eq. (5.1) is equivalent to exchanging indices 3 and 4, which obviously leaves the spectrum invariant.

While in the standard case the exponent is written as

*i*times the action, in our case we absorb this*i*in the definition of the corresponding parameters, masses and couplings, in order to simplify the notation. On the other hand, we assume convergence of the integrals which can be checked in the final results, but a more rigorous treatment would require a Wick rotation.Note that a different definition \({\textbf {B}}=\frac{m_1 {\textbf {z}}_1-m_2 {\textbf {z}}_2}{m_1-m_2}\) would decouple the oscillators but then the limit \(m_1\rightarrow m_2\) would become ill-defined. This is the reason why we choose the one in Eq. (A.24).

Because of the \(\delta \)-function in Eq. (5.9) one gets \(\textbf{q}_\beta =\textbf{q}_\alpha \), but due to the constrains, it needs to be evaluated with different sign.

## References

S. Jeon, U. Heinz, Introduction to hydrodynamics. Quark-Gluon Plasma

**5**, 131–187 (2016). https://doi.org/10.1142/9789814663717_0003P. Romatschke, U. Romatschke,

*Relativistic Fluid Dynamics In and Out of Equilibrium*, Cambridge Monographs on Mathematical Physics (Cambridge University Press, 2019). https://doi.org/10.1017/9781108651998. arXiv:1712.05815S. Schlichting, P. Tribedy, Collectivity in small collision systems: an initial-state perspective. Adv. High Energy Phys.

**2016**, 8460349 (2016). https://doi.org/10.1155/2016/8460349. arXiv:1611.00329B. Schenke, Origins of collectivity in small systems. Nucl. Phys. A

**967**, 105 (2017). https://doi.org/10.1016/j.nuclphysa.2017.05.017. arXiv:1704.03914C. Loizides, Experimental overview on small collision systems at the LHC. Nucl. Phys. A

**956**, 200 (2016). https://doi.org/10.1016/j.nuclphysa.2016.04.022. arXiv:1602.09138Z. Citron et al., Report from Working Group 5: Future physics opportunities for high-density QCD at the LHC with heavy-ion and proton beams. CERN Yellow Rep. Monogr.

**7**, 1159 (2019). https://doi.org/10.23731/CYRM-2019-007.1159. arXiv:1812.06772J.L. Nagle, W.A. Zajc, Small system collectivity in relativistic hadronic and nuclear collisions. Annu. Rev. Nucl. Part. Sci.

**68**, 211 (2018). https://doi.org/10.1146/annurev-nucl-101916-123209. arXiv:1801.03477F. Gelis, E. Iancu, J. Jalilian-Marian, R. Venugopalan, The color glass condensate. Annu. Rev. Nucl. Part. Sci.

**60**, 463 (2010). https://doi.org/10.1146/annurev.nucl.010909.083629. arXiv:1002.0333Y.V. Kovchegov, E. Levin,

*Quantum Chromodynamics at High Energy*, vol. 33 (Cambridge University Press, Cambridge, 2012). https://doi.org/10.1017/CBO9781139022187T. Lappi, L. McLerran, Some features of the glasma. Nucl. Phys. A

**772**, 200 (2006). https://doi.org/10.1016/j.nuclphysa.2006.04.001. arXiv:hep-ph/0602189T. Altinoluk, N. Armesto, Particle correlations from the initial state. Eur. Phys. J. A

**56**, 215 (2020). https://doi.org/10.1140/epja/s10050-020-00225-6. arXiv:2004.08185A. Dumitru, F. Gelis, L. McLerran, R. Venugopalan, Glasma flux tubes and the near side ridge phenomenon at RHIC. Nucl. Phys. A

**810**, 91 (2008). https://doi.org/10.1016/j.nuclphysa.2008.06.012. arXiv:0804.3858A. Dumitru, K. Dusling, F. Gelis, J. Jalilian-Marian, T. Lappi, R. Venugopalan, The ridge in proton–proton collisions at the LHC. Phys. Lett. B

**697**, 21 (2011). https://doi.org/10.1016/j.physletb.2011.01.024. arXiv:1009.5295Y.V. Kovchegov, D.E. Wertepny, Long-range rapidity correlations in heavy-light ion collisions. Nucl. Phys. A

**906**, 50 (2013). https://doi.org/10.1016/j.nuclphysa.2013.03.006. arXiv:1212.1195Y.V. Kovchegov, D.E. Wertepny, Two-gluon correlations in heavy-light ion collisions: energy and geometry dependence, IR divergences, and \(k_T\)-factorization. Nucl. Phys. A

**925**, 254 (2014). https://doi.org/10.1016/j.nuclphysa.2014.02.021. arXiv:1310.6701T. Altinoluk, N. Armesto, G. Beuf, A. Kovner, M. Lublinsky, Bose enhancement and the ridge. Phys. Lett. B

**751**, 448 (2015). https://doi.org/10.1016/j.physletb.2015.10.072. arXiv:1503.07126T. Altinoluk, N. Armesto, G. Beuf, A. Kovner, M. Lublinsky, Hanbury–Brown–Twiss measurements at large rapidity separations, or can we measure the proton radius in p-A collisions? Phys. Lett. B

**752**, 113 (2016). https://doi.org/10.1016/j.physletb.2015.11.033. arXiv:1509.03223T. Altinoluk, N. Armesto, G. Beuf, A. Kovner, M. Lublinsky, Quark correlations in the Color Glass Condensate: Pauli blocking and the ridge. Phys. Rev. D

**95**, 034025 (2017). https://doi.org/10.1103/PhysRevD.95.034025. arXiv:1610.03020M. Li, V.V. Skokov, First saturation correction in high energy proton-nucleus collisions. Part I. Time evolution of classical Yang-Mills fields beyond leading order. JHEP

**06**, 140 (2021). https://doi.org/10.1007/JHEP06(2021)140. arXiv:2102.01594M. Li, V.V. Skokov, First saturation correction in high energy proton-nucleus collisions. Part II. Single inclusive semi-hard gluon production. JHEP

**06**, 141 (2021). https://doi.org/10.1007/JHEP06(2021)141. arXiv:2104.01879M. Li, V.V. Skokov, First saturation correction in high energy proton–nucleus collisions. Part III. Ensemble averaging. JHEP

**01**, 160 (2022). https://doi.org/10.1007/JHEP01(2022)160. arXiv:2111.05304A. Kovner, V.V. Skokov, Does shape matter? \(v_2\) vs eccentricity in small x gluon production. Phys. Lett. B

**785**, 372 (2018). https://doi.org/10.1016/j.physletb.2018.09.001. arXiv:1805.09297T. Altinoluk, N. Armesto, A. Kovner, M. Lublinsky, V.V. Skokov, Angular correlations in pA collisions from CGC: multiplicity and mean transverse momentum dependence of \(v_2\). arXiv:2012.01810

A. Kovner, M. Lublinsky, V. Skokov, Exploring correlations in the CGC wave function: odd azimuthal anisotropy. Phys. Rev. D

**96**, 016010 (2017). https://doi.org/10.1103/PhysRevD.96.016010. arXiv:1612.07790Y.V. Kovchegov, V.V. Skokov, How classical gluon fields generate odd azimuthal harmonics for the two-gluon correlation function in high-energy collisions. Phys. Rev. D

**97**, 094021 (2018). https://doi.org/10.1103/PhysRevD.97.094021. arXiv:1802.08166A. Dumitru, V. Skokov, Anisotropy of the semiclassical gluon field of a large nucleus at high energy. Phys. Rev. D

**91**, 074006 (2015). https://doi.org/10.1103/PhysRevD.91.074006. arXiv:1411.6630A. Dumitru, A.V. Giannini, Initial state angular asymmetries in high energy p+A collisions: spontaneous breaking of rotational symmetry by a color electric field and C-odd fluctuations. Nucl. Phys. A

**933**, 212 (2015). https://doi.org/10.1016/j.nuclphysa.2014.10.037. arXiv:1406.5781P. Agostini, T. Altinoluk, N. Armesto, Non-eikonal corrections to multi-particle production in the Color Glass Condensate. Eur. Phys. J. C

**79**, 600 (2019). https://doi.org/10.1140/epjc/s10052-019-7097-5. arXiv:1902.04483P. Agostini, T. Altinoluk, N. Armesto, Effect of non-eikonal corrections on azimuthal asymmetries in the Color Glass Condensate. Eur. Phys. J. C

**79**, 790 (2019). https://doi.org/10.1140/epjc/s10052-019-7315-1. arXiv:1907.03668A. Accardi et al., Electron ion collider: the next QCD Frontier: understanding the glue that binds us all. Eur. Phys. J. A

**52**, 268 (2016). https://doi.org/10.1140/epja/i2016-16268-9. arXiv:1212.1701R. Abdul Khalek et al., Science Requirements and Detector Concepts for the Electron-Ion Collider: EIC Yellow Report. arXiv:2103.05419

T. Altinoluk, N. Armesto, G. Beuf, M. Martínez, C.A. Salgado, Next-to-eikonal corrections in the CGC: gluon production and spin asymmetries in pA collisions. JHEP

**07**, 068 (2014). https://doi.org/10.1007/JHEP07(2014)068. arXiv:1404.2219T. Altinoluk, N. Armesto, G. Beuf, A. Moscoso, Next-to-next-to-eikonal corrections in the CGC. JHEP

**01**, 114 (2016). https://doi.org/10.1007/JHEP01(2016)114. arXiv:1505.01400T. Altinoluk, A. Dumitru, Particle production in high-energy collisions beyond the shockwave limit. Phys. Rev. D

**94**, 074032 (2016). https://doi.org/10.1103/PhysRevD.94.074032. arXiv:1512.00279G.A. Chirilli, Sub-eikonal corrections to scattering amplitudes at high energy. JHEP

**01**, 118 (2019). https://doi.org/10.1007/JHEP01(2019)118. arXiv:1807.11435G.A. Chirilli, High-energy operator product expansion at sub-eikonal level. JHEP

**06**, 096 (2021). https://doi.org/10.1007/JHEP06(2021)096. arXiv:2101.12744T. Altinoluk, G. Beuf, A. Czajka, A. Tymowska, Quarks at next-to-eikonal accuracy in the CGC: forward quark-nucleus scattering. Phys. Rev. D

**104**, 014019 (2021). https://doi.org/10.1103/PhysRevD.104.014019. arXiv:2012.03886T. Altinoluk, G. Beuf, Quark and scalar propagators at next-to-eikonal accuracy in the CGC through a dynamical background gluon field. Phys. Rev. D

**105**, 074026 (2022). https://doi.org/10.1103/PhysRevD.105.074026. arXiv:2109.01620Y.V. Kovchegov, M.D. Sievert, Small-\(x\) helicity evolution: an operator treatment. Phys. Rev. D

**99**, 054032 (2019). https://doi.org/10.1103/PhysRevD.99.054032. arXiv:1808.09010F. Cougoulic, Y.V. Kovchegov, Helicity-dependent generalization of the JIMWLK evolution. Phys. Rev. D

**100**, 114020 (2019). https://doi.org/10.1103/PhysRevD.100.114020. arXiv:1910.04268Y.V. Kovchegov, D. Pitonyak, M.D. Sievert, Helicity evolution at small-x. JHEP

**01**, 072 (2016). https://doi.org/10.1007/JHEP01(2016)072. arXiv:1511.06737Y.V. Kovchegov, D. Pitonyak, M.D. Sievert, Helicity evolution at small \(x\): flavor singlet and non-singlet observables. Phys. Rev. D

**95**, 014033 (2017). https://doi.org/10.1103/PhysRevD.95.014033. arXiv:1610.06197Y.V. Kovchegov, M.G. Santiago, Quark sivers function at small \(x\): spin-dependent odderon and the sub-eikonal evolution. JHEP

**11**, 200 (2021). https://doi.org/10.1007/JHEP11(2021)200. arXiv:2108.03667F. Cougoulic, Y.V. Kovchegov, Helicity-dependent extension of the McLerran–Venugopalan model. Nucl. Phys. A

**1004**, 122051 (2020). https://doi.org/10.1016/j.nuclphysa.2020.122051. arXiv:2005.14688F. Cougoulic, Y.V. Kovchegov, A. Tarasov, Y. Tawabutr, Quark and gluon helicity evolution at small x: revised and updated. JHEP

**07**, 095 (2022). https://doi.org/10.1007/JHEP07(2022)095. arXiv:2204.11898shape Jefferson Lab Angular Momentum collaboration, First analysis of world polarized DIS data with small-x helicity evolution. Phys. Rev. D

**104**, L031501 (2021). https://doi.org/10.1103/PhysRevD.104.L031501. arXiv:2102.06159Y.V. Kovchegov, D. Pitonyak, M.D. Sievert, Small-\(x\) asymptotics of the quark helicity distribution: analytic results. Phys. Lett. B

**772**, 136 (2017). https://doi.org/10.1016/j.physletb.2017.06.032. arXiv:1703.05809J. Jalilian-Marian, Elastic scattering of a quark from a color field: longitudinal momentum exchange. Phys. Rev. D

**96**, 074020 (2017). https://doi.org/10.1103/PhysRevD.96.074020. arXiv:1708.07533J. Jalilian-Marian, Quark jets scattering from a gluon field: from saturation to high \(p_t\). Phys. Rev. D

**99**, 014043 (2019). https://doi.org/10.1103/PhysRevD.99.014043. arXiv:1809.04625J. Jalilian-Marian, Rapidity loss, spin, and angular asymmetries in the scattering of a quark from the color field of a proton or nucleus. Phys. Rev. D

**102**, 014008 (2020). https://doi.org/10.1103/PhysRevD.102.014008. arXiv:1912.08878A.V. Sadofyev, M.D. Sievert, I. Vitev, Ab initio coupling of jets to collective flow in the opacity expansion approach. Phys. Rev. D

**104**, 094044 (2021). https://doi.org/10.1103/PhysRevD.104.094044. arXiv:2104.09513C. Andres, F. Dominguez, A.V. Sadofyev, C.A. Salgado, Jet broadening in flowing matter—resummation. arXiv:2207.07141

J. Casalderrey-Solana, C.A. Salgado, Introductory lectures on jet quenching in heavy ion collisions. Acta Phys. Polon. B

**38**, 3731 (2007). arXiv:0712.3443Y. Mehtar-Tani, J.G. Milhano, K. Tywoniuk, Jet physics in heavy-ion collisions. Int. J. Mod. Phys. A

**28**, 1340013 (2013). https://doi.org/10.1142/S0217751X13400137. arXiv:1302.2579J.-P. Blaizot, Y. Mehtar-Tani, Jet structure in heavy ion collisions. Int. J. Mod. Phys. E

**24**, 1530012 (2015). https://doi.org/10.1142/S021830131530012X. arXiv:1503.05958P. Agostini, T. Altinoluk, N. Armesto, Multi-particle production in proton–nucleus collisions in the color glass condensate. Eur. Phys. J. C

**81**, 760 (2021). https://doi.org/10.1140/epjc/s10052-021-09475-0. arXiv:2103.08485A. Kovner, A.H. Rezaeian, Double parton scattering in the CGC: double quark production and effects of quantum statistics. Phys. Rev. D

**96**, 074018 (2017). https://doi.org/10.1103/PhysRevD.96.074018. arXiv:1707.06985A. Kovner, A.H. Rezaeian, Multiquark production in \(p+A\) collisions: quantum interference effects. Phys. Rev. D

**97**, 074008 (2018). https://doi.org/10.1103/PhysRevD.97.074008. arXiv:1801.04875T. Altinoluk, N. Armesto, A. Kovner, M. Lublinsky, Double and triple inclusive gluon production at mid rapidity: quantum interference in p-A scattering. Eur. Phys. J. C

**78**, 702 (2018). https://doi.org/10.1140/epjc/s10052-018-6186-1. arXiv:1805.07739Y. Mehtar-Tani, Relating the description of gluon production in pA collisions and parton energy loss in AA collisions. Phys. Rev. C

**75**, 034908 (2007). https://doi.org/10.1103/PhysRevC.75.034908. arXiv:hep-ph/0606236N. Armesto, L. McLerran, C. Pajares, Long range forward–backward correlations and the color glass condensate. Nucl. Phys. A

**781**, 201 (2007). https://doi.org/10.1016/j.nuclphysa.2006.10.074. arXiv:hep-ph/0607345C. Özonder, Triple-gluon and quadruple-gluon azimuthal correlations from glasma and higher-dimensional ridges. Phys. Rev. D

**91**, 034005 (2015). https://doi.org/10.1103/PhysRevD.91.034005. arXiv:1409.6347C. Özönder, Predictions on three-particle azimuthal correlations in proton-proton collisions. Turk. J. Phys.

**42**, 78 (2018). https://doi.org/10.3906/fiz-1710-6. arXiv:1712.05571T. Altinoluk, N. Armesto, D.E. Wertepny, Correlations and the ridge in the Color Glass Condensate beyond the glasma graph approximation. JHEP

**05**, 207 (2018). https://doi.org/10.1007/JHEP05(2018)207. arXiv:1804.02910L.D. McLerran, R. Venugopalan, Computing quark and gluon distribution functions for very large nuclei. Phys. Rev. D

**49**, 2233 (1994). https://doi.org/10.1103/PhysRevD.49.2233. arXiv:hep-ph/9309289L.D. McLerran, R. Venugopalan, Green’s functions in the color field of a large nucleus. Phys. Rev. D

**50**, 2225 (1994). https://doi.org/10.1103/PhysRevD.50.2225. arXiv:hep-ph/9402335K.J. Golec-Biernat, M. Wusthoff, Saturation effects in deep inelastic scattering at low Q**2 and its implications on diffraction. Phys. Rev. D

**59**, 014017 (1998). https://doi.org/10.1103/PhysRevD.59.014017. arXiv:hep-ph/9807513K.J. Golec-Biernat, M. Wusthoff, Saturation in diffractive deep inelastic scattering. Phys. Rev. D

**60**, 114023 (1999). https://doi.org/10.1103/PhysRevD.60.114023. arXiv:hep-ph/9903358B.G. Zakharov, Light cone path integral approach to the Landau–Pomeranchuk–Migdal effect. Phys. Atom. Nucl.

**61**, 838 (1998). arXiv:hep-ph/9807540J.-P. Blaizot, F. Dominguez, E. Iancu, Y. Mehtar-Tani, Medium-induced gluon branching. JHEP

**01**, 143 (2013). https://doi.org/10.1007/JHEP01(2013)143. arXiv:1209.4585L. Apolinário, N. Armesto, J.G. Milhano, C.A. Salgado, Medium-induced gluon radiation and colour decoherence beyond the soft approximation. JHEP

**02**, 119 (2015). https://doi.org/10.1007/JHEP02(2015)119. arXiv:1407.0599R. Baier, Y.L. Dokshitzer, A.H. Mueller, S. Peigne, D. Schiff, Radiative energy loss and p(T) broadening of high-energy partons in nuclei. Nucl. Phys. B

**484**, 265 (1997). https://doi.org/10.1016/S0550-3213(96)00581-0. arXiv:hep-ph/9608322R. Baier, Y.L. Dokshitzer, A.H. Mueller, S. Peigne, D. Schiff, Radiative energy loss of high-energy quarks and gluons in a finite volume quark-gluon plasma. Nucl. Phys. B

**483**, 291 (1997). https://doi.org/10.1016/S0550-3213(96)00553-6. arXiv:hep-ph/9607355B.G. Zakharov, Radiative energy loss of high-energy quarks in finite size nuclear matter and quark-gluon plasma. JETP Lett.

**65**, 615 (1997). https://doi.org/10.1134/1.567389. arXiv:hep-ph/9704255U.A. Wiedemann, Gluon radiation off hard quarks in a nuclear environment: opacity expansion. Nucl. Phys. B

**588**, 303 (2000). https://doi.org/10.1016/S0550-3213(00)00457-0. arXiv:hep-ph/0005129C. Grosche, F. Steiner,

*Handbook of Feynman Path Integrals*, vol. 145 (1998)

## Acknowledgements

PA, NA and FD have received financial support from Xunta de Galicia (Centro singular de investigación de Galicia accreditation 2019-2022), by European Union ERDF, by the “María de Maeztu” Units of Excellence program MDM-2016-0692, and by the Spanish Research State Agency under project PID2020-119632GB-I00. JGM acknowledges financial support from Fundaç\(\tilde{\textrm{a}}\)o para a Ci\(\hat{\textrm{e}}\)ncia e a Tecnologia (Portugal) under projects CERN/FIS-PAR/0024/2019 and CERN/FIS-PAR/0032/2021, and he gratefully acknowledges the hospitality of the CERN theory group. This work has been performed in the framework of the European Research Council project ERC-2018-ADG-835105 YoctoLHC and the MSCA RISE 823947 “Heavy ion collisions: collectivity and precision in saturation physics” (HIEIC), and has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No. 824093.

## Author information

### Authors and Affiliations

### Corresponding author

## Appendices

### A Harmonic oscillator path integrals

In this section we compute solve the path integrals that appear in Sect. 3. Before starting our discussion we define the following matrix:

with *x* a positive real scalar and \(i,j=1,\dots , N\). Therefore \(F_N\) is a \(N \times N\) tridiagonal matrix with *x* in the main diagonal and \(-1\) in the neighbour diagonals. Since this matrix will be a key ingredient in the calculation of this section we will study some of its properties.

By performing the Laplace expansion of this matrix it is easy to see that it follows the following recursive equation:

The solution of this equation is

where \(x_{\pm }=x \pm \sqrt{x^2-4}\).

The inverse of \(F_N(x)\) can be computed by using the expression of the inverse of a block matrix:

Noting that \(F_N(x)\) can be written as a block matrix with \(\textbf{A}=x\), \(\textbf{D}=F_{N-1}(x)\), \(\textbf{C}_i=-\delta _{i, 2}\) and \(\textbf{B}_j=-\delta _{j, 2}\), we can compute the inverse of its component (1, 1) by using the first equality of (A.4):

Analogously, we could use the second equality of (A.4) with \(\textbf{A}=F_{N-1}(x)\), \(\textbf{D}=x\), \(\textbf{C}_j=-\delta _{j, N-1}\) and \(\textbf{B}_i=-\delta _{i, N-1}\) to obtain the same recursive equation for \([F_N(x)^{-1}]_{N,N}\). Solving (A.5) we obtain

Inspecting (A.4), we realize that the rest of the components of \(F_N^{-1}(x)\) can be computed just as a function of its (1, 1) and (*N*, *N*) components. In this manuscript we will only need the (1, *N*) and (*N*, 1) components of this matrix which can be written

We now proceed to solve the path integrals.^{Footnote 16} The first path integral that we encounter in Sect. 3 is the one of a single harmonic oscillator of mass *m* and coupling \(\kappa \):

where, from (3.13), \(m=\frac{i k^+}{2}\) and \(\kappa =\frac{Q_s^2}{4 L^+}\). By discretizing the longitudinal direction into *N* slices we can write this equation as

where we have defined \(\Delta ^+=x^+-y^+\) and

In order to solve this integral it is convenient to write it in a matrix form. In order to do so we define the vector \(\textbf{x}=({\textbf {z}}_1,\dots ,{\textbf {z}}_{N-1})\) and we can make the following simplifications:

and

where the matrix \(F_N(x)\) is defined in (A.1).

Thus, by introducing the vector \(\textbf{x}\), we can write (A.9) as a multidimensional Gaussian integral:

where we have defined the vector

Since the solution of a multidimensional Gaussian integral is

where *M* is a constant matrix \(D \times D\) and *B* is a constant *D*-dimensional vector, we can write (A.13) as

Using (A.3) and taking the \(N\rightarrow \infty \) limit we obtain

where we define

and we have used the fact that

On the other hand, using (A.6) and (A.7) and taking the \(N\rightarrow \infty \) limit we get

where we have used that

Therefore we obtain the well known solution of the harmonic oscillator path integral [76]:

The second path integral that appears in Sect. 3 is the one analogous to the one for two coupled harmonic oscillators of masses \(m_1\) and \(-m_2\) and coupling \(\kappa \):

In order to solve this integral it is convenient to make the change of variables^{Footnote 17}

in such a way that the integral can be rewritten as

where we have introduced

In the discrete limit we can write the path integral as

with

The next step is to write, again, the path integral in a matrix form, i.e., as a multidimensional Gaussian integral. We therefore define \(\textbf{x}_{\textbf{r}}=({\textbf {r}}_1,\dots ,{\textbf {r}}_{N-1})\) and \(\textbf{x}_{\textbf{b}}=({\textbf {b}}_1,\dots ,{\textbf {b}}_{N-1})\) in order to make the following simplification:

Then the path integral can be simplified further by defining the matrix

the vector

and \(\textbf{x}=({\textbf {r}}_1,\dots ,{\textbf {r}}_N,{\textbf {b}}_0,\dots ,{\textbf {b}}_N)\). Thus we can write it as

where we have used (A.15) in order to solve the integral over \(\textbf{x}\).

The determinant of the matrix \({{{\mathcal {A}}}}\) can be evaluated by using the following property of block matrices:

By using this property we can write

where we have used the fact that \(F_N(x)+b F_N(y) =(1+b) F_N\left( \frac{x+by}{1+b}\right) \) and defined

Therefore, using (A.19) we can compute the limit \(N\rightarrow \infty \) in the prefactor in (A.32) to read

where

On the other hand, the limit \(N\rightarrow \infty \) of the argument of the exponential in (A.32) can be evaluated by using (A.4) and writing the inverse of the matrix \({{{\mathcal {A}}}}\) as

Therefore we can write

where we have introduced for simplicity

Finally, using (A.19) and (A.21) and the definitions of the coefficients \(a_i\) it is straightforward to see that

Thus we can write the solution of the double harmonic oscillator as

or, equivalently,

where \(\omega _{\pm }^2=\omega _1^2\pm \omega _2^2\). Note that this equation provides a well-defined limit for \(m_1\rightarrow m_2\equiv \omega _1\rightarrow \omega _2\).

### B Derivation of the terms in the target 2-point function

In this section we derive the explicit expressions of the functions appearing in Eq. (5.9). We start by using the definition of the reduced amplitude given in Eq. (2.2) and the definitions of the eikonal and non-eikonal dipole functions introduced in Eqs. (3.3), (3.2) and (3.1). Thus, the 2-point correlator given in Eq. (5.8) reads

where we have introduced the notation \(f_\alpha =y_\alpha ^+/L^+\) (and analogous expression for \(f_\beta \)) for convenience. Note that the factors \((-1)^\alpha \) account for whether the object \(\Lambda _\alpha \) is sitting on the right or on the left side of the cut. It is also straightforward to realize that after setting \(\alpha =1\) and \(\beta =2\), this equation reduces to the single inclusive case computed in Eq. (4.3).

In Eq. (B.1), integrals over the transverse coordinates can be performed using the definitions of the non-eikonal dipole functions given in Eqs. (3.4), (3.14) and (3.17), and also the solution of the Gaussian-like integral given in Eq. (4.5). In order to shorten the expressions, we use the notation \(\pm k_{a} \equiv (-1)^{a+1} k_a\) introduced in Sect. 5. Thus, upon performing the integrations over the transverse coordinates and introducing this notation we arrive at Eq. (5.9).

The first term in Eq. (5.9) corresponds to the case where both gluons are emitted after the interaction of the source with the target and it reads

Note that this term is independent of the sign of the transverse momenta since its argument does not carry the ± notation introduced in Eq. (5.9). Moreover, due to the appearance of the \(\delta \)-function in the first line of Eq. (5.9), this term is invariant under the change \(\alpha \leftrightarrow \beta \).

The second contribution in Eq. (5.9) corresponds to the case where the gluon represented by the index \(\alpha \), which we refer to as gluon \(\alpha \), is emitted after the interaction of the source with the target and the gluon \(\beta \) is emitted before. On the other hand, the inverse situation is given by the mirror term \(\alpha \leftrightarrow \beta \). Upon the performing the integrations, this term can be written as

where the \(k^+\) dependence is given implicitly in the definition of \(\epsilon \) (see Eq. (3.15)). This contribution carries dependence only on \(k_\beta ^+\) (and not on \(k_\alpha ^+\)) due to the fact that only the gluon \(\beta \) traverses the target.

The third contribution in Eq. (5.9) corresponds to the case in which both gluons are emitted before the target and it can be written as

where we have introduced the following notations for convenience:

and

This contribution to the non-eikonal multigluon spectra is somewhat peculiar. In a situation where the gluons have the same momenta but are sitting at different sides of the cut, which in our notation corresponds to \(\alpha =2m-1\) and \(\beta =2m\) (or vice-versa) for any integer *m*, but \(k_\beta ^+ = k_\alpha ^+\) and \(\textbf{k}_\beta =\textbf{k}_\alpha \) due to the constraint introduced in Eq. (5.5). However, since the momenta on the right side of the cut are evaluated with a different sign in Eq. (5.9), one would have \(\epsilon _+ \rightarrow 0\), \(\textbf{K} \rightarrow \textbf{k}_\alpha \) and^{Footnote 18}\(\textbf{Q} \rightarrow \textbf{q}_\alpha \), so that the non-eikonal corrections vanish and it leads to the same contribution as in the case of single inclusive production given in Eq. (4.10). On the other hand, Eq. (B.4) is symmetric under the exchange \(\alpha \leftrightarrow \beta \) as expected.

The fourth contribution in Eq. (5.9) corresponds to the case where the gluon \(\alpha \) is emitted after the interaction of the source with the target and the gluon \(\beta \) is emitted inside the target, with the opposite situation given by the mirror term \(\alpha \leftrightarrow \beta \). The expression of this contribution reads

where we have defined \({\tilde{f}}_\beta =1-f_\beta =\frac{L^+-y_\beta ^+}{L^+}\) that corresponds to the fraction of the medium traversed by gluon \(\beta \). This result is identical to the one evaluated for the single inclusive case given in Eq. (4.15), which can be explained as follows. The main difference between the multigluon and single gluon calculations relies on the form of the 2nd order NE dipole function. For the multigluon production presented in this section, the 2nd order NE dipole functions, in general, are evaluated at different longitudinal momenta whereas in the single inclusive case they are evaluated at the same longitudinal momenta. However, some of the contributions to the multigluon production, namely aft-aft, aft-bef and aft-in contributions (given in Eqs. (B.2), (B.3) and (B.7) respectively), contain only one longitudinal momentum and therefore their expressions are almost the same as the ones in the single inclusive case up to some overall factors.

The fifth contribution in Eq. (5.9) corresponds to the case where the gluon \(\alpha \) is emitted before the interaction of the source with the target and the gluon \(\beta \) is emitted inside the target, with the opposite situation given by its mirror term \(\alpha \leftrightarrow \beta \). We can write this contribution as

where we have introduced \(\Delta \textbf{k} = \textbf{k}_\alpha + \textbf{k}_\beta \). Since in this contribution we are including the non-eikonal effects that stem from the difference between the two longitudinal momenta, in the fifth term of Eq. (B.1) we performed the derivative with respect to \(\textbf{y}^{\lambda _\beta }\) and integrated over an extra transverse position \(\textbf{u}\). We would like to note that in the case of having the two momenta equal, i.e., when \(\Delta \textbf{k} \rightarrow 0\), \(\epsilon _+ \rightarrow 0\), \(\textbf{K} \rightarrow \textbf{k}_\alpha \) and \(\textbf{Q}\rightarrow \textbf{q}_\alpha \), this expression simplifies to the one presented for the single inclusive case, Eq. (4.17), up to some prefactors.

Finally, the last contribution to Eq. (5.9), \(I^{\lambda _\alpha \lambda _\beta }_{\mathrm{in-in}} \Big (\textbf{k}_\alpha ,\textbf{k}_\beta , \textbf{q}_\beta ; k_\alpha ^+, k_\beta ^+ ; f_\alpha , f_\beta \Big )\) corresponds to the case when both gluons are emitted inside the target. This expression involves an extra derivative with respect to \(\textbf{y}^{\lambda _\alpha }\) compared to the in-bef contribution given in Eq. (B.8). It can be written as

where \(\epsilon _{\pm }=\epsilon _\alpha ^2 \pm \epsilon _\beta ^2\). The five 2-dimensional integrals in this expression are of the form given in Eq. (4.5) and therefore they can be solved analytically. However, since its final form is too long we do not write it explicitly in this manuscript.

## Rights and permissions

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Funded by SCOAP^{3}. SCOAP^{3} supports the goals of the International Year of Basic Sciences for Sustainable Development.

## About this article

### Cite this article

Agostini, P., Altinoluk, T., Armesto, N. *et al.* Multiparticle production in proton–nucleus collisions beyond eikonal accuracy.
*Eur. Phys. J. C* **82**, 1001 (2022). https://doi.org/10.1140/epjc/s10052-022-10962-1

Received:

Accepted:

Published:

DOI: https://doi.org/10.1140/epjc/s10052-022-10962-1