NLO JIMWLK evolution with massive quarks

NLO evolution of the Jalilian-Marian-Iancu-McLerran-Weigert-Leonidov-Kovner (JIMWLK) equation with massless quarks was derived a few years ago. We make a step further to compute the evolution kernels focusing on the effects due to finite quark masses. To this goal, the light-cone wave function of a fast moving dilute hadronic projectile is computed up to ${\cal O}(g^3)$ in QCD coupling constant. Compared with the massless case, a new IR divergence emerges, which is eventually canceled by a mass dependent counter term. Our results extend the theoretical tools used in physics of gluon saturation and aim at improving precision in future phenomenological applications.


Introduction
High energy evolution of a hadronic observable O in a dilute-on-dense scattering (DIS-like processes) is governed by the JIMWLK equation of the forṁ direction over the last fifteen years or so. Particularly, NLO formulation of both the BK equation [10] and JIMWLK Hamiltonian became available [11,12] starting the era of NLO phenomenology [13,14]. It was discovered that at NLO, there emerge conceptual problems related to large energy independent new logarithmic corrections, threatening stability of the α s expansion. A lot of efforts have been invested over the years attempting to fix these problems [15][16][17][18][19][20][21][22]. While overcoming them is critical in order to gain full control over the calculations and achieve the desired NLO accuracy, these research directions are in a sense orthogonal to our goals below and hence would not be reviewed in any detail here. The NLO JIMWLK Hamiltonian has been derived for massless quarks only in previous literature, which at NLO contribute as loop corrections. Yet, it is known that quark masses could be significant phenomenologically, especially at realistic (not asymptotically large) energies. Say, at HERA, at relatively moderate Bjorken x ∼ 10 −2 − 10 −3 , the charm quarks contribute to F 2 about 10% [23]. Notice that this is precisely the kinematical regime of future EIC operation. Inclusion of quark mass effects into evolution is an important additional step towards improving precision of the phenomenological tools. By this publication we are filling this gap. Our work is complimentary to the most recent one [24], which computed the mass effects in the photon impact factors.
In the next Section, we review the basics of the JIMWLK, particularly the Light Cone Wave Function (LCWF) formalism, which is central to our approach. The derivation will closely follow Ref. [12]. Section 3 presents computation of the LCWF with massive quarks and derives the quark mass-modified kernels of the JIMWLK Hamiltonian. The last Section (Section 4) provides some concluding remarks. Several Appendices supplement our calculations.
2 From LCWF to JIMWLK -a brief review.
In this section, we revisit the key elements needed for the derivation of the NLO JIMWLK equation with massive quarks. We use the light cone quantization and the conventions used here closely follow those of Ref. [12].
In the LC gauge the LC Hamiltonian with massive quarks reads [25] for the kinetic part and the interaction part is [25]: where α i ≡ γ 0 γ i and β ≡ γ 0 . The mass dependent terms in the LC Hamiltonian include the fermion kinetic term in Eq. (2.1) and the three point (two fermions and one gluon) interaction terms in the last two lines of Eq. (2.2). Four-point interaction terms are mass independent. While dealing with massive quarks, it is convenient to operate with four component spinors for fermion fields (note that this convention is different from Ref. [12]). The final results are, of course, independent of the conventions. For 4-component spinors, we use the following notation which satisfy In the Dirac representation, the 4-component fermion spinors have the following simple projections: Quantization of the fermion fields proceeds as usual by introducing the Fourier mode expansion where α is the fermion color index, the 1/ √ 2 factor comes from (using Eq. (2.5)) Creation and annihilation operators satisfy the anti-commutation relations with all other anti-commutation relations vanishing.
The gluon fields Fourier expand as The commutation relations for creation and annihilation operators The following normalization conventions for free (anti-)quark and gluon states are used Conventions for Fourier transformations which also automatically fix the convention for free states in (transverse) coordinate space, In this work, when loop correction are being computed, dimensional regularization with modified minimal subtraction scheme (MS scheme) is used, i.e., whereμ is short for the subtraction scale µ MS in the MS scheme, and γ E is the Euler-Mascheroni constant.

