&&= - \frac{g^2N_c}{2\pi} \int^\infty_0 \frac{dk^-}{k^-} (y_\perp| \frac{  2x_BP^+k^- }{ p^2_\perp(2x_BP^+k^- + p^2_\perp)} |y_\perp) \Big( \int^{\infty}_{-\infty} dz^- e^{i x_BP^+ z^-} [\infty, z^-]^{an}_yF^n_{-j}(z^-, y_\perp)\Big)
\nonumber\\
&&+ \frac{ig^2N_c}{2\pi}\int^\infty_0 \frac{dk^-}{k^-} (y_\perp| \frac{p^s}{p^2_\perp} (2\delta^k_j \delta^m_s - g_{sj} g^{mk})\partial_k \Big( \int^{\infty}_{-\infty} dz^- e^{i x_BP^+ z^-} [\infty, z^-]^{an}F^n_{-m}(z^-)\Big)  \frac{ 1 }{ 2x_BP^+k^- + p^2_\perp} |y_\perp) \,,
\nonumber
\end{eqnarray}

which we can rewrite as,

\begin{eqnarray}
&&\langle \int^\infty_{-\infty} dy^- e^{i x_BP^+y^-} [\infty, y^-]^{an}_yF^{{\rm q+bg};n}_{-j}(y^-, y_\perp)\rangle_{B^{\rm bg}} = - \frac{ g^2N_c}{2\pi}  \int^\infty_0\frac{d k^-}{k^-} \int \dhd^2k_\perp \frac{ 2x_BP^+k^- }{k^2_\perp(2x_BP^+k^- + k^2_\perp)}
\nonumber\\
&&\times \int^{\infty}_{-\infty} dz^- e^{i x_BP^+ z^-} [\infty, z^-]^{an}_y F^n_{-j}(z^-, y_\perp) + \frac{g^2 N_c}{2\pi}  \int^\infty_0\frac{d k^-}{k^-} \int \dhd^2 p_\perp \int \dhd^2 k_\perp e^{i (p_\perp-k_\perp ) y_\perp } \frac{k^s}{k^2_\perp} (2 \delta^k_j \delta^m_s  - g_{js}g^{mk})
\nonumber\\
&&\times \frac{ (k - p)_k }{ 2x_BP^+k^- + p^2_\perp}
\int d^2z_\perp e^{i ( k_\perp - p_\perp ) z_\perp } \int^{\infty}_{-\infty} dz^- e^{i x_BP^+ z^-} [\infty, z^-]^{an}_z F^n_{-m}(z^-, z_\perp) \,.
\end{eqnarray}

A similar result can be obtained for the diagrams conjugated to the diagrams presented in Figs. <ref>a-d. Taking a sum of these virtual diagrams we obtain

\begin{eqnarray}
&&\int^\infty_{-\infty} dz^- e^{-i x_BP^+ z^-}   \langle p | \tilde{F}^m_{-i}(z^-, b_\perp) [z^-, \infty]^{ma}_b [\infty, 0^-]^{an}_0 F^n_{-j}(0^-, 0_\perp) |p\rangle^{\rm virt}
\nonumber\\
&&= -\frac{g^2N_c}{2\pi} \int^\infty_0 \frac{d k^-}{k^-}\Big\{ - \int \dhd^2p_\perp e^{ ip_\perp b_\perp} \int \dhd^2k_\perp e^{- ik_\perp b_\perp} \frac{k^s}{k^2_\perp} (2\delta^k_j\delta^m_s-g_{js}g^{mk}) \frac{(k-p)_k}{2x_B P^+ k^- + p^2_\perp}
\nonumber\\
&&\times \int d^2z_\perp e^{i (k_\perp - p_\perp ) z_\perp } \int^\infty_{-\infty} dz^- e^{-i x_BP^+ z^-}   \langle p | \tilde{F}^m_{-i}(z^-, z_\perp) [z^-, \infty]^{ma}_z [\infty, 0^-]^{an}_0 F^n_{-m}(0^-, 0_\perp) |p\rangle
\nonumber\\
&&+ \int \dhd^2k_\perp \frac{2x_B P^+ k^-}{k^2_\perp(2x_B P^+ k^- + k^2_\perp)} \int^\infty_{-\infty} dz^- e^{-i x_BP^+ z^-}   \langle p | \tilde{F}^m_{-i}(z^-, b_\perp) [z^-, \infty]^{ma}_b [\infty, 0^-]^{an}_0 F^n_{-j}(0^-, 0_\perp) |p\rangle \Big\}
\nonumber\\
&&-\frac{g^2N_c}{2\pi} \int^\infty_0 \frac{dk^-}{k^-} \Big\{ \int \dhd^2p_\perp e^{ip_\perp b_\perp} \int \dhd^2k_\perp e^{-ik_\perp b_\perp} \frac{(p-k)_k}{2x_B P^+ k^- + p^2_\perp} (2\delta^k_i\delta^l_s - g_{is}g^{kl}) \frac{k^s}{k^2_\perp}
\nonumber\\
&&\times \int d^2z_\perp e^{i (k_\perp-p_\perp ) z_\perp} \int^\infty_{-\infty} dz^- e^{-i x_BP^+ z^-}   \langle p | \tilde{F}^m_{-l}(z^-, z_\perp) [z^-, \infty]^{ma}_z [\infty, 0^-]^{an}_0 F^n_{-j}(0^-, 0_\perp) |p\rangle
\nonumber\\
&&+ \int \dhd^2k_\perp \frac{2x_B P^+ k^-}{k^2_\perp(2x_B P^+ k^- + k^2_\perp)}  \int^\infty_{-\infty} dz^- e^{-i x_BP^+ z^-}   \langle p | \tilde{F}^m_{-i}(z^-, b_\perp) [z^-, \infty]^{ma}_b [\infty, 0^-]^{an}_0 F^n_{-j}(0^-, 0_\perp) |p\rangle \Big\} \,.
\label{virt:kmin}
\end{eqnarray}

Introducing the variable $z$, see Eq.  (<ref>), we rewrite this equation in the following form

\begin{eqnarray}
&&\int^\infty_{-\infty} dz^- e^{-i x_BP^+ z^-}   \langle p | \tilde{F}^m_{-i}(z^-, b_\perp) [z^-, \infty]^{ma}_b [\infty, 0^-]^{an}_0 F^n_{-j}(0^-, 0_\perp) |p\rangle^{\rm virt} = \frac{g^2N_c}{2\pi} \int^1_0 \frac{dz}{z}\int \dhd^2p_\perp e^{ ip_\perp b_\perp}
\nonumber\\
&&\times \int \dhd^2k_\perp e^{- ik_\perp b_\perp} \Big[ \delta^l_i k^s (2\delta^k_j\delta^m_s-g_{js}g^{mk}) (k-p)_k + (k - p)_k (2\delta^k_i\delta^l_s - g_{is}g^{kl}) \delta^m_j k^s \Big] \frac{1}{k^2_\perp(z k^2_\perp + (1-z)p^2_\perp)}
\nonumber\\
&&\times   \int d^2z_\perp e^{i ( k_\perp - p_\perp) z_\perp} \int^\infty_{-\infty} dz^- e^{-i x_BP^+ z^-}   \langle p | \tilde{F}^m_{-l}(z^-, z_\perp) [z^-, \infty]^{ma}_z [\infty, 0^-]^{an}_0 F^n_{-m}(0^-, 0_\perp) |p\rangle
\nonumber\\
&& - \frac{g^2N_c}{\pi} \int^1_0 \frac{dz}{1 - z} \int \frac{ \dhd^2k_\perp }{k^2_\perp }  \int^\infty_{-\infty} dz^- e^{-i x_BP^+ z^-}   \langle p | \tilde{F}^m_{-i}(z^-, b_\perp) [z^-, \infty]^{ma}_b [\infty, 0^-]^{an}_0 F^n_{-j}(0^-, 0_\perp) |p\rangle \,.
\end{eqnarray}

The last line of this equation contains a scaleless integral over transverse momentum which is set to be zero in the dimensional regularization. However, as we discuss in the main text, the role of this contribution is non-trivial and we keep it for now. As a result, the equation reads

