High Energy QCD at NLO: from light-cone wave function to JIMWLK evolution

Soft components of the light cone wave-function of a fast moving projectile hadron is computed in perturbation theory to third order in QCD coupling constant. At this order, the Fock space of the soft modes consists of one-gluon, two-gluon, and a quark-antiquark states. The hard component of the wave-function acts as a non-Abelian background field for the soft modes and is represented by a valence charge distribution that accounts for non-linear density effects in the projectile. When scattered off a dense target, the diagonal element of the S-matrix reveals the Hamiltonian of high energy evolution, the JIMWLK Hamiltonian. This way we provide a new direct derivation of the JIMWLK Hamiltonian at the Next-to-Leading Order.


Introduction and Summary
The phenomenon of perturbative saturation is believed to occur in high energy hadronic collisions and is most pronounced in processes involving nuclei. The idea dates back to the seminal paper [1]. It has been long known theoretically [2] that at high energy gluonic density in a hadron grows exponentially. It was suggested in [1] that this growth does not continue indefinitely, but instead the gluon phase space density saturates when it reaches a critical value of order 1/α s .
In the years since the basic idea has been proposed the theoretical understanding of saturation has progressed significantly [5,6]. We now have a first principles QCD formulation of high energy evolution, including the saturation effects (see e.g. [3,4]). The basic framework was derived in [7] -the Balitsky's hierarchy, and in a series of papers [8] the Jalilian-Marian, Iancu Here the rapidity Y ∼ ln (energy). The JIMWLK Hamiltonian is applicable when one of the colliding particles is dilute (small parton number), such as in forward production in pA collisions at the LHC or in deep inelastic scattering (DIS). While in the dilute-dilute limit, the JIMWLK equation reduces to the linear BFKL equation [2] and its BKP extension [9], it is free from the famous infrared diffusion problem. Furthermore, the major success associated with the JIMWLK is in that it seemingly restores the s-channel unitarity of scattering amplitudes (see, however, [10] where the s-channel unitarity of the JIMWLK was argued to be incomplete).
Phenomenology based on complete LO JIMWLK has emerged recently [11][12][13][14] thanks to its reformulation in terms of stochastic Langevin equation [13]. Yet, for most phenomenological applications JIMWLK is commonly replaced by the Balitsky-Kovchegov (BK) equation [7,15,16], which is a mean field/large N c reduction of the JIMWLK. During the last decade, these developments have become a basis for phenomenological studies of saturation physics applied to high energy/low x collision data.
HERA has accumulated a vast amount of precise data in the low x region, which is particularly useful for calibration of gluon distributions used as an input for calculations of less inclusive processes at RHIC and LHC. Several generations of successful saturation-based fits to low x HERA data are available, from the first work [17], to the latest state of the art AAMQS phenomenology [18]. These provide comprehensive good quality description of the low x data, including the latest available combined H1-ZEUS set [19]. However, HERA was not convincing in experimentally confirming the saturation phenomenon.
Even though the main features of CGC are qualitatively understood, quantitative description is still at an unsatisfactory level as most of the studies have been performed to the lowest order (LO) in α s . It is, however, known that higher order corrections could be very significant and their inclusion crucial for qualitatively reliable phenomenology. Calculational control enables to discriminate between different approaches, some of which describe experimental data but do not involve saturation, is critically important and is the main challenge. Precision mandates inclusion of next-to-leading order (NLO) corrections to evolution and improvement of fixed order calculations for a variety of observables.
The era of saturation-based phenomenology with NLO precision just started. Among recent progress in fixed order NLO calculations it is important to mention the virtual photon impact factor computed in [20][21][22] and a series of papers on forward particle production in pA collisions [23,24] (see e.g [25] for a review). The full set of NLO corrections to BK equation was completed in [26] following on the earlier works [27,28]. When linearised, the NLO BK was shown to coincide with the NLO BFKL computed over a decade earlier [29]. Running coupling effects are known to be phenomenologically important. While they do not account for all the NLO effects in the evolution, they do carry important information about high order corrections. Those were analysed in [27,30,31]. Beyond the NLO accuracy, the BFKL was considered in [32,33].
The NLO extension of the JIMWLK Hamiltonian is imperative for calculation of more general amplitudes, beyond the quark dipole of the BK equation, which determine important experimental observables like single-and double-inclusive particle production. Moreover, precision also mandates going beyond the large N c approximation of the BK equation, which is realised by the JIMWLK. The NLO JIMWLK Hamiltonian was obtained in [34][35][36] building upon the calculations of [26,37], and in [38] by a direct calculation extending [26]. A conformal part of the Hamiltonian was also deduced in [39]. The NLO JIMWLK contains information both about the NLO BKP [40,41] and on α s corrections to the triple Pomeron vertex.
While improved equations are known, the key conceptual question remains to be stability of the α s expansion. In the framework of the linear BFKL equation the NLO corrections are big [29] and even may lead to negative cross sections. Numerical calculations [42] of forward particle production in pA collisions done in the hybrid formalism of [43] and including NLO corrections [23] displayed similar pathology. The recent study of the NLO BK [44] indicates that the problem is not cured by the non-linear effects. The origin of the problem is understood reasonably well: the LO resummation picks up enhanced logarithmic contributions from the region in the gluon emission phase space, where the underlying approximations break down. These spurious contributions are subtracted in higher orders, thus leading to large negative corrections. To stabilise the expansion, one is normally forced to perform additional collinear resummations introducing kinematical constraints [45][46][47][48][49] or more involved schemes compatible with the DGLAP equation [50]. It is clear that a good understanding of a systematic approach to rectifying this problem is urgently needed for reliable phenomenological applications of the NLO BK-JIMWLK framework.
In this paper, we report our contribution to the effort in establishing NLO accuracy in saturation physics. Our results are split into two parts. In the first part, the light cone wave function (LCWF) of a fast moving hadron is computed. More precisely, in a typical Born-Oppenheimer approximation, the Hilbert space of a fast moving projectile is split into a "valence" and a "soft" sectors. The valence gluons have longitudinal momenta greater than some "cutoff" implicitly defined by the collision energy. In the LCWF, these gluons are characterised by correlators of the colour charge densities ρ a (x) (x is a transverse coordinate).
Soft gluons with momenta smaller than the cutoff do not participate in scattering. Increasing the energy of the hadron increases longitudinal momenta of the soft gluons and they emerge from below the cutoff. Having done so they contribute to physical observables, such as scattering amplitude, particle production etc. This "soft" component of the LCWF at NLO schematically has the form |ψ = (1 − g 2 s κ 0 ρρ − g 4 s (δ 1 ρρ + δ 2 ρρρ + δ 3 ρρρρ) | no sof t gluons + + ( g s κ 1 ρ + g 3 s 1 ρ + g 3 s 2 ρρ + g 3 s 3 ρρρ) | one sof t gluon + g 2 s ( 4 ρ + 5 ρρ) | two sof t gluons + g 2 s 6 ρ | quark − antiquark . (1. 2) The coefficient κ 1 is nothing else but the LO Weizsacker Williams gluon emission vertex. Below we compute all the remaining coefficients in the third order perturbation theory in which the valence charges ρ a act as background fields. A major complexity related to this calculation is that these charges are non-commuting operators on the valence Hilbert space forming a local SU (N c ) algebra. Neither commute the matrix elements of perturbation operators averaged over the soft states. As a result, we had to rederive the perturbation theory while carefully controlling over the ordering between the matrix elements. Most of the coefficients in the first line in (1.2) are computed from the normalisation of the state. Yet, δ 2 turns out to be a tricky one: it vanishes in perturbation theory. In our calculation it emerges as a phase of the wavefunction, which could be determined from a constraint on applicability of the Born-Oppenheimer approximation. This phase is somewhat similar to the Berry phase, though its nature is apparently rooted in the non-commutativity of ρ rather than being sourced by any topology.
While our calculation and presentation of the results below are fully self-contained, admittedly some of them are not new and could be located in the literature. Particularly, the entire quark sector has been worked out in [30] in a formalism that is identical to ours. The coefficients in the second line in (1.2) are related to single inclusive gluon production and we believe that it might be possible to relate 2 with partial result of [51]. The comparison is however not straightforward and we have not pursued it. We also believe that some of our results could be extracted from the LCWF of a dense projectile [52], in which α s corrections enhanced by the density ρ were resummed.
There are many advantages in having the complete LCWF at NLO. First, it provides an alternative method for derivation of the NLO JIMWLK Hamiltonian. Given the complexity of the latter, it is critical to have an independent cross check of the results before moving forward with any phenomenological applications. This constitutes the second part of our paper. Second, the LCWF formalism is possibly a good starting point to reformulate the NLO JIMWLK as stochastic evolution suitable for numeric simulations. Finally, the formalism makes it possible to simultaneously address NLO corrections to the Hamiltonian and to semi-inclusive observables, such as single inclusive production.
As was mentioned above, with increase of collision energy, the soft modes emerge above the cutoff and contribute to scattering cross sections. Thus they are naturally identified with the source of high energy evolution [53]. Let's define the projectile averagedŜ-matrix The averaging is over the soft modes only leaving Σ[ρ] as operator on the valence Hilbert space. In fact, Σ is a rapidity evolution operator  [34][35][36]. There, we only assumed the structure (1.2) without computing any of the coefficients. This structure of the LCWF helped us to identify the most general form of H N LO JIM W LK . More constraints on the form of the Hamiltonian came from the symmetries of the theory. As discussed in detail in [54], the theory must have SU L (N ) × SU R (N ) symmetry, which in QCD terms is the gauge symmetry of |in and |out states and two discrete symmetries: the charge conjugation, and another Z 2 symmetry, which in [54] was identified with signature, and can be understood as the combination of charge conjugation and time reversal symmetry [55].
The Hamiltonian in [34] was parametrised by five kernels. They were initially fixed by demanding that the evolution equations for the quark dipole and SU(3) baryon generated by this Hamiltonian match the ones known from the literature [26,37]. The Hamiltonian thus constructed was suitable for evolution of gauge invariant operators only. It was later amended in [36] to also include non-gauge invariant pieces. Interestingly, the terms in the Hamiltonian that originate from the previously discussed δ 2 were found in [35] thanks to a constraint on conformal invariance of the JIMWLK Hamiltonian in N = 4. We quote these results for the Hamiltonian in Section 2.8. It will become a benchmark for comparison with our direct calculation.
Below we compute (1.3) with the LCWF (1.2). At α 2 s order, we find terms that are proportional to δY and terms that are proportional to δY 2 . The former are identified with the NLO JIMWLK Hamiltonian (1.6) and we find a complete agreement with our previous results. The latter are naturally identified with (H LO JIM W LK ) 2 , which is a non-trivial cross check on our calculation. This paper is organised as follows. Section 2 provides an introductory overview of the wavefunctional formalism. We start with presenting the QCD Hamiltonian in the light cone gauge, introduce field quantisation and then discuss the eikonal approximation for the LC Hamiltonian and scattering matrix. To prepare the stage, we derive the LO JIMWLK Hamiltonian and also quote the previous result for the NLO JIMWLK. The LCWF is computed in Section 3. The reader who is interested in the end result only may skip most of the section and proceed directly to the final result summarised in Section 3.5. Section 4 is devoted to the calculation of Σ and extraction of the NLO Hamiltonian. Several Appendices contain complementary materials and details of the calculations.

Basics of JIMWLK
In this section, we review the main elements of high energy QCD in the light cone Hamiltonian formalism.

Light Cone QCD Hamiltonian
Our starting point is QCD Hamiltonian in light-cone gauge [56][57][58][59][60]. In light-cone coordinates, four-vectors are for transverse components. We adopt mostly negative metric [56]. The light-cone gauge: Derivation of the Hamiltonian in the light cone gauge can be found in many reviews and original publications. For self-completeness of our presentation the derivation is briefly sketched in Appendix A: where the electric and magnetic pieces have the form: The interaction Hamiltonian H int reads The flavour index is suppressed in this section.
σ i are Pauli matrices, the matrices t a denote the SU (N c ) gauge group generators in fundamental representation (the gauge group generators in adjoint representation will be denoted by T a ), which obey the algebra t a , t b = if abc t c , where f abc are structure constants of the gauge group.

Field Quantisation
Quantisation of the fields is performed in usual manner introducing creation/annihilation operators and imposing commutation (anti-commutation) relations among them. For the gauge fields: The creation and annihilation operators obey the bosonic algebra: Transforming to coordinate space, the commutation relation becomes: For the quark fields: The polarisation vectors are: The anti-commutation relations: Transforming the fields to coordinate space a la (2.8): (2.14) Inserting the field expansions (2.6) and (2.10) into (2.4), the free Hamiltonian becomes: from which the dispersion relation for free quarks and gluons is E k = k 2 2k + . The interacting part of the Hamiltonian will be quantised in subsection 2.5.

Light Cone Wave-Function of a Fast Hadron
Consider a fast moving projectile hadron at some rapidity Y 0 . Its LCWF |Ψ Y 0 is an eigenfunction of H LC QCD . It is convenient to split the modes in the LCWF in accord to their longitudinal momenta, introducing a longitudinal cutoff Λ, which is implicitly related to the energy of the target hadron (see [61] for review of the formalism). We further assume that the soft modes in the LCWF whose longitudinal momenta k + is smaller than Λ are not energetic enough to contribute significantly to the scattering process (their contribution is power suppressed by the collision energy). The hard modes with k + > Λ are also referred as valence modes.
Total energy of the collision can be increased by boosting the projectile: Y = Y 0 + δY, whereas δY is a boost parameter. By boosting the projectile, all the longitudinal momenta of partons in the LCWF get shifted up ( Fig.1): As a result, some of the soft modes that lived below the cutoff Λ get lifted above Λ and start to contribute to the scattering process. After the boost, these modes occupy a window Λ < k + < Λ e δY . The valence modes get shifted above this window.
Ignoring the modes below Λ, the LCWF |Ψ Y 0 before the boost is dominated by the valence modes only, whose distribution is approximately independent of their longitudinal momenta. (2.17) Here |v represents a valence state with no soft gluons (vaccuum of soft gluons). We assume that this LCWF is known. The LCWF after the boost is where |ψ denotes the soft part of the wave function living in the window Λ < k + < Λ e δY . The factorisation (2.18) is a typical Born-Oppenheimer approximation. Our goal below will be to find |ψ by perturbatively diagonalising H LC QCD on the Hilbert space of the soft modes, assuming the valence sector is known and frozen.

Eigenstates of Free Hamiltonian
We denote by |0 the vacuum state of the free Hamiltonian H 0 on the soft Hilbert space. The energy of this state will be set to zero, E 0 = 0. The vacuum is defined as: for any Λ < k + < e δY Λ. The soft vacuum |0 should be understood as |0 ⊗|v and, because of its valence component, it differs from the QCD vacuum. As explained above, our prime objective in the next section will be to develop a perturbation theory for this soft vacuum state. Three eigenstates of H 0 will emerge in the NLO calculation below: • One gluon state: These is a normalised state with the normalisation This state has the energy E g (k) = k 2 2k + . • Two gluon state: |0 . (2.21) with the energy E gg (k, p) ≡ k 2 2k + + p 2 2p + .

Eikonal Approximation for QCD Hamiltonian
In section 2.3, we split the Hilbert space into soft and valence modes and introduced a Born-Oppenheimer approximation for the LCWF. Accordingly, we have to split the Hamiltonian H int as defined by (2.5) into "soft" and "valence". First, we introduce field decomposition: where the underlined fields contain only modes with longitudinal momenta in the interval k + ∈ Λ, e δY Λ (soft modes) while the fields with the bar contain modes with k + > e δY Λ (valence modes). Substitution of this decomposition into H int is done in Appendix B.
Our next step is to introduce an eikonal approximation for the interaction Hamiltonian. A typical longitudinal momentum of a valence mode is much larger than that of a soft mode. Consequently, all the terms in H int that are suppressed by the ratio of soft to valence longitudinal momenta can be systematically ignored. Keeping these terms would introduce non-eikonal power suppressed corrections to our calculation. Example of terms in the Hamiltonian which are neglected 1 The action of 1 ∂ + is equivalent to division by a very large valence momentum (i.e. 1 p + → 0 for any p + > e δY Λ). In Appendix B, we first split the modes between soft and valence in H int and then employ the eikonal approximation. The result is The first eight contributions defined by (B.5) − (B.12), account for soft-soft and soft-valence interactions, while H V involves valence-valence interactions only.
Here ρ a (−p) denotes the total valence current (charge density operator) which acts as a background for the soft fields, and is defined by: and The valence current ρ a (p) satisfies the SU (N c ) algebra: where ρ a (x) is a Fourier transform of ρ a (p): As a preparation for the next section, it is useful to express (B.5) − (B.9) in terms of the creation and annihilation operators. This is done by employing the expansions (2.6) and (2.10). Only the terms, which would yield a non-trivial contribution when inserted between the three states defined in section 2.3, and the soft vacuum, will be kept. We omit further details and write the final results, (2.32) (2.34)

Eikonal Scattering
At very high collision energies, a fast projectile parton moves along its straight-line classical trajectory and the only scattering effect is the eikonal phase factor acquired along its propagation path. We assume this approximation to be valid for all partons with longitudinal momenta above Λ. TheŜ matrix operator is diagonal in the coordinate space: S ab A (x) and S αβ (x) are Wilson lines in the adjoint and fundamental representations, (2.38) The valence part scattering:Ŝ|v = |Sv . The operator ρ gets affected by the scattering as well. As explained in [53], it is possible to express the action of ρ on the state |Sv as a functional Lie derivative: which act as right and left rotations on the Wilson lines S, .

(2.40)
In a similar manner we can also express ρ qq : (2.43) The following relations will be useful for the NLO calculations: (2.44) The algebra of the operators J a L/R (x) ≡ J a L/R, F (x) + J a L/R, adj (x) follows from that of ρ: while J a L (x), J b R (y) = 0. As mentioned above, J L/R act on S as left and right rotations: These algebra will be used in section 4.9, where we subtract the (LO) 2 effects.

LO JIMWLK Hamiltonian
To warm up, we will follow the approach of [61] and derive the LO JIMWLK Hamiltonian using the LCWF formalism. At leading order, the only relevant process is one soft gluon emission by valence current, which is displayed in Fig.2. The soft part of the wave-function is given by the first order perturbation theory: 47) where N LO is determined by the normalisation condition, ψ LO |ψ LO = 1: At order g, there is only one state, |g a i (k) that contributes to (2.47). The only non-vanishing matrix-element is that of H g , (2.29): Therefore, (2.47) becomes: We will proceed by calculating the Fourier transformation of ψ LO . The coordinate space representation will make it possible to compute the scattered wave function, since, as explained in the previous section, the S-matrix is diagonal in the x-space. The charge densities ρ a (−k) are Fourier transformed according to (2.28), and the transformation to the coordinates space is done using (C.21). The wave function becomes: (2.48): The scattered wave function reads: (2.53) By inserting the wave functions before and after the scattering process in (1.3) and using the relations (2.44), we arrive at: Following the definition (1.6), the LO Hamiltonian is H LO ≡ − δΣ LO δY , which gives (see Fig.3 for schematic illustration of the LO Hamiltonian) (2.55) The LO Hamiltonian is invariant under SU L (N c ) × SU R (N c ) rotations, which reflect gauge invariance of scattering amplitudes. In addition, the Hamiltonian is invariant under Z 2 transformation, S → S † together with J L → −J R , which in [54] was identified as signature, and the charge conjugation symmetry S → S * .