LCWF of a fast moving projectile
LCWF approach provides a simple way of deriving the JIMWLK equation (see Ref. [26] for an introduction to the formalism). Following the Born-Oppenheimer approximation, the longitudinal modes of a fast moving dilute hadron (projectile) are divided in longitudinal momentum into valence and soft modes, with a separation scale Λ, which implicitly depends on the collision target. Below the scale Λ, the modes are too soft to contribute to the Figure 1. The longitudinal momenta of the fast moving hadrons (projectiles) are divided into valence modes and soft modes with a separation scale Λ. After a differential boost of δY , the soft modes within the range δY right below the cutoff Λ will emerge above Λ.
scattering amplitude. Energy of the collision can be increased by boosting the projectile. How much a parton in a projectile is boosted is characterized by the boost parameter Y (or rapidity), which is defined in the light-cone momentum as where m is the mass of the parton, and k ⊥ is the transverse momentum of the parton (or k in our convention setup above). Since the rapidity is additive under a series of boosts, we have k + → k + = e δY k + under a boost of δY . After a boost of δY , the soft modes within the range of δY right below Λ emerge above Λ and start to contribute to the scattering amplitude. This is the essence of the rapidity evolution of the LCWF. Assume that the projectile has already been boosted to a rapidity Y 0 with the quantum state described by the valence modes only After the boost, the soft modes emerge above Λ, resulting in a new LCWF of the projectile, where |ψ refers to the soft modes of the state, whose longitudinal momenta lie in the window (Λ, e δY Λ), as shown in Fig. 1. Note that ⊗ here is symbolic since |ψ depends on the valence state |v (soft and valence modes are entangled in a complex way). Within the Born-Oppenheimer approximation, the valence part of the wave function in Eq. (2.20) is assumed to be frozen while the soft part, |ψ , can be calculated perturbatively, order by order in QCD coupling g. The valence modes manifest themselves as a non-Abelian (noncommutative) background charge distribution, which is intrinsically non-perturbative. The quark and gluon fields are decomposed into valence and soft modes where ψ and A are the soft modes, and ψ and A denote the valence modes. Substituting this decomposition into the Hamiltonian (2.2), one obtains an effective Hamiltonian for the soft-modes, including their self-interaction, and interaction of the soft modes with the valence [12]. The latter are treated as classical background fields. The effective soft Hamiltonian can be further simplified due to eikonal approximation: power suppressed terms (or higher order terms) will be ignored, where the power counting parameter will be a small parameter of the order λ ∼ p + soft /p + valence . For instance, where the former is counted as λ 0 and the latter as λ −1 . For the former, 1 ∂ + is counted as 1/p + valence and ∂ + A c j as p + valence , so the term counts as λ 0 . For the latter, 1 ∂ + is counted as p + soft (the sum of two valence momenta could be soft) and ∂ + A c j as p + valence , so the term is counted as λ −1 . Below we will only focus on the terms in the effective Hamiltonian which are relevant for the calculation of the quark mass effects. These are which are in fact mass independent and identical to those used in Ref. [12]. The first term is the leading soft gluon-valence interaction, responsible for the LO JIMWLK evolution. The second term is the instantaneous interaction of the quarks with the valence sector. It becomes important for quark bubble resummation. The only quark mass dependent interaction is the gluon-quark-quark term (2.24) From Eq. (2.24), one obtains H gqq as a second quantized operator on the soft Hilbert space At NLO, a pair of soft quarks is emitted from a soft gluon with the mass dependent transition matrix element Additional matrix elements, which are involved in the calculation below, are all mass independent and are induced by the interactions (2.22) and (2.23) [12]: and (2.32)