\begin{eqnarray}
&&\int^\infty_{-\infty} dz^- e^{-i x_BP^+ z^-}   \langle p | \tilde{F}^m_{-i}(z^-, b_\perp) [z^-, \infty]^{ma}_b [\infty, 0^-]^{an}_0 F^n_{-j}(0^-, 0_\perp) |p\rangle^{\rm virt} = \frac{g^2N_c}{2\pi} \int^1_0 \frac{dz}{z}\int \dhd^2p_\perp e^{ ip_\perp b_\perp}
\nonumber\\
&&\times \int \dhd^2k_\perp e^{- ik_\perp b_\perp} \Big[ \delta^l_i k^s (2\delta^k_j\delta^m_s-g_{js}g^{mk}) (k-p)_k + (k - p)_k (2\delta^k_i\delta^l_s - g_{is}g^{kl}) \delta^m_j k^s \Big] \frac{1}{k^2_\perp(z k^2_\perp + (1-z)p^2_\perp)}
\nonumber\\
&&\times   \int d^2z_\perp e^{i ( k_\perp - p_\perp) z_\perp} \int^\infty_{-\infty} dz^- e^{-i x_BP^+ z^-}   \langle p | \tilde{F}^m_{-l}(z^-, z_\perp) [z^-, \infty]^{ma}_z [\infty, 0^-]^{an}_0 F^n_{-m}(0^-, 0_\perp) |p\rangle
\nonumber\\
&&- \frac{g^2N_c}{\pi} \int^1_0 \frac{dz}{1 - z} \int \frac{ \dhd^2k_\perp }{k^2_\perp }  \int^\infty_{-\infty} dz^- e^{-i x_BP^+ z^-}   \langle p | \tilde{F}^m_{-i}(z^-, b_\perp) [z^-, \infty]^{ma}_b [\infty, 0^-]^{an}_0 F^n_{-j}(0^-, 0_\perp) |p\rangle \,.
\label{virt:z}
\end{eqnarray}

§ CALCULATION USING BACKGROUND-FEYNMAN GAUGE

In this section, we shall repeat the above calculation with the choice of “quantum" field gauge to be,

\begin{eqnarray}
&&D^{\mu}B^{\rm q}_{\mu}(x^-,x_{\perp})=0\,,
\end{eqnarray}

where the covariant derivative is constructed from the background field $B^{\rm bg}_\mu$. This gauge is called the background-Feynman gauge. For this gauge choice, the propagator takes the form,

\begin{eqnarray}
&&i\langle B^{{\rm q} a}_{\mu}(x) B^{{\rm q} b}_{\nu}(y) \rangle_{B^{\rm bg}}=(x|\frac{1}{ P^2 g^{\mu\nu}+ 2i F^{\mu\nu}+i\epsilon}|y)^{ab}\,,
\end{eqnarray}

where $P_\mu\:=p_\mu + g B^{\rm bg}_\mu$. To simplify the calculation, in this section we'll also assume that the background field has only $B^{\rm bg}_-(x^-,x_{\perp})$ non-zero component.

Up to two gluon accuracy, the propagator can be written in the following form,

\begin{eqnarray}
&&i\langle B_{\mu}^{{\rm q}a}(x)B_{\nu}^{{\rm q} b}(y)\rangle_{B^{\rm bg}}\Big|_{x^- > y^-}
\label{feyn:propagator}\\
&&=-\frac{1}{2\pi}\int_0^{\infty} \frac{dp^-}{2p^-}e^{-ip^- (x-y)^+}(x_{\perp}|e^{-i \frac{p_{\perp}^2 }{ 2p^-}x^-}\Big(i g_{\mu\nu}-gg_{\mu\nu}\int_{y^-}^{x^-} dz^- e^{i\frac{p_{\perp}^2 }{ 2p^-} z^-} B^{\rm bg}_-(z^-)e^{-i\frac{p_{\perp}^2 }{ 2p^-} z^-}
\nonumber\\
&&- \frac{ig}{p^-}\int_{y^-}^{x^-} dz^- e^{i\frac{p_{\perp}^2 }{ 2p^-} z^-} F_{\mu\nu}(z^-) e^{-i\frac{p_{\perp}^2 }{ 2p^-} z^-} + \frac{gg_{\mu+}g_{\nu+} }{ 2 (p^-)^2}\int_{y^-}^{x^-} dz^- e^{i \frac{p_{\perp}^2 }{ 2 p^-}z^-} \partial^k F_{-k}e^{-i \frac{p_{\perp}^2 }{ 2 p^-}z^-}\Big) e^{i\frac{p_{\perp}^2 }{ 2p^-}y^-}|y_{\perp})^{ab}\,,
\nonumber
\end{eqnarray}

where the third term with $g_{\mu +}\:g_{\nu +}$ is the contribution from the quark background fields $\partial^k F^a_{-k} = g\bar{\psi}\gamma^+ t^a \psi$. The propagator which shall be used for real emissions takes the form,

\begin{eqnarray}
\label{feyn:propagator momentum}
&&i\lim_{k^2\rightarrow 0} k^2\langle B_{\mu}^{{\rm q}a}(k) B_{\nu}^{{\rm q} b}(y)\rangle = -ie^{ik^- y^+}(k_{\perp}|\Big(i g_{\mu\nu}-gg_{\mu\nu}\int_{y^-}^{\infty} dz^- e^{i\frac{p_{\perp}^2}{2k^-} z^-}B^{\rm bg}_-(z^-)\:e^{-i\frac{p_{\perp}^2}{2k^-} z^-}
\\
&&- \frac{ig}{k^-}\int_{y^-}^{\infty} dz^- e^{i\frac{p_{\perp}^2}{2k^-} z^-} F_{\mu\nu}(z^-) e^{-i\frac{p_{\perp}^2}{ 2k^-} z^-} + \frac{gg_{\mu+}g_{\nu+}}{ 2 (k^-)^2} \int_{y^-}^{\infty} dz^- e^{i \frac{p_{\perp}^2 }{ 2 k^-}z^-} \partial^k F_{-k}e^{-i \frac{p_{\perp}^2 }{ 2 k^-}z^-}\Big) e^{i\frac{p_{\perp}^2 }{ 2k^-}y^-}|y_{\perp})^{ab} \,.
\nonumber
\end{eqnarray}

First, we will compute the real emission vertex, see Fig. <ref>, using the above gauge and match it with (<ref>). Apart for the gauge condition (<ref>), we will also fix $B^{\rm bg}_i=0$, which means that in this calculation we only track the $\partial_i B^{\rm bg}_-$ component of the $F_{-i}$ strength tensor, which is sufficient to obtain the final result (<ref>).

Writing (<ref>) slightly differently, we have

\begin{eqnarray}
&& L_{\mu j}^{ab}(k,y_{\perp},x_B)=\int_{-\infty}^{\infty}dy^- e^{ix_B P^+ y^-}\Big(-i x_B P^+ [\infty,y^-]^{bd}i\lim_{k^2\rightarrow 0}k^2\langle B_{\mu}^{{\rm q}a}(k)B_j^{{\rm q}d}(y^-,y_{\perp})\rangle_{B^{\rm bg}}
\label{feyn:initial}
\\
&&- [\infty,y^-]^{bd} i\lim_{k^2\rightarrow 0}k^2\langle B^{{\rm q} a}_{\mu}(k) \partial_j B^{{\rm q} d}_-(y^-,y_{\perp})\rangle_{B^{\rm bg}}
\nonumber\\
&&+ ig \int_{y^-}^{\infty} dz^- \left([\infty,z^-] T^d [z^-,y^-]\right)^{bm} F^m_{-j}(y^-,y_{\perp}) i\lim_{k^2\rightarrow 0}k^2\langle B^{{\rm q}a}_{\mu}(k)B_-^{{\rm q}d}(z^-,y_{\perp})\rangle \Big)\,,
\nonumber
\end{eqnarray}

where the last propagator is free, without any background fields.

Substituting (<ref>) into the expression (<ref>), we have

