A unified description of DGLAP, CSS, and BFKL: TMD factorization bridging large and small x

The concept of factorization can be elegantly introduced within the background field method. Starting with a matrix element of an arbitrary operator $\mathcal{O}(\hat{C})$ which depends on quark and gluon fields,²²2Here we understand $C$ as a general notation for the quark and gluon fields. one can write it as a functional integral over those fields

$\displaystyle\langle P_{1}|\mathcal{O}(\hat{C})|P_{2}\rangle=\int\mathcal{D}C~{}\Psi^{\ast}_{P_{1}}(\vec{C}(t_{f}))\mathcal{O}(C)\Psi_{P_{2}}(\vec{C}(t_{i}))e^{iS_{\rm QCD}(C)}\,,$

(9)

where $\Psi_{P_{2}}$ and $\Psi_{P_{1}}$ are initial and final state wave functions at $t_{i}\to-\infty$ and $t_{f}\to\infty$.

Evaluation of this functional integral is hard in general. Hence, we split the field into different components, and rewrite the integral as independent functional integrals over those components. This is the logic of the background field method. This separation is arbitrary in general, but in the context of the QCD factorization, we shall assume that there is a set of factorization scales $\sigma$ that separate the components from each other, see Fig. 2a.

For example, in the case of deep inelastic scattering (DIS) in the Bjorken limit, it is sufficient to split the field into two field modes $A$ and $B$,

$\displaystyle C_{\mu}\to A_{\mu}+B_{\mu}\,.$

(10)

The fields are separated by the value of the transverse momenta, i.e. one can introduce a cut-off for the transverse momenta $\sigma=\mu$, such that the fields $A$ are defined as having $k^{2}_{\perp}>\mu^{2}$, and fields $B$ are characterized with $k^{2}_{\perp}<\mu^{2}$. This is called the collinear factorization scheme which is defined by a strong ordering of emission in transverse momenta. Another example where the factorized QCD medium can be described by splitting the field into two modes is the high energy rapidity factorization in DIS. In this case, the field mode $A$ has longitudinal momentum fraction $k^{-}>\sigma$ and the mode $B$ has $k^{-}<\sigma$. The factorization scale $\sigma$ in this case is the rigid cut-off in $k^{-}$.

For brevity, let’s assume that the factorization can be described with only two modes as in (10). And a scale $\sigma$ is a factorization scale separating components $A$ and $B$. In general, the structure of dynamical field modes can be more complicated and depend on the nuances of the factorization in the process. For example, as we will see in the next subsection, to describe the TMD factorization one has to introduce four modes: one hard mode (“quantum”), two collinear modes (background), and a soft mode for an intersection between the collinear modes.

We will call the $A$ fields “quantum” to indicate that we aim to integrate over this mode. The fields $B$ are usually called the background fields. The functional integral over this mode contains the non-perturbative part of the scattering process. One can think about the background fields as fields associated with the target, and the “quantum” fields as fields of the hard scattering.

After splitting (10) we can rewrite the functional integral in Eq. (9) as

$\displaystyle\langle P_{1}|\mathcal{O}(\hat{C})|P_{2}\rangle=\int\mathcal{D}B~{}\Psi^{\ast}_{P_{1}}(\vec{B}(t_{f}))\tilde{\mathcal{O}}(B,\sigma)\Psi_{P_{2}}(\vec{B}(t_{i}))e^{iS_{\rm QCD}(B)}\,,$

(11)

where the integral over “quantum” mode $A$ is,

$\displaystyle\tilde{\mathcal{O}}(B,\sigma)=\int\mathcal{D}A~{}\mathcal{O}(A+B)e^{iS_{\rm bQCD}(A,B)}\,.$

(12)

The integration over $A$ should be performed with the QCD action in the background field,

$\displaystyle S_{\rm bQCD}(A,B)=S_{\rm QCD}(A+B)-S_{\rm QCD}(B)\,.$

(13)

In Eq. (11) we assume that due to the factorization regime in the high-energy scattering, the wave functions depend only on the $B$ fields, which, as we mentioned above, should be associated with the fields of the target.

The integral over “quantum” fields $A$ in Eq. (12), i.e. fields of the projectile, can be evaluated perturbatively. However, it usually requires some sort of expansion in factorization scales, e.g. twist expansion in the collinear factorization, etc. The result of integration has a general form

$\displaystyle\tilde{\mathcal{O}}(B,\sigma)=\sum_{i}H_{i}(\sigma)\otimes\mathcal{V}_{i}(B,\sigma)\,.$

(14)