Deriving the NLO JIMWLK Hamiltonian from the LCWF
The JIMWLK Hamiltonian is derived from the diagonal matrix element of the second quantizedŜ-matrix operator in the LCWF, From Σ one reads the Hamiltonian: At high energies, eikonal approximation for the scattering is appropriate: the effect of a dense target on the projectile is only to color rotate the partons in the projectile, Here S αβ is a Wilson line in the fundamental representation of the gauge fields along the target light-cone direction, while S ab A is the one in the adjoint representation. The action ofŜ on the valence charge distribution ρ a (x) can be expressed succinctly as Lie derivatives with respect to the Wilson lines [27] ρ a g (x)Ŝ|v = J a R,adj (x)|Ŝv ,Ŝρ a g (x)|v = J a L,adj (x)|Ŝv , and for quarks, (2.40) The Wilson lines S parametrize the target. No additional information about the target or any explicit form of S is needed for the derivation of the Hamiltonian. The algorithm for the derivation is as follows. One first computes the LCWF |ψ in perturbation theory. At LO, it is entirely determined by the matrix element (2.28). The LCWF at LO has, in addition to the vacuum, a soft gluon component, where E g (k) ≡ k 2 2k + is the gluon LC energy. At a next step, one calculates Σ in Eq. (2.33), from which the JIMWLK Hamiltonian is read off as explained above. For massless quarks, the complete and detailed derivation along this line can be found in Ref. [12]. Below we quote the final result for the sake of completeness of our presentation and easy comparison. In the next section, we will revise the calculation of [12] by extending it to include finite quark masses.
The NLO JIMWLK Hamiltonian with massless quarks are [11,12]: Quark mass corrections will affect, as will be shown in the next sections, the kernels K qq and K JSJ only, which we will denote as K m qq and K m JSJ (the upperscript m indicates dependence on mass). The other kernels have been presented for completeness only and are irrelevant for the rest of the paper.

NLO JIMWLK with massive quarks
Similarly to the massless case [12], there are two terms in the LCWF with quark-antiquark state. Those are depicted in Fig. 2. Let |ψ 1 qq denote the wave function corresponding to Fig. 2 (A), The energy denominator including quark masses is The relevant matrix elements are given in Eq. (2.27) and (2.28). Plugging in these elements, |ψ 1 qq becomes Let |ψ 2 qq be the wave function component corresponding to the instantaneous quark pair emission diagram (B) in Fig. 2, where the minus sign in front is due to perturbative expansion with one interaction insertion (see Eq. (3.1) of [12], for example). Plugging in Eq. (2.29) and the massive quark energy denominator, we have The sum of the two diagrams in Fig. 2 (A) and (B), |ψq qρ = ψ 1 qq + ψ 2 qq , where ξ ≡ p + /k + . Notice that the instantaneous interaction in Fig. 2 (B) is canceled exactly by one of the terms from diagram (A), which is shown in the second equality of Eq. (3.6).
With the wave function Eq. (3.6) in hand, one can extract the kernel K m qq of the JIMWLK Hamiltonian, as explained in the previous section. Denote the forward scattering amplitude for the quark anti-quark pair as Substituting Eq. (3.6) into Eq. (3.7) and using the equality, Here again, ξ ≡ p + /k + and θ ≡ v + /u + . After a few algebraic steps, detailed in Appendix A, we obtain The multiple momentum integral k,p,u,ṽ in Eq. (3.9) can be explicitly carried out with the help of the integrals collected in Appendix D, and the final result is The functions H i are defined below in terms of the Modified Bessel functions of the second kind K 0 and K 1 :
Eq. (3.10) can be massaged a little bit further to take the following form From the expression above one reads off the quark mass-modified kernel K m qq : We have not succeeded to carry out the ξ-integral analytically. Hence, the final result cannot be presented in terms of elementary functions. Yet, the ξ-integral should be relatively easy for numerical evaluation. The result has been written for a single flavor only. Summation over different flavors/masses is obviously implied.

3.2
K m JSJ with massive quarks There are two mass-dependent contributions to the kernel K m JSJ . The first one is direct contribution of the quark loop to the one gluon component of the LCWF (see Fig. 3). Another one originates from the subtraction of the K m qq kernel, to be discussed below. The quark loop correction to the single gluon component of the LCWF is displayed in Fig. 3. In the massless limit, the diagram is UV divergent, giving rise to running of the coupling (the N f term in the β-function), after standard dimensional regularization is applied [12]. While renormalization is straightforward in the massless case, the massive case turns out to be much more involved. Most of this subsection is devoted to carefully carrying out the renormalization procedure, which requires introduction of a mass-dependent counter term.