\begin{eqnarray}
&& L_{\mu j}^{ab}(k,y_{\perp},x_B)= -i x_B P^+ \int_{-\infty}^{\infty}dy^- e^{i x_B P^+ y^-}(k_{\perp}|\Big( g_{\mu j} + i g g_{\mu j} \int_{y^-}^{\infty} dz^- e^{i \frac{p_{\perp}^2}{ 2 k^-} z^-} B^{\rm bg}_{-}(z^-) e^{-i\frac{p_{\perp}^2 }{ 2k^-}z^-}
\\
&&- \frac{ g }{ k^-} \int_{y^-}^{\infty} dz^- e^{i \frac{p_{\perp}^2 }{ 2 k^-}z^-} F_{\mu j}(z^-) e^{-i \frac{p_{\perp}^2 }{ 2k^-}z^-}\Big) e^{i \frac{p_{\perp}^2 }{ 2k^-}y^-} |y_{\perp})^{ab}
+ g g_{\mu j} x_B P^+ \int_{-\infty}^{\infty}dy^- e^{i x_B P^+ y^-}
\nonumber\\
&&\times\int_{y^-}^{\infty} dz^- B_{-}^{{\rm bg}ba}(z^-) e^{i \frac{k_{\perp}^2 }{ 2k^-}y^-} e^{-i k_{\perp}y_{\perp}}
- \int_{-\infty}^{\infty}dy^- e^{i x_B P^+ y^-}(k_{\perp}|\Big(g_{\mu-}+i g g_{\mu-}\int_{y^-}^{\infty} dz^- e^{i \frac{p_{\perp}^2 }{ 2k^-}z^-} B^{\rm bg}_-(z^-)e^{-i \frac{p_{\perp}^2 }{ 2k^-}z^-}
\nonumber\\
&&- \frac{g }{ k^-}\int_{y^-}^{\infty} dz^- e^{i \frac{p_{\perp}^2 }{ 2 k^-}z^-} F_{\mu-}(z^-)e^{-i \frac{p_{\perp}^2 }{ 2 k^-}z^-}
- \frac{i g g_{\mu+}}{2 (k^-)^2} \int_{-\infty}^{\infty}dy^- e^{i x_B P^+ y^-}\int_{y^-}^{\infty} dz^- e^{i \frac{p_{\perp}^2 }{ 2 k^-}z^-} \partial^k F_{-k}(z^-) e^{-i \frac{p_{\perp}^2 }{ 2k^-}z^-}\Big)
\nonumber\\
&&\times e^{i \frac{p_{\perp}^2 }{ 2k^-}y^-}(ip_j) |y_{\perp})^{ab} - ig g_{\mu-} \int_{-\infty}^{\infty}dy^- e^{i x_B P^+ y^-} \int_{y^-}^{\infty} dz^- B^{{\rm bg}ab}_-(z^-) e^{i \frac{k_{\perp}^2 }{ 2 k^-}y^-} e^{-i k_{\perp}y_{\perp}}  (ik_j)
\nonumber\\
&&+ ig g_{\mu -}\int_{-\infty}^{\infty}dy^- e^{i x_B P^+ y^-} \int_{y^-}^{\infty} dz^- F_{-j}^{ab}(y^-,y_{\perp}) e^{i \frac{k_{\perp}^2 }{ 2 k^-} z^-} e^{-i k_{\perp}y_{\perp}}\,.
\nonumber
\end{eqnarray}

Carrying out the integral over $y^-$ and simplifying, we have

\begin{eqnarray}
&& L_{\mu j}^{ab}(k,y_{\perp},x_B) = -ig g_{\mu j} \int_{-\infty}^{\infty} dz^- e^{i(x_BP^+ + \frac{k_{\perp}^2 }{ 2k^-})z^-}(k_{\perp}|B^{\rm bg}_-(z^-)\frac{2 x_B P^+ k^-}{ 2x_B P^+k^- + p_{\perp}^2}|y_{\perp})^{ab} + \frac{g }{ k^-}
\\
&&\times \int_{-\infty}^{\infty} dz^- e^{i(x_BP^+ + \frac{k_{\perp}^2 }{ 2k^-})z^-} (k_{\perp}|F_{\mu j}(z^-) \frac{2 x_B P^+ k^-}{ 2 x_B P^+ k^- + p_{\perp}^2}|y_{\perp})^{ab}
-ig g_{\mu j} \int_{-\infty}^{\infty} dz^- e^{i(x_BP^+ + \frac{k_{\perp}^2 }{ 2k^-})z^-} B_-^{{\rm bg}ab}(z^-)
\nonumber\\
&&\times \frac{2x_B P^+ k^-}{ 2 x_B P^+ k^-+ k_{\perp}^2} e^{-ik_{\perp}y_{\perp}} -ig g_{\mu -} \int_{-\infty}^{\infty} dz^- e^{i(x_B P^+ + \frac{k_{\perp}^2 }{ 2k^-})z^-}(k_{\perp}|B^{\rm bg}_-(z^-) \frac{2k^- p_j }{ 2  x_B P^+ k^- + p_{\perp}^2}|y_{\perp})^{ab}
\nonumber\\
&&+\frac{g }{ k^-} \int_{-\infty}^{\infty} dz^- e^{i(x_B P^+ + \frac{k_{\perp}^2 }{ 2k^-})z^-}(k_{\perp}|F_{\mu-}(z^-) \frac{2p_j k^- }{ 2x_B P^+ k^-+ p_{\perp}^2}|y_{\perp})^{ab}
+ \frac{ig g_{\mu+} }{ 2 (k^-)^2} \int_{-\infty}^{\infty} dz^- e^{i(x_B P^+ + \frac{k_{\perp}^2 }{ 2k^-})z^-}
\nonumber\\
&&\times(k_{\perp}|\partial_k F_{-k}(z^-)\frac{2p_jk^- }{ 2 x_B P^+ k^-+ p_{\perp}^2}|y_{\perp})^{ab} +ig g_{\mu-} \int_{-\infty}^{\infty} dz^- e^{i(x_B P^+ + \frac{k_{\perp}^2 }{ 2k^-})z^-} B_{-}^{{\rm bg}ab}(z^-)\frac{2 k^- k_j }{ 2 x_B P^+ k^- + k_{\perp}^2} e^{-i k_{\perp}y_{\perp}}
\nonumber\\
&& - \frac{2g g_{\mu-}k^- }{ k_{\perp}^2}\int_{-\infty}^{\infty} dz^- e^{i(x_B P^+ + \frac{k_{\perp}^2 }{ 2k^-})z^-} F_{-j}^{ab}(z^-)e^{-i k_{\perp}y_{\perp}}\,.
\nonumber
\end{eqnarray}

Using $p_i \mathcal{O}\:=\:i\partial_i \mathcal{O}+\mathcal{O}p_i$ to simplify further,

\begin{eqnarray}
&&L_{\mu j}^{ab}(k,y_{\perp},x_B) = g \int_{-\infty}^{\infty} dz^- e^{i(x_B P^+ + \frac{k_{\perp}^2 }{ 2k^-})z^-}(k_{\perp}|\frac{ 2 x_B P^+k^- g_{\mu j} + 2g_{\mu -}k^- k_j }{ 2 x_B P^+ k^- + k_{\perp}^2}\partial_k B^{\rm bg}_-(z^-)\frac{ p_k + k_k }{ 2x_B P^+k^- + p_{\perp}^2}|y_{\perp})^{ab}
\nonumber\\
&&+ \frac{g }{ k^-} \int_{-\infty}^{\infty} dz^- e^{i(x_BP^+ + \frac{k_{\perp}^2 }{ 2k^-})z^-}(k_{\perp}|F_{\mu j}(z^-)\frac{2 x_B P^+ k^- }{ 2 x_B P^+ k^-+ p_{\perp}^2}|y_{\perp})^{ab} -g g_{\mu-} (k_{\perp}|\int_{-\infty}^{\infty} dz^- e^{i(x_BP^+ + \frac{k_{\perp}^2 }{ 2k^-})z^-}
\nonumber\\
&&\times \partial_j B^{\rm bg}_{-}(z^-)\frac{2k^- }{ 2x_B P^+k^- + p_{\perp}^2}|y_{\perp})^{ab} + \frac{gg_{\mu+} }{ 2 (k^-)^2} \int_{-\infty}^{\infty} dz^- e^{i(x_B P^+ + \frac{k_{\perp}^2 }{ 2k^-})z^-} (k_{\perp}|(k^k F_{-k}(z^-) - F_{-k}(z^-)p^k)
\nonumber\\
&&\times \frac{2p_j k^- }{ 2 x_B P^+ k^- + p_{\perp}^2}|y_{\perp})^{ab}
+ \frac{g }{ k^-} \int_{-\infty}^{\infty} dz^- e^{i(x_B P^+ + \frac{k_{\perp}^2 }{ 2k^-})z^-}(k_{\perp}|F_{\mu-}(z^-) \frac{2p_jk^- }{ 2x_B P^+ k^- + p_{\perp}^2}|y_{\perp})^{ab}
\nonumber\\
&&- \frac{2g g_{\mu-}k^- }{ k_{\perp}^2}(k_{\perp}|\int_{-\infty}^{\infty} dz^- e^{i(x_B P^+ + \frac{k_{\perp}^2 }{ 2k^-})z^-} F_{-j}(z^-)|y_{\perp})^{ab}\,.
\end{eqnarray}

Now, replacing $\partial_j B^{\rm bg}_{-} \to -F_{-j}$ and restoring the Wilson lines, we obtain