where $\mathcal{V}_{i}$ are operators constructed from background mode $B$. Since the background fields $B$ depend on the factorization scale $\sigma$, so do the operators, as we explicitly indicate. $H_{i}(\sigma)$ is a hard function that we obtain after integration over the “quantum” fields $A$.

After this, Eq. (11) reads

$\displaystyle\langle P_{1}|\mathcal{O}(\hat{C})|P_{2}\rangle=\sum_{i}H_{i}(\sigma)\otimes\langle P_{1}|\mathcal{V}_{i}(\sigma)|P_{2}\rangle\,,$

(15)

where the matrix elements of operators $\mathcal{V}_{i}$ are generated by the functional integral over background fields $B$,

$\displaystyle\langle P_{1}|\mathcal{V}_{i}(\sigma)|P_{2}\rangle\equiv\int\mathcal{D}B~{}\Psi^{\ast}_{P_{1}}(\vec{B}(t_{f}))\mathcal{V}_{i}(B,\sigma)\Psi_{P_{2}}(\vec{B}(t_{i})).$

(16)

Here, the index $i$ enumerates different types of operators. These matrix elements can be further parameterized in terms of distribution functions. For example, in the collinear factorization, those are the well-known collinear PDFs.

Eq.(15) is a factorization formula that we obtain by assuming the QCD medium develops different dynamical modes, so the separation (10) can be done and the integration over “quantum” modes $A$ can be considered independently in a fixed configuration of background fields $B$.

Note that the structure of factorization depends on the scattering process and kinematics. However, due to a wide separation of scales, each type of factorization is characterized by a set of large logarithms. These logarithms should be resummed. This is usually done by solving the RGEs which are differential equations with respect to the factorization scales, see for example Eqs. (4) and (5) for the TMD factorization.

One can use the following approach to extract the dependence on the factorization scales in the background field method. Starting with a matrix element (16), we introduce a new set of factorization scales $\sigma^{\prime}$ and split the background field $B$ as $B\to B^{\rm q}+B^{\rm bg}$. Here the fields $B^{\rm q}$ correspond to dynamical modes between scales $\sigma$ and $\sigma^{\prime}$, and $B^{\rm bg}$ is a new background field corresponding to fields below the scale $\sigma^{\prime}$, see Fig. 2b.

Since the cut-offs $\sigma$ act as upper cut-offs for the field modes, we will refer to them as the UV scales. Similarly, since the cut-offs $\sigma^{\prime}$ act as lower cut-offs for the $B^{\rm q}$ modes, we will refer to them as the IR scales.

To study the dependence on the set of UV factorization scales $\sigma$, we need to integrate over $B^{\rm q}$ fields. This integration is done perturbatively and leads to the answer

$\displaystyle\langle P_{1}|\mathcal{V}_{i}(\sigma)|P_{2}\rangle=\sum_{j}C_{ij}(\sigma,\sigma^{\prime})\otimes\langle P_{1}|\mathcal{V}_{j}(\sigma^{\prime})|P_{2}\rangle\,,$

(17)

where coefficients $C_{ij}(\sigma,\sigma^{\prime})$ are obtained by integrating over the $B^{\rm q}$ modes and describe dynamics between the UV scales $\sigma$ and the IR scales $\sigma^{\prime}$. One can now differentiate with respect to the UV scales $\sigma$ and obtain a system of evolution equations that define the dependence of the matrix element (16) and the corresponding distribution functions on the set of factorization scales $\sigma$, e.g. Eqs. (4) and (5) in the TMD factorization.

At the same time, the coefficients $C_{ij}(\sigma,\sigma^{\prime})$ also describe the perturbative component of the IR structure of the matrix element through the dependence on the IR cut-off $\sigma^{\prime}$. The non-perturbative part of the IR structure is encoded in the matrix element of the operators $\mathcal{V}_{j}(\sigma^{\prime})$.

This paper aims to calculate the $C_{ij}(\sigma,\sigma^{\prime})$ coefficients for the gluon TMDPDFs. The role of the UV scales $\sigma$ is played by the cut-offs $(\mu^{2}_{\rm UV},\nu)$, while the IR scales $\sigma^{\prime}$ correspond to the cut-offs $(\mu^{2}_{\rm IR},\rho)$, see Fig. 1. More details can be found in sections I and III.

Another crucial but technical point is how we can actually introduce different field modes in practice. We can use rigid cut-offs, for example, as it is done in the high-energy rapidity factorization approach. These cut-offs subsequently appear in the perturbative calculation of functional integrals, e.g. in Eq. (12), as cut-offs in loop integrals over momenta. These cut-offs regulate all divergences of the loop integrals. However, these rigid cut-offs are rarely used due to technical difficulties.

