1 Introduction

Nowadays, it is an ordinary fact that most processes occurring at high energies do not involve just one large energy scale. Consequently, the standard collinear factorization often is not sufficient or even does not apply. Example of a wide class of such processes are those involving large measurable internal transverse momenta of partons. A consistent theoretical treatment of such processes was initiated by Diakonov, Dokshitzer and Troyan [1] and turned into the Transverse Momentum Dependent (TMD) factorization [2, 3] (for a recent review see [4]). This concept, not only resumes the large logarithms, but also defines, within the QCD theory, more general (and interesting) objects than usual collinear parton distribution functions (PDFs) – the TMD parton distribution functions. The high energy factorization (or \(k_{T}\)-factorization) [5,6,7,8] and Color Glass Condensate (CGC) effective theory [9,10,11,12] also address transverse momentum dependence of gluons in a hadron, although in somewhat different kinematic regime, namely in the so-called small-x limit, where the gluonic degrees of freedom dominate due to large logarithms of energy that have to be resummed. The collinear-factorization-based Monte Carlo event generators like Pythia [13] or Herwig [14] are also capable of simulating the semi-hard processes by constructing explicit parton branching mechanisms, based on the Dokshitzer–Gribov–Lipatov–Altarelli–Parisi (DGLAP) evolution kernel, i.e. the Sudakov form factors (although they have to implement several model-dependent mechanisms to maintain the momentum conservation, regulate singularities, etc...). The Cascade Monte Carlo event generator [15] attacks similar problems employing the evolution equations which resum also the logarithms enhanced at small x.

Although the applicability of the strict TMD factorization theorems is limited to few processes only (like Drell–Yan or semi-inclusive DIS), the basic objects appearing in the formalism – the TMD parton distributions – can be studied in the broader context. They are defined as the Fourier transforms of the hadronic matrix elements of bilocal field operators with non-light-like separation. To ensure the gauge invariance the Wilson links connecting the two space-time points must be inserted. For the gauge invariance itself the shape of the links is not relevant. In the TMD factorization however, the shape of the links is determined by the hard process accompanying the TMD parton distribution. This happens because the collinear gluons (to the incoming hadron), which couple to various components of the hard process have to be considered as a part of the nonperturbative wave function. They can be resummed into the Wilson links attached to each external leg by means of the Ward identity. Since the external legs are connected by certain color matrix, so are the pieces of Wilson links and this is how the process dependence enters (see [16, 17] for details). For simple processes like Drell–Yan pairs production, the color flow in the hard process is rather simple because of only two colored partons. Consequently, the resulting TMD parton distribution has also simple structure. On the contrary, for processes with several colored partons one gets multiple nonequivalent structures (including Wilson loops) which cannot be eliminated by a gauge choice.

Although, as mentioned, the strict all-order factorization theorems fail for more than two colored partons participating in the hard collision [18], in the nonlinear small-x regime the lowest order TMDs are of great phenomenological importance. In Ref. [19] a leading power limit of the expressions for dijet production in pA collisions within the CGC was studied. They found that the correlators of Wilson lines averaged over color sources according to the CGC theory correspond exactly to the TMD gluon distributions for \(2\rightarrow 2\) processes, provided the hadronic matrix elements are traded for the color source averages. Not only the correlators agree, but also the hard factors. Although it is not known whether this correspondence survives beyond the leading order, it opened new phenomenological opportunities to study with better theoretical control semi-hard jets in the gluon saturation domain, see [20,21,22,23,24]. In particular, in Ref. [20] a beyond-leading-power extension of the TMD factorization for forward dijets in pA collisions was proposed, such that it coincides with the leading power of CGC in the dense nucleus regime, and with the all-power high energy factorization in the dilute nucleus limit. One should understand the notion ’factorization’ here in the following sense. First, the overall kinematic conditions justify so called hybrid approach [25], i.e. where the projectile proton is treated as a dilute state so that an average parton coming from it is a large-x parton modeled from ordinary collinear parton distribution functions. The target nucleus is probed in the dense state, so it is modeled basing on the small-x dynamics (note that the operator definitions of the TMD distributions formally are valid also at small x). Second, it is a generalized factorization, i.e. the formulae involve several TMD gluon distributions for nucleus. In the formal leading-power TMD factorization, even the generalized factorization breaks, because one is unable to define the separate correlators whilst more than two colored partons are present [18]. In the small-x approach for dilute-dense collisions described above, however, there is only one correlator with transverse separation. Therefore the complications leading to the lack of possibility to separate Wilson links into TMD operators, formally do not appear here. Outside the small-x limit for dilute-dense collisions these results might also be useful: for example to access the factorization breaking effects.

In the TMD factorization formalism the TMD parton distributions have operator definitions and evolve according to the renormalization group equations [3]. In the small x regime with gluon saturation playing significant role, which is of main interest in the context of this work, the evolution equations are nonlinear and thus more complicated than the equations at moderate x [26,27,28]. In addition, the program of obtaining the renormalization group evolution equations for all possible TMD operators is nowhere near the end. Hopefully, the correspondence of the small-x TMD gluon distributions and CGC correlators [19] allows for a treatment of evolution in the strict small x limit using the Balitsky–Jalilian–Marian–Iancu–McLerran–Weigert–Leonidov–Kovner (B-JIMWLK) equations [29,30,31,32,33,34,35] following Ref. [22]. At small x, but in the linear regime, where the saturation scale is much smaller than the typical scale of the internal transverse momenta, it seems that the various TMD gluon distributions converge to one universal distribution, which may be identified with the so-called unintegrated gluon distribution [21, 22, 24]. This object is much better understood and constrained from data. There are several approaches to their evolution. First, there are extensions of the original Balitsky–Fadin–Kuraev–Lipatov (BFKL) equation (see e.g. [36] for a review) like the Catani-Ciafaloni-Fiorani-Marchesini (CCFM) equation [37,38,39,40], Kwieciński–Martin–Staśto (KMS) equation [41] or the Kimber–Martin–Ryskin (KMR) approach [42]. As the linear evolution can be solved through the explicit branching process, it allows for a natural determination of the unintegrated PDFs from Monte Carlo simulations [43,44,45]. The complete set of evolution equations in the linear regime can be also derived through the projector method [46], see [47, 48] for a recent approach. There have been many calculations attacking various processes where the usage of unintegrated parton distributions is important, see for example [49,50,51,52,53,54,55,56,57,58,59,60] where mostly forward jet observables in hadroproduction were addressed. These calculations, however, use universal unintegrated gluon distributions, the same for any color flow. While in the linear regime or in certain phase space regions, this is a good approximation, it is definitely not the case in the region where the gluon saturation may dominate [19, 61].

Motivated by the phenomenological usability of the non-universal TMD gluon distributions discussed above (and demonstrated in [21]), we will present explicit results for the operator structures for all five and six colored parton processes. Instead of working with particular Feynman diagrams and calculating the corresponding operator structure, we choose to work with color decomposition of amplitudes (see e.g. [62]). This is motivated simply by the way the amplitudes are calculated at present in practice. Such a procedure for the operator structures in the TMD gluon distributions was for the first time used in [20] for four parton processes.

The paper is organized as follows. We will start with definitions of the TMD parton distribution functions and summary of color decomposition of scattering amplitudes (Sect. 2). Next, in Sect. 3, we will introduce the color flow diagrams for the operators appearing in the definitions of the TMD distributions. Basing on these rules, we list all structures appearing in arbitrary process in Sect. 4. The explicit results for four, five and six parton processes will be given in Sect. 5. We will summarize our work in Sect. 7.

2 Preliminaries

We start off by providing necessary definitions and conventions for the TMD gluon distributions and color decomposition.

We adopt the light cone basis defined using two null four vectors \(n=\left( 1,0,0,-1\right) /\sqrt{2}\) and \({\tilde{n}}=\left( 1,0,0,1\right) /\sqrt{2}\). They define the ’plus’ and ’minus’ components of a four vector v: \(v^{+}=n\cdot v\), \(v^{-}={\tilde{n}}\cdot v\), so that the four vector has a decomposition

$$\begin{aligned} v^{\mu }=v^{+}\,{\tilde{n}}^{\mu }+v^{-}\,n^{\mu }+v_{T}^{\mu }. \end{aligned}$$
(1)

The light-cone coordinates are \(\left( v^{+},v^{-},\vec {v}_{T}\right) \), where the Euclidean transverse vector is defined (in canonical coordinates) as \(v_{T}^{\mu }=\left( 0,\vec {v}_{T},0\right) \).

We consider n-parton processes with a gluon in the initial state

$$\begin{aligned} g\left( k_{1}\right) +b_{n}\left( k_{n}\right) \rightarrow b_{2}\left( k_{2}\right) +\dots +b_{n-1}\left( k_{n-1}\right) , \end{aligned}$$
(2)

where the partons \(b_{i}\) can be quarks or gluons (restricted by the flavor number conservation of course). The initial state gluon with momentum \(k_{1}\) carries an x fraction of the parent hadron with momentum \(P^{\mu }=P^{+}{\tilde{n}}^{\mu }\):

$$\begin{aligned} k_{1}^{\mu }\simeq xP^{\mu }+k_{T}^{\mu }. \end{aligned}$$
(3)

Above, the minus component is suppressed as it is neglected in the hard part. The transverse component is also neglected within the leading twist collinear and TMD factorization. In the more general case, the gluon may be off-shell and a suitable redefinition of the hard process is required to maintain the gauge invariance (see e.g. [63,64,65,66,67]). Even then, at least formally, the principles to obtain the TMD distributions still hold, therefore we shall not distinguish these situations here.

2.1 TMD gluon distributions

In the present work we will be concerned with the gluon TMD distributions as explained in the Introduction. A generic distribution is defined as the following matrix element:

$$\begin{aligned} {\mathcal {F}}\left( x,k_{T}\right)= & {} 2\int \frac{d\xi ^{-}d^{2}\xi _{T}}{\left( 2\pi \right) ^{3}P^{+}}\,e^{ixP^{+}\xi ^{-}-i\vec {k}_{T}\cdot \vec {\xi }_{T}}\,\nonumber \\&\times \left\langle P\right| \mathrm {Tr}\left\{ {\hat{F}}^{i+}\left( 0\right) {\mathcal {U}}_{C_{1}}{\hat{F}}^{i+}\left( \xi ^{+}=0,\xi ^{-},\vec {\xi }_{T}\right) {\mathcal {U}}_{C_{2}}\right\} \left| P\right\rangle ,\nonumber \\ \end{aligned}$$
(4)

where \(\left| P\right\rangle \) is a hadron state, \({\hat{F}}^{\mu \nu }\left( x\right) =F_{a}^{\mu \nu }\left( x\right) t^{a}\) is the \(\mathrm {SU}\left( N_{c}\right) \) algebra-valued field strength tensor (we use \(\mathrm {Tr}\left( t^{a}t^{b}\right) =T_{F}\delta ^{ab}\), \(T_{F}=1/2\) convention for the generators),  \({\mathcal {U}}_{C_{1}}\), \({\mathcal {U}}_{C_{2}}\) are certain fundamental representation Wilson lines joining space-time points \(\left( \xi ^{+}=0,\xi ^{-}=0,\vec {\xi }_{T}=\vec {0}_{T}\right) \) and \(\left( \xi ^{+}=0,\xi ^{-},\vec {\xi }_{T}\right) \), multiplied by possible traces of Wilson loops. The exact shape of Wilson lines will depend on the hard process coupled to the TMD. Their calculation for multiple partons is the main goal of the present work.

The above generic definition represents a bare TMD distribution. In QCD there are divergences, in particular the rapidity divergence that have to be regulated. In the present work we will not be considering the renormalization of these operators. Recent studies of that matter in the small-x limit which mainly motivates the present work are given in [26,27,28].

A generic Wilson line joining x and y through a path C is defined as

$$\begin{aligned} {\mathcal {U}}_{C}={\mathcal {P}}\exp \left\{ -ig\int _{C}dz_{\mu }{\hat{A}}^{\mu }\left( z\right) \right\} . \end{aligned}$$
(5)

The Wilson line can be defined also in the adjoint representation, by replacing generators \(t^{a}\) by \(\left( T^{a}\right) _{bc}=-if^{abc}\). In the case where the path is a straight line segment, we will use the following notation

$$\begin{aligned} {\mathcal {U}}_{C}\equiv \left[ x,y\right] . \end{aligned}$$
(6)

2.2 Color decomposition

The calculation of the operator structure entering the TMD distributions is nicely systematized not by considering a particular diagrams, but rather by considering various color flows in the amplitude (squared) under consideration. Such systematization is achieved by using gauge invariant decomposition of amplitudes into so-called color-ordered amplitudes (called also partial, or dual amplitudes). Here we are presenting only necessary definitions and properties, see e.g. [62] for a complete review.

Let us start with pure gluonic tree-level amplitudes. For the sake of this section we assume that all the momenta are outgoing (later, it will become necessary to distinguish incoming and outgoing legs). The most standard decomposition reads