NLO JIMWLK Hamiltonian
The NLO JIMWLK Hamiltonian deduced in [34] is: It is important to stress that in (2.56) all the rotation operators J L and J R are assumed not to act on S in the Hamiltonian itself. We now list the kernels: where: (2.60) The terms which depend on one variable only, either x or y, will not contribute when the Hamiltonian is taken to act on gauge invariant operators. (2.61) Explicitly: The following useful equality holds: (2.63) The kernel K JSJ reads 2 : Here µ 2 M S is the normalisation point in the M S scheme and b is the first coefficient of the An alternative representation for K JSJ : The quark sector: The kernels K JJSSJ and K JJSJ are anti-symmetric under x and y replacements, while the kernels K JSJ , K JSSJ , and K qq are symmetric. Additional properties of the kernels and relations among them can be found in Appendix D.

The Light Cone Wave Function at NLO
In this section, we compute the LCWF |ψ perturbatively in H int , up to order g 3 . Normalisation of the LCWF is done to order g 4 . Apart of the usual technical aspects, such as divergencies, related to loop calculations, there are two unfamiliar subtleties. The first one originates from the fact that the matrix elements of H int (2.5) over the soft states defined in section 2.4 are non-commuting operators on the valence Hilbert space. We thus have to revise the standard perturbation theory in order to account for non-commutativity of the matrix elements, carefully tracing their ordering. This is done in Appendix E. The second subtlety is related to the Born-Oppenheimer approximation and its applicability. It leads to a non-trivial condition on the phase of the LCWF, which otherwise cannot be determined from the perturbation theory. Summary of our results is presented in section 3.5.