We use a different method in this paper, known as the renormalization approach. In this scheme, each divergent integral is regulated with a regulator accompanied by a renormalization scale. This scale plays the role of a cut-off separating different modes. As mentioned before, we shall use dimensional regularization to separate the modes in the transverse momentum. To separate the modes in rapidity, we will use the regularization of the longitudinal momentum fraction integrals as in Eq. (8).³³3There is a number of alternative methods: deviation from the light-cone direction [1, 7], $\delta$-regulator [101, 102], exponential regulator [35], and analytic regulator [103, 104, 105]. Note that in the renormalization approach, we will find some unphysical poles which is a consequence of the separation of different dynamical field modes. These poles need to be removed by a proper renormalization of operators.

In the previous section, we showed how the background field method can be used to construct a factorization scheme. While our discussion was general, different factorization schemes can be substantially distinct. This can be most easily seen in the structure of the dynamical modes, e.g., compare the structure of the high-energy rapidity factorization in DIS, see Fig. 2, and the TMD factorization presented in Fig. 1.

In this section, we will give an overview of the transverse momentum dependent (TMD) factorization. In particular, we will discuss the structure of modes in Fig. 1 using the framework of the background field approach and give a formal definition of the corresponding operators.

The TMD factorization is applied to high-energy scattering reactions when a final state with a small transverse momentum $q_{\perp}$ is measured. The scattering reaction includes two energetic hadrons moving along $n^{\mu}=(1,0,0,-1)$ and $\bar{n}^{\mu}=(1,0,0,1)$ directions, e.g. two colliding hadrons in the Drell-Yan scattering with momenta $\bar{P}^{\mu}=P^{-}n^{\mu}+M^{2}\bar{n}^{\mu}/P^{-}$ and $P^{\mu}=P^{+}\bar{n}^{\mu}+M^{2}n^{\mu}/P^{+}$ in the center of mass frame, where $M^{2}$ is hadron’s mass.

Let’s now introduce the TMD factorization scheme in the background field approach. Following the logic of the previous section, we start with the background field $B$. Since the scattering process involves two hadrons, we need to introduce two types of background field modes. We will refer to the modes associated with the $P^{\mu}$ hadron as collinear modes and modes of the $\bar{P}^{\mu}$ hadron as anti-collinear. Both of these modes describe the dynamics of the corresponding hadron.

The hard scattering process is described by the “quantum” mode $A$. This mode is separated from the background modes by a factorization scale $\mu^{2}_{\rm UV}$, see Fig. 1. The integration over fields of this mode yields a hard function in the factorization formula, c.f. Eqs. (1) and (15). In integration over “quantum” modes collinear and anti-collinear background modes yield two TMD operators, i.e. operators $\mathcal{V}_{i}$ in the notation of Eq. (15). The matrix elements of these operators can be parameterized in terms of the TMDPDFs, see Eq. (1). For brevity, we refer to these matrix elements as TMDPDFs.

The resulting matrix element of the TMD operator constructed from the gluon collinear modes has the form

$\displaystyle\mathcal{B}_{ij}(x_{B},b_{\perp})=\int^{\infty}_{-\infty}dz^{-}e^{-ix_{B}P^{+}z^{-}}\langle P,S|\bar{T}\{F^{m}_{-i}(z^{-},b_{\perp})[z^{-},\infty]^{ma}_{b}\}T\{[\infty,0^{-}]^{an}_{0}F^{n}_{-j}(0^{-},0_{\perp})\}|P,S\rangle\,.$

(18)

Although there is a similar quark TMD operator and our approach can be applied to it as well in this paper we focus only on the gluon case. Moreover, one can also write a similar set of operators for the anti-collinear background modes. Since the analysis is identical, we do not explicitly present it in this paper. Finally, we will only consider the unpolarized matrix elements to simplify the discussion.

In Eq. (18), $i$ and $j$ are Lorentz indices. The adjoint Wilson lines are along the light-cone direction

$\displaystyle[x^{-},y^{-}]_{z_{\perp}}=\mathcal{P}\exp[ig\int^{x^{-}}_{y^{-}}dz^{-}A_{-}(z^{-},z_{\perp})\Big{]}\,.$

(19)

We assume that the Wilson lines are connected at infinity with a transverse gauge link, which makes the operator gauge invariant. We choose future-point Wilson lines which correspond to SIDIS. The $\tilde{T}$ notation is to indicate that the corresponding operators should be inverse-time ordered. The Fourier transformation of the matrix element is defined as

$\displaystyle\mathcal{B}_{ij}(x_{B},p_{\perp})=\int d^{2}b_{\perp}e^{ip_{\perp}b_{\perp}}\mathcal{B}_{ij}(x_{B},b_{\perp})\,.$