\begin{eqnarray}
&&L_{\mu j}^{ab}(k,y_{\perp},x_B)
\\
&&=-g (k_{\perp}|\frac{2 x_B P^+k^- g_{\mu j} + 2g_{\mu -}k^- k_j }{ 2x_B P^+k^- + k_{\perp}^2}\Big( \int^\infty_{-\infty} dz^- e^{i (x_BP^+ + \frac{k^2_\perp}{2k^-})z^-} [\infty, z^-]^{ae}F^{eb}_{-k}(z^-) \Big)\frac{ p_k + k_k }{ 2x_B P^+k^- + p_{\perp}^2}|y_{\perp})
\nonumber\\
&&+ \frac{gg_{\mu +}}{ k^-}(k_{\perp}|\Big( \int^\infty_{-\infty} dz^- e^{i (x_BP^+ + \frac{k^2_\perp}{2k^-})z^-} [\infty, z^-]^{ae}F^{eb}_{-j}(z^-) \Big)\frac{2k^- x_B P^+ }{ 2x_B P^+k^-  + p_{\perp}^2}|y_{\perp})
\nonumber\\
&&+g g_{\mu-}(k_{\perp}|\Big( \int^\infty_{-\infty} dz^- e^{i (x_BP^+ + \frac{k^2_\perp}{2k^-})z^-} [\infty, z^-]^{ae}F^{eb}_{-j}(z^-) \Big)\frac{2k^- }{ 2x_B P^+k^- + p_{\perp}^2}|y_{\perp})^{ab}
\nonumber\\
&&+ \frac{gg_{\mu+} }{ 2 (k^-)^2}(k_{\perp}|\Big( \int^\infty_{-\infty} dz^- e^{i (x_BP^+ + \frac{k^2_\perp}{2k^-})z^-} [\infty, z^-]^{ae}F^{eb}_{-k}(z^-) \Big)\frac{2p_jk^- (k^k - p^k) }{ 2k^- x_B P^+ + p_{\perp}^2}|y_{\perp})
\nonumber\\
&&- \frac{g }{ k^-}(k_{\perp}|\Big( \int^\infty_{-\infty} dz^- e^{i (x_BP^+ + \frac{k^2_\perp}{2k^-})z^-} [\infty, z^-]^{ae}F^{eb}_{-\mu}(z^-) \Big)\frac{2p_jk^- }{ 2k^-x_B P^+ + p_{\perp}^2}|y_{\perp})
\nonumber\\
&&-\frac{2gg_{\mu-}k^- }{ k_{\perp}^2}(k_{\perp}|\Big( \int^\infty_{-\infty} dz^- e^{i (x_BP^+ + \frac{k^2_\perp}{2k^-})z^-} [\infty, z^-]^{ae}F^{eb}_{-j}(z^-) \Big)|y_{\perp}) \,.
\nonumber
\end{eqnarray}

The above expression can be checked to be gauge covariant. As a result, since we aim to calculate a product (<ref>), we can make a replacement

\begin{eqnarray}
&&g_{\mu -}k^- \to g_{\mu -}k^- - k_{\mu}\:=\:-k^+g_{\mu+}-k^\perp_ \mu\,,
\end{eqnarray}

since one can add a term proportional to $k^{\mu}$ without any cost. We also notice that the terms proportional to $g_{\mu+}$ do not contribute to the square of two emission vertex (<ref>). Removing such terms, we get