Third Order Perturbation Theory
General perturbative expression for the wave function before normalisation, correct up to order g 3 , is computed in Appendix E. After simplifying notations ( n (i) −→ |i ) and also normalising the state when dividing by its norm: is a contribution to the norm from the first order perturbation theory. Up to phase, N N LO is determined from normalisation of the state: Summation over repeated indices i, j and k is over the eigenstates of free Hamiltonian H 0 , excluding the vacuum state |0 . The ones that are relevant for the NLO calculation are defined in (2.20), (2.21), and (2.22). E i denotes the energy of the respective eigenstate, and H int is defined by (2.5). (3.1) is correct even when the matrix elements are operator valued and do not commute with each other. Notice that the overall phase of the wave function is not determined by the expansion (3.2). We will discuss below how this phase can be nevertheless uniquely fixed.
Our objective is to compute ψ N LO . It is convenient to split the wave function according to its soft component content: where the subscripts g, gg and qq correspond to the soft particle content in each state. ψ LO g ρ is the LO contribution (see section 2.7): The following states (integration over the momenta and summation over the quantum numbers of each particle are assumed) contribute and will be computed one by one below. Order g 2 states: Order g 3 states: The following contributions vanish:

Matrix Elements
We write down all the matrix elements that are relevant for the NLO calculation. These expressions are computed based on (2.29) − (2.34). At NLO there are two matrix elements for the soft-soft interactions and four matrix elements for the soft-valence interactions: • Gluon splits into quark and anti-quark pair • Valence adds an extra gluon to an existing one • Valence instantaneously creates a quark and anti-quark pair • Valence instantaneously creates two gluons • Valence instantaneously interacts with a soft gluon