Quark loop effects in |ψ g
By |ψ 1 g we denote the renormalized component. Un-renormalized component |ψ 1 g is expressed in perturbation theory, The matrix element of H gqq is given in Eq. (2.27). The detailed calculation of Eq. (3.20) can be found in Appendix B. The final, dimensionally regularized result is where we had to introduce δ as an IR regulator, modifying the 1/k 4 pole in the second line of Eq. (3.21) to 1/(k 2 + δ 2 ) 2 . K 0 is logarithmically divergent as δ → 0, and The IR divergence which we encounter is a furious divergence due to use of the LC quantization. This divergence gets canceled by introduction of a mass dependent counter term (for a more detailed discussion, consult Refs. [28][29][30]). The mass dependent counter term is introduced in order to ensure that the gluon remains massless. In other words, the renormalization condition is absence of any mass corrections to the gluon mass, to all orders in perturbation theory. The dispersion relation on the light-cone coordinate is The gluon's pole mass m 2 has to vanish. At fixed k + and k ⊥ , the mass correction δm 2 is proportional to the LC energy correction δE g (k), So, vanishing of δm 2 is equivalent to vanishing of δE g (k). The LC energy (self-energy) correction at one loop (order O(g 2 )) shown in Fig. 4 (a) is Derivation of Eq. (3.25) is nearly identical to that of Eq. (3.20) (see Appendix B) except for the difference in the energy denominators. The denominator in δE g (k) where, again, ξ ≡ p + /k + andp ≡ p − ξk. In order to impose δm 2 = 0, the following mass dependent counter term ( Fig. 4 (b)) is added to the effective LC Hamiltonian While it is not obvious from Eq. (3.28), ∆ is proportional to m 2 . In dimensional regularization, The counter term (3.27) contributes to the one gluon component of the wave function Using the matrix element it is clear that the energy denominator 1/(E g (r)E g (k)) give rise to 1/k 4 pole, which cancels the similar pole in Eq. (3.21). Introducing the IR regulator δ, as was done in Eq. (3.21), one immediately obtains (with the help of Eq. (D.8)) Combining the perturbative result with the counter term In fact, introduction of the IR regulator is not necessary, if Eq. (3.28) is used for the counter term in its integral form, instead of the expression (3.29). Adding the counter term (3.30) to the one gluon component (3.20) prior to momentum integration, we obtain The second line of Eq. (3.34) is . (3.36) And then the inverse Fourier Transform of k-integral in Eq. (3.34) back to x, z is straightforward using the integrals collected in Appendix D. The result reads 38) The final result in (3.37) is IR finite and has a well defined massless limit. Using small argument expansion of the Bessel functions, we obtain The ln m 2 in the expansion of g 1 (m|X|) is canceled by the explicit ln(m 2 /μ 2 ) term (3.37). This way, the massless result [12] is recovered.