(20)

In the TMD factorization scheme, the collinear and anti-collinear modes are only separated by their rapidity $y=1/2\ln(k^{-}/k^{+})$. There might be an intersection region depending on how these modes are separated. In the TMD factorization framework, this is called the soft region. From the point of view of the background field approach, the fields of this mode should still be considered as the background fields.

To resolve the potential double counting in the intersection region, one must introduce a rapidity cut-off $\nu$ and assume that the collinear fields in Eq. (18) are constructed only from the modes below this cut-off. In Fig. 1, these fields correspond to the blue and red regions. Similarly, the anti-collinear modes correspond to the yellow and purple sectors.

To take into account the contribution of the soft region, we need to multiply (18) with a soft function constructed from the soft modes. Note that the structure of the soft factor depends on the details of separation in rapidity between the fields of the collinear and anti-colliner modes.⁴⁴4In principle, one can construct the collinear and anti-collinear modes to avoid the intersection and thus render the soft factor trivial, see Refs. [106, 107]. We use a renormalization approach and define the collinear fields using the $\eta$ regulator, see Eq. (8). In this case the soft factor $\sqrt{\mathcal{S}(b_{\perp})}$ is given by a vacuum matrix element

$\displaystyle\mathcal{S}(b_{\perp})=\frac{1}{N^{2}_{c}-1}\langle 0|{\rm Tr}[S^{\dagger}_{\bar{n}}(b_{\perp})S_{n}(b_{\perp})S^{\dagger}_{n}(0_{\perp})S_{\bar{n}}(0_{\perp})]|0\rangle\,,$

(21)

where the Wilson lines,

$\displaystyle S_{n}(b_{\perp})=P\exp[ig\int^{\infty}_{0}d(x\cdot\bar{n})n\cdot A(x\cdot\bar{n},b_{\perp})\Big{]}$

(22)

are constructed from soft modes defined by the regularization

$\displaystyle\int\frac{dk^{+}dk^{-}}{k^{+}k^{-}}=\nu^{\eta}2^{-\eta/2}\int\frac{dk^{+}dk^{-}}{k^{+}k^{-}}|k^{+}-k^{-}|^{-\eta}\,.$

(23)

The fields of the soft mode correspond to the green sector in Fig. 1.

The full matrix element of the background modes is defined as

$\displaystyle f_{ij}(x_{B},b_{\perp})=\sqrt{\mathcal{S}(b_{\perp})}\mathcal{B}_{ij}(x_{B},b_{\perp})\,.$

(24)

To study the dependence of this matrix element on the factorization scales $\sigma=(\mu^{2}_{\rm UV},\nu)$, we implement the approach discussed in the previous section. We introduce a new set of factorization scales $\sigma^{\prime}=(\mu^{2}_{\rm IR},\rho)$ to split the background collinear mode into the $B^{\rm q}$ mode which lies between the scales $\sigma$ and $\sigma^{\prime}$, see the blue region in Fig. 1, and the $B^{\rm bg}$ mode which lies below the scales $\sigma^{\prime}$, i.e. the red sector in Fig. 1.

We aim to integrate over $B^{\rm q}$ and soft modes to obtain a hard function containing the perturbative structure of the matrix element (24) (see Eq. (17)). The non-perturbative structure is encoded in the matrix element constructed from the fields of the $B^{\rm bg}$ mode. The details and results of this calculation are presented in the next section.

Before moving to the main part of the paper, we would like to point out that our definition of these modes differs from the Soft Collinear Effective Theory (SCET). The structure of modes in SCET is summarized in Fig. 3.

The collinear modes are defined to have the following scaling of momenta: $p_{n}=(p_{n}^{+},p_{n}^{-},\vec{p}_{n\perp})\sim Q(\lambda^{2},1,\lambda)$ and $p_{\bar{n}}\sim Q(1,\lambda^{2},\lambda)$, where $\lambda\sim q_{\perp}/Q^{2}\ll 1$. Both collinear fields are approximately on the mass-shell $p^{2}_{n}\sim p^{2}_{\bar{n}}\sim q^{2}_{\perp}$, so they lie on the same hyperbola. The hard fields are defined as having $p^{2}\geq Q^{2}$. The soft emission is defined as having $p_{s}\sim Q(\lambda,\lambda,\lambda)$. The main difference with the background field approach is that in the latter the field modes are not necessarily on the mass-shell; thus in principle, we consider more general kinematics. Also, to define the modes we do not use any scaling of momenta but rather use a set of cut-offs to delineate regions in momentum space.