Technical Aspects of the Calculation
Before diving into the calculation of the LCWF (3.1) in the next subsection, let us first introduce several important technical aspect, which are used throughout this calculation.

• The Decomposition Procedure
As explained in the Introduction, the result of the calculation of the LCWF will have the form (1.2). Many of the terms we are to compute will have a product of more than one ρ operator. These operators do not commute and their mutual ordering is very important. (3.1) will be later used to compute the NLO JIMWLK Hamiltonian. In order to compare with section 2.8, we would have to represent our result as in (2.56), including the symmetry properties for the kernels (the ones mentioned at the end of section 2.8). These symmetry properties do not come automatically in our calculation and we would have to manipulate our final results. It turns out that in order to achieve this, it suffices to symmetrise any product of two ρ operators in the LCWF via the following identity: The replacement (3.27) will be referred to as "decomposition procedure". The physics behind this decomposition is that it explicitly separates between contributions from two emitters being in different transverse points in coordinate space from that when both emitters are in the same point. The latter is then equivalent to a contribution from a single emitter, in accordance to the colour algebra. After implementing the decomposition procedure where applicable, each component of the LCWF will receive an extra label in accord to the number of ρ operators it contains. For example, in section 3.4.1 we show that ψ 1 qq contains only one ρ operator, and this component will be denoted thereafter by ψ 1 qq ρ . •

Regularisation of Divergent Integrals
There are two types of divergencies that appear throughout the calculation. These are the usual divergencies in loop calculations. First, we will face divergencies that come from the longitudinal phase space integrals of the type dk + /k + , which correspond to emission of soft partons, just as in the case of LO calculation. These divergencies are regularised by confining the soft modes to live within the window Λ < k + < Λe δY . When computing the Hamiltonian in the next section, we will also encounter this longitudinal divergence "squared". Physically, such double divergence corresponds to independent (uncorrelated) emission of two soft gluons, which is accounted for by iteration of the LO Hamiltonian.
Λ acts as an IR cutoff which will be set to zero (and δY to infinity) whenever it is safe to do so, without stating it explicitly.
UV divergences in transverse space integrals will be treated via dimensional regularisation, switching to d ≡ 2 − dimensions, (3.28) In addition, we will be frequently using the identity: A couple of relevant integrals in dimensional regularisation are quoted in Appendix C. As usual, (C.33) will be used for expansion. We will retain the divergent 1 terms at the level of the wave function, as by itself the wave function is not a physical observable. Coupling constant renormalisation is done for physical observables, such as Σ (or the JIMWLK Hamiltonian). This will be implemented in the M S scheme as a default action without any further explanations. The calculation also involves overlapping singularities when both the longitudinal and transverse integrals diverge. Those are treated similarly as above.
In addition to the above discussed UV divergencies, there are also potential IR divergencies at small transverse momenta. Those, however, do not emerge in our calculation explicitly because they always get regularised by the charge density operators ρ. Yet, this depends on a non-specified IR behaviour of the valence degrees of freedom and a priori is not guaranteed to lead to a safe IR. Particularly the IR divergencies resurface when the resulting JIMWLK Hamiltonian is taken to act on a gauge non-invariant operator, such as in the problem of gluon reggeization. The gluon trajectory is then found to be IR divergent even at LO.
• Phase of the LCWF The perturbation theory cannot fix the phase of the LCWF, that is the phase of N N LO . Normally, the phase is unobservable as it cancels in physical observables. However, this is not the case at hand, as we are not computing the full wave function but only its soft component. The normalisation coefficient N N LO is operator valued on the valence Hilbert space and thus its phase can make a non-trivial contribution to the observables.
How can we find this phase which is clearly beyond the perturbation theory? We recall that our formalism is based on the Born-Oppenheimer adiabatic approximation 3 . Obviously, there is some freedom in the phase, as, in principle, it can be always absorbed into the valence part of the LCWF. Yet, the factorisation (2.18) assumes that where H V is a Hamiltonian for the valence degrees of freedom. The second expectation value is taken over the soft modes only. The condition (3.30) means that the dynamics of the soft modes does not significantly affect that of the valence. This assumption is implicit in the way the high energy factorisation is set, eliminating the freedom in the phase.
While the explicit form of H V is not known, the condition (3.30) presumably should be valid for any operator on the valence Hilbert space. Particularly, H V can be thought as been constructed from some powers of ρ, with (3.30) implying that they commute with |ψ on the valence Hilbert space. Commuting one ρ with |ψ is equivalent to taking derivative δ/δρ. We thus require the following condition to hold Notice that (3.31) is the Berry connection, which we require to vanish. It is easy to check that at LO the condition is fulfilled by |ψ LO automatically. At NLO, this condition becomes a non-trivial constraint on the phase, which is worked out in Appendix F.