\begin{eqnarray}
&&L_{k j}^{ab}(k,y_{\perp},x_B) = -2g(k_{\perp}|\Big( \int^\infty_{-\infty} dz^- e^{i (x_BP^+ + \frac{k^2_\perp}{2k^-})z^-} [\infty, z^-]^{ae}F^{eb}_{-m}(z^-) \Big)
\\
&&\times\Big(\frac{  x_B P^+ k^- g_{k j} - k_k k_j}{ 2 x_B P^+k^- + k_{\perp}^2}\frac{ p_m + k_m }{ 2x_B P^+k^- + p_{\perp}^2} - \frac{k_k g_{mj} + p_j g_{mk} }{ 2x_B P^+k^- + p_{\perp}^2} + \frac{k_k g_{mj} }{ k_{\perp}^2} |y_{\perp})
\nonumber
\end{eqnarray}

which is same as (<ref>).Now we calculate the virtual correction. To do that, we start with

\begin{eqnarray}
&&\langle \int^\infty_{-\infty} dy^- e^{i x_BP^+y^-} [\infty, y^-]^{an}_yF^{{\rm q+bg};n}_{-j}(y^-, y_\perp)\rangle_{B^{\rm bg}}
\nonumber\\
&&= gx_B P^+\int_{-\infty}^{\infty} dy^-e^{ix_B P^+ y^-} \int_{y^-}^{\infty} dz^- ([\infty,z^-]T^e[z^-,y^-])^{an} \langle B_{-}^{{\rm q}e}(z^-) B_j^{{\rm q}n}(y^-,y_{\perp})\rangle_{B^{\rm bg}}
\\
&&-ig \int_{-\infty}^{\infty} dy^- e^{ix_B P^+ y^-} \int_{y^-}^{\infty} dz^- ([\infty,z^-]T^e[z^-,y^-])^{an} \langle B_-^{{\rm q}e}(z^-)\partial_j B_-^{{\rm q}n}(y^-,y_{\perp}) \rangle_{B^{\rm bg}}\,.
\nonumber
\end{eqnarray}

Note, that a diagram in Fig. <ref>c doesn't contribute in the background-Feynman gauge.

The explicit propagators (<ref>) which contribute are

\begin{eqnarray}
&&i\langle B_-^{{\rm q}e}(z^-,y_{\perp}) B_j^{{\rm q}n}(y^-,y_{\perp}) \rangle_{B^{\rm bg}}
\\
&&= \frac{ig}{ 2\pi} \int_0^{\infty} \frac{dp^- }{ 2 (p^-)^2} (y_{\perp}|e^{-i\frac{p_{\perp}^2 }{ 2p^-}z^-}
\Big(\int_{y^-}^{z^-} dx^-e^{i\frac{p_{\perp}^2 }{ 2p^-}x^-} F_{-j}(x^-) e^{-i\frac{p_{\perp}^2 }{ 2p^-}x^-}\Big) e^{i\frac{p_{\perp}^2 }{ 2p^-}y^-} |y_{\perp})^{en} \,,
\nonumber
\end{eqnarray}

\begin{eqnarray}
&&i \langle B_-^{{\rm q}e}(z^-,y_{\perp}) \partial_j B_-^{{\rm q}n}(y^-,y_{\perp}) \rangle_{B^{\rm bg}}
\\
&&= -\frac{g }{ 2\pi} \int_0^{\infty} \frac{dp^- }{ 4 (p^-)^3} (y_{\perp}|e^{-i \frac{p_{\perp}^2 }{ 2p^-}z^-}
\Big(\int_{y^-}^{z^-} dx^-e^{i\frac{p_{\perp}^2 }{ 2p^-}x^-} \partial^kF_{-k}(x^-) e^{-i \frac{p_{\perp}^2 }{ 2p^-}x^-} \Big)e^{i \frac{p_{\perp}^2 }{ 2p^-}y^-} ip_j|y_{\perp})^{en}\,.
\nonumber
\end{eqnarray}

Integrating over $y^-$ and $z^-$, we get

\begin{eqnarray}
&&\langle \int^\infty_{-\infty} dy^- e^{i x_BP^+y^-} [\infty, y^-]^{an}_yF^{{\rm q+bg};n}_{-j}(y^-, y_\perp)\rangle_{B^{\rm bg}}
\\
&&= - \frac{g^2 N_c}{ 2\pi} \int_0^{\infty} \frac{dp^-}{ p^-} (y_{\perp}|\frac{1}{ p_{\perp}^2}\int_{-\infty}^{\infty} dx^- e^{i x_B P^+ x^-} F_{-j}^n(x^-) \frac{2x_B P^+p^- }{ 2x_B P^+ p^- + p_{\perp}^2}|y_{\perp})
\nonumber\\
&&- \frac{i g^2 N_c }{ 2\pi} \int_0^{\infty} \frac{dp^- }{ p^-} (y_{\perp}|\frac{1}{ p_{\perp}^2}\int_{-\infty}^{\infty} dx^- e^{ix_B P^+ x^-} \partial^k F^n_{-k}(x^-) \frac{p_j }{ 2x_B P^+p^- + p_{\perp}^2}|y_{\perp}) \,.
\end{eqnarray}

Again using the operator relation $p_i \mathcal{O}\:=\:i\partial_i \mathcal{O}+\mathcal{O}p_i$, and restoring the Wilson lines, we rewrite the equation in the following form

\begin{eqnarray}
&&\langle \int^\infty_{-\infty} dy^- e^{i x_BP^+y^-} [\infty, y^-]^{an}_yF^{{\rm q+bg};n}_{-j}(y^-, y_\perp)\rangle_{B^{\rm bg}}
\\
&&= - \frac{g^2 N_c }{ 2\pi} \int_0^{\infty} \frac{dp^-}{ p^-}  \Big( \int^\infty_{-\infty} dy^- e^{i x_BP^+y^-} [\infty, y^-]^{an}_yF^n_{-j}(y^-, y_\perp) \Big)(y_{\perp}|\frac{2x_B P^+ p^- }{ p_{\perp}^2 (2x_B P^+ p^- + p_{\perp}^2)}|y_{\perp})
\nonumber\\
&&+ \frac{ig^2 N_c }{ 2\pi} \int_0^{\infty} \frac{dp^-}{ p^-}(y_{\perp}|\frac{p^s(2 \delta^m_s \delta^k_j - g_{js} g^{km}) }{ p_{\perp}^2} \partial_k \Big( \int^\infty_{-\infty} dy^- e^{i x_BP^+y^-} [\infty, y^-]^{an}_yF^n_{-m}(y^-, y_\perp) \Big)\frac{1}{ 2x_B P^+ p^- + p_{\perp}^2}|y_{\perp})
\nonumber
\end{eqnarray}

which is same as (<ref>). Note that in this derivation we ignore the fully transverse strength tensor.

[1]
John Collins.
Foundations of perturbative QCD, volume 32.
Cambridge University Press, 11 2013.

[2]
John C. Collins and Davison E. Soper.
Back-To-Back Jets in QCD.
Nucl. Phys. B, 193:381, 1981.
[Erratum: Nucl.Phys.B 213, 545 (1983)].

[3]
John C. Collins, Davison E. Soper, and George F. Sterman.
Transverse Momentum Distribution in Drell-Yan Pair and W and Z Boson
Nucl. Phys. B, 250:199–224, 1985.

[4]
John C. Collins and Davison E. Soper.
The Theorems of Perturbative QCD.
Ann. Rev. Nucl. Part. Sci., 37:383–409, 1987.

[5]
John C. Collins, Davison E. Soper, and George F. Sterman.
Factorization of Hard Processes in QCD.
Adv. Ser. Direct. High Energy Phys., 5:1–91, 1989.

[6]
Ruibin Meng, Fredrick I. Olness, and Davison E. Soper.
Semiinclusive deeply inelastic scattering at small q(T).
Phys. Rev. D, 54:1919–1935, 1996.

[7]
Xiang-dong Ji, Jian-ping Ma, and Feng Yuan.
QCD factorization for semi-inclusive deep-inelastic scattering at
low transverse momentum.
Phys. Rev. D, 71:034005, 2005.

[8]
Xiang-dong Ji, Jian-Ping Ma, and Feng Yuan.
QCD factorization for spin-dependent cross sections in DIS and
Drell-Yan processes at low transverse momentum.
Phys. Lett. B, 597:299–308, 2004.

[9]
Renaud Boussarie et al.
TMD Handbook.
4 2023.

[10]
Miguel G. Echevarria, Ahmad Idilbi, Zhong-Bo Kang, and Ivan Vitev.
QCD Evolution of the Sivers Asymmetry.
Phys. Rev. D, 89:074013, 2014.

[11]
Zhong-Bo Kang, Alexei Prokudin, Peng Sun, and Feng Yuan.
Extraction of Quark Transversity Distribution and Collins
Fragmentation Functions with QCD Evolution.
Phys. Rev. D, 93(1):014009, 2016.

[12]
Alessandro Bacchetta, Filippo Delcarro, Cristian Pisano, Marco Radici, and
Andrea Signori.
Extraction of partonic transverse momentum distributions from
semi-inclusive deep-inelastic scattering, Drell-Yan and Z-boson production.
JHEP, 06:081, 2017.
[Erratum: JHEP 06, 051 (2019)].

[13]
Alessandro Bacchetta, Valerio Bertone, Chiara Bissolotti, Giuseppe Bozzi,
Filippo Delcarro, Fulvio Piacenza, and Marco Radici.
Transverse-momentum-dependent parton distributions up to N$^{3}$LL
from Drell-Yan data.
JHEP, 07:117, 2020.

[14]
Ignazio Scimemi and Alexey Vladimirov.
Non-perturbative structure of semi-inclusive deep-inelastic and
Drell-Yan scattering at small transverse momentum.
JHEP, 06:137, 2020.

[15]
Valerio Bertone, Ignazio Scimemi, and Alexey Vladimirov.
Extraction of unpolarized quark transverse momentum dependent parton
distributions from Drell-Yan/Z-boson production.
JHEP, 06:028, 2019.

[16]
Alessandro Bacchetta, Filippo Delcarro, Cristian Pisano, and Marco Radici.
The 3-dimensional distribution of quarks in momentum space.
Phys. Lett. B, 827:136961, 2022.

[17]
Miguel G. Echevarria, Zhong-Bo Kang, and John Terry.
Global analysis of the Sivers functions at NLO+NNLL in QCD.
JHEP, 01:126, 2021.

[18]
Zhong-Bo Kang, Jared Reiten, Ding Yu Shao, and John Terry.
QCD evolution of the gluon Sivers function in heavy flavor dijet
production at the Electron-Ion Collider.
JHEP, 05:286, 2021.

[19]
Justin Cammarota, Leonard Gamberg, Zhong-Bo Kang, Joshua A. Miller, Daniel
Pitonyak, Alexei Prokudin, Ted C. Rogers, and Nobuo Sato.
Origin of single transverse-spin asymmetries in high-energy
Phys. Rev. D, 102(5):054002, 2020.

[20]
Marcin Bury, Alexei Prokudin, and Alexey Vladimirov.
Extraction of the Sivers function from SIDIS, Drell-Yan, and
$W^\pm/Z$ boson production data with TMD evolution.
JHEP, 05:151, 2021.

[21]
I. Balitsky and A. Tarasov.
Higher-twist corrections to gluon TMD factorization.
JHEP, 07:095, 2017.

[22]
I. Balitsky and A. Tarasov.
Power corrections to TMD factorization for Z-boson production.
JHEP, 05:150, 2018.

[23]
Markus A. Ebert, Ian Moult, Iain W. Stewart, Frank J. Tackmann, Gherardo Vita,
and Hua Xing Zhu.
Subleading power rapidity divergences and power corrections for
JHEP, 04:123, 2019.

[24]
Ian Balitsky.
Gauge-invariant TMD factorization for Drell-Yan hadronic tensor at
small x.
JHEP, 05:046, 2021.

[25]
Alexey Vladimirov.
Kinematic power corrections in TMD factorization theorem.
7 2023.

[26]
John C. Collins and Davison E. Soper.
Back-To-Back Jets: Fourier Transform from B to K-Transverse.
Nucl. Phys. B, 197:446–476, 1982.

[27]
Alexander M. Polyakov.
Gauge Fields as Rings of Glue.
Nucl. Phys. B, 164:171–188, 1980.

[28]
G. P. Korchemsky and A. V. Radyushkin.
Loop Space Formalism and Renormalization Group for the Infrared
Asymptotics of QCD.
Phys. Lett. B, 171:459–467, 1986.

[29]
S. Moch, J. A. M. Vermaseren, and A. Vogt.
The Three loop splitting functions in QCD: The Nonsinglet case.
Nucl. Phys. B, 688:101–134, 2004.

[30]
P. A. Baikov, K. G. Chetyrkin, A. V. Smirnov, V. A. Smirnov, and
M. Steinhauser.
Quark and gluon form factors to three loops.
Phys. Rev. Lett., 102:212002, 2009.

[31]
Ye Li, Andreas von Manteuffel, Robert M. Schabinger, and Hua Xing Zhu.
Soft-virtual corrections to Higgs production at N$^3$LO.
Phys. Rev. D, 91:036008, 2015.

[32]
T. Gehrmann, E. W. N. Glover, T. Huber, N. Ikizlerli, and C. Studerus.
Calculation of the quark and gluon form factors to three loops in
JHEP, 06:094, 2010.

[33]
R. N. Lee, A. V. Smirnov, and V. A. Smirnov.
Analytic Results for Massless Three-Loop Form Factors.
JHEP, 04:020, 2010.

[34]
Ye Li, Duff Neill, and Hua Xing Zhu.
An exponential regulator for rapidity divergences.
Nucl. Phys. B, 960:115193, 2020.

[35]
Ye Li and Hua Xing Zhu.
Bootstrapping Rapidity Anomalous Dimensions for Transverse-Momentum
Phys. Rev. Lett., 118(2):022004, 2017.

[36]
Alexey A. Vladimirov.
Correspondence between Soft and Rapidity Anomalous Dimensions.
Phys. Rev. Lett., 118(6):062001, 2017.

[37]
Ignazio Scimemi and Alexey Vladimirov.
Systematic analysis of double-scale evolution.
JHEP, 08:003, 2018.

[38]
John C. Collins and Davison E. Soper.
Parton Distribution and Decay Functions.
Nucl. Phys. B, 194:445–492, 1982.

[39]
Zhong-Bo Kang and Jian-Wei Qiu.
QCD evolution of naive-time-reversal-odd parton distribution
Phys. Lett. B, 713:273–276, 2012.

[40]
Peng Sun and Feng Yuan.
Transverse momentum dependent evolution: Matching semi-inclusive
deep inelastic scattering processes to Drell-Yan and W/Z boson production.
Phys. Rev. D, 88(11):114012, 2013.

[41]
Ling-Yun Dai, Zhong-Bo Kang, Alexei Prokudin, and Ivan Vitev.
Next-to-leading order transverse momentum-weighted Sivers asymmetry
in semi-inclusive deep inelastic scattering: the role of the three-gluon
Phys. Rev. D, 92(11):114024, 2015.

[42]
V. M. Braun, A. N. Manashov, and B. Pirnay.
Scale dependence of twist-three contributions to single spin
Phys. Rev. D, 80:114002, 2009.
[Erratum: Phys.Rev.D 86, 119902 (2012)].

[43]
Miguel G. Echevarria, Ignazio Scimemi, and Alexey Vladimirov.
Universal transverse momentum dependent soft function at NNLO.
Phys. Rev. D, 93(5):054004, 2016.

[44]
Ignazio Scimemi, Andrey Tarasov, and Alexey Vladimirov.
Collinear matching for Sivers function at next-to-leading order.
JHEP, 05:125, 2019.

[45]
S. Catani and M. Grazzini.
Higgs Boson Production at Hadron Colliders: Hard-Collinear
Coefficients at the NNLO.
Eur. Phys. J. C, 72:2013, 2012.
[Erratum: Eur.Phys.J.C 72, 2132 (2012)].

[46]
Stefano Catani, Leandro Cieri, Daniel de Florian, Giancarlo Ferrera, and
Massimiliano Grazzini.
Vector boson production at hadron colliders: hard-collinear
coefficients at the NNLO.
Eur. Phys. J. C, 72:2195, 2012.

[47]
Thomas Gehrmann, Thomas Luebbert, and Li Lin Yang.
Calculation of the transverse parton distribution functions at
next-to-next-to-leading order.
JHEP, 06:155, 2014.

[48]
Thomas Lübbert, Joel Oredsson, and Maximilian Stahlhofen.
Rapidity renormalized TMD soft and beam functions at two loops.
JHEP, 03:168, 2016.

[49]
Miguel G. Echevarria, Ignazio Scimemi, and Alexey Vladimirov.
Unpolarized Transverse Momentum Dependent Parton Distribution and
Fragmentation Functions at next-to-next-to-leading order.
JHEP, 09:004, 2016.

[50]
Yuri L. Dokshitzer.
Calculation of the Structure Functions for Deep Inelastic Scattering
and e+ e- Annihilation by Perturbation Theory in Quantum Chromodynamics.
Sov. Phys. JETP, 46:641–653, 1977.

[51]
V. N. Gribov and L. N. Lipatov.
Deep inelastic e p scattering in perturbation theory.
Sov. J. Nucl. Phys., 15:438–450, 1972.

[52]
Guido Altarelli and G. Parisi.
Asymptotic Freedom in Parton Language.
Nucl. Phys. B, 126:298–318, 1977.

[53]
Ignazio Scimemi and Alexey Vladimirov.
Analysis of vector boson production within TMD factorization.
Eur. Phys. J. C, 78(2):89, 2018.

[54]
F. Landry, R. Brock, Pavel M. Nadolsky, and C. P. Yuan.
Tevatron Run-1 $Z$ boson data and Collins-Soper-Sterman resummation
Phys. Rev. D, 67:073016, 2003.

[55]
Anton V. Konychev and Pavel M. Nadolsky.
Universality of the Collins-Soper-Sterman nonperturbative function
in gauge boson production.
Phys. Lett. B, 633:710–714, 2006.

[56]
Thomas Becher, Matthias Neubert, and Daniel Wilhelm.
Electroweak Gauge-Boson Production at Small $q_T$: Infrared Safety
from the Collinear Anomaly.
JHEP, 02:124, 2012.

[57]
Umberto D'Alesio, Miguel G. Echevarria, Stefano Melis, and Ignazio Scimemi.
Non-perturbative QCD effects in $q_{T}$ spectra of Drell-Yan and
Z-boson production.
JHEP, 11:098, 2014.

[58]
Fabio Dominguez, Bo-Wen Xiao, and Feng Yuan.
$k_t$-factorization for Hard Processes in Nuclei.
Phys. Rev. Lett., 106:022301, 2011.

[59]
Fabio Dominguez, Cyrille Marquet, Bo-Wen Xiao, and Feng Yuan.
Universality of Unintegrated Gluon Distributions at small x.
Phys. Rev. D, 83:105005, 2011.

[60]
Bo-Wen Xiao, Feng Yuan, and Jian Zhou.
Transverse Momentum Dependent Parton Distributions at Small-x.
Nucl. Phys. B, 921:104–126, 2017.

[61]
Tolga Altinoluk, Néstor Armesto, Guillaume Beuf, Mauricio Martínez, and
Carlos A. Salgado.
Next-to-eikonal corrections in the CGC: gluon production and spin
asymmetries in pA collisions.
JHEP, 07:068, 2014.

[62]
Tolga Altinoluk, Néstor Armesto, Guillaume Beuf, and Alexis Moscoso.
Next-to-next-to-eikonal corrections in the CGC.
JHEP, 01:114, 2016.

[63]
Yuri V. Kovchegov, Daniel Pitonyak, and Matthew D. Sievert.
Helicity Evolution at Small-x.
JHEP, 01:072, 2016.
[Erratum: JHEP 10, 148 (2016)].

[64]
Tolga Altinoluk and Adrian Dumitru.
Particle production in high-energy collisions beyond the shockwave
Phys. Rev. D, 94(7):074032, 2016.

[65]
I. Balitsky and A. Tarasov.
Rapidity evolution of gluon TMD from low to moderate x.
JHEP, 10:017, 2015.

[66]
Yuri V. Kovchegov, Daniel Pitonyak, and Matthew D. Sievert.
Helicity Evolution at Small $x$: Flavor Singlet and Non-Singlet
Phys. Rev. D, 95(1):014033, 2017.

[67]
I. Balitsky and A. Tarasov.
Gluon TMD in particle production from low to moderate x.
JHEP, 06:164, 2016.

[68]
Pedro Agostini, Tolga Altinoluk, and Néstor Armesto.
Non-eikonal corrections to multi-particle production in the Color
Glass Condensate.
Eur. Phys. J. C, 79(7):600, 2019.

[69]
Pedro Agostini, Tolga Altinoluk, and Néstor Armesto.
Effect of non-eikonal corrections on azimuthal asymmetries in the
Color Glass Condensate.
Eur. Phys. J. C, 79(9):790, 2019.

[70]
Florian Cougoulic and Yuri V. Kovchegov.
Helicity-dependent generalization of the JIMWLK evolution.
Phys. Rev. D, 100(11):114020, 2019.

[71]
Tolga Altinoluk, Guillaume Beuf, Alina Czajka, and Arantxa Tymowska.
Quarks at next-to-eikonal accuracy in the CGC: Forward quark-nucleus
Phys. Rev. D, 104(1):014019, 2021.

[72]
Tolga Altinoluk and Guillaume Beuf.
Quark and scalar propagators at next-to-eikonal accuracy in the CGC
through a dynamical background gluon field.
Phys. Rev. D, 105(7):074026, 2022.

[73]
Yuri V. Kovchegov and M. Gabriel Santiago.
Quark sivers function at small $x$: spin-dependent odderon and the
sub-eikonal evolution.
JHEP, 11:200, 2021.
[Erratum: JHEP 09, 186 (2022)].

[74]
Giovanni Antonio Chirilli.
High-energy operator product expansion at sub-eikonal level.
JHEP, 06:096, 2021.

[75]
Florian Cougoulic, Yuri V. Kovchegov, Andrey Tarasov, and Yossathorn Tawabutr.
Quark and gluon helicity evolution at small x: revised and updated.
JHEP, 07:095, 2022.

[76]
Pedro Agostini, Tolga Altinoluk, and Néstor Armesto.
Finite width effects on the azimuthal asymmetry in proton-nucleus
collisions in the Color Glass Condensate.
Phys. Lett. B, 840:137892, 2023.

[77]
Victor S. Fadin, E. A. Kuraev, and L. N. Lipatov.
On the Pomeranchuk Singularity in Asymptotically Free Theories.
Phys. Lett. B, 60:50–52, 1975.

[78]
E. A. Kuraev, L. N. Lipatov, and Victor S. Fadin.
Multi - Reggeon Processes in the Yang-Mills Theory.
Sov. Phys. JETP, 44:443–450, 1976.

[79]
E. A. Kuraev, L. N. Lipatov, and Victor S. Fadin.
The Pomeranchuk Singularity in Nonabelian Gauge Theories.
Sov. Phys. JETP, 45:199–204, 1977.

[80]
I. I. Balitsky and L. N. Lipatov.
The Pomeranchuk Singularity in Quantum Chromodynamics.
Sov. J. Nucl. Phys., 28:822–829, 1978.

[81]
I. Balitsky.
Operator expansion for high-energy scattering.
Nucl. Phys. B, 463:99–160, 1996.

[82]
Ian Balitsky.
Operator expansion for diffractive high-energy scattering.
AIP Conf. Proc., 407(1):953, 1997.

[83]
Yuri V. Kovchegov.
Small x F(2) structure function of a nucleus including multiple
pomeron exchanges.
Phys. Rev. D, 60:034008, 1999.

[84]
L. F. Abbott.
The Background Field Method Beyond One Loop.
Nucl. Phys. B, 185:189–203, 1981.

[85]
L. F. Abbott.
Introduction to the Background Field Method.
Acta Phys. Polon. B, 13:33, 1982.

[86]
Jui-yu Chiu, Ambar Jain, Duff Neill, and Ira Z. Rothstein.
The Rapidity Renormalization Group.
Phys. Rev. Lett., 108:151601, 2012.

[87]
Jui-Yu Chiu, Ambar Jain, Duff Neill, and Ira Z. Rothstein.
A Formalism for the Systematic Treatment of Rapidity Logarithms in
Quantum Field Theory.
JHEP, 05:084, 2012.

[88]
Sean Fleming.
The role of Glauber exchange in soft collinear effective theory and
the Balitsky–Fadin–Kuraev–Lipatov
Phys. Lett. B, 735:266–271, 2014.

[89]
Ira Z. Rothstein and Iain W. Stewart.
An Effective Field Theory for Forward Scattering and Factorization
JHEP, 08:025, 2016.

[90]
Zhong-Bo Kang and Xiaohui Liu.
Power Counting the Small-$x$ Observables.
10 2019.

[91]
Hao-yu Liu, Kexin Xie, Zhongbo Kang, and Xiaohui Liu.
Single inclusive jet production in pA collisions at NLO in the
small-x regime.
JHEP, 07:041, 2022.

[92]
Miguel G. Echevarria, Ahmad Idilbi, and Ignazio Scimemi.
Factorization Theorem For Drell-Yan At Low $q_T$ And Transverse
Momentum Distributions On-The-Light-Cone.
JHEP, 07:002, 2012.

[93]
Jui-yu Chiu, Andreas Fuhrer, Andre H. Hoang, Randall Kelley, and Aneesh V.
Soft-Collinear Factorization and Zero-Bin Subtractions.
Phys. Rev. D, 79:053007, 2009.

[94]
M. Beneke and T. Feldmann.
Factorization of heavy to light form-factors in soft collinear
effective theory.
Nucl. Phys. B, 685:249–296, 2004.

[95]
Jui-yu Chiu, Frank Golf, Randall Kelley, and Aneesh V. Manohar.
Electroweak Sudakov corrections using effective field theory.
Phys. Rev. Lett., 100:021802, 2008.

[96]
Thomas Becher and Guido Bell.
Analytic Regularization in Soft-Collinear Effective Theory.
Phys. Lett. B, 713:41–46, 2012.

[97]
Thomas Becher and Matthias Neubert.
Drell-Yan Production at Small $q_T$, Transverse Parton Distributions
and the Collinear Anomaly.
Eur. Phys. J. C, 71:1665, 2011.

[98]
Thomas Becher, Matthias Neubert, and Daniel Wilhelm.
Higgs-Boson Production at Small Transverse Momentum.
JHEP, 05:110, 2013.

[99]
George Leibbrandt.
The Light Cone Gauge in Yang-Mills Theory.
Phys. Rev. D, 29:1699, 1984.

[100]
Tolga Altinoluk, Guillaume Beuf, and Jamal Jalilian-Marian.
Renormalization of the gluon distribution function in the background
field formalism.
5 2023.

[101]
A. H. Mueller, Bo-Wen Xiao, and Feng Yuan.
Sudakov Resummation in Small-$x$ Saturation Formalism.
Phys. Rev. Lett., 110(8):082301, 2013.

[102]
A. H. Mueller, Bo-Wen Xiao, and Feng Yuan.
Sudakov double logarithms resummation in hard processes in the
small-x saturation formalism.
Phys. Rev. D, 88(11):114010, 2013.

[103]
Paul Caucal, Farid Salazar, and Raju Venugopalan.
Dijet impact factor in DIS at next-to-leading order in the Color
Glass Condensate.
JHEP, 11:222, 2021.

[104]
Paul Caucal, Farid Salazar, Björn Schenke, and Raju Venugopalan.
Back-to-back inclusive dijets in DIS at small x: Sudakov suppression
and gluon saturation at NLO.
JHEP, 11:169, 2022.

[105]
Paul Caucal, Farid Salazar, Björn Schenke, Tomasz Stebel, and Raju
Back-to-back inclusive dijets in DIS at small x: gluon
Weizsäcker-Williams distribution at NLO.
JHEP, 08:062, 2023.

[106]
Paul Caucal, Farid Salazar, Björn Schenke, Tomasz Stebel, and Raju
Back-to-back inclusive dijets in DIS at small $x$: Complete NLO
results and predictions.
7 2023.

[107]
Anna Stasto, Shu-Yi Wei, Bo-Wen Xiao, and Feng Yuan.
On the Dihadron Angular Correlations in Forward $pA$ collisions.
Phys. Lett. B, 784:301–306, 2018.

[108]
Ian Balitsky and Giovanni A. Chirilli.
Rapidity evolution of TMDs with running coupling.
Phys. Rev. D, 106(3):034007, 2022.

[109]
Ian Balitsky.
Rapidity-only TMD factorization at one loop.
JHEP, 03:029, 2023.

[110]
Jian Zhou.
The evolution of the small x gluon TMD.
JHEP, 06:151, 2016.

[111]
Martin Hentschinski, Krzysztof Kutak, and Andreas van Hameren.
Forward Higgs production within high energy factorization in the
heavy quark limit at next-to-leading order accuracy.
Eur. Phys. J. C, 81(2):112, 2021.
[Erratum: Eur.Phys.J.C 81, 262 (2021)].

[112]
Martin Hentschinski.
Transverse momentum dependent gluon distribution within high energy
factorization at next-to-leading order.
Phys. Rev. D, 104(5):054014, 2021.

[113]
Duff Neill, Aditya Pathak, and Iain W. Stewart.
Small-x factorization from effective field theory.
JHEP, 09:089, 2023.

[114]
Iain Stewart and Varun Vaidya.
Power Counting to Saturation.
5 2023.

[115]
S. Catani, M. Ciafaloni, and F. Hautmann.
High-energy factorization and small x heavy flavor production.
Nucl. Phys. B, 366:135–188, 1991.

[116]
S. Catani, M. Ciafaloni, and F. Hautmann.
GLUON CONTRIBUTIONS TO SMALL x HEAVY FLAVOR PRODUCTION.
Phys. Lett. B, 242:97–102, 1990.

[117]
S. Catani and F. Hautmann.
High-energy factorization and small x deep inelastic scattering
beyond leading order.
Nucl. Phys. B, 427:475–524, 1994.

[118]
G. P. Salam.
A Resummation of large subleading corrections at small x.
JHEP, 07:019, 1998.

[119]
Guido Altarelli, Richard D. Ball, and Stefano Forte.
Resummation of singlet parton evolution at small x.
Nucl. Phys. B, 575:313–329, 2000.

[120]
M. Ciafaloni, D. Colferai, and G. P. Salam.
A collinear model for small x physics.
JHEP, 10:017, 1999.

[121]
M. Ciafaloni, D. Colferai, and G. P. Salam.
Renormalization group improved small x equation.
Phys. Rev. D, 60:114036, 1999.

[122]
Guido Altarelli, Richard D. Ball, and Stefano Forte.
Factorization and resummation of small x scaling violations with
running coupling.
Nucl. Phys. B, 621:359–387, 2002.

[123]
Guido Altarelli, Richard D. Ball, and Stefano Forte.
An Anomalous dimension for small x evolution.
Nucl. Phys. B, 674:459–483, 2003.

[124]
M. Ciafaloni, D. Colferai, G. P. Salam, and A. M. Stasto.
Renormalization group improved small x Green's function.
Phys. Rev. D, 68:114003, 2003.

[125]
M. Ciafaloni, D. Colferai, D. Colferai, G. P. Salam, and A. M. Stasto.
Extending QCD perturbation theory to higher energies.
Phys. Lett. B, 576:143–151, 2003.

[126]
Marcello Ciafaloni, Dimitri Colferai, Gavin P. Salam, and Anna M. Stasto.
The Gluon splitting function at moderately small x.
Phys. Lett. B, 587:87–94, 2004.

[127]
Simone Marzani.
Combining $Q_T$ and small-$x$ resummations.
Phys. Rev. D, 93(5):054047, 2016.

[128]
Dimitri Colferai, Wanchen Li, and Anna M. Stasto.
Renormalization group improved photon impact factors and the high
energy virtual photon scattering.
11 2023.

[129]
Xiangdong Ji, Peng Sun, Xiaonu Xiong, and Feng Yuan.
Soft factor subtraction and transverse momentum dependent parton
distributions on the lattice.
Phys. Rev. D, 91:074009, 2015.

[130]
Xiangdong Ji, Lu-Chang Jin, Feng Yuan, Jian-Hui Zhang, and Yong Zhao.
Transverse momentum dependent parton quasidistributions.
Phys. Rev. D, 99(11):114006, 2019.

[131]
Markus A. Ebert, Iain W. Stewart, and Yong Zhao.
Determining the Nonperturbative Collins-Soper Kernel From Lattice
Phys. Rev. D, 99(3):034505, 2019.

[132]
Markus A. Ebert, Iain W. Stewart, and Yong Zhao.
Towards Quasi-Transverse Momentum Dependent PDFs Computable on the
JHEP, 09:037, 2019.

[133]
Markus A. Ebert, Iain W. Stewart, and Yong Zhao.
Renormalization and Matching for the Collins-Soper Kernel from
Lattice QCD.
JHEP, 03:099, 2020.

[134]
Xiangdong Ji, Yizhuang Liu, and Yu-Sheng Liu.
TMD soft function from large-momentum effective theory.
Nucl. Phys. B, 955:115054, 2020.

[135]
Xiangdong Ji, Yizhuang Liu, and Yu-Sheng Liu.
Transverse-momentum-dependent parton distribution functions from
large-momentum effective theory.
Phys. Lett. B, 811:135946, 2020.

[136]
Markus A. Ebert, Stella T. Schindler, Iain W. Stewart, and Yong Zhao.
One-loop Matching for Spin-Dependent Quasi-TMDs.
JHEP, 09:099, 2020.

[137]
Xiangdong Ji, Yizhuang Liu, Andreas Schäfer, and Feng Yuan.
Single Transverse-Spin Asymmetry and Sivers Function in Large
Momentum Effective Theory.
Phys. Rev. D, 103(7):074005, 2021.

[138]
Xiangdong Ji and Yizhuang Liu.
Computing light-front wave functions without light-front
quantization: A large-momentum effective theory approach.
Phys. Rev. D, 105(7):076014, 2022.

[139]
Markus A. Ebert, Stella T. Schindler, Iain W. Stewart, and Yong Zhao.
Factorization connecting continuum & lattice TMDs.
JHEP, 04:178, 2022.

[140]
Daniel Adamiak, Nicholas Baldonado, Yuri V. Kovchegov, W. Melnitchouk, Daniel
Pitonyak, Nobuo Sato, Matthew D. Sievert, Andrey Tarasov, and Yossathorn
Global analysis of polarized DIS & SIDIS data with improved
small-$x$ helicity evolution.
8 2023.

[141]
Julian S. Schwinger.
On gauge invariance and vacuum polarization.
Phys. Rev., 82:664–679, 1951.

§.§ Background field method

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,[Here we understant $C$ as a general notation for the quark and gluon fields.] one can write it as a functional integral over those fields

\begin{eqnarray}
&&\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)}
\label{int:mel-form1}
\end{eqnarray}

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.

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$,