$$\begin{aligned}&{\mathcal {M}}^{a_{1}\dots a_{n}}\left( k_{1},\dots ,k_{n}\right) \nonumber \\&\quad =\sum _{\pi \in S_{n}/Z_{n}}\mathrm {Tr}\left( t^{a_{\pi \left( 1\right) }}\dots t^{a_{\pi \left( n\right) }}\right) \,{\mathcal {A}}\left( \pi \left( 1\right) ,\dots ,\pi \left( n\right) \right) ,\nonumber \\ \end{aligned}$$
(7)

where the sum runs over all noncyclic permutations \(\pi \) of an n-element set. Three important properties of the above decomposition are: i) the partial amplitudes \({\mathcal {A}}\) are gauge invariant, ii) the partial amplitudes contain only planar diagrams; consequently the full amplitude squared satisfies \(\left| {\mathcal {M}}\right| ^{2}=C\sum _{S_{n-1}}\left| {\mathcal {A}}\left( 1,\mathcal {\pi }\left( 2\right) ,\dots ,\pi \left( n\right) \right) \right| ^{2}+{\mathcal {O}}\left( 1/N_{c}^{2}\right) \), with C being a color factor, iii) the amplitudes \({\mathcal {A}}\) satisfy so-called Ward identities: \({\mathcal {A}}\left( 1,\dots ,n\right) +{\mathcal {A}}\left( 1,\dots ,n,n-1\right) +\dots +{\mathcal {A}}\left( 1,n,2,\dots \right) =0\) (and similar for other partial amplitudes). Because of the last property, sometimes more desirable is a decomposition which utilizes only \(\left( n-2\right) !\) independent partial amplitudes, instead of \(\left( n-1\right) !\) as in the fundamental-representation (7). Such decomposition uses the adjoint generators [68]:

$$\begin{aligned} {\mathcal {M}}^{a_{1}\dots a_{n}}\left( k_{1},\dots ,k_{n}\right)= & {} \sum _{\pi \in S_{n-2}}\left( T^{a_{\pi \left( 2\right) }}\dots T^{a_{\pi \left( n-1\right) }}\right) _{a_{1}a_{n}}\nonumber \\&\times {\mathcal {A}}\left( 1,\pi \left( 2\right) ,\dots ,\pi \left( n-1\right) ,n\right) ,\nonumber \\ \end{aligned}$$
(8)

with \(\left( T^{a}\right) _{bc}=-if^{abc}\). The partial amplitudes above are the same as in the fundamental-representation decomposition.

Finally, let us recall the so-called color flow decomposition [69]. It will be useful especially for processes with quarks as it treats gluons and quarks on equal footing. The basic idea is to work with the gluon fields as the elements of the \(\mathrm {SU}\left( N_{c}\right) \) algebra, i.e. matrices \({\hat{A}}_{j}^{i}\equiv A_{a}\left( t^{a}\right) _{j}^{i}\). That is, a gluon is characterized by a pair of fundamental and anti-fundamental representation indices \(i,j=\left\{ 1,\dots ,N_{c}\right\} \). In this representation, the amplitude can be decomposed as

$$\begin{aligned} {\mathcal {M}}_{j_{1}\dots j_{n}}^{i_{1}\dots i_{n}}\left( k_{1},\dots ,k_{n}\right)= & {} 2^{-n/2}\sum _{\pi \in S_{n-1}}\delta _{j_{\pi \left( 2\right) }}^{i_{1}}\delta _{j_{\pi \left( 3\right) }}^{i_{\pi \left( 2\right) }}\delta _{j_{\pi \left( 4\right) }}^{i_{\pi \left( 3\right) }}\dots \delta _{j_{1}}^{i_{\pi \left( n\right) }}\nonumber \\&\quad \times {\mathcal {A}}\left( 1,\pi \left( 2\right) ,\dots ,\pi \left( n\right) \right) , \end{aligned}$$
(9)

again with exactly the same partial amplitudes as in the other two representations.

For processes with quarks, we use the color flow decomposition as it treats the quarks and gluons uniformly, and is best for easy calculation of the TMD operator structures. The decomposition for a process with one quark–anti-quark pair,

$$\begin{aligned} g\left( k_{1}\right) q\left( k_{2}\right) g\left( k_{3}\right) \dots g\left( k_{n-1}\right) {\bar{q}}\left( k_{n}\right) \rightarrow \emptyset , \end{aligned}$$

reads:

$$\begin{aligned}&{\mathcal {M}}_{j_{1}j_{3}\dots j_{n-1}j_{n}^{{\bar{q}}}}^{i_{1}i_{2}^{q}i_{3}\dots i_{n-1}}\left( k_{1},\dots ,k_{n}\right) \nonumber \\&\quad =2^{-(n-2)/2}\sum _{\pi \in S_{n-2}}\delta _{j_{\pi \left( 1\right) }}^{i_{2}^{q}}\delta _{j_{\pi \left( 3\right) }}^{i_{\pi \left( 1\right) }}\delta _{j_{\pi \left( 4\right) }}^{i_{\pi \left( 3\right) }}\dots \delta _{j_{n}^{{\bar{q}}}}^{i_{\pi \left( n-1\right) }}\nonumber \\&\qquad \times {\mathcal {A}}\left( 2^{q},\pi \left( 1\right) ,\pi \left( 3\right) ,\dots ,\pi \left( n-1\right) ,n^{{\bar{q}}}\right) . \end{aligned}$$
(10)

Above we have put superscripts q, \({\bar{q}}\) to remind which indices belong to a quark (anti-quark). The decomposition for a process with two quark–anti-quark pairs,

$$\begin{aligned} g\left( k_{1}\right) q\left( k_{2}\right) {\bar{q}}\left( k_{3}\right) q\left( k_{4}\right) g\left( k_{5}\right) \dots g\left( k_{n-1}\right) {\bar{q}}\left( k_{n}\right) \rightarrow \emptyset , \end{aligned}$$

reads:

$$\begin{aligned}&{\mathcal {M}}_{j_{1}j_{3}^{{\bar{q}}}\dots j_{n-1}j_{n}^{{\bar{q}}}}^{i_{1}i_{2}^{q}i_{4}^{q}i_{5}\dots i_{n-1}}\left( k_{1},\dots ,k_{n}\right) \nonumber \\&\quad =2^{-(n-4)/2}\sum _{\pi \in S_{n-3}}\delta _{j_{\{3}^{{\bar{q}}}}^{i_{2}^{q}}\delta _{j_{\pi \left( 1\right) }}^{i_{4\}}^{q}}\delta _{j_{\pi \left( 5\right) }}^{i_{\pi \left( 1\right) }}\delta _{j_{\pi \left( 6\right) }}^{i_{\pi \left( 5\right) }}\dots \delta _{j_{n}^{{\bar{q}}}}^{i_{\pi \left( n-1\right) }}\nonumber \\&\qquad \times {\mathcal {A}}\left( 2^{q},\pi \left( 3^{{\bar{q}}},4^{q}\right) ,\pi \left( 1\right) ,\pi \left( 5\right) ,\dots ,\pi \left( n-1\right) ,n^{{\bar{q}}}\right) \nonumber \\&\qquad -\frac{1}{N_{c}}\sum _{r\in \left\{ 1,5,\dots ,n-1\right\} }\sum _{\pi \in S_{n-4}}\left( \delta _{j_{\pi \left( 1\right) }}^{i_{2}^{q}}\delta _{j_{\pi \left( 5\right) }}^{i_{\pi \left( 1\right) }}\dots \delta _{j_{n}^{{\bar{q}}}}^{i_{\pi \left( r\right) }}\right) \nonumber \\&\qquad \times \left( \delta _{j_{\pi \left( r+1\right) }}^{i_{4}^{q}}\delta _{j_{\pi \left( r+2\right) }}^{i_{\pi \left( r+1\right) }}\dots \delta _{j_{3}^{{\bar{q}}}}^{i_{\pi \left( n-1\right) }}\right) \nonumber \\&\qquad \times {\mathcal {A}}\left( 2^{q},\pi \left( 1\right) ,\dots ,\pi \left( r_{i}\right) ,n^{{\bar{q}}},4^{q},\pi \left( r_{i+1}\right) ,\dots ,\pi \left( n-1\right) ,3^{{\bar{q}}}\right) .\nonumber \\ \end{aligned}$$
(11)

In the decomposition above, the first sum runs over all permutations of the \(n-4\) gluons and a quark–anti-quark pair (the curly brackets in deltas denote that the enclosed indices should be permuted together, according to the permutation \(\pi \)), while the second sum runs over various partitions of the two quark–anti-quark pairs with gluon insertions. The second sum is genuinely suppressed by \(1/N_{c}\) in case of distinct quark–anti-quark pairs; for identical pairs subleading terms will contribute to both sums in the partial amplitudes. In the present work, we shall explicitly consider processes with up to 6 partons, thus we do not give decomposition for more quark–anti-quark pairs.

In a sense, there is a price for the simplicity of the color flow decomposition. Namely, to each final state gluon we have to apply the projector

$$\begin{aligned} {\mathcal {P}}_{jj'}^{ii'}=\delta ^{ii'}\delta _{jj'}-\frac{1}{N_{c}}\delta _{j}^{i}\delta _{j'}^{i'}, \end{aligned}$$
(12)

which removes the redundant degrees of freedom from the sum over colors. For pure gluon amplitude they are actually not needed, but must be applied to the quark amplitudes.

3 Color flow Feynman rules for TMD operators

The color flow Feynman rules (see e.g. [69]) are useful for calculating color factors. It turns out that they are also very useful in the context of calculation of the structure of the TMD operators in (4), especially, when quarks are involved. We shall supplement the standard color flow rules for color-ordered diagrams (see Table 1) with a set of additional rules which are simple color flow representations of the rules derived in [16] for calculation of a TMD operator structure in an arbitrary process.

Table 1 Standard color flow Feynman rules for partial amplitudes. All momenta are outgoing. In the middle column we show the color part only

The original procedure effectively leads to the following recipe. For each final state we assign the gauge link \({\mathcal {U}}^{\left[ +\right] }\), which joins the points 0 and \(\xi \) (see Sect. 2.1) through the point in \(+\infty \):

$$\begin{aligned} {\mathcal {U}}^{\left[ +\right] }= & {} \left[ \left( 0^{+},0^{-},\vec {0}_{T}\right) ,\left( 0^{+},\infty ^{-},\vec {0}_{T}\right) \right] \nonumber \\&\times \left[ \left( 0^{+},\infty ^{-},\vec {0}_{T}\right) ,\left( 0^{+},\infty ^{-},\vec {\xi }_{T}\right) \right] \nonumber \\&\times \left[ \left( 0^{+},\infty ^{-},\vec {\xi }_{T}\right) ,\left( 0^{+},\xi ^{-},\vec {\xi }_{T}\right) \right] . \end{aligned}$$
(13)

In case of gluons the gauge link is to be defined in adjoint representation. The Wilson link replaces the deltas for color summation when the amplitude is squared: \(\delta _{i'i}\rightarrow \left( {\mathcal {U}}^{\left[ +\right] }\right) _{i'i}\) for quarks, \(\delta ^{jj'}\rightarrow \left( {\mathcal {U}}^{\left[ +\right] \dagger }\right) ^{jj'}\) for anti-quarks and \(\delta _{a'a}\rightarrow \left( {\mathcal {U}}^{\left[ +\right] }\right) _{a'a}\) for gluons (here and in what follows \(i,j,k,\dots \) are fundamental color indices, while \(a,b,c,\dots \) are adjoint). For the initial state (not connected to the TMD gluon distribution), the resummation of the initial state interactions leads to the Wilson line extending to \(-\infty \):

$$\begin{aligned} {\mathcal {U}}^{\left[ -\right] }= & {} \left[ \left( 0^{+},0^{-},\vec {0}_{T}\right) ,\left( 0^{+},-\infty ^{-},\vec {0}_{T}\right) \right] \nonumber \\&\times \left[ \left( 0^{+},-\infty ^{-},\vec {0}_{T}\right) ,\left( 0^{+},-\infty ^{-},\vec {\xi }_{T}\right) \right] \nonumber \\&\times \left[ \left( 0^{+},-\infty ^{-},\vec {\xi }_{T}\right) ,\left( 0^{+},\xi ^{-},\vec {\xi }_{T}\right) \right] . \end{aligned}$$
(14)