Subtraction of K m qq
In the massless case, the kernel K m qq (2.51) has a non-integrable UV divergence when z → z , even though the full Hamiltonian is free of UV divergencies. The same divergence appears in (3.19), since introducing masses for quarks does not modify the UV behavior.
In the Hamiltonian, the UV divergence is regularized by subtracting the term (see the third line in Eq. (2.42)) The very same term is then added back in the first line in Eq. (2.42) (the real emission term) which is then absorbed into the definition of K JSJ (x, y, z). This contribution is given by z K qq x, y, z, z . (3.42) We postpone the discussion of K m JSJ until the next subsection. Meanwhile we focus on the subtraction. Hence, our goal is to calculate the z -integral of Eq. (3.19). It turns out, however, that it is difficult to carry out the integration directly. Instead, we start from Eq. (A.5) and perform the z -integral first, before doing the Fourier transform (more details on the derivation are given in Appendix C): The z -integral is now trivial and contributes a delta-function, which is then used to eliminate another integral, the v-integral. 1 Since the phase factor of the Fourier transform does not involve p, we first carry out the p-integral: where |ψ LO g is the LO one-gluon emission wave function given in Eq. (2.41), and g 1 (x) is defined in Eq. (3.40). Then, discarding the 1/ term, the evolution kernel reads where 1/128 comes from 1/(64 × 2) due to the convention that K JSJ corresponds to the coefficient of 2J a L S ab J b R . In the massless limit of g 1 (m|X|) and g 1 (m|Y |), Eq. (3.47) reproduces the corresponding massless limit of Ref. [12] (Eq. (4.36) therein). As discussed previously, the final result K m JSJ is obtained after subtraction with the subtraction contribution given in Eq. (3.44). We note in passing that part of our result can be identified with the mass effects in the running coupling, since the coupling running is related to the resummation of the quark loop diagram in Fig. 4. We, however, don't see it instructive to separate these effects from the rest of the kernel. Finally, we would like to remark on the virtual terms of the type J a L J a L . They are expected to come with the same kernel K m JSJ , though we have not explicitly demonstrated this here. These terms appear from the normalization of the wave function. It is straightforward to see [12] that the kernel in front of the virtual terms is indeed the same as the one in front of the real term, including the mass corrections.

CWZ subtraction scheme
We would like to comment more on the subtraction scheme used in our work. 2 So far, we have been using MS subtraction scheme for the massive quarks, which means a quadratic gluon counter terms are introduced to cancel the 1/ UV divergence in Eq. (3.37). Note that the 1/ canceling counter term (which exists also for diagrams with massless quarks) is in addition to the mass dependent counter term in Figure 4 (b). When quark masses are relevant in a physical process, it is also useful to use the Collins-Wilczek-Zee (CWZ) subtraction scheme [34][35][36][37], which makes the decoupling of heavy quark manifest in the large mass limit. In this subtraction scheme, the ordinary MS subtraction is carried out for diagrams with active quarks (light quarks) and gluons, while for inactive quarks (heavy quarks) zero momentum subtractions are implemented. Figure 5 shows how zero momentum subtraction is carried out for the gluon vacuum polarization: the loop integral with external momentum k → 0 is subtracted from the original loop integral. This vacuum polarization diagram is a sub-diagram for the loop correction diagram as shown in Figure 3. So a zero momentum subtraction should be Figure 5. Zero momentum subtraction of the sub-diagram (meaning the gluon self-energy diagram) implemented to calculated the wave function due to loop diagram in Figure 3.
implemented for Eq. (3.36). By setting k → 0, one obtains for the wave function which is to be subtracted from Eq. (3.34). The zero momentum subtracted wave function is then simply obtained by replacingμ by m and removing 2 , as is usually the case for CWZ subtraction in ordinary equal time quantization. Note that since the loop integral (i.e., Eq. (3.36) and Figure 5) is independent of k + , i.e., the k + component does not flow into the loop, the zero momentum subtraction simply means zero k ⊥ momentum subtraction. Similar argument shows that the CWZ implementation for dz K m qq (x, y, z, z ) is equivalent to replacingμ → m. After replacingμ → m in Eq. ) have a modified Bessel function in every term (either K 0 [m|X|] or K 1 [m|Y |] for instance), which in the heavy mass limit are exponentially suppressed, so that the decoupling of heavy quark is manifest in the m → ∞ limit.
In the above, we have done the zero momentum subtraction for massive quarks in all external momentum region (i.e., |k| ∈ (0, ∞) in Eq. (3.36)), which shows the explicit decoupling of massive quarks in the large mass limit, and which is also a good cross check of our computations. In practice, however, we should carry out the CWZ subtraction of Eq. (3.36) step by step when multiple quarks with different masses are involved (similar to the running coupling that should be matched step by step in the CWZ scheme). So with multiple massive quarks involved, after CWZ subtraction, Eq. (3.36) should be replaced with a sum of step functions as follows which is then plugged back in Eq. (3.34). Similar manipulations should be done for the subtraction contribution of K m JSJ , i.e., z K m qq (x, y, z, z ). However, the Fourier Transforms back to coordinate space (the counter parts of Eq. (3.37) and Eq. (3.44)) become complex, and we do not discuss them further in the current work.

