Gluon TMD in particle production from low to moderate x

We study the rapidity evolution of gluon transverse momentum dependent distributions appearing in processes of particle production and show how this evolution changes from small to moderate Bjorken x.


Introduction
The TMDs [1][2][3] (also called unintegrated parton distributions) are widely used in the analysis of various scattering processes like SIDIS or Drell-Yan. The TMD generalizes the usual concept of parton density by allowing PDFs to depend on intrinsic transverse momenta in addition to the usual longitudinal momentum fraction variable. At low energies the relevant quantities are quark TMDs and there is a vast literature on the application JHEP06(2016)164 of quark TMDs for analysis of cross sections of processes measured at JLab and elsewhere (see e.g. refs. [4][5][6][7][8][9][10][11][12][13][14][15][16][17][18], for review see also refs. [19][20][21]). However, since at the future EIC accelerator the majority of the produced particles will be gluons one needs to study also the evolution of gluon TMDs. Moreover, the EIC energies may be in the intermediate region between hard physics described by linear CSS evolution [22] and low-x physics described by non-linear BK/JIMWLK evolution [23][24][25][26][27][28][29][30][31][32][33] so one needs to study the transition of the evolution of gluon TMDs between these two regimes. 1 The gluon TMD (unintegrated gluon distribution) is defined as [36] D(x B , k ⊥ , η) = d 2 z ⊥ e i(k,z) ⊥ D(x B , z ⊥ , η), where |P is an unpolarized target with momentum p (typically proton) and n is a light-like vector. Hereafter we use the notation F a ξ (z ⊥ + un) ≡ n µ gF m µξ (un + z ⊥ )[un + z ⊥ , −∞n + z ⊥ ] ma (1.2) where [x, y] denotes straight-line gauge link connecting points x and y: There are more involved definitions with eq. (1.1) multiplied by some Wilson-line factors [3,37] following from CSS factorization [22] but we will discuss the "primordial" TMD (1.1). It is well known, however, that gluon TMDs are not universal in a sense that the direction of gauge links providing gauge invariance depends on the type of processes under consideration, see ref. [38]. For example, TMDs entering the description of processes particle production have light-like gauge links starting at minus infinity as in eq. (1.2), but TMDs which appear in the analysis of semi-inclusive processes have gauge links stretching to plus infinity (so the corresponding expression for TMDs is obtained by replacement −∞ ↔ ∞ in eq. (1.2)). For a more complicated processes the structure of gauge links may be even more involved, see e.g. ref. [39].
In our recent paper [40] we have obtained the leading-order evolution equation for gluon TMDs for semi-inclusive processes like semi-inclusive deep inelastic scattering (SIDIS). The obtained equation describes the rapidity evolution of gluon TMDs in the whole region of small to moderate Bjorken x B and for any transverse momentum. It interpolates between the linear DGLAP and Sudakov evolution equations at moderate x B and the non-linear BK equation for small x B . In this paper we extend our analysis to the case of gluon TMDs appearing in particle production processes with gauge links extending to minus infinity in the light-cone (LC) time direction. The analysis is very close to the study of our paper [40] so we will streamline the presentation of technical details paying attention to differences between these two cases (with links going to plus or minus infinity). The final evolution equations are similar (but in general not the same!) to TMDs with gauge links extending to plus infinity.

JHEP06(2016)164
The paper is organized as follows. In section 2 we remind the logic of rapidity factorization for the inclusive particle production and rapidity evolution. In section 3 we discuss rapidity evolution of gluon TMDs and calculate the leading-order kernel of the evolution equation. We present the final form of the evolution equation in section 4 and discuss BK, Sudakov and DGLAP limits in section 5 and linearized equation in section 6. Section 7 contains conclusions and outlook. The necessary formulas for propagators near the light cone and in the shock-wave background can be found in appendices.