Computation of the NLO Wave Function
In this section we compute the wave function (3.1). For each of the contributions (3.5) -(3.17), a separate subsection is devoted. Lengthy computations are available in Appendix G. The reader who is not interested in these details can proceed directly to the final results in section 3.5.

Quark anti−Quark State
Computation of ψ 1 qq ψ 1 qq ρ is defined in (3.5) (Fig. 4a). . (3.32) The physics content of (3.32) becomes clear by reading the matrix elements from the right to left: the first one corresponds to initial emission of a soft gluon from the valence current; then this soft gluon subsequently splits to a quark and anti-quark pair. The momentum conservation δ (3) (q − k + p), which is explicit in the second matrix element (3.21), removes q integral and effectively turns it into k − p. The process is illustrated in Fig. 4a.
After inserting the relevant matrix elements, (2.49) and (3.21), and integrating over q: The variable ξ is defined as: This definition will be used repeatedly throughout this paper.
Here 1 2 is to avoid double-counting in two-gluon states. Inserting the matrix elements, (2.49) and (3.23): The physics content of (3.37) is quite obvious and corresponds to subsequent emissions of two soft gluons. The momentum conservation imposed via the delta functions turns the momentum q into either k or k − p, which is illustrated in Fig. 5a. After integration over q, (3.38) becomes: (3.39) Here ξ ≡ p + k + as already defined in (3.34). The result (3.39) can be decomposed according to (3.27), and represented as contributions with one and two ρ operators: and Computation of ψ 2 gg ψ 2 gg is defined in (3.8) (Fig. 5b). .
The condition q = k − p is enforced via momentum conservation, which is reflected on Fig.  5b. By inserting the relevant matrix elements, (2.49) and (3.22), and integrating over q: .
After inserting the relevant matrix element (3.24) and changing variables according to (3.34) and (3.45), we can write the last result as: While focusing on two-parton production, we have not encountered any divergencies. This is obviously because we have not yet ran into any loop integrations. At this stage we can only anticipate their appearance in the next section, where the overlaps of two wavefunctions are computed. Particularly, it is pretty obvious that ψ 1 qq will lead to the logarithmic singularity in longitudinal momenta integration squared, which corresponds to second iteration of the LO JIMWLK Hamiltonian. Somewhat less obvious, but ψ 3 qq also contributes to square of the LO. This is due to the kinematical region where one of the two gluons is much softer than another. In the next subsection, we will encounter both longitudinal and UV divergencies in the LCWF itself. The divergencies will be present in most of the loop corrections to the single gluon state.
. The physics of (3.49) is transparent from Fig. 6a. The figure reflects the momentum conservation, which eventually imposes q = k − p and r = k. The rest of the calculation is available in Appendix G.1. At the end of calculation there, we arrive at the following result: where we have introduced µ 2 M S ≡ 4πe −γ µ 2 . 2/ reflects dimensional regularisation used to regulate the usual UV divergence of the quark loop integral. The renormalisation procedure will be applied in the next section, where we compute the physical Hamiltonian: the 2/ singularity will be then absorbed into redefinition of the coupling constant.
Just as in the previous calculation and as demonstrated on Fig. 6b, both q and r become q = k − p and r = k after realisation of the momentum conservation imposing delta functions. The rest of the calculation is available in Appendix G.2. The result is: Apart the terms which are easily recognised as the gluon loop contribution to the βfunction, we notice additional terms and particularly the ln Λ/k + , which are potentially divergent in the Λ → 0 limit. These contributions arise from the phase space region where one of the gluons in the loop is much softer than another. This is precisely the kinematical domain where the triple gluon vertex H ggg admits eikonal approximation. It thus becomes clear that these logarithmic contributions correspond to iteration of the LO emission process. Indeed, in the next section we will demonstrate that all the ln Λ/k + terms eventually contribute δY 2 to the evolution identified with second iteration of the LO JIMWLK Hamiltonian. . (3.53) Due to the instantaneous interaction, the gluon that was originally emitted can either gain or lose longitudinal momentum. The longitudinal integration over p + in (3.53) can be split into two intervals: Accordingly, we also write ψ 3 g as a sum of two contributions ψ 3d g and ψ 3u g , ψ 3 and illustrated in Figure 7a, where our standard relation p + = ξ k + is also used. The rest of the calculation is available in Appendix G.3. The end result integrated over ξ with the aid of (C.7) reads: (3.56) We seem to again encounter the ln(k + /Λ) type divergence. The latter is, however, fictitious and does not represent any physics. The divergence will cancel when we collect all similar contributions.
Next is ψ 3u g . We find it convenient to exchange the dummy momenta k and p, k ↔ p. The expression for ψ 3u g is then given by . (3.58) Introducing ξ for compactness the last result is: The contributions ψ 4 g and ψ 5 g have very similar integrand structures. This is the reason we founnd it convenient to combine the two. ψ 4 g is defined in (3.13) (Fig. 7c), . (3.60) The matrix element g d l (p) |H g | g b j (q) g a i (k − q) contains a momentum conserving delta function which eventually sets q = p (or q = k − p). Fig. 7c illustrates the process. ψ 4 g is computed in Appendix G.4. ψ 5 g is defined in (3.14) (Fig. 8a), . (3.61) The momentum conservation in ψ 5 g imposes q = p and r = k − p (or q = k − p and r = p), which is the case displayed in Fig. 8a. The calculation of ψ 5 g can be found in Appendix G.4.
We now introduce sum of ψ 4 g and ψ 5 g , Next, ψ 4+5 g is split according to (3.27) into two parts involving one and two ρ operators: All the algebra is performed in Appendix G.4. Below, we quote the final results only.
• One ρ part (3.64) One of the longitudinal integrals, either ξ or k + , in (3.64) could be computed, but we have again chosen to leave the expressions as is until the final result is assembled. Yet, we see that the longitudinal integral would result in a logarithmic divergence, although as usual regularised by Λ. The origin of this divergence can be traced to ψ 5 g and particularly to the phase space region where H ggg admits eikonal approximation. Thus ψ 5 g has a piece, which can be interpreted as a soft gluon emission from a LO evolved LCWF. In the next section, such contributions will be identified with (LO) 2 terms.
Furthermore, 2/ points to a regularised UV divergence. Several diagrams below will also contribute such terms. Most, however, will cancel against each other. Some UV divergence will nevertheless survive and that is a feature of the LCWF at NLO. It will not, however, contribute any UV divergence at the level of the Hamiltonian. There they will cancel against a subtraction term made of two produced gluons, which cross the shock wave at the same transverse point.
• Two ρ part From (G.26) (after the substitution p ↔ k): A convenient way to represent (3.65): The advantage of the last representation is obvious: in the next subsection we will add up ψ 4+5 g ρρ with ψ 3u g ρρ and the mutual cancelation becomes transparent via representation (3.66).
Computation of ψ 6 g ψ 6 g is defined in (3.15) (Fig. 8b). . (3.67) Just as in the previous calculations, the momentum conservation imposes q = k and r = k − p. This is how this process is illustrated in Fig. 8b with our usual definition of ξ. The rest of the calculation is available in Appendix G.5, from which we quote: (3.68) The contributions with one and two ρ operators can be isolated according to (3.27): (3.70) • Two ρ part .
(3.72) ψ 7 g accounts for two distinct processes depicted, after the dummy momenta are integrated over, in Fig.9a and Fig.9b. They correspond to two different light cone orderings in gluon emission/absorption process. Intuitively, we expect ψ 7 g to mostly contribute to (LO) 2 effects. Indeed the structure of singularities, see below, is consistent with our intuition. Particularly, we are to see a logarithmic divergence in the longitudinal integration, which can be identified as a virtual contribution at the LO, followed by emission of another soft gluon. The calculation is available in Appendix G.6, from which we quote the result: which can be written as a sum: (3.75) • Two ρ part This part can be directly deduced from (3.73): (3.76) Figure 9. One gluon production diagrams with contribution to three ρ. (a,b)ψ 7 g , (c)ψ 8 g .