In this section, we present an analysis of the real emission diagrams. Some typical diagrams are given in Fig. 4. As per the background field approach, the loop is constructed from the “quantum” fields $B^{\rm q}$, and the calculation is done in the background field $B^{\rm bg}$.

It is convenient to start the derivation with a calculation of the emission vertex

$\displaystyle L^{ab}_{\mu j}(k,y_{\perp},x_{B})\equiv i\lim_{k^{2}\to 0}k^{2}\langle B^{{\rm q}a}_{\mu}(k)\int^{\infty}_{-\infty}dy^{-}e^{ix_{B}P^{+}y^{-}}[\infty,y^{-}]^{bd}_{y}F^{{\rm q+bg};d}_{-j}(y^{-},y_{\perp})\rangle_{B^{\rm bg}}\,,$

(25)

which corresponds to a sum of diagrams in Fig. 5. Once this vertex is calculated, the sum of real emission diagrams can be easily obtained by taking a product of two vertices as

$\displaystyle\langle\tilde{\mathcal{F}}^{a}_{i}(x_{B},x_{\perp})\mathcal{F}^{a}_{j}(x_{B},y_{\perp})\rangle^{\rm real}=-\int\frac{{\textstyle d}\lower 0.12915pt\hbox{\kern-3.80005pt${}^{\scriptstyle-}$}\kern-0.50003pt{}k^{-}}{2k^{-}}\int{\textstyle d}\lower 0.12915pt\hbox{\kern-3.80005pt${}^{\scriptstyle-}$}\kern-0.50003pt{}^{2}k_{\perp}\tilde{L}^{\ \mu ba}_{i}(k,x_{\perp},x_{B})L^{ab}_{\mu j}(k,y_{\perp},x_{B})\,.$

(26)

The details of the calculation of the emission vertex can be found in Appendix A. Taking the product, Eq.(26), of two emission vertices Eq.(112) and Eq.(113), it is easy to obtain an expression for the matrix element in Eq.(18),

	$\displaystyle\mathcal{B}^{\rm q(1)+bg;real}_{ij}(x_{B},b_{\perp})$		(27)
	$\displaystyle=-16\pi\alpha_{s}N_{c}\int{\textstyle d}\lower 0.12915pt\hbox{\kern-3.80005pt${}^{\scriptstyle-}$}\kern-0.50003pt{}^{2}p_{\perp}e^{ip_{\perp}b_{\perp}}\int\frac{{\textstyle d}\lower 0.12915pt\hbox{\kern-3.80005pt${}^{\scriptstyle-}$}\kern-0.50003pt{}k^{-}}{2k^{-}}\int{\textstyle d}\lower 0.12915pt\hbox{\kern-3.80005pt${}^{\scriptstyle-}$}\kern-0.50003pt{}^{2}k_{\perp}\Big{[}\Big{(}\frac{(p+k)_{l}}{2x_{B}P^{+}k^{-}+k^{2}_{\perp}}\frac{x_{B}P^{+}k^{-}\delta^{k}_{i}-k^{k}k_{i}}{2x_{B}P^{+}k^{-}+p^{2}_{\perp}}-\frac{\delta^{k}_{l}p_{i}+g_{li}k^{k}}{2x_{B}P^{+}k^{-}+p^{2}_{\perp}}\Big{)}$
	$\displaystyle\times\Big{(}\frac{(p+k)_{m}}{2x_{B}P^{+}k^{-}+k^{2}_{\perp}}\frac{x_{B}P^{+}k^{-}g_{kj}-k_{k}k_{j}}{2x_{B}P^{+}k^{-}+p^{2}_{\perp}}-\frac{g_{mk}p_{j}+g_{mj}k_{k}}{2x_{B}P^{+}k^{-}+p^{2}_{\perp}}\Big{)}+\frac{g_{il}}{k^{2}_{\perp}}\Big{(}\frac{(p+k)_{m}}{2x_{B}P^{+}k^{-}+k^{2}_{\perp}}\frac{x_{B}P^{+}k^{-}k_{j}+k^{2}_{\perp}k_{j}}{2x_{B}P^{+}k^{-}+p^{2}_{\perp}}$
	$\displaystyle-\frac{k_{m}p_{j}-g_{mj}k^{2}_{\perp}}{2x_{B}P^{+}k^{-}+p^{2}_{\perp}}\Big{)}+\Big{(}\frac{(p+k)_{l}}{2x_{B}P^{+}k^{-}+k^{2}_{\perp}}\frac{x_{B}P^{+}k^{-}k_{i}+k^{2}_{\perp}k_{i}}{2x_{B}P^{+}k^{-}+p^{2}_{\perp}}-\frac{k_{l}p_{i}-g_{li}k^{2}_{\perp}}{2x_{B}P^{+}k^{-}+p^{2}_{\perp}}\Big{)}\frac{g_{mj}}{k^{2}_{\perp}}-\frac{g_{il}g_{mj}}{k^{2}_{\perp}}\Big{]}$
	$\displaystyle\times\int d^{2}z_{\perp}e^{i(k-p)_{\perp}z_{\perp}}\mathcal{B}^{\rm bg}_{lm}(x_{B}+\frac{k^{2}_{\perp}}{2P^{+}k^{-}},z_{\perp})\,,$