TMDs in particle production
To simplify the description of particle production, let us consider the model where a (colorless) scalar particle can be produced by gluon-gluon fusion through the vertex coming from the Lagrangian One may consider this as a model of Higgs production by gluon fusion in the region where transverse momentum of produced Higgs boson is smaller than the mass of the top quark. Let us consider the production of this Φ-boson in the high-energy scattering of a virtual photon with virtuality ∼ few GeV off the hadron target. As demonstrated in the appendix A the cross section of Φ-boson production can be represented by a double functional integral × Ψ * p ( Ã (t i ),ψ(t i ))e −iS QCD (Ã,ψ) e iS QCD (A,ψ)j µ (w)F 2 (x)F 2 (y)j ν (0)Ψ p ( A(t i ), ψ(t i )) where Ψ p are proton wave functionals at the initial time t i → −∞. (The boundary condi-tionÃ( x, t f → ∞) = A( x, t f → ∞) and similar condition for quark fields reflects the sum over all intermediate states X). We will analyze the energy dependence of this cross section using the high-energy OPE in Wilson lines. To this end, we integrate over rapidities greater than the rapidity of the produced Φ-boson Y > η φ and leave the fields with Y < η φ to be integrated over later. The result of the integration over Y > η φ is the coefficient function (called "impact factor") in front of the Wilson-line operator(s) made of gluons (and quarks) with rapidities Y < η φ . (Strictly speaking, we integrate over rapidities Y > η φ − so the vertex of Φ-boson production is included into the impact factor). To make connections with parton model we will have in mind the frame where target's velocity is large and call the small α fields by the name "fast fields" and large α fields by "slow" fields. Of course, "fast" vs "slow" depends on frame but we will stick to naming fields as they appear in the projectile's frame. (Note that in refs. [23,24] the terminology is opposite, as appears in the target's frame). As discussed in refs. [23,24], the interaction of "slow" gluons of large Y with "fast" fields of small Y is described by eikonal gauge factors and the integration over slow fields results in Feynman diagrams in the background of fast fields which form a thin shock wave JHEP06(2016)164 due to Lorentz contraction. 2 In the spirit of high-energy OPE, the rapidity of the gluons is restricted from above by the "rapidity divide" η separating the impact factor and the matrix element so the proper definition of U x is where the Sudakov variable α is defined as usual, k = αp 1 +βp 2 +k ⊥ . We define the light-like vectors p 1 and p 2 close to projectile and target's momenta q and p so that q = p 1 + q 2 s p 2 and p = p 2 + m 2 s p 1 . We use metric g µν = (1, −1, −1, −1) so that p · q = (α p β q + α q β p ) s 2 − (p, q) ⊥ . For the coordinates we use the notations x • ≡ x µ p µ 1 and x * ≡ x µ p µ 2 for dimensionless lightcone coordinates (x * = s 2 x + and x • = s 2 x − ). In accordance with general background-field formalism we separate the gluon field into the "classical" background part and "quantum" part where the "classical" fields are fast (α < σ = e η ) and "quantum" fields are slow (α > σ = e η ). It should be emphasized that our "classical" field does not satisfy the equation D µ F cl µν = 0; rather, (D µ F cl µν ) a = −gψγ ν t a ψ where ψ are the "classical" (i.e. fast) quark fields. The first-order term in the expansion of the operator F m •i (y * , y ⊥ )[y * , −∞] ma y in quantum fields has the form (to save space, we omit the label cl from classical fields).
In the leading order the impact factor is given by the diagram shown in figure 1. The quark propagator in the external field has the form ψ (x)ψ(y) (2.5) x * ,y * <0 where σ = e η is the lower rapidity cutoff for the impact factor (and upper cutoff for α's in Wilson lines). Hereafter we use Schwinger's notations 2 An exceptional case discussed later is when the transverse momenta of the external field are much smaller than the characteristic transverse momenta in the impact factor. In this case the "shock wave" is no longer narrow and one needs the light-cone approximation rather than the shock-wave one. However, if the virtuality of the photon is ∼ few GeV the characterisic transverse momenta of the impact factor and of the fast "external fields" are of the same order of magnitude so the shock-wave approximation is applicable. Figure 1. Rapidity factorization for particle production. The dashed lines denote gauge links.
Note that unlike the case of total cross section, here we consider particle production so the gluon lines in figure 1 terminate at the Φ-boson emission point leading to gluon TMDs rather than proper Wilson lines (stretching from minus to plus infinity in LC-time direction). Indeed, the gluon propagator with one point in the shock wave has the form of the free propagator multiplied by the gauge link going from point y to −∞ in the p 1 direction [40]: Since the propagators (2.5) and (2.7) have simple structure one can calculate the integrals in figure 1 and the result has the form where tr{. . .} is the color trace in the fundamental representation and I ij µν (z 1 ⊥ , z 2 ⊥ , x ⊥ , y ⊥ ; σ) is the impact factor with the lower rapidity cutoff η = ln σ. 3 Hereafter we use the short-

JHEP06(2016)164
hand notations for gauge links and As discussed in refs. [23,24], the fast fields at light-cone time ±∞ are pure gauge so the precise form of the contour in eq. (2.10) is irrelevant. The calculation of the impact factor I ij µν (z 1 ⊥ , z 2 ⊥ , x ⊥ , y ⊥ ; η) is similar to the calculation of the NLO photon impact factor for the DIS structure functions carried out in refs. [41,42]. Since the explicit form of I ij µν is irrelevant for our purpose of finding the evolution of gluon TMDs and since in the real life the contribution of the diagram shown in figure 1 is a tiny correction to the total cross section of Higgs production in DIS we did not attempt to calculate this impact factor. In the case of proton-proton scattering the impact factor should be given by another gluon TMD made of Wilson lines stretched in p 2 direction. We intend to discuss the obtained factorization in a separate publication.
As demonstrated in appendix A (see eq. (A.8)), the double functional integral (2.8) represents the matrix element Note that all the gluon operators in the r.h.s. of this equation are separated either by space-like or by light-like distances. In both cases, the operators commute 4 so one can eraseT and T signs and get the matrix element Moreover, as we mentioned above, for the fast gluons the precise form of gauge link at infinity does not matter so we can connect points x ⊥ and y ⊥ by a straight-line gauge link proportional to gluon TMD (1.1). Note, however, that forward matrix element of this operator has an unbounded integration over x * − y * . It is convenient to introduce the notation

JHEP06(2016)164
p|O|p for the forward matrix element of the operator O stripped of this integration With this notation the unintegrated gluon TMD (1.1) can be represented as Returning to eq. (2.13), since the dependence on z i ⊥ is gone from the matrix element, we can integrate the impact factor over z 1 ⊥ and z 2 ⊥ and get the cross section as a convolution of the new impact factor I µν (x ⊥ , y ⊥ ; η) with the gluon TMD Note that the Wilson-line operators U † z and U z in eq. (2.11) cancel only when we take a sum over all intermediate states. If we are interested in, say, production of another particle (at lower rapidity), we need to consider the full double functional integral (2.8).

Rapidity factorization and evolution of TMDs in the leading order
We will study the rapidity evolution of the operator Matrix elements of this operator between unpolarized hadrons can be parametrized as [36] where m is the mass of the target hadron (typically proton). The reason we study the evolution of the operator (3.1) with non-convoluted indices i and j is that, as we shall see below, the rapidity evolution mixes functions D and H. It should be also noted that our final equation for the evolution of the operator (3.1) is applicable for polarized targets as well. In the spirit of rapidity factorization, in order to find the evolution of TMD with respect to rapidity cutoff η (see eq. (2.3)) one should integrate in the matrix element (3.3) over gluons and quarks with rapidities η > Y > η and temporarily "freeze" fields with Y < η to be integrated over later. (For a review, see refs. [44,45].) In this case, we obtain functional integral of eq. (A.8) type over fields with η > Y > η in the "external" fields with Y < η . In terms of Sudakov variables we integrate over gluons with α between σ = e η and σ = e η and, in the leading order, only the diagrams with gluon emissions are relevant -the quark diagrams will enter as loops at the next-to-leading (NLO) level.
To calculate diagrams, one needs to return to a double functional integral representation of gluon TMD (3.3): Now, in accordance with general background-field formalism we separate the gluon field into the "classical" background part with Y < η and "quantum" part with η > Y > η and integrate over quantum fields. In the leading order there are two types of diagrams: with and without gluon production, see figure 2 (we assume that there are no gluons with η > Y > η in the proton wave function).

Production part of the LO kernel
The first-order term in the expansion of the operator F m •i (y * , y ⊥ )[y * , −∞] ma y in quantum fields has the form (to save space, we omit the label cl from classical fields). As it was proved in ref. [40], to find the evolution kernel in the leading order in α s it is sufficient to consider the classical

JHEP06(2016)164
background field of the form where the absence of x • in the argument corresponds to α = 0. Using the gluon propagator (B.23) from section B.3 we obtain the result for the diagram in figure 2a in the form where O denotes the expectation value of operator O in the external field. Note that in this paper we perform calculations of diagrams in the background field (3.6) in the light-like gauge p µ 2 A µ (x) = 0 (3.8) We will make necessary comparisons with the background-Feynman gauge calculations of ref. [40] in appendix D. Let us consider now the remaining integral over "classical" fields with Y < η . It has the form where Tr(. . .) is the trace in the adjoint representation. As discussed in appendix A, the double functional integral (3.9) represents the matrix element

JHEP06(2016)164
As we mentioned above, all operators in the r.h.s. of eq. (3.10) commute since they are separated either by space-like or by light-like distance. In addition, from eq. (B.6) we see that Substituting eq. (3.11) in eq. (3.7) we get At this point we compare (3.12) to the evolution equation for p|F m Repeating steps which lead us from eq. (3.7) to eq. (3.12) we obtain We see that the production part of the evolution equation (3.12) can be obtained from eq. (3.13) by formally replacing +∞ by −∞ everywhere. Consequently, the final expression for the production part of the evolution equation for the matrix element (3.3) can be obtained from eq. (4.28) from ref. [40] by replacement ∞ ↔ −∞. 6

Virtual part of the evolution kernel
The virtual part of the kernel comes from the diagrams of the figure 2b type. The secondorder term in the expansion of the operator F m •i (y * , y ⊥ )[y * , −∞] ma y in quantum fields has 6 In the appendix D we show that the eq. (3.13), obtained in the light-like gauge, agrees with the calculations in ref. [40] performed in the background-Feynman gauge.

JHEP06(2016)164
the form (cf. eq. (3.5)) Using gluon propagator (B.18) we get Let us start with the last term in the r.h.s. of the above equation. We will prove that (y ⊥ , y * |p j 1 in our approximation. Indeed, in the "light-cone" case (when the characteristic transverse momenta of background field l ⊥ are much smaller than the momenta of the "quantum" fields p ⊥ ) it is evident since and terms ∼ O(F •j ) exceed our accuracy. (The second term in the l.h.s. of eq. (3.16) is proportional to δ(y * − y * ) and [y * , y * ] y = 1 is introduced for convenience.) In the "shock-wave" case when l ⊥ ∼ p ⊥ , if the points y and y are outside of the shock wave, the formula is trivial (y and y can only be both to the right of the shock wave since y lies inside). If y or both of them are inside the shock wave, one can again use the light-cone expansion (see the discussion in ref. [40]) and get the result (3.17). Thus, in both cases we can use eq. (3.16) so 4i It is convenient to change α ↔ −α and β ↔ −β (which is equivalent to changing y * ↔ y * ) and get Thus, we obtain the result for the last term in the r.h.s. of eq. (3.15) in the form Next we turn our attention to the first term in the r.h.s. of eq. (3.15) and start with the light-cone case l ⊥ p ⊥ : where O y * α is defined in eq. (B.9). The first term in the r.h.s. of this equation yields As to the second term in the r.h.s. of eq. (3.21), it vanishes in our approximation. Thus, the first term in r.h.s. of eq. (3.15) in the light-cone case has the form Let us now consider the shock-wave case. It is convenient to start with the representation of this term by the second line in eq. (3.15) Using eq. (B.5) for Feynman propagator one obtains If the point y * is outside the shock wave this gives

JHEP06(2016)164
If the point y is inside the shock wave we can again use the light-cone expansion and get eq. (3.25). It is easy to see that in both cases we can approximate the first term in eq. (3.15) by with our accuracy. Adding the contribution (3.20) of the second term in r.h.s. of eq. (3.15) we finally obtain the second-order virtual correction in the form where we put upper and lower cutoffs for the rapidity integrals, see the discussion following eq. (3.3). After Fourier transformation eq. (3.30) turns to Note that this equation can be obtained from eq. (4.56) from ref. [40] by reversing the sign of β B . In doing so one should go around the singularity at αβ B s = p 2 ⊥ according to Feynman rules since it corresponds to the diagram in figure 2b with cut gluon propagator.
The virtual part in the complex conjugate amplitude can be similarly obtained from eq. (4.60) from ref. [40] by replacement β B → −β B . The singular denominators should look like 1 αβ B s−p 2 ⊥ −i as appropriate for the complex conjugate amplitude.

Evolution equation for gluon TMDs
Now we are in a position to assemble all leading-order contributions to the rapidity evolution of gluon TMDs. As we discussed, in the production part of the evolution equation for the matrix element (3.3) can be obtained from eq. (4.28) from ref. [40] by replacement ∞ ↔ −∞. Adding the virtual correction to the amplitude (3.31) and its complex conjugate

JHEP06(2016)164
we obtain the evolution equation for gluon TMD operator (3.1) in the form: Here the operatorsF i (β) and F j (β) are defined as Again, this equation can be reconstructed from eq. (5.2) from ref. [40]. It should be emphasized that the reconstruction is by no means trivial: one should change ∞p 1 ↔ −∞p 1 in the production part of the amplitude and change ∞p 1 ↔ −∞p 1 and β B ↔ −β B in the virtual part. 7 7 The difference between the changes in the real and virtual part of the kernel comes from the fact that in the production part we insert the full set of out-states and use double functional integral (3.9) afterwards. The "total" replacement of lightcone time ∞ ↔ −∞ would imply also the insertion of the full set of instates. In this case the real part of the kernel will also undergo the replacement βB ↔ −βB leading to singularities 1 αβ B s−p 2 ⊥ in the production part of the amplitude. In addition, there will be diagrams with both F•i and F j • on one side of the cut which will probably cancel these singularities. In any case, the good way to avoid these complications is to insert full set of out-states but use "group law" (3.11) for O operators to set the endpoints of gauge links to −∞.

JHEP06(2016)164
The evolution equation (4.1) can be rewritten in the form where cancellation of IR and UV divergencies is evident The evolution equation (4.3) is one of the main results of this paper. It describes the rapidity evolution of the operator at any Bjorken x B ≡ β B and any transverse momenta.
When we consider the evolution of gluon TMD (1.1) given by the matrix element (3.3) of the operator we need to take into account the kinematical constraint k 2 ⊥ ≤ α(1 − β B )s in the production part of the amplitude coming from the fact that matrix element p|F i β B + This equation describes the rapidity evolution of gluon TMD (3.3) with rapidity cutoff (2.3) in the whole range of β B = x B and k ⊥ (∼ |x − y| −1 ⊥ ). In the next section we will consider some specific cases.

Small-x case: BK evolution of the Weizsacker-Williams distribution
First, let us consider the evolution of Weizsacker-Williams (WW) unintegrated gluon distribution which can be obtained from eq. (4.4) by setting β B = 0. Moreover, in the small-x regime it is assumed that the energy is much higher than anything else so the characteristic transverse σs and F i (β B ) can be replaced by U † i∂ i U (and similarly forF i ). After some algebra one obtains (cf. eq. (6.1) from ref. [40]) which agrees with ref. [46]. This equation can be rewritten as where all indices are 2-dimensional and Tr stands for the trace in the adjoint representation. Note that the expression in the square brackets is actually the BK kernel [23][24][25][26]. One should also mention that eq. (5.3) coincides with eq. (12) from ref. [47] after some algebra.
Similarly to +∞ case, the eq.

Large transverse momenta and the light-cone limit
Now let us discuss the case when β B = x B ∼ 1 and (x − y) −2 ⊥ ∼ s. At the start of the evolution (at σ ∼ 1) the cutoff in p 2 ⊥ in the integrals eq. (4.4) is ∼ (x − y) −2 ⊥ . However, as the evolution in rapidity (∼ ln σ) progresses the characteristic p 2 ⊥ become smaller due to the kinematical constraint p 2 ⊥ < σ(1 − β B )s. Due to this kinematical constraint evolution in σ is correlated with the evolution in p 2 ⊥ : if σ σ the corresponding transverse momenta of background fields p ⊥ 2 are much smaller than p 2 ⊥ in quantum loops. This means that during the evolution we are always in the light-cone case considered in section 3 and therefore the evolution equation for β B = x B ∼ 1 and (x − y) −2 ⊥ ∼ s takes the form

JHEP06(2016)164
which reduces to the system of evolution equations for gluon TMDs D(β B , |z ⊥ |, ln σ) and H(β B , |z ⊥ |, ln σ) in the case of unpolarized hadron. The evolution equation (5.4) can be rewritten as a system of evolution equations for D and H functions (z ≡ σsβ B . 8 The above equation is our final result for the rapidity evolution of gluon TMDs (1.1) in the near-light-cone case.
If we take the light-cone limit x ⊥ = y ⊥ (⇔ z ⊥ = 0) we get the (one-loop) DGLAP equation: One immediately recognizes the expression in the square brackets as gluon-gluon DGLAP kernel.
There is a subtle point in comparison of our rapidity evolution of light-ray operators to the conventional µ 2 evolution described by renorm-group equations: the self-energy diagrams are not regulated by our rapidity cutoff so the δ-function terms in our version of the DGLAP equations are absent. 9 Indeed, in our analysis we do not change the UV treatment of the theory, we just define the Wilson-line (or light-ray) operators by the requirement that gluons emitted by those operators have rapidity cutoff (2.3). The UV divergences in self-energy and other internal loop diagrams appearing in higher-order calculations are absorbed in the usual Z-factors. So, in a way, we will have two evolution equations for our operators: the trivial µ 2 evolution described by anomalous dimensions of corresponding gluon (or quark) fields and the rapidity evolution. Combined together, the two should describe the Q 2 evolution of DIS structure functions. Presumably, the argument of coupling constant in LO equation (5.6) (which is µ 2 by default) will be replaced by σβ B s in 8 Careful analysis shows that virtual correction ∼ V.p. σβ B s k 2 ⊥ (σβ B s−k 2 ⊥ ) leads to the same (. . .)+ prescription as the virtual correction ∼ σβ B s k 2 ⊥ (σβ B s+k 2 ⊥ ) for the operator F•i[y * , +∞] so the eq. (5.5) coincides with eq. (3.29) from ref. [40] 9 For Eq. (5.6) the absence of these terms is accidental, due to an extra αs in the definition (2.15).

JHEP06(2016)164
accordance with common lore that this argument is determined by characteristic transverse momenta. 10 We plan to return to this point in the future NLO analysis.

Sudakov logarithms
Finally, let us consider the evolution of D(x B , k ⊥ , η = ln σ) in the region where x B ≡ β B ∼ 1 and k 2 ⊥ ∼ (x−y) −2 ⊥ ∼ few GeV 2 . In this case the integrals over p 2 ⊥ in the production part of the kernel (4.4) are ∼ (x − y) −2 ⊥ ∼ k 2 ⊥ so that p 2 ⊥ σβ B s for the whole range of evolution 1 > σ > As to the virtual part the two last lines can be omitted. To prove this we follow the logic of ref. [40] and consider two cases: the "light-cone case" l 2 ⊥ p 2 ⊥ and the "shock-wave" situation when l 2 ⊥ ∼ p 2 ⊥ . It is easy to see that in the light-cone case the two last terms in the r.h.s. of eq. (5.8) reduce to the operators of higher collinear twist. In the shock-wave case we need to consider two sub-cases: if p 2 ⊥ σβ B s and p 2 ⊥ ∼ σβ B s. In the first (sub)case the two last terms in the r.h.s. of eq. (5.8) are again trivially negligible in comparison to the first term in the r.h.s. of that equation. In the second (sub)case (when p 2 ⊥ ∼ σβ B s) one can expand the operator

The first term in the r.h.s. of this equation is obviously zero while the second is
σβ B s in comparison to the leading first term in the r.h.s. of eq. (5.8) (the transverse momenta inside the hadron target are ∼ m N ∼ 1GeV).

JHEP06(2016)164
Thus, we obtain the following rapidity evolution equation in the Sudakov region: Similarly to ref. [40], there is a double-log region where 1 σ In that region only the second term in the r.h.s. of eq. (5.9) survives so the evolution equation reduces to which can be rewritten for the TMD (1.1) as leading to the usual Sudakov double-log result It is worth noting that the coefficient in front of ln 2 σs k 2 ⊥ is determined by the cusp anomalous dimension of two light-like Wilson lines going from point y to ∞p 1 and ∞p 2 directions (with our cutoff α < σ), see the discussion in ref. [40].

Rapidity evolution of unintegrated gluon distribution in linear approximation
It is instructive to present the evolution kernel (4.4) in the linear (two-gluon) approximation. Since in the r.h.s. of eq. (4.4) we already haveF k and F l (and each of them has at least one gluon) all factors U andŨ in the r.h.s. of eq. (4.4) can be omitted and we get (η ≡ ln σ) where we performed Fourier transformation to the momentum space. Also, the forward matrix element Eliminating this factor and rewriting in terms of R ij (see eq. (3.2)) we obtain (η ≡ ln σ) As we demonstrated in ref. [40] in the low-x limit β B → 0 the above equation reduces to the BFKL equation and the evolution of is governed by the DGLAP equation (5.6). It would be interesting to find how the linear evolution equation (6.2) is connected with other results on the combined small-and large-x resummation, see refs. [48][49][50][51][52][53].

JHEP06(2016)164 7 Conclusions
We have described the rapidity evolution of gluon TMD (1.1) with Wilson lines going to −∞ in the whole range of Bjorken x B and the whole range of transverse momentum k ⊥ . It should be emphasized that with our definition of rapidity cutoff (2.3) the leadingorder matrix elements of TMD operators are UV-finite so the rapidity evolution is the only evolution and it describes all the dynamics of gluon TMDs (1.1) in the leading-log approximation. In the next-to-leading order one should expect usual renorm-group on the top of rapidity evolution so the coupling constant α s in our equation will become running coupling, presumably dependent on some transverse momenta distances as in the NLO BK equation [54][55][56][57][58].
It should be emphasized that rapidity evolution equations for gauge links to + and − infinity are not identical: the virtual part of the kernel in eq. (4.4) is different from eq. (5.5) from our previous paper [40] (the real part is the same). However, this difference disappears in all interesting limits: DGLAP, Sudakov and small-x, so we think that it will be important only for transition between linear and non-linear evolution where one should take into account the whole eq. (4.4).
For completeness, let us present the description of various cases of linear vs nonlinear evolution repeating the discussion in ref. [40].
The rapidity evolution of gluon TMD (1.1) with rapidity cutoff (2.3) is given by (4.4) and, in general, is non-linear. Nevertheless, for some specific cases the equation x B ≡ β B k 2 ⊥ s , we get an interplay of linear and non-linear evolutions. If we change the rapidity cutoff σ (↔ rapidity) from 1 to k 2 ⊥ s first there will be Sudakov-type double logarithmic evolution (5.11) from σ = 1 to σ = k 2 ⊥ β B s , then there will be a transitional region at σ ∼ k 2 ⊥ β B s , and after that we will have the non-linear evolution (5.3) at s (the interplay of the non-linear evolution and Sudakov double logarithms in this region was studied in refs. [59][60][61] at the NLO level). Needless to say, the transition between the linear evolution (5.11) and the non-linear one (5.3) should be described by the full equation

JHEP06(2016)164
In conclusion, let us mention that an obvious outlook project is to present the "impact factor for the photon" in eq. (2.11) for the cross section as another TMD with gauge links aligned along the proton's momentum. The hope is to get k T -factorization in the form of product of the two TMDs in the whole range of Bjorken x and make the connection between k T -factorization and collinear factorization. This study is in progress.
The authors are grateful to G.A. Chirilli, J.C. Collins, Yu. Kovchegov, A. Prokudin, A.V. Radyushkin, T. Rogers, and F. Yuan for valuable discussions. This work was supported by contract DE-AC05-06OR23177 under which the Jefferson Science Associates, LLC operate the Thomas Jefferson National Accelerator Facility, and by the grant DE-FG02-97ER41028.

A Inclusive particle production as double functional integral
In this section we will prove that the amplitude of inclusive particle production is given by the double functional integral (2.2).
The cross section of the production of Φ-meson in deep inelastic scattering is given by where X denotes the sum over full set of "out" states. Using standard LSZ formula we reduce eq. (A.1) to where F 2 ≡ F m αβ F mαβ for brevity. Now, |X and |p may be considered as eigenstates of the full QCD HamiltonianĤ |X = E X |X ,Ĥ|p = E p |p so one can rewrite X|T {F 2 (y)j ν (0)}|p as where t i → −∞ is the initial time and t f → ∞ is the final time.

JHEP06(2016)164
so the cross section (A.2) takes the form At this point it is convenient to switch to the sum over all states in the "coordinate representation" where Ψ p ( A(t i ), ψ(t i )) is the proton wave function at the initial time t i . In the same way one can demonstrate that a general matrix element can be represented by a double functional integral: with the boundary conditionÃ( x, t = ∞) = A( x, t = ∞) (and similarly for quark fields) reflecting the sum over all intermediate states X.

JHEP06(2016)164 B Propagators in fast background fields
In this section we will obtain propagators for the double functional integral (3.4) in external low-α fields. As we proved in ref. [40], it is sufficient to consider the external field of the type A • (x * , x ⊥ ) (and quark fields p 1 ψ(x * , x ⊥ )) with all other components being zero. 11 Indeed, if the characteristic transverse momenta of fast fields (l ⊥ ) and slow fields (k ⊥ ) are comparable, the usual rescaling of refs. [23,24] applies so only A • (x * , x ⊥ ) of the type of shock wave survives. Conversely, if k ⊥ l ⊥ the fast fields do not necessarily shrink to a shock wave but we can apply the light-cone expansion of propagators. The parameter of the light-cone expansion is the twist of the operator and we will expand up to operators of leading collinear twist two. Such operators are built of two gluon operators ∼ F •i F •j or quark onesψ p 1 ψ and gauge links. To get coefficients in front of these operators it is sufficient to consider the external gluon field of the type A • (z * , z ⊥ ) with A i = A * = 0.

B.1 Scalar Feynman propagator
For simplicity we will first perform the calculation for "scalar propagator" (x| 1 P 2 +i |y). As we mentioned above, we assume that the only nonzero component of the external field is A • and it does not depend on z • so the operator α = i ∂ ∂z• commutes with all background fields. The propagator in the external field A • (z * , z ⊥ ) has the form The Pexp in the r.h.s. of eq. (B.1) can be transformed to Now we expand This is an expansion around the light cone z ⊥ + 2 s z * p 1 . We are keeping the first three terms of the expansion which is sufficient in both shock-wave case l ⊥ ∼ k ⊥ and "lightcone" case l ⊥ k ⊥ . In the shock-wave case it is obvious since the parameter of the so the the scalar propagator in the fast external field takes the form Note that O α (x * , y * ) trivially satisfies the group property For future use we present also two equivalent expressions with derivative operators to the right and to the left of the field operators:

JHEP06(2016)164
Here we display right or left p ⊥ in the notation for O to indicate whether we use representation (B.7) or (B.8).
To finish the proof of eq. (B.5) we need to demonstrate that it is correct in the lightcone case. We will need the general formula In the light-cone case one expands the external field either around the light cone y ⊥ + 2 s z * p 1 or x ⊥ + 2 s z * p 1 . Let us consider the first case (the second is equivalent). The Pexp in the r.h.s. of eq. (B.1) can be transformed to This is effectively expansion around the light ray y ⊥ + 2 s y * p 1 with the parameter of the expansion ∼ |l ⊥ | |p ⊥ | 1. As we mentioned, we expand up to the operators of twist two.
Similarly, one can demonstrate that the propagator in the complex conjugate amplitude has the form ⊥ αs x * = O x * α (x * , y * ; p ⊥ ) and rewriting according to eq. (B.8) this equation coincides with eq. (A.12) from ref. [40].

B.2 Scalar propagator of Wightman type
The scalar propagator from point x to the left of the cut to point y to the right of the cut reads 14) It is convenient to represent this equation as an integral of product of two amplitudes of particle emission found in ref. [40]: In the shock-wave case l ⊥ ∼ k ⊥ these formulas coincide with eqs. (B.18) and (B.20) from ref. [40]; in the light-cone case one needs to rewrite them as after which they coincide with eqs. (A.14) and (A.16) from ref. [40]. Using eq. (B.15) one easily obtains whereÕ is built of theÃ fields in the left functional integral in eq. (A.8).

B.3 Gluon propagator in the light-like gauge
The general expression for Feynman gluon propagator in the light-like gauge p µ 2 A µ = 0 in the background field (3.6) has the form Using the expression (B.5) for 1 P 2 +i we get For the complex conjugate amplitude one obtains in a similar way where we used eq. (B.13) for 1 P 2 −i . The "cut" propagator in the background field (3.6) is given by eq. (C.4) Using eq. (B.17) for scalar propagator we obtain where, as usual,Õ is built of theÃ fields in the left functional integral in eq. (A.8).

C Feynman diagrams for the gluon propagator in the light-like gauge
The formulas (B.18) and (B.20) can be easily obtained from general formula for the propagator in the light-like gauge in refs. [23,24]. However, the expression (B.22) for Wightman gluon propagator needs derivation and the easiest way is to analyze Feynman diagrams in the background field (3.6) (cf. ref. [65]). Let us consider a typical diagram shown in figure 3. The perturbative gluon propagators in the light-like p µ 2 A µ = 0 gauge has the form

JHEP06(2016)164
First, we prove that only one term in the three-gluon vertex survives. Indeed, consider a typical 3-gluon vertex It is easy to see that the two last terms do not contribute since the vertex is multiplied by d αµ (k) and d νβ (k + q) so we are left with the first term which is a vertex of emission of the gluon by scalar propagator multiplied by g µν . Second, let us consider the product of numerators of gluon propagators in figure 3 d αµ 1 (k)d µ 1 µ 2 (k + q 1 )d µ 2 µ 3 (k + q 1 + q 2 ) . . . d µnβ (k + q 1 + · · · + q n ) (C.2) It is clear that for all d µν 's, except the first and the last ones, we can replace d µν (k) by g ⊥ µν since terms ∼ p 2µ vanish. For the same reason, only two terms in the first and in the last d µν 's survive: Thus, the gluon propagator in the background field (3.6) in the light-like p µ 2 A µ = 0 gauge differs from the scalar propagator in the same background field (B.14) only by two JHEP06(2016)164

D.1 Light-cone case
Let us start with the "light-cone case" when the characteristic transverse momenta of background field l ⊥ are much smaller than the momenta of the "quantum" fields p ⊥ . As

D.2 Shock-wave case
If the characteristic transverse momenta of background field l ⊥ are of the same order of magnitude as the momenta of the "quantum" fields p ⊥ we have a "shock-wave case" when longitudinal size of background fields σ * ∼ σs l 2 ⊥ is much smaller than typical distances in quantum Feynman diagrams ∼ αs l 2 ⊥ (recall that α σ). As in ref. [40], we must consider separately two cases: y * inside and outside of the shock wave. The first case is simple: To prove eq. (D.9) we first notice that at y * > 0 (and outside of the shock wave) the eq. (D.9) vanishes since F •i (z * ) = 0. Second, if y * < 0 the integral