\begin{eqnarray}
&&C_\mu \to A_\mu + B_\mu
\label{int:Bj-spt}
\end{eqnarray}

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\gg \mu^2$, and fields $B$ are characterized with $k_\perp \ll \mu^2$. This is called the collinear factorization scheme which is defined by a strong ordering of emission in transverse momenta.

Another example when the factorized QCD medium can described by spliting into two modes is the high energy rapidity factorization in DIS. In this case the field mode $A$ has longitudinal momentum fraction $k^- \gg \sigma$ and the mode $B$ has $k^- \ll \sigma$. The factorization scale $\sigma$ in this case in the rigid cut-off in $k^-$, see Fig. <ref>a.

The high -energy rapidity factorization for DIS with a rigid cut-off $\sigma$ in $k^-$. a) The "quantum" fields $A$ are above the $\sigma$ cut-off. The background fields $B$ (grey) are below the cut-off. b) To study dependence on the factorization scale we split the background fields into two components $B \to B^{\rm q} + B^{\rm bg}$. We integrate over $B^{\rm q}$ (blue) fields in a fixed background of $B^{\rm bg}$ (red) fields.

For brevity, let's assume that there the factorization can be described with only two modes as in (<ref>). 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'll see in the next subsection, to describe the TMD factorization one has to introduce four modes: one hard modes ("quantum"), two collinear modes (background) and a soft modes for an intersection between the collinear modes.