Conclusions
In this paper, we generalized the work of Ref. [12], exploring high energy evolution with non-vanishing quark masses. We computed the hadronic light-cone wave function of a fast moving projectile, focusing on the effects originating from non-vanishing quark masses. As an application of the wave function, we derived the NLO evolution kernels for the JIMWLK Hamiltonian, particularly the kernels K m qq and K m JSJ . Both enter the dipole evolution at NLO [10] through the form presented in [11].
There are only two types of diagrams that involve massive quarks other than those that have already been calculated in Ref. [12]: the quark splitting diagram (Fig. 2) and the virtual quark loop diagram for one gluon component of the wave function (Fig. 3). At every step, we carefully checked that the massless limit of Ref. [12] is recovered.
While a priori inclusion of mass looks like a very straightforward generalization of the previous calculations, we have encountered several complications. First of all, on the technical level, some integrals became too complex to be performed analytically. Furthermore, for the quark loop diagram, a new IR divergence (divergence coming from small k ⊥ ) emerges. This is somewhat surprising since masses generally cure IR divergences and not the other way around. In fact, this IR divergence is, to some extent, a known "artifact" of the LC quantization. It turns out that the cancellation of this divergence is realized by introducing a mass dependent gauge field counter term, which is required for the gluon to remain massless to all orders in perturbation theory [28][29][30].
Unsurprisingly, we have not been able to obtain simple fully analytical expressions for the kernels K m qq and K m JSJ analogous to the ones of the massless limit. Our results for the kernels are expressed as one dimensional integrals. We are hopeful that they are nevertheless sufficiently simple and could be implemented in future numerical studies and phenomenology. Obviously, a complete NLO calculation is already very challenging, particularly due to multiple coordinate integrations. In principle, one might expect that considering various quark mass limits (light/heavy limits) might lead to more tractable analytical expressions. We have considered both limits but have not observed any simplifications worth reporting. This is largely due to the absence of any clear scale separation between the mass and momenta which are integrated over in the kernels.

A Computation of Σq q at NLO
In this Appendix, details on derivation of Eq. (3.10) are provided. We start from Eq. (3.6) and the definition of Σ NLŌ qq in Eq. (3.7): where θ ≡ v + /u + . After Fourier-transforming ρ b (u) and ρ a (−k) back to coordinate space (using Eq. (2.32)), plugging in Eq. (3.8), and also converting where Z, X , Y are defined in Eq. (2.44), and λ χ λ χ † λ = Λ + is used to obtain the trace in the last equality. After some Dirac algebra, the trace part becomes which for massless quarks reproduces Eq. (4.12) in Ref. [12]. Eq. (3.9) is obtained after the following change of variables is implemented, where the matrix elements (2.27) and (2.28) are used. The polarization indices λ 1 and λ 2 are summed over. The following part of the integrand can be converted (after summing over λ 1 and λ 2 ) to has been used. After some Dirac algebra, the trace expression (B.3) becomes where ξ ≡ p + /k + , andp ≡ p − ξk. Substituting (B.5) back into Eq. (B.1) (and changing variables to ξ andp), where tr[t a t b ] = δ ab /2 is used. The d 2p integral is UV divergent. To regularize it, we use dimensional regularization with the dimension 2 → d = 2 − : Now we can Fourier transform ρ a (−k) and |g a i (k) (employing the conventions (2.32) and (2.15)) In contrast to the massless case, however, the k integral is IR divergent. In order to regularize the divergence, a "gluon mass" δ can be introduced, i.e., we make the following replacement 1 (B.9) With this IR regularization, d 2 k and dξ integrals can be performed (using the Fourier integral Eq. (D.9)), and we obtain the final expression in Eq. (3.21). As discussed in the main text, the IR regulator is eventually removed by proper renormalization, which preserves masslessness of the gluon.

D Integrals
In addition to some integrals tabulated in Ref. [12], we also used the following integrals. x · y x 2 y 2 ln (x + y) 2 y 2 .