where we indicate that the matrix element on the right has operators constructed from $B^{\rm bg}$ fields, but the one on the left has contributions from both $B^{\rm bg}$ and $B^{\rm q}$. Superscript $q(1)$ indicates that “quantum” fields are evaluated at order $\alpha_{s}$.

To separate the rapidity divergences appearing in the integral over $k^{-}$ from the divergences in the transverse integral over $k_{\perp}$, it is convenient to introduce the variable

$\displaystyle z\equiv\frac{x_{B}}{x_{B}+\frac{k^{2}_{\perp}}{2P^{+}k^{-}}}\,.$

(28)

In terms of this variable, the equation can be re-written in a compact form,

	$\displaystyle\mathcal{B}^{\rm q(1)+bg;real}_{ij}(x_{B},b_{\perp})=-4\alpha_{s}N_{c}\int{\textstyle d}\lower 0.12915pt\hbox{\kern-3.80005pt${}^{\scriptstyle-}$}\kern-0.50003pt{}^{2}p_{\perp}e^{ip_{\perp}b_{\perp}}\int^{1}_{0}\frac{dz}{z(1-z)}$		(29)
	$\displaystyle\times\int{\textstyle d}\lower 0.12915pt\hbox{\kern-3.80005pt${}^{\scriptstyle-}$}\kern-0.50003pt{}^{2}k_{\perp}\Big{[}\mathcal{R}^{a}_{ij;lm}(z,p_{\perp},k_{\perp})+\mathcal{R}^{b}_{ij;lm}(z,p_{\perp},k_{\perp})+\mathcal{R}^{c}_{ij;lm}(k_{\perp})\Big{]}\int d^{2}z_{\perp}e^{i(k-p)_{\perp}z_{\perp}}\mathcal{B}^{\rm bg}_{lm}(\frac{x_{B}}{z},z_{\perp})\,,$

where the real emission kernels are coming from Eq. (27),

	$\displaystyle\mathcal{R}^{a}_{ij;lm}(z,p_{\perp},k_{\perp})\equiv(1-z)^{2}\Big{(}\frac{1}{2}\frac{(p+k)_{l}}{k^{2}_{\perp}}\frac{zk^{2}_{\perp}\delta^{k}_{i}-2(1-z)k^{k}k_{i}}{zk^{2}_{\perp}+(1-z)p^{2}_{\perp}}-\frac{\delta^{k}_{l}p_{i}+g_{li}k^{k}}{zk^{2}_{\perp}+(1-z)p^{2}_{\perp}}\Big{)}$
	$\displaystyle\times\Big{(}\frac{1}{2}\frac{(p+k)_{m}}{k^{2}_{\perp}}\frac{zk^{2}_{\perp}g_{kj}-2(1-z)k_{k}k_{j}}{zk^{2}_{\perp}+(1-z)p^{2}_{\perp}}-\frac{g_{mk}p_{j}+g_{mj}k_{k}}{zk^{2}_{\perp}+(1-z)p^{2}_{\perp}}\Big{)},$		(30)

	$\displaystyle\mathcal{R}^{b}_{ij;lm}(z,p_{\perp},k_{\perp})\equiv(1-z)\frac{g_{il}}{k^{2}_{\perp}}\Big{(}\frac{(p+k)_{m}}{2}\frac{zk_{j}+2(1-z)k_{j}}{zk^{2}_{\perp}+(1-z)p^{2}_{\perp}}-\frac{k_{m}p_{j}-g_{mj}k^{2}_{\perp}}{zk^{2}_{\perp}+(1-z)p^{2}_{\perp}}\Big{)}$
	$\displaystyle+(1-z)\Big{(}\frac{(p+k)_{l}}{2}\frac{zk_{i}+2(1-z)k_{i}}{zk^{2}_{\perp}+(1-z)p^{2}_{\perp}}-\frac{k_{l}p_{i}-g_{li}k^{2}_{\perp}}{zk^{2}_{\perp}+(1-z)p^{2}_{\perp}}\Big{)}\frac{g_{mj}}{k^{2}_{\perp}},$		(31)