(3.77)
Another useful representation is obtained by changing p → k − p in the second term, (3.78) • Three ρ part This part can be directly deduced from (3.73): After integration over ξ: (3.80) From the general structure of the NLO JIMWLK Hamiltonian, we know that ψ 7 g ρρρ is not expected to contribute to it. At the level of Hamiltonian, ψ 7 g ρρρ gives rise to ρ 4 terms, which are not present in the NLO JIMWLK Hamiltonian. Consistency demands that all such terms are expected to be absorbed into (LO) 2 Hamiltonian. After a trivial cancelation against ψ 8 g , which is to be computed next, a combined contribution is indeed found of the right form.

The Final Result
In this section we will assemble together all the different contributions that were computed in section 3.4. It turns out to be useful to replace the representation of the LCFW as in (??) by a new representation schematically outlined in (1.2), in which the classification of the contributions will be made according to the soft particle content (qq, gg or g) and the number of valence current operators (ρ, ρρ or ρρρ). This information will be indicated by a subscript. For example, the state |ψ gg ρρ contains two soft gluons and two ρ operators. The NLO wave function is decomposed as follows: ψ N LO = N N LO |0 + ψ LO g ρ + |ψ qq ρ + |ψ gg ρ + |ψ gg ρρ + |ψ g ρ + |ψ g ρρ + |ψ g ρρρ . (3.84) ψ LO g ρ was defined in (3.4). Additional states are:

(3.95)
After transformation to coordinate space using (C.21): where the parameter b was defined in (2.65) and µ 2 M S ≡ 4πe −γ µ 2 . The transformation to coordinate space is possible using (C.21) and (C.22): In section 4, we will use this form of |ψ g ρ . However, after defining a new variablek + = ξk + , one of the longitudinal integrals can be evaluated. For completeness of the presentation, we also quote this result: |ψ g ρ appears to be one of the most non-trivial parts of our calculation, which to our knowledge has not been computed before. It has a mixture of various terms, including singular ones. It is hard to give interpretation to every term at the level of the LCWF. We only comment here that |ψ g ρ contributes both to NLO JIMWLK Hamiltonian (the kernel K JSJ ) and second iteration of the LO one. |ψ g ρ is also a source of the β-function of the running coupling.
• The one gluon state with two ρ By inserting (3.56), (3.59), (3.66), (3.71), and (3.76) to (3.88) we arrive at: Again, despite the possibility to evaluate one of the integrals, in section 4, we will use this form of |ψ g ρρ . For completeness of presentation, we also quote the final result: We did not manage to compute a coordinate representation for this component. Together with |ψ g ρ , |ψ g ρρ constitutes one of the main new non-trivial result. In the next section, we will demonstrate how |ψ g ρρ contributes both to NLO JIMWLK Hamiltonian and the LO iteration. Particularly |ψ g ρρ will give rise to the NLO kernel K JJSJ . The kernel K JJSJ is related to the kernel K JJSSJ through the equation (D.4). In its turn, the kernel K JJSSJ will be obtained from |ψ gg ρρ (3.42). We notice that while there is a relation between the kernels, it does not seem to be possible to establish any interesting relation at the level of the LCWF components.
• The one gluon state with three ρ Substituting (3.80) and (3.83) in (3.89), we arrive at: Inserting the identity d 2 q p·q q 2 δ (2) (q−p) = 1 in (3.102), the latter can be straightforwardly written in coordinate space: (3.103) The |ψ g ρρρ component of the LCWF is pretty transparent. It is already proportional to δY because of the virtual soft gluon, which we have integrated over. One can recognise in the integrand of (3.103) the kernel of the LO JIMWLK Hamiltonian. The extra soft gluon g(k) in |ψ g ρρρ corresponds to another iteration of the kernel as will be demonstrated in the next section.