Similar to final states, one needs to replace the color deltas for initial states by the matrix elements of \(\,{\mathcal {U}}^{\left[ -\right] }\). The remaining initial state (connected to the TMD) is attached to \(F_{a}^{i+}\left( \xi \right) \) in the amplitude and to \(F_{a'}^{i+}\left( 0\right) \) in the conjugate amplitude. The rest of the procedure is similar to calculating color factors: one extracts the color structure of the pertinent amplitude and makes all the contractions (here with Wilson lines and field strength tensors instead of deltas). In the end one needs to divide-out the color factor for a process without gauge links.

Table 2 Color flow Feynman rules for the gauge links. The diagrams correspond to the cut lines, as denoted by the vertical dotted line. The routing in the color loops is clock-wise

Passing to the color flow representation is straightforward. Nothing really is to be done for quarks and anti-quarks. For gluons, we first need to make a connection of the adjoint Wilson line with the trace of fundamental-representation instances of the same Wilson line, and next project it onto the fundamental color quantum numbers with the help of the Fierz identity. All rules with graphical representation are collected in Table 2. The procedure of calculating the TMD operator structures is now reduced to considering all possible color flows and applying the rules. Although, in principle, we could consider all standard Feynman diagrams, draw them in the color flow representation and calculate TMD operator structures, fortunately, we do not need to do this. Instead we can just use the color flow decomposition described in Sect. 2.2. This will also ensure, that we work with gauge invariant sets from the start.

When constructing the TMD operators, the initial and final states are treated differently, i.e. they are assigned different gauge links. Therefore, we have to adjust the color flow decomposition (9)–(11) to take into account the fact, that two legs are incoming (recall, that these decomposition are within the standard convention of all outgoing partons). This is fixed by making the replacement \(i_{1}\longleftrightarrow j_{1}\), \(i_{n}\longleftrightarrow j_{n}\), as in our convention always the first and the last partons are incoming.

Below, we present some examples to better illustrate the procedure.

3.1 Examples

Let us first illustrate the usage of color flow Feynman rules to calculate the structure of the TMD operator for the following diagram:

(15)

This diagram contributes to the process \(g\left( k_{1}\right) q\left( k_{4}\right) \rightarrow q\left( k_{2}\right) g\left( k_{3}\right) \) and represents the diagram squared and summed over final/initial colors (except insertions of the field operators). The arrows indicate whether the line is incoming/outgoing. Let us stress, that considering particular diagrams is not the way we will ultimately proceed; instead we will consider various color flows as defined in Eqs. (9)–(11). The structure of the TMD operator for this diagram was calculated in [16]. In the color flow representation we have to consider two diagrams:

(16)

The diagram with dashed line represents an exchange of the \(U\left( 1\right) \) gluon (a colorless gluon). To calculate the diagrams we simply look for the closed quark loops and make the trace of the objects appearing in the loop. The direction of the trace is clockwise. The dashed lines carry no color, thus they do not make any traces (they also always accompany \(1/N_{c}\) factors). Note, we calculate only color part (with possible \(\mathrm {SU}\left( N_{c}\right) \) matrix insertions) – we are not concerned with any kinematic factors. For the first diagram, we have

$$\begin{aligned} \mathrm {Tr}\left\{ F\left( \xi \right) {\mathcal {U}}^{\left[ +\right] \dagger }F\left( 0\right) {\mathcal {U}}^{\left[ +\right] }\right\} \mathrm {Tr}\left\{ {\mathcal {U}}^{\left[ \square \right] }\right\} , \end{aligned}$$
(17)

where the first trace corresponds to the bottom loop, the second to the top loop. Above, we defined the Wilson loop [16]

$$\begin{aligned} {\mathcal {U}}^{\left[ \square \right] }={\mathcal {U}}^{\left[ -\right] \dagger }{\mathcal {U}}^{\left[ +\right] }. \end{aligned}$$
(18)

We also use shorthand notation \(F\left( \xi \right) \equiv {\hat{F}}^{i+}\Big (\xi ^{+}=0,\xi ^{-},\vec {\xi }_{T}\Big )\). The second diagram reads

$$\begin{aligned} -\frac{1}{N_{c}}\,\mathrm {Tr}\left( F\left( \xi \right) {\mathcal {U}}^{\left[ -\right] \dagger }F\left( 0\right) {\mathcal {U}}^{\left[ +\right] }\right) . \end{aligned}$$
(19)

To get the final result, the sum of the two contributions must be divided by the sum of the color factors (without the Wilson lines), with open indices where the field operators are attached:

(20)

The part multiplying the open indices reads

$$\begin{aligned} N_{c}-\frac{1}{N_{c}}=\frac{N_{c}^{2}-1}{N_{c}}. \end{aligned}$$
(21)

Thus the TMD operator reads

$$\begin{aligned}&\mathrm {Tr}\left\{ F\left( \xi \right) \left[ \frac{N_{c}^{2}}{N_{c}^{2}-1}\,\frac{\mathrm {Tr}{\mathcal {U}}^{\left[ \square \right] }}{N_{c}}\,{\mathcal {U}}^{\left[ +\right] \dagger }\right. \right. \nonumber \\&\qquad \left. \phantom {\left\{ F\left( \xi \right) \left[ \frac{N_{c}^{2}}{N_{c}^{2}-1}\,\frac{\mathrm {Tr}{\mathcal {U}}^{\left[ \square \right] }}{N_{c}}\,{\mathcal {U}}^{\left[ +\right] \dagger }\right. \right. }\left. \phantom {\left\{ F\left( \xi \right) \left[ \frac{N_{c}^{2}}{N_{c}^{2}-1}\,\frac{\mathrm {Tr}{\mathcal {U}}^{\left[ \square \right] }}{N_{c}}\,{\mathcal {U}}^{\left[ +\right] \dagger }\right. \right. }\quad -\frac{1}{N_{c}^{2}-1}{\mathcal {U}}^{\left[ -\right] \dagger }\right] F\left( 0\right) {\mathcal {U}}^{\left[ +\right] }\right\} , \end{aligned}$$
(22)

which exactly agrees with the result quoted in [16].

As an illustration of a more complicated structure, let us consider an example contribution to the process \(gg\rightarrow q{\bar{q}}gg\):

(23)

Applying the color flow rules gives immediately the operator structure for the leading color flow displayed on the r.h.s.:

$$\begin{aligned} N_{c}\mathrm {Tr}\left\{ F\left( \xi \right) {\mathcal {U}}^{\left[ +\right] \dagger }F\left( 0\right) {\mathcal {U}}^{\left[ +\right] }\right\} \,\mathrm {Tr}{\mathcal {U}}^{\left[ \square \right] }\,\mathrm {Tr}{\mathcal {U}}^{\left[ \square \right] \dagger }. \end{aligned}$$
(24)

Above, the \(N_{c}\) factor comes from the second loop from the bottom, \(\mathrm {Tr}\left\{ {\mathcal {U}}^{\left[ +\right] \dagger }{\mathcal {U}}^{\left[ +\right] }\right\} =\mathrm {Tr}{\mathbf {1}}=N_{c}\).

To close this section, let us stress, that the problem of proliferation of color flow diagrams compared to ordinary diagrams, will not concern us at all. As mentioned, we shall use the color flow decomposition, which sets the color flow without need to consider particular diagrams.

4 The operator basis for arbitrary TMD gluon distribution

Using the color flow Feynman rules from the previous section we can easily determine all possible ’basis’ operators, from which a TMD gluon distribution for arbitrary process can be constructed. Alternatively, one can think about ’basis’ TMD gluon distributions.

Plenty of different operators already appear for processes with four colored partons considered in [16]. In order to find all of them, we use the following facts. First, there are at most two \({\mathcal {U}}^{\left[ -\right] }\) Wilson lines. This is the case for initial state gluons where \({\mathcal {U}}^{\left[ -\right] }\) and \({\mathcal {U}}^{\left[ -\right] \dagger }\) appear. Thus, we can built at most two Wilson loops (18), when they are looped with \({\mathcal {U}}^{\left[ +\right] }\) or \({\mathcal {U}}^{\left[ +\right] \dagger }\) (see the last example in Sect. 3). Second, any color flow loop will contribute trace of at most first power of \({\mathcal {U}}^{\left[ \pm \right] }\), \({\mathcal {U}}^{\left[ \pm \right] \dagger }\) (and \(F\left( \xi \right) \), \(F\left( 0\right) \), or both), in addition to mentioned Wilson loops (at most \({\mathcal {U}}^{\left[ \square \right] }\) and \({\mathcal {U}}^{\left[ \square \right] \dagger }\)). This is because for a color flow loop with many Wilson lines (contributed by many final states), most of the Wilson lines will collapse to unity, \({\mathcal {U}}^{\left[ +\right] \dagger }{\mathcal {U}}^{\left[ +\right] }={\mathbf {1}}\), leaving only at most single instances of \({\mathcal {U}}^{\left[ \pm \right] }\), \({\mathcal {U}}^{\left[ \pm \right] \dagger }\), \({\mathcal {U}}^{\left[ \square \right] }\), \({\mathcal {U}}^{\left[ \square \right] \dagger }\).

Basing on the above, below we list all ’basis’ TMD gluon distributions, from which an arbitrary TMD is given as a linear combination. We assume here, that the correlators are real valued functions.

$$\begin{aligned} {\mathcal {F}}_{qg}^{(1)}\left( x,k_{T}\right)= & {} 2\int \frac{d\xi ^{-}d^{2}\xi _{T}}{\left( 2\pi \right) ^{3}P^{+}}\,e^{ixP^{+}\xi ^{-}-i\vec {k}_{T}\cdot \vec {\xi }_{T}}\, \nonumber \\&\times \left<\mathrm {Tr}\left[ {\hat{F}}^{i+}\left( \xi \right) {\mathcal {U}}^{[-]\dagger }{\hat{F}}^{i+}\left( 0\right) {\mathcal {U}}^{[+]}\right] \right>\,\nonumber \\&=2\int \frac{d\xi ^{-}d^{2}\xi _{T}}{\left( 2\pi \right) ^{3}P^{+}}\,e^{ixP^{+}\xi ^{-}-i\vec {k}_{T}\cdot \vec {\xi }_{T}}\,\nonumber \\&\times \left<\mathrm {Tr}\left[ {\hat{F}}^{i+}\left( \xi \right) {\mathcal {U}}^{[+]\dagger }{\hat{F}}^{i+}\left( 0\right) {\mathcal {U}}^{[-]}\right] \right>, \end{aligned}$$
(25)
$$\begin{aligned} {\mathcal {F}}_{qg}^{(2)}\left( x,k_{T}\right)= & {} 2\int \frac{d\xi ^{-}d^{2}\xi _{T}}{\left( 2\pi \right) ^{3}P^{+}}\,e^{ixP^{+}\xi ^{-}-i\vec {k}_{T}\cdot \vec {\xi }_{T}}\,\nonumber \\&\times \left<\frac{\mathrm {Tr}\left[ {\mathcal {U}}^{[\square ]}\right] }{N_{c}}\mathrm {Tr}\left[ {\hat{F}}^{i+}\left( \xi \right) {\mathcal {U}}^{[+]\dagger }{\hat{F}}^{i+}\left( 0\right) {\mathcal {U}}^{[+]}\right] \right>\,\nonumber \\= & {} 2\int \frac{d\xi ^{-}d^{2}\xi _{T}}{\left( 2\pi \right) ^{3}P^{+}}\,e^{ixP^{+}\xi ^{-}-i\vec {k}_{T}\cdot \vec {\xi }_{T}}\,\nonumber \\&\times \left<\frac{\mathrm {Tr}\left[ {\mathcal {U}}^{[\square ]\dagger }\right] }{N_{c}}\mathrm {Tr}\left[ {\hat{F}}^{i+}\left( \xi \right) {\mathcal {U}}^{[+]\dagger }{\hat{F}}^{i+}\left( 0\right) {\mathcal {U}}^{[+]}\right] \right>, \end{aligned}$$
(26)
$$\begin{aligned} {\mathcal {F}}_{qg}^{(3)}\left( x,k_{T}\right)= & {} 2\int \frac{d\xi ^{-}d^{2}\xi _{T}}{\left( 2\pi \right) ^{3}P^{+}}\,e^{ixP^{+}\xi ^{-}-i\vec {k}_{T}\cdot \vec {\xi }_{T}}\, \nonumber \\&\times \left<\mathrm {Tr}\left[ {\hat{F}}^{i+}\left( \xi \right) {\mathcal {U}}^{[+]\dagger }{\hat{F}}^{i+}\left( 0\right) {\mathcal {U}}^{[\square ]}{\mathcal {U}}^{[+]}\right] \right>\,\nonumber \\= & {} 2\int \frac{d\xi ^{-}d^{2}\xi _{T}}{\left( 2\pi \right) ^{3}P^{+}}\,e^{ixP^{+}\xi ^{-}-i\vec {k}_{T}\cdot \vec {\xi }_{T}}\,\nonumber \\&\times \left<\mathrm {Tr}\left[ {\hat{F}}^{i+}\left( \xi \right) {\mathcal {U}}^{[\square ]\dagger }{\mathcal {U}}^{[+]\dagger }{\hat{F}}^{i+}\left( 0\right) {\mathcal {U}}^{[+]}\right] \right>, \end{aligned}$$
(27)
$$\begin{aligned} {\mathcal {F}}_{gg}^{(1)}\left( x,k_{T}\right)= & {} 2\int \frac{d\xi ^{-}d^{2}\xi _{T}}{\left( 2\pi \right) ^{3}P^{+}}\,e^{ixP^{+}\xi ^{-}-i\vec {k}_{T}\cdot \vec {\xi }_{T}}\,\nonumber \\&\times \left<\frac{\mathrm {Tr}\left[ {\mathcal {U}}^{[\square ]\dagger }\right] }{N_{c}}\mathrm {Tr}\left[ {\hat{F}}^{i+}\left( \xi \right) {\mathcal {U}}^{[-]\dagger }{\hat{F}}^{i+}\left( 0\right) {\mathcal {U}}^{[+]}\right] \right>\,\nonumber \\= & {} 2\int \frac{d\xi ^{-}d^{2}\xi _{T}}{\left( 2\pi \right) ^{3}P^{+}}\,e^{ixP^{+}\xi ^{-}-i\vec {k}_{T}\cdot \vec {\xi }_{T}}\,\nonumber \\&\left<\frac{\mathrm {Tr}\left[ {\mathcal {U}}^{[\square ]}\right] }{N_{c}}\mathrm {Tr}\left[ {\hat{F}}^{i+}\left( \xi \right) {\mathcal {U}}^{[+]\dagger }{\hat{F}}^{i+}\left( 0\right) {\mathcal {U}}^{[-]}\right] \right>, \end{aligned}$$
(28)
$$\begin{aligned} {\mathcal {F}}_{gg}^{(2)}\left( x,k_{T}\right)= & {} 2\int \frac{d\xi ^{-}d^{2}\xi _{T}}{\left( 2\pi \right) ^{3}P^{+}}\,e^{ixP^{+}\xi ^{-}-i\vec {k}_{T}\cdot \vec {\xi }_{T}}\,\frac{1}{N_{c}}\nonumber \\&\times \left<\mathrm {Tr}\left[ {\hat{F}}^{i+}\left( \xi \right) {\mathcal {U}}^{[\square ]\dagger }\right] \mathrm {Tr}\left[ {\hat{F}}^{i+}\left( 0\right) {\mathcal {U}}^{[\square ]}\right] \right>\,\nonumber \\= & {} 2\int \frac{d\xi ^{-}d^{2}\xi _{T}}{\left( 2\pi \right) ^{3}P^{+}}\,e^{ixP^{+}\xi ^{-}-i\vec {k}_{T}\cdot \vec {\xi }_{T}}\,\frac{1}{N_{c}}\nonumber \\&\times \left<\mathrm {Tr}\left[ {\hat{F}}^{i+}\left( \xi \right) {\mathcal {U}}^{[\square ]}\right] \mathrm {Tr}\left[ {\hat{F}}^{i+}\left( 0\right) {\mathcal {U}}^{[\square ]\dagger }\right] \right>, \end{aligned}$$
(29)
$$\begin{aligned} {\mathcal {F}}_{gg}^{(3)}\left( x,k_{T}\right)= & {} 2\int \frac{d\xi ^{-}d^{2}\xi _{T}}{\left( 2\pi \right) ^{3}P^{+}}\,e^{ixP^{+}\xi ^{-}-i\vec {k}_{T}\cdot \vec {\xi }_{T}}\, \nonumber \\&\times \left<\mathrm {Tr}\left[ {\hat{F}}^{i+}\left( \xi \right) {\mathcal {U}}^{[+]\dagger }{\hat{F}}^{i+}\left( 0\right) {\mathcal {U}}^{[+]}\right] \right>, \end{aligned}$$
(30)
$$\begin{aligned} {\mathcal {F}}_{gg}^{(4)}\left( x,k_{T}\right)= & {} 2\int \frac{d\xi ^{-}d^{2}\xi _{T}}{\left( 2\pi \right) ^{3}P^{+}}\,e^{ixP^{+}\xi ^{-}-i\vec {k}_{T}\cdot \vec {\xi }_{T}}\,\nonumber \\&\times \left<\mathrm {Tr}\left[ {\hat{F}}^{i+}\left( \xi \right) {\mathcal {U}}^{[-]\dagger }{\hat{F}}^{i+}\left( 0\right) {\mathcal {U}}^{[-]}\right] \right>, \end{aligned}$$
(31)
$$\begin{aligned} {\mathcal {F}}_{gg}^{(5)}\left( x,k_{T}\right)= & {} 2\int \frac{d\xi ^{-}d^{2}\xi _{T}}{\left( 2\pi \right) ^{3}P^{+}}\,e^{ixP^{+}\xi ^{-}-i\vec {k}_{T}\cdot \vec {\xi }_{T}}\,\nonumber \\&\times \left<\mathrm {Tr}\left[ {\hat{F}}^{i+}\left( \xi \right) {\mathcal {U}}^{[\square ]\dagger }{\mathcal {U}}^{[+]\dagger }{\hat{F}}^{i+}\left( 0\right) {\mathcal {U}}^{[\square ]}{\mathcal {U}}^{[+]}\right] \right>, \end{aligned}$$
(32)
$$\begin{aligned} {\mathcal {F}}_{gg}^{(6)}\left( x,k_{T}\right)= & {} 2\int \frac{d\xi ^{-}d^{2}\xi _{T}}{\left( 2\pi \right) ^{3}P^{+}}\,e^{ixP^{+}\xi ^{-}-i\vec {k}_{T}\cdot \vec {\xi }_{T}}\,\nonumber \\&\times \left<\frac{\mathrm {Tr}\left[ {\mathcal {U}}^{[\square ]}\right] }{N_{c}}\frac{\mathrm {Tr}\left[ {\mathcal {U}}^{[\square ]\dagger }\right] }{N_{c}}\mathrm {Tr}\left[ {\hat{F}}^{i+}\left( \xi \right) {\mathcal {U}}^{[+]\dagger }{\hat{F}}^{i+}\left( 0\right) {\mathcal {U}}^{[+]}\right] \right>, \end{aligned}$$
(33)
$$\begin{aligned} {\mathcal {F}}_{gg}^{(7)}\left( x,k_{T}\right)= & {} 2\int \frac{d\xi ^{-}d^{2}\xi _{T}}{\left( 2\pi \right) ^{3}P^{+}}\,e^{ixP^{+}\xi ^{-}-i\vec {k}_{T}\cdot \vec {\xi }_{T}}\,\nonumber \\&\times \left<\frac{\mathrm {Tr}\left[ {\mathcal {U}}^{[\square ]}\right] }{N_{c}}\mathrm {Tr}\left[ {\hat{F}}^{i+}\left( \xi \right) {\mathcal {U}}^{[\square ]\dagger }{\mathcal {U}}^{[+]\dagger }{\hat{F}}^{i+}\left( 0\right) {\mathcal {U}}^{[+]}\right] \right>\,\nonumber \\= & {} 2\int \frac{d\xi ^{-}d^{2}\xi _{T}}{\left( 2\pi \right) ^{3}P^{+}}\,e^{ixP^{+}\xi ^{-}-i\vec {k}_{T}\cdot \vec {\xi }_{T}}\,\nonumber \\&\times \left<\frac{\mathrm {Tr}\left[ {\mathcal {U}}^{[\square ]\dagger }\right] }{N_{c}}\mathrm {Tr}\left[ {\hat{F}}^{i+}\left( \xi \right) {\mathcal {U}}^{[+]\dagger }{\hat{F}}^{i+}\left( 0\right) {\mathcal {U}}^{[\square ]}{\mathcal {U}}^{[+]}\right] \right>. \end{aligned}$$
(34)

In the definitions above, the average should be understand as the hadronic matrix elements, cf. Eq. (4). Two new structures appear in addition to those known in the literature: \({\mathcal {F}}_{qg}^{\left( 3\right) }\) and \({\mathcal {F}}_{gg}^{\left( 7\right) }\).

Here the subscripts refer to a partonic process to which a given TMD distribution belongs – whether this is a pure gluonic process or a process with quarks.Footnote 1 This notation was first introduced in [19] in the context of the small-x limit and we stick to that notation in the present work.

The above set of basic TMD gluon distribution constitutes the basis for any TMD gluon distribution to be convoluted with a hard process, assuming there are no other TMD operators involved. As discussed in the Introduction, this assumption is meant to be used at small-x where it can be justified from the CGC effective theory. It is important to note, that it is the complete basis within the rules of [17] – it does not represent a basis for a gluon correlator with arbitrary gauge link structure. At least formally, the basis structures are independent; however, in the large \(k_T\) limit they start to be degenerate (or vanish), as discussed in the Introduction.

5 Operator structures for 5 and 6 colored partons

In this section we present the main result of the present work, i.e. the form of the TMD gluon distributions for the processes with 5 and 6 colored partons, together with their large \(N_{c}\) limit, which might be useful phenomenologically in short run. We will start with a derivation of 4 parton TMD operators, to demonstrate the procedure utilizing the color decomposition and, more importantly, to introduce the general notation we shall use for more complicated processes (the operator structures for 4 parton processes were first obtained in [16], and in [20] using the color decomposition).

5.1 Outline of the method using 4 parton example

As the color decomposition is most straightforward for pure gluonic amplitude, let us start with the process

$$\begin{aligned} g\left( k_{1}\right) g\left( k_{4}\right) \rightarrow g\left( k_{2}\right) g\left( k_{3}\right) . \end{aligned}$$
(35)

For gluons, three color decompositions can be used: the fundamental (7), the color flow (9), and the adjoint (8). First two involve 6 partial amplitudes, while the last one only two. As mentioned in Sect. 2.2 the 6 partial amplitudes are not independent, but their squares give the leading contribution in the large \(N_{c}\) limit – a property which we will use in Sect. 6. Here, we are interested in the full answer, thus we use the adjoint color decomposition (for processes with quarks we will use exclusively color flow decomposition). It reads

$$\begin{aligned}&{\mathcal {M}}^{a_{1}a_{2}a_{3}a_{4}}\left( k_{1},k_{2},k_{3},k_{4}\right) \nonumber \\&\quad =\left( T^{a_{2}}T^{a_{3}}\right) _{a_{1}a_{4}}\,{\mathcal {A}}(1,2,3,4)\nonumber \\&\qquad +\left( T^{a_{3}}T^{a_{2}}\right) _{a_{1}a_{4}}\,{\mathcal {A}}(1,3,2,4). \end{aligned}$$
(36)

The square of the amplitude, summed over colors, can in general be written in a matrix form

$$\begin{aligned} \left| {\mathcal {M}}\right| ^{2}=\vec {{\mathcal {A}}}^{\,\dagger }\,{\mathbf {C}}\,\vec {{\mathcal {A}}}, \end{aligned}$$
(37)

where \({\mathbf {C}}\) is the color matrix and \(\vec {{\mathcal {A}}}\) is a column vector constructed from the partial amplitudes. For the present simple case

$$\begin{aligned} {\mathbf {C}}=N_{c}^{2}N_{A}\left( \begin{array}{cc} 1 &{} \frac{1}{2}\\ \frac{1}{2} &{} 1 \end{array}\right) , \end{aligned}$$
(38)

and \(\vec {{\mathcal {A}}}\) is given in Table 3.

Table 3 Definitions of the vector of partial amplitudes \(\vec {{\mathcal {A}}}\) for all four-parton processes. The subscripts in the sub-process indication correspond to the momenta enumeration

In order to calculate the TMD operator structure, we need to insert the appropriate gauge links instead of deltas summing over colors, as reviewed in Sect. 3

$$\begin{aligned}&{\mathcal {M}}^{a_{1}a_{2}a_{3}a_{4}}{\mathcal {M}}_{a'_{1}a'_{2}a'_{3}a'_{4}}^{\dagger }\left( {\mathcal {U}}^{[+]}\right) ^{a'_{2}a_{2}}\left( {\mathcal {U}}^{[+]}\right) ^{a'_{3}a_{3}}\nonumber \\&\quad \left( {\mathcal {U}}^{[-]\dagger }\right) ^{a_{4}a'_{4}}F_{a_{1}}^{i+}\left( \xi \right) F_{a'_{1}}^{i+}\left( 0\right) . \end{aligned}$$
(39)

The structure of TMD operators is most conveniently (and inevitably) expressed in the fundamental representation. Thus the Wilson lines are transformed to the fundamental representation using

$$\begin{aligned} \left( {\mathcal {U}}^{[\pm ]}\right) _{ab}=\frac{1}{T_{F}}\mathrm {Tr}\left[ t^{a}{\mathcal {U}}^{[\pm ]}t^{b}{\mathcal {U}}^{[\pm ]\dagger }\right] . \end{aligned}$$
(40)

Next, the decomposition (36) is used to represent the above expression in the following general form

$$\begin{aligned} \vec {{\mathcal {A}}}^{\,\dagger }\,{\varvec{F}}\,\vec {{\mathcal {A}}}, \end{aligned}$$
(41)

where \({\varvec{F}}\) is the matrix of the TMD operators containing implicitly the color factors of the hard process. In most cases, it is reasonable to keep these color factors together with the hard matrix elements. Thus, to avoid double counting, we divide the elements of \({\varvec{F}}\) by the corresponding color factors of the square of the amplitude, but without the summation of indices where the field operators are attached (this corresponds to the elements of the matrix \({\mathbf {C}}\) (38) divided by \(N_{A}\)). This leads to the following definition of the TMD distribution matrix

$$\begin{aligned} {{\varvec{\Phi }}}=2\int \frac{d\xi ^{-}d^{2}\xi _{T}}{\left( 2\pi \right) ^{3}P^{+}}\,e^{ixP^{+}\xi ^{-}-i\vec {k}_{T}\cdot \vec {\xi }_{T}}\,\left\langle P \left| {\varvec{F}}\oslash \left( \frac{1}{N_{A}}{\mathbf {C}}\right) \right| P\right\rangle ,\nonumber \\ \end{aligned}$$
(42)

where the symbol \(\oslash \) represents the Hadamard division, i.e. the element-wise division: \(\left( {\varvec{A}}\oslash {\varvec{B}}\right) _{ij}={\varvec{A}}_{ij}/{\varvec{B}}_{ij}\). It may happen, for certain multiparticle processes, that some elements of the color matrix \({\mathbf {C}}\) vanish, but the corresponding elements of \({\varvec{F}}\) are non-zero. In that case, we need to modify the above prescription. We shall come back to this point when discussing processes where this happens. An additional motivation to divide out the color factors from the TMD operators is that one could in principle use the results with matrix elements not represented in the color-ordered form.

With the above definitions, the cross section for a collinear parton a to scatter off a gluon with some internal transverse momentum and producing certain number of colored partons, can be generically written as

$$\begin{aligned} d\sigma _{ag\rightarrow X}=\int \,\vec {{\mathcal {A}}}^{\,\dagger }\left( {\varvec{C}}\circ {{{\varvec{\Phi }}}}_{ag\rightarrow X}\right) \vec {{\mathcal {A}}}\,\,d\varGamma , \end{aligned}$$
(43)

where \(d\varGamma \) represents all pre-factors, phase space, and convolution in x and \(k_{T}\). The symbol \(\circ \) is the Hadamard (element-wise) multiplication, \(\left( {\varvec{A}}\circ {\varvec{B}}\right) _{ij}={\varvec{A}}_{ij}{\varvec{B}}_{ij}\).

In the present example of four gluons, the TMD gluon distribution matrix reads

$$\begin{aligned} {{{\varvec{\Phi }}}}_{gg\rightarrow gg}=\left( \begin{array}{cc} {\varPhi }_{1} &{} {\varPhi }_{2}\\ {\varPhi }_{2} &{} {\varPhi }_{1} \end{array}\right) , \end{aligned}$$
(44)

with two independent TMD gluon distributions expressed in terms of the basis distributions:

$$\begin{aligned} {\varPhi }_{1}= & {} \frac{1}{2N_{c}^{2}}\left( N_{c}^{2}{\mathcal {F}}_{\text {gg}}^{(1)}-2{\mathcal {F}}_{\text {gg}}^{(3)}+{\mathcal {F}}_{\text {gg}}^{(4)}+{\mathcal {F}}_{\text {gg}}^{(5)}+N_{c}^{2}{\mathcal {F}}_{\text {gg}}^{(6)}\right) , \nonumber \\\end{aligned}$$
(45)
$$\begin{aligned} {\varPhi }_{2}= & {} \frac{1}{N_{c}^{2}}\left( N_{c}^{2}{\mathcal {F}}_{\text {gg}}^{(2)}-2{\mathcal {F}}_{\text {gg}}^{(3)}+{\mathcal {F}}_{\text {gg}}^{(4)}+{\mathcal {F}}_{\text {gg}}^{(5)}+N_{c}^{2}{\mathcal {F}}_{\text {gg}}^{(6)}\right) .\nonumber \\ \end{aligned}$$
(46)

For more complicated processes with gluons it is useful to write the above equations in matrix form:

$$\begin{aligned} \left( \begin{array}{c} {\varPhi }_{1}\\ \vdots \\ {\varPhi }_{k} \end{array}\right) ={\mathbf {M}}\left( \begin{array}{c} {\mathcal {F}}_{\text {gg}}^{(1)}\\ {\mathcal {F}}_{\text {gg}}^{(2)}\\ \vdots \\ {\mathcal {F}}_{\text {gg}}^{(7)} \end{array}\right) , \end{aligned}$$
(47)

where \({\mathbf {M}}\) is a matrix with k rows and 7 columns. For the present case, this matrix reads

$$\begin{aligned} {\mathbf {M}}_{gg\rightarrow gg}=\left( \begin{array}{ccccccc} \frac{1}{2} &{} 0 &{} -\frac{1}{N_{c}^{2}} &{} \frac{1}{2N_{c}^{2}} &{} \frac{1}{2N_{c}^{2}} &{} \frac{1}{2} &{} 0\\ 0 &{} 1 &{} -\frac{2}{N_{c}^{2}} &{} \frac{1}{N_{c}^{2}} &{} \frac{1}{N_{c}^{2}} &{} 1 &{} 0 \end{array}\right) . \end{aligned}$$
(48)

In a similar fashion, one can derive the matrices \({{{\varvec{\Phi }}}}\) and \({\mathbf {M}}\) for other 4 parton channels. The only difference is that for processes with quarks, we always use the color flow color decomposition of an amplitude. For the channel

$$\begin{aligned} g\left( k_{1}\right) q\left( k_{4}\right) \rightarrow g\left( k_{2}\right) q\left( k_{3}\right) , \end{aligned}$$
(49)

we obtain

$$\begin{aligned} {{{\varvec{\Phi }}}}_{gq\rightarrow gq}=\left( \begin{array}{cc} {\varPhi }_{2} &{} {\varPhi }_{1}\\ {\varPhi }_{1} &{} {\varPhi }_{1} \end{array}\right) , \end{aligned}$$
(50)

with the \({\varPhi }_{i}\) given in Table 4. For a similar process with an anti-quark we get

$$\begin{aligned} {{{\varvec{\Phi }}}}_{g{\bar{q}}\rightarrow g{\bar{q}}}=\left( \begin{array}{cc} {\varPhi }_{1} &{} {\varPhi }_{1}\\ {\varPhi }_{1} &{} {\varPhi }_{2} \end{array}\right) . \end{aligned}$$
(51)

Finally, for

$$\begin{aligned} g\left( k_{1}\right) g\left( k_{4}\right) \rightarrow q\left( k_{2}\right) {\bar{q}}\left( k_{3}\right) \end{aligned}$$
(52)

we have

$$\begin{aligned} {{{\varvec{\Phi }}}}_{gg\rightarrow q{\bar{q}}}=\left( \begin{array}{cc} {\varPhi }_{1} &{} {\varPhi }_{2}\\ {\varPhi }_{2} &{} {\varPhi }_{1} \end{array}\right) . \end{aligned}$$
(53)

The partial amplitude vectors \(\vec {{\mathcal {A}}}\) for the above cases are listed in Table 3.

Table 4 Matrices \({\mathbf {M}}\) of structures appearing in four-parton processes. The subscripts in the sub-process indication correspond to the momenta enumeration

5.2 Five partons

The calculation of the TMD gluon distributions with 5 colored partons proceeds in the same fashion, but is technically more complicated. Also, a new feature appears. Certain color factors, building up the matrix \({\mathbf {C}}\), vanish for some processes. However, some of the corresponding TMD operators do not vanish (more precisely, we mean here corresponding elements of the \({\mathbf {F}}\) matrix). It is a special property of the TMD factorization: certain color flows would not contribute in the collinear factorization (where only the matrix \({\mathbf {C}}\) appears), but they do contribute if the TMD gluon distributions are considered. Thus we need to modify the definition of the TMD gluon distribution matrix \({{{\varvec{\Phi }}}}\) (42) and the Eq. (43) for such processes. In both formulas, instead of the matrix \({\mathbf {C}}\) (which has zeros), we use the matrix \({\mathbf {C}}'\) with elements

$$\begin{aligned} {\mathbf {C}}'_{ij}={\left\{ \begin{array}{ll} {\mathbf {C}}_{ij} &{}\quad \text {if }{\mathbf {C}}_{ij}\ne 0\\ 1 &{} \quad \text {if }{\mathbf {C}}_{ij}=0 \end{array}\right. }. \end{aligned}$$
(54)

This is a simple way to extract the hard matrix element color factors only from those TMD operators, for which the color factor is nonzero. For reader’s convenience, the color factors for 5 parton processes in the color-ordered-amplitude representation are collected in “Appendix D” (they were cross-checked with [68, 71]).

Below, we present the TMD gluon distribution matrices \({{{\varvec{\Phi }}}}\) for various channels. The vectors \(\vec {{\mathcal {A}}}\) of partial amplitudes, corresponding to the entries of the matrices \({{{\varvec{\Phi }}}}\), are given in Table 5 in “Appendix A”. The TMD gluon distributions \({\varPhi }_{i}\) building up these matrices, are expressed through the ’basis’ distributions (25)–(34), as given by the \({\mathbf {M}}\) matrices listed in Table 6 (“Appendix A”). The \({\mathbf {M}}\) matrices for processes in which an incoming and outgoing quarks are replaced by incoming and outgoing anti-quarks are the same.

For the pure gluonic process,

$$\begin{aligned} g\left( k_{1}\right) g\left( k_{5}\right) \rightarrow g\left( k_{2}\right) g\left( k_{3}\right) g\left( k_{4}\right) , \end{aligned}$$
(55)

we obtain

$$\begin{aligned} {{{\varvec{\Phi }}}}_{gg\rightarrow ggg}=\left( \begin{array}{cccccc} {\varPhi }_{1} &{} {\varPhi }_{2} &{} {\varPhi }_{2} &{} {\varPhi }_{3} &{} {\varPhi }_{3} &{} {\varPhi }_{4}^{*}\\ {\varPhi }_{2} &{} {\varPhi }_{1} &{} {\varPhi }_{3} &{} {\varPhi }_{4}^{*} &{} {\varPhi }_{2} &{} {\varPhi }_{3}\\ {\varPhi }_{2} &{} {\varPhi }_{3} &{} {\varPhi }_{1} &{} {\varPhi }_{2} &{} {\varPhi }_{4}^{*} &{} {\varPhi }_{3}\\ {\varPhi }_{3} &{} {\varPhi }_{4}^{*} &{} {\varPhi }_{2} &{} {\varPhi }_{1} &{} {\varPhi }_{3} &{} {\varPhi }_{2}\\ {\varPhi }_{3} &{} {\varPhi }_{2} &{} {\varPhi }_{4}^{*} &{} {\varPhi }_{3} &{} {\varPhi }_{1} &{} {\varPhi }_{2}\\ {\varPhi }_{4}^{*} &{} {\varPhi }_{3} &{} {\varPhi }_{3} &{} {\varPhi }_{2} &{} {\varPhi }_{2} &{} {\varPhi }_{1} \end{array}\right) . \end{aligned}$$
(56)

The \({\varPhi }_{i}\) gluon distributions are listen in the first row of Table 6. This process has the property mentioned in the beginning of this section. The entries for which the color factors are zero are marked with the asterix \(^{*}\).

For

$$\begin{aligned} g\left( k_{1}\right) g\left( k_{5}\right) \rightarrow q\left( k_{2}\right) {\bar{q}}\left( k_{3}\right) g\left( k_{4}\right) , \end{aligned}$$
(57)

we get

$$\begin{aligned} {{{\varvec{\Phi }}}}_{gg\rightarrow q{\bar{q}}g}=\left( \begin{array}{cccccc} {\varPhi }_{1} &{} {\varPhi }_{2} &{} {\varPhi }_{2} &{} {\varPhi }_{3} &{} {\varPhi }_{3} &{} {\varPhi }_{4}\\ {\varPhi }_{2} &{} {\varPhi }_{2} &{} {\varPhi }_{5} &{} {\varPhi }_{6} &{} {\varPhi }_{3} &{} {\varPhi }_{3}\\ {\varPhi }_{2} &{} {\varPhi }_{5} &{} {\varPhi }_{2} &{} {\varPhi }_{3} &{} {\varPhi }_{6} &{} {\varPhi }_{3}\\ {\varPhi }_{3} &{} {\varPhi }_{6} &{} {\varPhi }_{3} &{} {\varPhi }_{2} &{} {\varPhi }_{5} &{} {\varPhi }_{2}\\ {\varPhi }_{3} &{} {\varPhi }_{3} &{} {\varPhi }_{6} &{} {\varPhi }_{5} &{} {\varPhi }_{2} &{} {\varPhi }_{2}\\ {\varPhi }_{4} &{} {\varPhi }_{3} &{} {\varPhi }_{3} &{} {\varPhi }_{2} &{} {\varPhi }_{2} &{} {\varPhi }_{1} \end{array}\right) , \end{aligned}$$
(58)

with the TMD gluon distributions given in the second row of Table 6.

For the process with initial state quark

$$\begin{aligned} g\left( k_{1}\right) q\left( k_{5}\right) \rightarrow g\left( k_{2}\right) g\left( k_{3}\right) q\left( k_{4}\right) , \end{aligned}$$
(59)

or anti-quark, we obtain, respectively

$$\begin{aligned} {{{\varvec{\Phi }}}}_{gq\rightarrow ggq}=\left( \begin{array}{cccccc} {\varPhi }_{1} &{} {\varPhi }_{2} &{} {\varPhi }_{3} &{} {\varPhi }_{4} &{} {\varPhi }_{5} &{} {\varPhi }_{4}\\ {\varPhi }_{2} &{} {\varPhi }_{1} &{} {\varPhi }_{5} &{} {\varPhi }_{4} &{} {\varPhi }_{3} &{} {\varPhi }_{4}\\ {\varPhi }_{3} &{} {\varPhi }_{5} &{} {\varPhi }_{3} &{} {\varPhi }_{4} &{} {\varPhi }_{6} &{} {\varPhi }_{4}\\ {\varPhi }_{4} &{} {\varPhi }_{4} &{} {\varPhi }_{4} &{} {\varPhi }_{4} &{} {\varPhi }_{4} &{} {\varPhi }_{4}\\ {\varPhi }_{5} &{} {\varPhi }_{3} &{} {\varPhi }_{6} &{} {\varPhi }_{4} &{} {\varPhi }_{3} &{} {\varPhi }_{4}\\ {\varPhi }_{4} &{} {\varPhi }_{4} &{} {\varPhi }_{4} &{} {\varPhi }_{4} &{} {\varPhi }_{4} &{} {\varPhi }_{4} \end{array}\right) , \end{aligned}$$
(60)

and

$$\begin{aligned} {{{\varvec{\Phi }}}}_{g{\bar{q}}\rightarrow gg{\bar{q}}}=\left( \begin{array}{cccccc} {\varPhi }_{4} &{} {\varPhi }_{4} &{} {\varPhi }_{4} &{} {\varPhi }_{4} &{} {\varPhi }_{4} &{} {\varPhi }_{4}\\ {\varPhi }_{4} &{} {\varPhi }_{4} &{} {\varPhi }_{4} &{} {\varPhi }_{4} &{} {\varPhi }_{4} &{} {\varPhi }_{4}\\ {\varPhi }_{4} &{} {\varPhi }_{4} &{} {\varPhi }_{3} &{} {\varPhi }_{3} &{} {\varPhi }_{6} &{} {\varPhi }_{5}\\ {\varPhi }_{4} &{} {\varPhi }_{4} &{} {\varPhi }_{3} &{} {\varPhi }_{1} &{} {\varPhi }_{5} &{} {\varPhi }_{2}\\ {\varPhi }_{4} &{} {\varPhi }_{4} &{} {\varPhi }_{6} &{} {\varPhi }_{5} &{} {\varPhi }_{3} &{} {\varPhi }_{3}\\ {\varPhi }_{4} &{} {\varPhi }_{4} &{} {\varPhi }_{5} &{} {\varPhi }_{2} &{} {\varPhi }_{3} &{} {\varPhi }_{1} \end{array}\right) , \end{aligned}$$
(61)

with the TMD gluon distributions given in the third row of Table 6. These matrices differ by the permutations of the entries, which has its origin in a slightly different color decomposition for quarks and anti-quarks. Namely the order of quark–anti-quark lines (with the outgoing-momenta convention) is reversed in one case with respect to the other.

Finally, the processes with two quark–anti-quark pairs, with incoming quark

$$\begin{aligned} g\left( k_{1}\right) q\left( k_{5}\right) \rightarrow q\left( k_{2}\right) {\bar{q}}\left( k_{3}\right) q\left( k_{4}\right) , \end{aligned}$$
(62)

or anti-quark, involve respectively

$$\begin{aligned} {{{\varvec{\Phi }}}}_{gq\rightarrow q{\bar{q}}q}=\left( \begin{array}{cccc} {\varPhi }_{1} &{} 0 &{} {\varPhi }_{1} &{} {\varPhi }_{1}\\ 0 &{} {\varPhi }_{2} &{} {\varPhi }_{3} &{} {\varPhi }_{1}\\ {\varPhi }_{1} &{} {\varPhi }_{3} &{} {\varPhi }_{2} &{} 0\\ {\varPhi }_{1} &{} {\varPhi }_{1} &{} 0 &{} {\varPhi }_{1} \end{array}\right) , \end{aligned}$$
(63)

and

$$\begin{aligned} {{{\varvec{\Phi }}}}_{g{\bar{q}}\rightarrow q{\bar{q}}{\bar{q}}}=\left( \begin{array}{cccc} {\varPhi }_{1} &{} 0 &{} {\varPhi }_{1} &{} {\varPhi }_{1}\\ 0 &{} {\varPhi }_{2} &{} {\varPhi }_{1} &{} {\varPhi }_{3}\\ {\varPhi }_{1} &{} {\varPhi }_{1} &{} {\varPhi }_{1} &{} 0\\ {\varPhi }_{1} &{} {\varPhi }_{3} &{} 0 &{} {\varPhi }_{2} \end{array}\right) . \end{aligned}$$
(64)

The TMD gluon distributions appearing in these matrices are listed in the fourth row of Table 6. Interestingly, for this case, not only some of the color factors vanish, but also the corresponding TMDs.

As the potential phenomenological application of the results (in short run) concerns rather the large \(N_{c}\) limit, we present the relevant matrices in this limit in “Appendix C” (Table 12).

5.3 Six partons

Six parton processes do not involve new features, except more channels and more involved calculations. The vectors \(\vec {{\mathcal {A}}}\) of the partial amplitudes, and the \({\mathbf {M}}\) matrices are given in Tables 7 and 8, 9, 10, 11 in “Appendix A”. The \({\mathbf {M}}\) matrices for processes in which an incoming and outgoing quarks are replaced by incoming and outgoing anti-quarks are the same. Below, we present results for the TMD gluon distribution matrices \({{{\varvec{\Phi }}}}\) for all channels. The number of partial amplitudes necessitates the use of block matrices to compactify the notation.

For the six-gluon process,

$$\begin{aligned} g\left( k_{1}\right) g\left( k_{6}\right) \rightarrow g\left( k_{2}\right) g\left( k_{3}\right) g\left( k_{4}\right) g\left( k_{5}\right) , \end{aligned}$$
(65)

the \({{{\varvec{\Phi }}}}\) matrix is

$$\begin{aligned} {{{\varvec{\Phi }}}}_{gg\rightarrow gggg}=\left( \begin{array}{cccc} T_{1} &{} T_{2} &{} T_{3} &{} T_{4}\\ T_{2} &{} T_{1} &{} T_{5} &{} T_{6}\\ T_{3}^{\intercal } &{} T_{5} &{} T_{1} &{} T_{7}\\ T_{4}^{\intercal } &{} T_{6}^{\intercal } &{} T_{7} &{} T_{1} \end{array}\right) , \end{aligned}$$
(66)

where \(T_{i}\) are 6\(\times \)6 block matrices given by Eqs. (B.1)–(B.4) in “Appendix B.1”.

In the present case we have two TMD operators, for which the color factor vanishes – \({\varPhi }_{4}^{*}\) and \({\varPhi }_{8}^{*}\) (we remind, that we mark these matrix elements with an asterix). The full list of the TMD gluon distributions is given in the Table 8.

Next consider the process

$$\begin{aligned} g\left( k_{1}\right) g\left( k_{6}\right) \rightarrow q\left( k_{2}\right) {\bar{q}}\left( k_{3}\right) g\left( k_{4}\right) g\left( k_{5}\right) . \end{aligned}$$
(67)

The TMD matrix reads

$$\begin{aligned} {{{\varvec{\Phi }}}}_{gg\rightarrow q{\bar{q}}gg}=\left( \begin{array}{cccc} T_{1} &{} T_{2} &{} T_{3} &{} T_{4}\\ T_{2}^{\intercal } &{} T_{5} &{} T_{6} &{} T_{7}\\ T_{3}^{\intercal } &{} T_{6} &{} T_{5} &{} T_{8}\\ T_{4}^{\intercal } &{} T_{7}^{\intercal } &{} T_{8}^{\intercal } &{} T_{9} \end{array}\right) , \end{aligned}$$
(68)

where the blocks gathered in Eqs. (B.8)–(B.12) in “Appendix B.2”. The TMD gluon distributions are given in the Table 9.

For the process

$$\begin{aligned} g\left( k_{1}\right) q\left( k_{6}\right) \rightarrow g\left( k_{2}\right) g\left( k_{3}\right) g\left( k_{4}\right) q\left( k_{5}\right) , \end{aligned}$$
(69)

the TMD matrix reads

$$\begin{aligned} {{{\varvec{\Phi }}}}_{gq\rightarrow gggq}=\left( \begin{array}{cccc} T_{1} &{} T_{2} &{} T_{3} &{} T_{4}\\ T_{2}^{\intercal } &{} T_{5} &{} T_{6} &{} T_{7}\\ T_{3}^{\intercal } &{} T_{6} &{} T_{5} &{} T_{8}\\ T_{4}^{\intercal } &{} T_{7}^{\intercal } &{} T_{8} &{} T_{5} \end{array}\right) , \end{aligned}$$
(70)

with the blocks expressed by Eqs. (B.13)–(B.16) in “Appendix B.3”. The TMD distributions are in Table 10.

Similarly, for the process with the anti-quark

$$\begin{aligned} g\left( k_{1}\right) {\bar{q}}\left( k_{6}\right) \rightarrow g\left( k_{2}\right) g\left( k_{3}\right) g\left( k_{4}\right) {\bar{q}}\left( k_{5}\right) , \end{aligned}$$
(71)

we get

$$\begin{aligned} {{{\varvec{\Phi }}}}_{g{\bar{q}}\rightarrow ggg{\bar{q}}}=\left( \begin{array}{cccc} T_{1} &{} T_{1} &{} T_{1} &{} T_{1}\\ T_{1} &{} T_{2} &{} T_{3} &{} T_{4}\\ T_{1} &{} T_{3} &{} T_{2} &{} T_{5}\\ T_{1} &{} T_{4}^{\intercal } &{} T_{5} &{} T_{2} \end{array}\right) , \end{aligned}$$
(72)

with the blocks given by Eqs. (B.17)–(B.19) in “Appendix B.4”. The TMD distributions are in Table 10.

Processes with two quark–anti-quark pairs have smaller number of partial amplitudes. For the process

$$\begin{aligned} g\left( k_{1}\right) g\left( k_{6}\right) \rightarrow q\left( k_{2}\right) {\bar{q}}\left( k_{3}\right) q\left( k_{4}\right) {\bar{q}}\left( k_{5}\right) , \end{aligned}$$
(73)

we obtain

$$\begin{aligned} {{{\varvec{\Phi }}}}_{gg\rightarrow q{\bar{q}}q{\bar{q}}}=\left( \begin{array}{cc} T_{1} &{} T_{2}\\ T_{2} &{} T_{1} \end{array}\right) , \end{aligned}$$
(74)

with only two blocks:

$$\begin{aligned} T_{1}= & {} \left( \begin{array}{cccccc} {\varPhi }_{1} &{} 0 &{} {\varPhi }_{2} &{} {\varPhi }_{3} &{} 0 &{} {\varPhi }_{2}\\ 0 &{} {\varPhi }_{4} &{} 0 &{} 0 &{} {\varPhi }_{5} &{} 0\\ {\varPhi }_{2} &{} 0 &{} {\varPhi }_{1} &{} {\varPhi }_{2} &{} 0 &{} {\varPhi }_{3}\\ {\varPhi }_{3} &{} 0 &{} {\varPhi }_{2} &{} {\varPhi }_{1} &{} 0 &{} {\varPhi }_{2}\\ 0 &{} {\varPhi }_{5} &{} 0 &{} 0 &{} {\varPhi }_{4} &{} 0\\ {\varPhi }_{2} &{} 0 &{} {\varPhi }_{3} &{} {\varPhi }_{2} &{} 0 &{} {\varPhi }_{1} \end{array}\right) ,\nonumber \\ T_{2}= & {} \left( \begin{array}{cccccc} {\varPhi }_{1} &{} {\varPhi }_{3} &{} {\varPhi }_{1} &{} {\varPhi }_{3} &{} {\varPhi }_{1} &{} {\varPhi }_{3}\\ {\varPhi }_{3} &{} {\varPhi }_{6} &{} {\varPhi }_{1} &{} {\varPhi }_{1} &{} {\varPhi }_{6} &{} {\varPhi }_{3}\\ {\varPhi }_{1} &{} {\varPhi }_{1} &{} {\varPhi }_{1} &{} {\varPhi }_{3} &{} {\varPhi }_{3} &{} {\varPhi }_{3}\\ {\varPhi }_{3} &{} {\varPhi }_{1} &{} {\varPhi }_{3} &{} {\varPhi }_{1} &{} {\varPhi }_{3} &{} {\varPhi }_{1}\\ {\varPhi }_{1} &{} {\varPhi }_{6} &{} {\varPhi }_{3} &{} {\varPhi }_{3} &{} {\varPhi }_{6} &{} {\varPhi }_{1}\\ {\varPhi }_{3} &{} {\varPhi }_{3} &{} {\varPhi }_{3} &{} {\varPhi }_{1} &{} {\varPhi }_{1} &{} {\varPhi }_{1} \end{array}\right) . \end{aligned}$$
(75)

The TMD distributions are given in the Table 11.

For the process

$$\begin{aligned} g\left( k_{1}\right) q\left( k_{6}\right) \rightarrow g\left( k_{2}\right) q\left( k_{3}\right) {\bar{q}}\left( k_{4}\right) q\left( k_{5}\right) , \end{aligned}$$
(76)

we have

$$\begin{aligned} {{{\varvec{\Phi }}}}_{gq\rightarrow gq{\bar{q}}q}=\left( \begin{array}{cc} T_{1} &{} T_{2}\\ T_{2}^{\intercal } &{} T_{3} \end{array}\right) , \end{aligned}$$
(77)

with three different blocks

$$\begin{aligned} T_{1}= & {} \left( \begin{array}{cccccc} {\varPhi }_{1} &{} 0 &{} {\varPhi }_{2} &{} {\varPhi }_{3} &{} 0 &{} {\varPhi }_{2}\\ 0 &{} {\varPhi }_{4} &{} {\varPhi }_{5}^{*} &{} 0 &{} {\varPhi }_{3} &{} {\varPhi }_{6}^{*}\\ {\varPhi }_{2} &{} {\varPhi }_{5}^{*} &{} {\varPhi }_{4} &{} {\varPhi }_{3} &{} 0 &{} {\varPhi }_{4}\\ {\varPhi }_{3} &{} 0 &{} {\varPhi }_{3} &{} {\varPhi }_{3} &{} 0 &{} {\varPhi }_{3}\\ 0 &{} {\varPhi }_{3} &{} 0 &{} 0 &{} {\varPhi }_{3} &{} 0\\ {\varPhi }_{2} &{} {\varPhi }_{6}^{*} &{} {\varPhi }_{4} &{} {\varPhi }_{3} &{} 0 &{} {\varPhi }_{4} \end{array}\right) ,\nonumber \\ T_{2}= & {} \left( \begin{array}{cccccc} {\varPhi }_{1} &{} {\varPhi }_{3} &{} {\varPhi }_{1} &{} {\varPhi }_{3} &{} {\varPhi }_{1} &{} {\varPhi }_{3}\\ {\varPhi }_{2} &{} {\varPhi }_{3} &{} {\varPhi }_{1} &{} {\varPhi }_{2} &{} {\varPhi }_{7} &{} {\varPhi }_{3}\\ {\varPhi }_{2} &{} {\varPhi }_{3} &{} {\varPhi }_{1} &{} {\varPhi }_{2} &{} {\varPhi }_{2} &{} {\varPhi }_{3}\\ {\varPhi }_{3} &{} {\varPhi }_{3} &{} {\varPhi }_{3} &{} {\varPhi }_{3} &{} {\varPhi }_{3} &{} {\varPhi }_{3}\\ {\varPhi }_{3} &{} {\varPhi }_{3} &{} {\varPhi }_{3} &{} {\varPhi }_{3} &{} {\varPhi }_{3} &{} {\varPhi }_{3}\\ {\varPhi }_{2} &{} {\varPhi }_{3} &{} {\varPhi }_{3} &{} {\varPhi }_{7} &{} {\varPhi }_{2} &{} {\varPhi }_{3} \end{array}\right) , \end{aligned}$$
(78)
$$\begin{aligned} T_{3}= & {} \left( \begin{array}{cccccc} {\varPhi }_{4} &{} 0 &{} {\varPhi }_{2} &{} {\varPhi }_{4} &{} {\varPhi }_{8}^{*} &{} {\varPhi }_{3}\\ 0 &{} {\varPhi }_{3} &{} 0 &{} 0 &{} {\varPhi }_{3} &{} 0\\ {\varPhi }_{2} &{} 0 &{} {\varPhi }_{1} &{} {\varPhi }_{2} &{} 0 &{} {\varPhi }_{3}\\ {\varPhi }_{4} &{} 0 &{} {\varPhi }_{2} &{} {\varPhi }_{4} &{} {\varPhi }_{9}^{*} &{} {\varPhi }_{3}\\ {\varPhi }_{8}^{*} &{} {\varPhi }_{3} &{} 0 &{} {\varPhi }_{9}^{*} &{} {\varPhi }_{4} &{} 0\\ {\varPhi }_{3} &{} 0 &{} {\varPhi }_{3} &{} {\varPhi }_{3} &{} 0 &{} {\varPhi }_{3} \end{array}\right) . \end{aligned}$$
(79)

Note, that in this process there appear both the vanishing structures for vanishing color factors and non-vanishing structures for vanishing color factors. The list of the TMD distributions is given in the Table 11. Similarly for the process with an anti-quark, we get:

$$\begin{aligned} {{{\varvec{\Phi }}}}_{g{\bar{q}}\rightarrow gq{\bar{q}}{\bar{q}}}=\left( \begin{array}{cc} T_{1} &{} T_{2}\\ T_{2} &{} T_{3} \end{array}\right) , \end{aligned}$$
(80)

with

$$\begin{aligned} T_{1}= & {} \left( \begin{array}{cccccc} {\varPhi }_{3} &{} 0 &{} {\varPhi }_{3} &{} {\varPhi }_{3} &{} 0 &{} {\varPhi }_{3}\\ 0 &{} {\varPhi }_{4} &{} {\varPhi }_{6}^{*} &{} 0 &{} {\varPhi }_{3} &{} {\varPhi }_{5}^{*}\\ {\varPhi }_{3} &{} {\varPhi }_{6}^{*} &{} {\varPhi }_{4} &{} {\varPhi }_{2} &{} 0 &{} {\varPhi }_{4}\\ {\varPhi }_{3} &{} 0 &{} {\varPhi }_{2} &{} {\varPhi }_{1} &{} 0 &{} {\varPhi }_{2}\\ 0 &{} {\varPhi }_{3} &{} 0 &{} 0 &{} {\varPhi }_{3} &{} 0\\ {\varPhi }_{3} &{} {\varPhi }_{5}^{*} &{} {\varPhi }_{4} &{} {\varPhi }_{2} &{} 0 &{} {\varPhi }_{4} \end{array}\right) ,\nonumber \\ T_{2}= & {} \left( \begin{array}{cccccc} {\varPhi }_{3} &{} {\varPhi }_{3} &{} {\varPhi }_{3} &{} {\varPhi }_{3} &{} {\varPhi }_{3} &{} {\varPhi }_{3}\\ {\varPhi }_{3} &{} {\varPhi }_{7} &{} {\varPhi }_{2} &{} {\varPhi }_{1} &{} {\varPhi }_{3} &{} {\varPhi }_{2}\\ {\varPhi }_{3} &{} {\varPhi }_{2} &{} {\varPhi }_{7} &{} {\varPhi }_{3} &{} {\varPhi }_{3} &{} {\varPhi }_{2}\\ {\varPhi }_{3} &{} {\varPhi }_{1} &{} {\varPhi }_{3} &{} {\varPhi }_{1} &{} {\varPhi }_{3} &{} {\varPhi }_{1}\\ {\varPhi }_{3} &{} {\varPhi }_{3} &{} {\varPhi }_{3} &{} {\varPhi }_{3} &{} {\varPhi }_{3} &{} {\varPhi }_{3}\\ {\varPhi }_{3} &{} {\varPhi }_{2} &{} {\varPhi }_{2} &{} {\varPhi }_{1} &{} {\varPhi }_{3} &{} {\varPhi }_{2} \end{array}\right) , \end{aligned}$$
(81)
$$\begin{aligned} T_{3}= & {} \left( \begin{array}{cccccc} {\varPhi }_{3} &{} 0 &{} {\varPhi }_{3} &{} {\varPhi }_{3} &{} 0 &{} {\varPhi }_{3}\\ 0 &{} {\varPhi }_{4} &{} {\varPhi }_{9}^{*} &{} 0 &{} {\varPhi }_{3} &{} {\varPhi }_{8}^{*}\\ {\varPhi }_{3} &{} {\varPhi }_{9}^{*} &{} {\varPhi }_{4} &{} {\varPhi }_{2} &{} 0 &{} {\varPhi }_{4}\\ {\varPhi }_{3} &{} 0 &{} {\varPhi }_{2} &{} {\varPhi }_{1} &{} 0 &{} {\varPhi }_{2}\\ 0 &{} {\varPhi }_{3} &{} 0 &{} 0 &{} {\varPhi }_{3} &{} 0\\ {\varPhi }_{3} &{} {\varPhi }_{8}^{*} &{} {\varPhi }_{4} &{} {\varPhi }_{2} &{} 0 &{} {\varPhi }_{4} \end{array}\right) . \end{aligned}$$
(82)

The large \(N_{c}\) limits of gluon distributions for 6 parton processes were gathered in Tables 13, 14, 15, 16 in “Appendix C”. Additionally, we collect the color factors for all processes in “Appendix D”.

6 Large \(N_{c}\) analysis for arbitrary number of gluons

In this section, we shall utilize the color flow method to give the large \(N_{c}\) results for a process with n gluons

$$\begin{aligned} g\left( k_{1}\right) g\left( k_{n}\right) \rightarrow g\left( k_{2}\right) \dots g\left( k_{n-1}\right) . \end{aligned}$$
(83)

We shall use the fact that the color flow decomposition (9) involves all \(\left( n-1\right) !\) partial amplitudes which are the same as in the fundamental decomposition (7). Therefore, the leading \(N_{c}\) contribution is given by the partial amplitudes squared (the interference terms are subleading) [72]

$$\begin{aligned} \left| {\mathcal {M}}\right| ^{2}={\mathcal {C}}\sum _{\pi \in S_{n}/Z_{n}}\left\{ \left| {\mathcal {A}}\left( \pi \left( 1\right) ,\dots ,\pi \left( n\right) \right) \right| ^{2}+{\mathcal {O}}\left( \frac{1}{N_{c}^{2}}\right) \right\} ,\nonumber \\ \end{aligned}$$
(84)

with \({\mathcal {C}}\) being some color coefficient. Note, that if we used the adjoint color decomposition to reduce the number of partial amplitudes only to the linearly independent ones, as we did in the previous section, we would not be able to claim (84). Consequently, the general analysis of large \(N_{c}\) would be very difficult. Therefore, there is a trade off: switching to a general argumentation requires giving up the advantage of using minimal number of amplitudes. In practice, however, any partial amplitude can be easily calculated numerically, so the real loss is not so big.

Based on the above, the idea is to calculate first the diagonal elements of matrix \({{{\varvec{\Phi }}}}\), as they will definitely contribute in the large \(N_{c}\) limit. This would be the final answer, if there is no enhancement of powers of \(N_{c}^{2}\) for some of the non-diagonal elements. In fact, as we shall see, the enhancement indeed occurs, but still the TMD gluon distribution appearing off the diagonal is numerically small.

Let us start with calculating the diagonal elements of the TMD gluon distribution matrix \({{{\varvec{\Phi }}}}\). It is sufficient to consider only the following diagrams:

(85)

The first diagram from the left corresponds to the partial amplitude squared \(\left| {\mathcal {A}}\left( 1,2,\dots ,n\right) \right| ^{2}\) and the TMD operator reads (after dividing by the corresponding color factor)

$$\begin{aligned} \frac{N_{c}^{n-3}}{N_{c}^{n-2}}\mathrm {Tr}\left\{ F\left( \xi \right) {\mathcal {U}}^{\left[ -\right] \dagger }F\left( 0\right) {\mathcal {U}}^{\left[ +\right] }\right\} \mathrm {Tr}{\mathcal {U}}^{\left[ \square \right] \dagger }\,\rightsquigarrow \,{\mathcal {F}}_{gg}^{\left( 1\right) }, \end{aligned}$$
(86)

i.e. it corresponds to the TMD \({\mathcal {F}}_{gg}^{\left( 1\right) }\), Eq. (28). However, any permutation of the following \(\left( n-2\right) \) final state legs will give the same contribution, thus, the set

$$\begin{aligned} \left\{ \left. \left| {\mathcal {A}}\left( 1,\pi \left( 2\right) ,\pi \left( 3\right) ,\dots ,\pi \left( n-1\right) ,n\right) \right| ^{2}\right| \pi \in S_{n-2}\right\} \rightsquigarrow \,{\mathcal {F}}_{gg}^{\left( 1\right) }.\nonumber \\ \end{aligned}$$
(87)

The second diagram, corresponding to \(|{\mathcal {A}}\big (1,2,\dots ,n,n-1\big )|^{2}\), gives

$$\begin{aligned} \frac{N_{c}^{n-4}}{N_{c}^{n-2}}\mathrm {Tr}\left\{ F\left( \xi \right) {\mathcal {U}}^{\left[ +\right] \dagger }F\left( 0\right) {\mathcal {U}}^{\left[ +\right] }\right\} \mathrm {Tr}{\mathcal {U}}^{\left[ \square \right] \dagger }\mathrm {Tr}{\mathcal {U}}^{\left[ \square \right] }\rightsquigarrow \,{\mathcal {F}}_{gg}^{\left( 6\right) }.\nonumber \\ \end{aligned}$$
(88)

Not only any permutation of final states will give the same result, but also any diagram with leg \(k_{n}\) permuted with \(\left\{ 3,\dots ,n-2\right\} \). Thus

$$\begin{aligned}&\left\{ \left. \left| {\mathcal {A}}\left( 1,\pi \left( 2\right) ,\dots ,\pi \left( n-2\right) ,n,\pi \left( n-1\right) \right) \right| ^{2}\right| \pi \in S_{n-2}\right\} \nonumber \\&\quad \cup \left\{ \left. \left| {\mathcal {A}}\left( 1,\pi \left( 2\right) ,\dots ,\pi \left( n-3\right) ,n,\pi \left( n-2\right) ,\pi \left( n-1\right) \right) \right| ^{2}\right| \pi \in S_{n-2}\right\} \nonumber \\&\qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \dots \nonumber \\&\quad \cup \left\{ \left. \left| {\mathcal {A}}\left( 1,\pi \left( 2\right) ,n,\pi \left( 3\right) ,\dots ,\pi \left( n-2\right) ,\pi \left( n-1\right) \right) \right| ^{2}\right| \pi \in S_{n-2}\right\} \nonumber \\&\quad \rightsquigarrow \,{\mathcal {F}}_{gg}^{\left( 6\right) }. \end{aligned}$$
(89)

Finally, the third diagram, gives complex conjugate of the operator in (86), thus also \({\mathcal {F}}_{gg}^{\left( 1\right) }\), because of our assumption of the reality of the correlators. We get therefore

$$\begin{aligned} \left\{ \left. \left| {\mathcal {A}}\left( 1,n,\pi \left( 2\right) ,\pi \left( 3\right) ,\dots ,\pi \left( n-1\right) \right) \right| ^{2}\right| \pi \in S_{n-2}\right\} \rightsquigarrow \,{\mathcal {F}}_{gg}^{\left( 1\right) }. \end{aligned}$$

Now let us put together the above results, using the matrix notation as in Sect. 5. Let us define the partial amplitude vector so that it preserves the block structure emerging above:

$$\begin{aligned} \vec {{\mathcal {A}}}=\left( \begin{array}{c} {\mathcal {A}}({\widehat{1}},2,\dots ,n-2,n-1,{\widehat{n}})\\ \vdots \\ {\mathcal {A}}({\widehat{1}},2,3,\dots ,n-2,{\widehat{n}},n-1)\\ \vdots \\ {\mathcal {A}}({\widehat{1}},2,3,\dots ,{\widehat{n}},n-2,n-1)\\ \vdots \\ {\mathcal {A}}({\widehat{1}},2,{\widehat{n}},3,\dots ,n-2,n-1)\\ \vdots \\ {\mathcal {A}}({\widehat{1}},{\widehat{n}},3,\dots ,n-2,n-1)\\ \vdots \end{array}\right) , \end{aligned}$$
(90)

where we have used hats to denote momenta with fixed position in a given group (the actual ordering in each group doesn’t matter). For this choice of the vector \(\vec {{\mathcal {A}}}\), the diagonal contribution to the matrix \({{{\varvec{\Phi }}}}\) at large \(N_{c}\) reads

$$\begin{aligned} {{{\varvec{\Phi }}}}_{\mathrm {diag}}=\left( \begin{array}{ccc} T_{1} &{} 0 &{} 0\\ 0 &{} T_{2} &{} 0\\ 0 &{} 0 &{} T_{1} \end{array}\right) , \end{aligned}$$
(91)

where

$$\begin{aligned} T_{1}={\mathcal {F}}_{gg}^{\left( 1\right) }{\varvec{1}}_{\left( n-2\right) !},\,\,\,\,\,\,T_{2}={\mathcal {F}}_{gg}^{\left( 6\right) }{\varvec{1}}_{(n-3)\left( n-2\right) !}. \end{aligned}$$
(92)

For example, for \(n=4\), we have explicitly

$$\begin{aligned} {{{\varvec{\Phi }}}}_{\mathrm {diag}}=\left( \begin{array}{cccccc} {\mathcal {F}}_{gg}^{\left( 1\right) } &{} 0 &{} 0 &{} 0 &{} 0 &{} 0\\ 0 &{} {\mathcal {F}}_{gg}^{\left( 1\right) } &{} 0 &{} 0 &{} 0 &{} 0\\ 0 &{} 0 &{} {\mathcal {F}}_{gg}^{\left( 6\right) } &{} 0 &{} 0 &{} 0\\ 0 &{} 0 &{} 0 &{} {\mathcal {F}}_{gg}^{\left( 6\right) } &{} 0 &{} 0\\ 0 &{} 0 &{} 0 &{} 0 &{} {\mathcal {F}}_{gg}^{\left( 1\right) } &{} 0\\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} {\mathcal {F}}_{gg}^{\left( 1\right) } \end{array}\right) . \end{aligned}$$
(93)

Now, let us consider the nondiagonal elements. As said above, these elements will be convoluted with partial amplitudes (interference terms) whose color factors are suppressed by at least \(1/N_{c}^{2}\) (to say it differently, the non-diagonal elements of the color matrix \({\mathbf {C}}\) in (43), if it is calculated in fundamental color decomposition, are subleading of at least \(1/N_{c}^{2}\)). Therefore, they do not contribute in large \(N_{c}\), unless some off-diagonal TMD gluon distribution is enhanced by at least \(N_{c}^{2}\). This still might not be enough, but is a sign that a careful analysis has to be carried.

The most suspicious non-diagonal elements are those, which correspond to color flow diagrams with least number of loops. This is slightly counter-intuitive, but we have to keep in mind that, by definition, we divide the color factors out of the TMD (there are no vanishing color factors for gluons in the color flow representation, unlike for the adjoint representation). Thus, the enhancement may happen if the diagrams with Wilson lines have much more loops than the pure color factor diagrams. It is best to illustrate this by an explicit example. Consider a 4 gluon process and the following interference term:

$$\begin{aligned} {\mathcal {A}}\left( 1,2,3,4\right) {\mathcal {A}}^{*}\left( 1,4,2,3\right) . \end{aligned}$$
(94)

We have the following two leading diagrams for the color factor:

(95)

The \(U\left( 1\right) \) colorless propagator for the \(k_{1}\) leg stems from the projectors that have to be inserted for the final state gluons, cf. (12). Recall, that in general these color factor diagrams have to be divided by \(N_{A}\) to get the color factor with ’open indices’. The first, Möbius-loop-like diagram, cancels with the second:

$$\begin{aligned} \frac{1}{N_{A}}\left( N_{c}^{2}-\frac{1}{N_{c}}N_{c}^{3}\right) =0. \end{aligned}$$
(96)

(We called the first diagram ’Möbius-loop-like diagram’ because one of the internal loops shares its border with the external loop.) Similar cancellation happens for the diagrams, where the \(U\left( 1\right) \) gluon appears for legs \(k_{2}\) and \(k_{3}\). The sub-leading diagrams are those where \(U\left( 1\right) \) colorless gluon is \(k_{4}\), i.e. it crosses the other legs:

(97)

In this case, we get

$$\begin{aligned} \frac{1}{N_{A}}\left( -\frac{1}{N_{c}}N_{c}^{3}+\frac{1}{N_{c}^{2}}N_{c}^{2}\right) =-1. \end{aligned}$$
(98)

Now, let us look at the leading diagram for the TMD operator:

(99)

It reads

$$\begin{aligned} N_{c}\mathrm {Tr}\left\{ F\left( \xi \right) {\mathcal {U}}^{\left[ \square \right] }\right\} \mathrm {Tr}\left\{ F\left( 0\right) {\mathcal {U}}^{\left[ \square \right] \dagger }\right\} =N_{c}^{2}{\mathcal {F}}_{gg}^{\left( 2\right) }. \end{aligned}$$
(100)

Dividing by the leading color factor (divided by \(N_{A}\)) we get finally the following non-diagonal element of the \({{{\varvec{\Phi }}}}\) matrix in the large \(N_{c}\) limit:

$$\begin{aligned} \left( {{{\varvec{\Phi }}}}\right) _{15}=-N_{c}^{2}{\mathcal {F}}_{gg}^{\left( 2\right) }. \end{aligned}$$
(101)

As the color factor for \({\mathcal {A}}\left( 1,2,3,4\right) {\mathcal {A}}^{*}\left( 1,4,3,2\right) \) is suppressed with respect to diagonal elements by \(1/N_{c}^{2}\), the TMD gluon distribution \(-{\mathcal {F}}_{gg}^{\left( 2\right) }\) indeed contributes in the large \(N_{c}\) limit. In a similar manner, but considering much more complicated diagrams with maximal number of crossed lines, one can deduce, that this will be always the case for some non-diagonal elements for any multi gluon process. For example, after a similar but tedious calculation for 5 gluon process, we find that the dominant non-diagonal element is \(-N_{c}^{4}{\mathcal {F}}_{gg}^{\left( 2\right) }/4\). We will always get the \({\mathcal {F}}_{gg}^{\left( 2\right) }\) TMD gluon distribution, because of the Möbius-loop-like structure, which gives the two traces appearing in the definition (29).

While perhaps it is possible to derive the answer for the non-diagonal leading \(N_{c}\) elements for any n, let us note that \({\mathcal {F}}_{gg}^{\left( 2\right) }\) gives numerically rather small contribution to the cross section, compared to the other gluon distributions [21, 22]. Indeed, it vanishes very quickly with \(k_{T}\), so that it is small for transverse momenta around the saturation scale. Moreover, it does not survive the collinear limit. Therefore, in possible phenomenological studies of multigluon production, it is safe to set

$$\begin{aligned} {{{\varvec{\Phi }}}}_{gg\rightarrow g\dots g}={{{\varvec{\Phi }}}}_{\mathrm {diag}}. \end{aligned}$$
(102)

The study of large \(N_{c}\) limit for multiparton processes with quarks, and for gluons without the approximation described above, is left for a separate work.

7 Summary

In the present paper we have faced the task of calculating the TMD gluon distributions for processes with five and six colored partons, following the procedure of [17]. So far, in the literature, processes with four partons were considered. Although it is known that within the formal TMD factorization the generalized factorization fails for processes with more than two colored partons, it was argued from the CGC theory that at small-x for dilute-dense collisions exactly such structures appear. At leading order, our results are sufficient to calculate three and four jet production in the gluon saturation regime in dilute-dense collisions, provided the two new basic TMD gluon distributions, \({\mathcal {F}}_{qg}^{\left( 3\right) }\) and \({\mathcal {F}}_{gg}^{\left( 7\right) }\), defined respectively in (27) and (34), are determined. This can be done using the B-JIMWLK equation, as in [22]. At tree level, the hard matrix elements can be easily obtained from available software for automatic calculations, for example KaTie [73], which can deal with on-shell and off-shell initial states.

Instead of calculating the structure of the operators for particular Feynman diagrams, we have used the color decomposition, which is the most efficient way of dealing with the multi-particle QCD amplitudes. In particular, for processes with quarks, we have used the color flow decomposition, which treats quarks and gluons on equal footing. In addition, we formulated straightforward color flow Feynman rules for the gauge links, which allow immediate derivation of the TMD operator for a given color flow.

The color flow Feynman rules are particularly convenient for large \(N_{c}\) analysis. In the present work, as a first step towards this goal, we attacked multigluon processes with arbitrary number of legs. We find a general answer, but in a certain approximation, motivated by known numerical studies of a particular TMD gluon distribution. In the large \(N_{c}\) limit, we find that only two structures contribute, for any number of legs. This is similar to the conclusion made in [74], where the universality at large \(N_{c}\) was found in the multiparton production in the Color Glass Condensate: only dipoles and quadrupoles contribute.

Finally, it would be very interesting to compare the TMD factorization formulae for three jet production (by factorization we mean the approach described in detail in the Introduction) with the leading power of the corresponding CGC result, which was recently derived in [75].