$\displaystyle\mathcal{R}^{c}_{ij;lm}(k_{\perp})\equiv-\frac{g_{il}g_{mj}}{k^{2}_{\perp}}.$

(32)

The above kernels are generated by different real emission diagrams. Specifically, the functions $\mathcal{R}^{a}$ and $\mathcal{R}^{b}$ represent diagrams with the same topology as the first two diagrams in Fig. 4.

It is straightforward to see that these terms do not contain any singularities. In particular, a potential rapidity divergence at $z\to 0$ is regulated by a non-zero value of $x_{B}$. Therefore, these terms do not provide any logarithmic contribution and yield finite terms in our final result.

Yet, in Sec. IV, we will find that although these terms are finite for a generic kinematics, in the collinear limit (when $k_{\perp}-p_{\perp}\to 0$), they become large and diverge logarithmically. These terms ultimately generate a part of the DGLAP splitting function.

At the same time, we find that the last term in Eq. (29), i.e. $\mathcal{R}^{c}$, contains divergences. This term is given by the last diagram in Fig. 4.⁵⁵5Note that we perform the calculation in the axial gauge. In other gauges, the distribution of terms between diagrams can differ, but the sum of diagrams should not change. In particular, we also calculated the background Feynman gauge with the same result, see Appendix C. The singularities are regulated by the factorization cut-offs. As discussed in the section I, we introduce these cut-offs using the renormalization approach. Now, let’s understand how this procedure works.

The contribution of the last term in Eq. (29) has a simple form,

$\displaystyle\mathcal{B}^{\rm q(1)+bg;real,c}_{ij}(x_{B},b_{\perp})=4\alpha_{s}N_{c}\int{\textstyle d}\lower 0.12915pt\hbox{\kern-3.80005pt${}^{\scriptstyle-}$}\kern-0.50003pt{}^{2}p_{\perp}e^{ip_{\perp}b_{\perp}}\int^{1}_{0}\frac{dz}{z(1-z)}\int\frac{{\textstyle d}\lower 0.12915pt\hbox{\kern-3.80005pt${}^{\scriptstyle-}$}\kern-0.50003pt{}^{2}k_{\perp}}{k^{2}_{\perp}}\int d^{2}z_{\perp}e^{i(k-p)_{\perp}z_{\perp}}\mathcal{B}^{\rm bg}_{ij}(\frac{x_{B}}{z},z_{\perp})\,.$

(33)

Consider with the rapidity divergence first. While the contribution does not contain a divergence as $z\to 0$, which is regulated by a non-zero value of $x_{B}$, there is a divergence as $z\to 1$. In terms of the $k^{-}$ variable, this divergence corresponds to $k^{-}\to\infty$. To separate this divergence, we rewrite the contribution as

	$\displaystyle\mathcal{B}^{\rm q(1)+bg;real,c}_{ij}\left(x_{B},b_{\perp}\right)$		(34)
	$\displaystyle=4\alpha_{s}N_{c}\int\frac{{\textstyle d}\lower 0.12915pt\hbox{\kern-3.80005pt${}^{\scriptstyle-}$}\kern-0.50003pt{}^{2}k_{\perp}}{k^{2}_{\perp}}e^{ik_{\perp}b_{\perp}}\Big{(}\int^{1}_{0}\frac{dz}{(1-z)_{+}}\mathcal{B}^{\rm bg}_{ij}(\frac{x_{B}}{z},b_{\perp})+\int^{1}_{0}\frac{dz}{z}\mathcal{B}^{\rm bg}_{ij}(\frac{x_{B}}{z},b_{\perp})+\int^{1}_{0}\frac{dz}{1-z}\mathcal{B}^{\rm bg}_{ij}(x_{B},b_{\perp})\Big{)}\,,$

where the plus distribution is defined in the usual way,

$\displaystyle(f(z))_{+}\equiv f(z)-\delta(1-z)\int^{1}_{0}dz^{\prime}f(z^{\prime}).$

(35)

Now, the rapidity divergence is in the last term of Eq. (34), which we regulate by making the replacement (8). In terms of the $z$ variable, this corresponds to

$\displaystyle\int^{1}_{0}\frac{dz}{1-z}\to\Big{(}\frac{\nu}{x_{B}P^{+}}\Big{)}^{\eta}\int^{1}_{0}\frac{dz}{1-z}\Big{(}\frac{1-z}{z}\Big{)}^{-\eta}\,,$

(36)