• Normalisation
The norm of the wave function can be represented as: From normalisation condition on the expansion (3.84), The computation relevant for N N LO is postponed until section 4.7. Here we quote the final result: (3.106) where K JSJ (x, y, z) is defined in (2.66). The first non-trivial term is a LO contribution, which is of order α s . The second term is an α 2 s correction and, because of its linear dependence on δY , it will be contributing to virtual terms in the NLO JIMWLK Hamiltonian. The last term will contribute to subtraction of (LO) 2 terms in the evolution.
The phase of the LCWF is computed in Appendix F from the condition (3.31): where K JJSSJ is defined in (2.57). Up to order α 2 s ,

The Next to Leading Order JIMWLK Hamiltonian
The aim of this section is to compute Σ (1.3), that is theŜ-matrix expectation value in the LCWF (3.84), which was calculated in the previous section. Each soft quark or gluon in the LCWF will contribute one eikonal factor S defined respectively in (2.38) and (2.37). The ρ factors in the ket/bra states will be converted into left/right rotation operators J L/R according to (2.39). There are nine overlaps that emerge in this computation, which we classify according to the number of J operators and the Wilson lines S that they contain: Σ LO was defined in (2.54). The other entrees (all of order α 2 s ) are defined below: Σ JJSSJ ≡ ψ gg ρ |Ŝ |ψ gg ρρ + ψ gg ρρ |Ŝ |ψ gg ρ , (4.3) For each of these contributions a dedicated subsection will be devoted. Some details of the calculations are provided in Appendix H.
As anticipated, we will find not only contributions linear in δY, but also contributions of order (δY) 2 . To make separation between these very different types of contributions crystal clear, we will introduce an additional label (upper script) − NLO or (δY) 2 respectively. Each Σ ... defined in (4.2) − (4.9) is decomposed according to: At the end, in section 4.8, we will assemble together all the contributions at NLO and write down the NLO JIMWLK Hamiltonian defined in (1.6).
Finally, in section 4.9 we will collect all the (δY) 2 contributions, and show that they correspond to contributions which are generated from the LO Hamiltonian applied twice. Based on the expansion (1.4), we will show that the anticipated relation: holds indeed. This relation serves as a strong consistency check on our calculation.

Computation of Σ qq
Σ qq is defined in (4.2) as an overlap between the incoming state |ψ qq ρ (3.90), before and after passing through the shockwave (Fig.10a).
Representing the operators ρ in terms of the operators J L and J R as explained in section 2.6, and replacing the products of Pauli matrices using the identity σ i σ j = iε ijk σ k + δ ij I, We have used the fact that the scattering matrix is diagonal in coordinate space: Here ξ ≡ p + k + and ϑ ≡ v + u + . Equation ln X 2 (X ) 2 (4.14) The result contains terms of order δY only. To indicate this, we add the upper script N LO, so that Σ N LO qq ≡ Σ qq . (4.14) can be equivalently written as: where K qq (x, y, z, z ) is defined in (2.67).

Computation of Σ JSSJ
Σ JSSJ is defined in (4.4) as an overlap of |ψ gg ρ , (3.94), before and after passing through the shockwave. The relevant diagrams appears in Fig.10b.
where the following matrix element was used: The remaining computation is moved to Appendix H.2. The result is: where K JSSJ (x, y, z, z ) is defined in (2.59), and (4.20) Figure 11. The diagrams for Σ JJSSJ . Effective two gluon emission vertex corresponds to |ψ gg ρ .

Computation of Σ JJSSJ
Σ JJSSJ is defined in (4.3) as an S-matrix overlap between |ψ gg ρ (3.94) and |ψ gg ρρ (3.96). The relevant diagrams appear in Fig.11. With the aid of (4.17), w, x, y, z, z Notice that the minus sign between the terms in the last line in (4.21) is due to complex conjugation. The integration over ξ is done using (C.2), (C.5) and (C.6): w, x, y, z, z It is important to notice that J L and J R in (4.22) by construction do not act on S A inside (4.22). As anticipated, the result contains both terms proportional to δY and (δY) 2 , which are split as: JJSSJ emerges from the term ln Λ k + , which becomes proportional to (δY) 2 after integration over k + . Σ N LO JJSSJ reads: where K JJSSJ (w, x, y, z, z ) is defined in (2.57). (4.25)
After the ξ integral is evaluated using (C. 18) and (C.19): (4.27) As usual, this expression can be split according to: The remaining computation is available in Appendix H.3. The results: (4.29) With K JJSJ (w, x, y, z) defined in (2.58). In addition:

Computation of Σ JSJ
Σ JSJ is defined in (4.6) as an S-matrix element between the state |ψ g ρ , (3.98), and the state ψ LO g ρ , (3.4). The relevant diagrams appears in Fig.13. Σ JSJ = − g 4 64π 5 x, y, z Using (C.11), (C.13) and (C.15) we can integrate over ξ: and obtain: x, y, z Now, after integration over k + , Σ JSJ is split according to: with where The result for K JSJ can be identified (up to 2γ − ln 4) in equation (87) of [26]. The relation between K JSJ (x, y, z) and K JSJ (x, y, z) is given below, in (4.58). Also, which can be equivalently written as: (4.38)

Computation of Σ JJSSJJ and Σ JJJSJ
Σ JJSSJJ is defined in (4.7). This contribution is driven by |ψ gg ρρ , (3.96), interfering with itself before and after passing though the shockwave, as shown in Fig.14.   whereȲ = y −z andV = v −z . The rest of this calculation is available in Appendix H.4. The result after integration over k + is: Now let us compute Σ JJJSJ as defined in (4.8). The contribution comes from the S-matrix element between the state |ψ g ρρρ (3.103) and the state ψ LO g ρ (3.4), as shown in Fig.15. Which we write as:

Virtual Contributions
Σ virtual is defined in (4.9). This contribution comes from the scattering of the vacuum component N N LO |0 of the LCWF. To save space, we will not illustrate this contribution, but, quite obviously, these are all the graphs with no soft partons crossing the shockwave, similarly to Fig. (3a,3b). With (3.108) and (3.105) in mind, and counting powers of g, To organise the computation, we introduce the following decomposition for Σ virtual , where, after inserting ||N N LO ||, as defined in equation (3.105), into (4.43), (4.45) The contribution with two J operators The contribution with three J operators 0| ( ψ gg ρρ |ψ gg ρ + ψ gg ρ |ψ gg ρρ )Ŝ |0 − 1 2 0|Ŝ ( ψ gg ρρ |ψ gg ρ + ψ gg ρ |ψ gg ρρ ) |0 , (4.52) The contribution with four J operators (4.53) Finally, In the rest of this section we present the results for all the virtual terms defined in (4.47)−(4.54). In fact, it is straightforward to deduce most of them by recycling our previous computations. We notice that most of the overlaps that appear in (4.47)−(4.53) could be read off from the corresponding real terms, when the limit S −→ 1 is taken. In addition, this limit has to be taken with understanding that all the factors ρ in (4.47)−(4.53) (apart of the last line in (4.53)) are placed either on the left or right ofŜ. This basically implies a replacement in the real terms J L −→ J R (J R −→ J L ) in the caseŜ is placed to the right (left) of the matrix elements. One should, however, be careful with the ordering of the operators. The correct procedure involves first writing all the J R operators to the left of J L operators. Second is to convert the J operators back into ρ in accord to (2.44). Finally, the string of ρ gets converted again to a string of J operators.