We will call the $A$ fields as "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 contain 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 (<ref>) we can rewrite the functional integral in Eq. (<ref>) as

\begin{eqnarray}
&&\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)}
\label{int:bBqA}
\end{eqnarray}

where the integral over "qunatum" mode $A$ is

\begin{eqnarray}
&&\tilde{\mathcal{O}}(B, \sigma) = \int \mathcal{D}A ~\mathcal{O}(A + B) e^{iS_{bQCD}(A, B)}
\label{int:AbInt}
\end{eqnarray}

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

\begin{eqnarray}
&&S_{bQCD}(A, B) = S_{\rm QCD}(A+B) - S_{\rm QCD}(B)
\end{eqnarray}

In Eq. (<ref>) 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. (<ref>), i.e. fields of the projectile, can be evaluated perturbatively. Though 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

\begin{eqnarray}
&&\tilde{\mathcal{O}}(B, \sigma) = \sum_i H_i(\sigma)\otimes\mathcal{V}_i(B, \sigma)
\label{int:aAi}
\end{eqnarray}

where $\mathcal{V}_i$ are operators constructed from background mode $B$. Since the background fields $B$ depend on the factorization scale $\sigma$, the operators also depend on it, which we indicate explicitly. $H_i(\sigma)$ is a hard function which we obtain after integration over "quantum" fields $A$.