where the factorization scale $\nu$ should be understood as the upper cut-off in rapidity. Performing integration over $z$ we obtain

	$\displaystyle\mathcal{B}^{\rm q(1)+bg;real,c}_{ij}(x_{B},b_{\perp})=4\alpha_{s}N_{c}\int\frac{{\textstyle d}\lower 0.12915pt\hbox{\kern-3.80005pt${}^{\scriptstyle-}$}\kern-0.50003pt{}^{2}k_{\perp}}{k^{2}_{\perp}}e^{ik_{\perp}b_{\perp}}$		(37)
	$\displaystyle\times\Big{(}\int^{1}_{0}\frac{dz}{(1-z)_{+}}\mathcal{B}^{\rm bg}_{ij}(\frac{x_{B}}{z},b_{\perp})+\int^{1}_{0}\frac{dz}{z}\mathcal{B}^{\rm bg}_{ij}(\frac{x_{B}}{z},b_{\perp})-\Big{(}\frac{\nu}{x_{B}P^{+}}\Big{)}^{\eta}\Gamma(1-\eta)\Gamma(\eta)\mathcal{B}^{\rm bg}_{ij}(x_{B},b_{\perp})\Big{)}\,.$

Now we notice that apart from the rapidity divergence, there is a divergence in the integral over transverse momenta. We regulate it using the dimensional regularization,

	$\displaystyle\mathcal{B}^{\rm q(1)+bg;real,c}_{ij}(x_{B},b_{\perp})=\frac{\alpha_{s}N_{c}}{\pi}\Big{(}\frac{\mu^{2}b^{2}_{\perp}}{4e^{-2\gamma_{E}}}\Big{)}^{\epsilon}e^{-\epsilon\gamma_{E}}\Gamma(-\epsilon)$		(38)
	$\displaystyle\times\Big{(}\int^{1}_{0}\frac{dz}{(1-z)_{+}}\mathcal{B}^{\rm bg}_{ij}(\frac{x_{B}}{z},b_{\perp})+\int^{1}_{0}\frac{dz}{z}\mathcal{B}^{\rm bg}_{ij}(\frac{x_{B}}{z},b_{\perp})-\Big{(}\frac{\nu}{x_{B}P^{+}}\Big{)}^{\eta}\Gamma(1-\eta)\Gamma(\eta)\mathcal{B}^{\rm bg}_{ij}(x_{B},b_{\perp})\Big{)}\,,$

where $\mu$ is the $\overline{\rm{MS}}$ scale,

$\displaystyle\mu^{2}\equiv\frac{4\pi}{e^{\gamma_{E}}}\mu^{2}_{0}\,.$

(39)

Expanding in powers of $\epsilon$ and introducing the notation,

$\displaystyle L^{\mu}_{b}\equiv\ln\Big(\frac{b^{2}_{\perp}\mu^{2}}{4e^{-2\gamma_{E}}}\Big{missing})\,,$

(40)

one can easily obtain

	$\displaystyle\mathcal{B}^{\rm q(1)+bg;real,c}_{ij}(x_{B},b_{\perp})=-\frac{\alpha_{s}N_{c}}{\pi}\Big{(}\frac{1}{\epsilon}+L^{\mu}_{b}\Big{)}$		(41)
	$\displaystyle\times\Big{(}\int^{1}_{0}\frac{dz}{(1-z)_{+}}\mathcal{B}^{\rm bg}_{ij}(\frac{x_{B}}{z},b_{\perp})+\int^{1}_{0}\frac{dz}{z}\mathcal{B}^{\rm bg}_{ij}(\frac{x_{B}}{z},b_{\perp})-\Big{(}\frac{\nu}{x_{B}P^{+}}\Big{)}^{\eta}\Gamma(1-\eta)\Gamma(\eta)\mathcal{B}^{\rm bg}_{ij}(x_{B},b_{\perp})\Big{)}\,.$

Also expanding in $\eta$, the equation reads

	$\displaystyle\mathcal{B}^{\rm q(1)+bg;real,c}_{ij}(x_{B},b_{\perp})=-\frac{\alpha_{s}N_{c}}{\pi}\Big{(}\frac{1}{\epsilon}+L^{\mu}_{b}\Big{)}$		(42)
	$\displaystyle\times\Big{(}\int^{1}_{0}\frac{dz}{(1-z)_{+}}\mathcal{B}^{\rm bg}_{ij}(\frac{x_{B}}{z},b_{\perp})+\int^{1}_{0}\frac{dz}{z}\mathcal{B}^{\rm bg}_{ij}(\frac{x_{B}}{z},b_{\perp})-\left(\frac{1}{\eta}+\ln(\frac{\nu}{x_{B}P^{+}})\right)\mathcal{B}^{\rm bg}_{ij}(x_{B},b_{\perp})\Big{)}\,.$