The NLO JIMWLK Hamiltonian Assembled
We define Σ N LO as a sum of all the contributions of order δY:  [26], we subtract (and later add) terms which correspond to the case in which both gluons or the quark and anti-quark pair pass through the shockwave at the same transverse point. This corresponds to single gluon scattering, thanks to the identities S ab A (z) = 2tr S † (z)t a S(z)t b and f abc f def S be A (z) S cf A (z) = N c S ad A (z). Accordingly, we introduce the subtraction terms As was mentioned in the previous subsection, proper dimensional regularisation of z integrals is needed, and the results for the integrations are quoted in (D.5) and (D.6). With the subtraction terms added, Following (4.58), The UV divergence is canceled in K JSJ . The NLO JIMWLK Hamiltonian is obtained from (4.73) via the definition (1.6). The result is identical to the one quoted in (2.56).

Reduction of the (LO) 2 Contribution
In addition to the contributions to the NLO JIMWLK Hamiltonian we found terms proportional to (δY) 2 : Our objective in this section is to show that all these contributions are generated by second iteration of the LO Hamiltonian H LO JIM W LK . Let's first summarise all the (δY) 2 contributions that were found in subsections 4.3 -4.7. From (4.69), (4.76) Adding together the contribution with one S factor, (4.30), (4.38) and (4.42): Adding together the contributions which contain two S factors, (4.25), (4.20), and (4.40): Now let us write down the result that is obtained by applying H LO twice:   The QCD Lagrangian: The tensor F µν defined as F µν ≡ t a F a µν = − i g [D µ , D v ], satisfies the following equations of motion: The canonical momenta: Neglecting the quark masses, the Dirac equation becomes (after multiplication by γ 0 ): Two-coulmn vector ψ + = ψ 2 ψ 3 and ψ − = ψ 1 ψ 4 are defined by keeping only the nonvanishing components. The equations of motion for ψ ± can be found after multiplying equation (A.4) by the projections Λ ± : from which ψ − is determined in terms of ψ + . The following constraints are obtained from the + component of (A.2): Eliminating ψ − and A + , the LC Hamiltonian is: (2.4) and (2.5) are obtained using the following relations:

B Deriving the Eikonal Approximation for LC Hamiltonian
In order to introduce the eikonal approximation for the LC QCD Hamiltonian, we split the gluon and quark fields into the soft and valence modes: Similarly for the quark field: to H int defined in (2.5), we define the following components of the interaction Hamiltonian: (B.12) The eikonal approximation means that all contributions involving division by large (valence) longitudinal momenta are neglected, such as 1 Furthermore, we only keep terms which are relevant for our NLO calculation, the ones that contain either interaction between soft modes or interaction between soft and valence modes. Some soft interaction terms that would contribute starting from NNLO only, such as instantaneous gluon-quark interaction term, as well as the valence Hamiltonian are not written down explicitly.

C Integrals and Fourier Transformations
In the following Appendix we assemble all the integrals that were used throughout our calculations in this paper.
The integrals below are evaluated under the assumption Λ → 0 and δY → ∞. While explicitly divergent terms will be retain, we will take the limit in all the non-divergent terms without further notification. (C.10)

D Properties of the NLO JIMWLK Kernels
Where n (i) are orthogonal states, corresponding to i-th order in perturbation theory. By inserting (E.2) in (E.1), expanding and equating terms in the same order of g, we arrive at:

F The Phase of the Wave Function
The condition (3.31) imposed in order to find the phase reads It is trivial to check that the LO LCWF satisfies (F.1). This happens automatically thanks to cancelation between contributions of normalisation against diagonal contribution from one gluon component. Similar cancelations happen when ψ N LO (3.84) is inserted in (F.1).
It is easy to see that the following ansatz for iφ N LO solves (F. 3), G Supplementary for Section 3 G.1 Supplement for Computation of ψ 1 g ψ 1 g is defined in (3.49). The calculation in this section is identical to [30]. After inserting the relevant matrix elements, (2.49) and (3.21): By using After change of variables according to (3.34) and (3.45): The UV divergent p-integral is regularised via dimensional regularisation according to (3.28). Then, using (3.29) together with the following identity: we arrive at: Integration over p is done with the aid of (C.31): (G.7) After expanding with (C.33) and taking → 0 limit, the result becomes: Finally, integrating over ξ according to (C.1): which we can write as in (3.50).
G.2 Supplement for Computation of ψ 2 g ψ 2 g is defined in (3.51). After inserting the matrix elements, (2.49) and (3.22), By changing variables according to (3.34) and (3.45) we can write the last result as: (G.11) After some algebra: By replacing the measure according to (3.28): We can now compute the integration over p with the aid of the integral (C.31): After expanding with the aid of (C.33) taking the → 0 limit and, the last result becomes: The relevant integral for the integration over ξ is (C.10), the result after the integration is: G.3 Supplement for Computation of ψ 3 g ψ 3 g is defined in (3.55). After inserting the matrix elements, (2.49) and (3.26): (G.17) By changing variables according to (3.34), we arrive at:

(G.19)
The relevant integral for the integration over ξ is (C.7), the result after the integration appears in (3.56).
• One ρ part Let us now focus on the contribution which involves only one ρ operator: We rewrite the denominator of ψ 4+5 g by using Feynman parameter. We introduce the following definitions: (G.28) After integration over ξ with the aid of (C. 13) and (C.14), the last result becomes: , ρ a (−k + p) (G. 45) Notice that under the change p → k − p the second summand inside the rectangled brackets becomes exactly the same as the fourth summand. We write the result which is obtained after this change in (3.71).
After simplification of the last expression we arrive at: (G.47) By using the algebra of ρ operators, as in (2.27), we can write the last expression as: We can now perform Fourier transformation by using (C.23) − (C.24). It is possible to simplify the denominators using: Then, we arrive at:

(H.3)
Rewriting the scalar products: we arrive at: It is possible to integrate over ξ by using the integrals (C.2), (C.3), and (C.4), the result is shown in (4.15).

H.4 Supplement for Computation of Σ JJSSJJ
From (4.39) we deduce that: (H.20) Integrating over ξ and reordering the J operators: (H.21) The result after integration over k + is shown in equation (4.40).