After this, Eq. (<ref>) reads

\begin{eqnarray}
&&\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
\label{int:fact}
\end{eqnarray}

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

\begin{eqnarray}
&&\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)) e^{iS_{\rm QCD}(B)}
\label{int:meB}
\end{eqnarray}

where 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 well known collinear PDFs.

Eq. (<ref>) is a factorization formula which we obtain in assumption that the QCD medium develops different dynamical modes, so the separation (<ref>) 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. This 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. (<ref>) and (<ref>) for the TMD factorization.

To extract the dependence on the factorization scales in the background field method one can use the following approach. Starting with a matrix element (<ref>), we introduce a new set of factorization scales $\sigma'$ 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'$, and $B^{\rm bg}$ is a new background field corresponding to fields below the scale $\sigma'$. For the high-energy rapidity factorization scheme this separation is shown in Fig. <ref>b.

Since $\sigma$ is an upper cut-off for field modes, we will often refer to it as a UV scale (or a set of scales). Similarly, since the $\sigma'$ scale is a lower scale for the $B^{\rm q}$ modes, we will refer to it as an IR scale.

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 following answer

\begin{eqnarray}
&&\langle P_1|\mathcal{V}_i( \sigma)|P_2\rangle = \sum_j C_{ij}(\sigma, \sigma')\otimes \langle P_1|\mathcal{V}_j(\sigma')|P_2\rangle
\label{mel:ptb}
\end{eqnarray}

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

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

In this paper we aim to calculate the $C_{ij}(\sigma, \sigma')$ 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'$ coresponds to the cut-offs $(\mu^2_{\rm IR}, \rho)$, see Fig. <ref>. More details can be found in the Introduction and the next section.

Another important, though rather technical aspect, is the way of how different field modes can be introduced in practice. To do that one can for example use rigid cut-offs as it is done in the high-energy rapidity factorization approach. This cut-offs subsequently appear in the perturbative calculation of functional integrals, e.g. in Eq. (<ref>), as cut-offs in loop integrals over momenta. These cut-offs regulate all divergences of the loop integrals. However, due to technical difficulties this rigid cut-off are rarely used.

An alternative scheme, which we use in this paper, is the renormalization approach. In this approach each divergent integral is regulated with a regulator which is accompanied with a renormalization scale, which plays a role of a cut-off separating different modes. As we mentioned before we will use the 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 faraction integrals as in Eq. (<ref>).[There is a number of alternative methods: deviation from the light-cone direction [...], $\delta$-regulator [...], exponential regulator [...], and analytic regulator [...].] Note that in the renormalization approach we will find some unphysical poles which is a subsequence of the separation of different dynamical field modes. This poles as usual should be removed by a proper renormalization of operators.