Matching generalised transverse-momentum-dependent distributions onto generalised parton distributions at one loop

Abstract

The operator definition of generalised transverse momentum-dependent (GTMD) distributions is exploited to compute for the first time the full set of one-loop corrections to the off-forward matching functions. These functions allow one to obtain GTMDs in the perturbative regime in terms of generalised parton distributions (GPDs). In the unpolarised case, non-perturbative corrections can be incorporated using recent determinations of transverse-momentum-dependent (TMD) distributions. Evolution effects for GTMDs closely follow those for TMDs and can thus be easily accounted for up to next-to-next-to-leading logarithmic accuracy. As a by-product, the relevant one-loop anomalous dimensions are derived, confirming previous results. As a practical application, numerical results for a specific kind of GTMD are presented, highlighting some salient features.

1 Introduction

At present, the study of the hadronic structure is a particularly lively field of research. On the one hand, this is a very interesting topic in itself, and on the other hand, it is instrumental to precision physics at present and future high-energy colliders.

The case of the unpolarised collinear parton-distribution functions (PDFs) of the proton is emblematic of the huge effort that is being put into the study of the hadronic structure. Driven by the experimental activity of colliders such as HERA, Tevatron, and the Large Hadron Collider (LHC), and by a steady methodological and theoretical progress, PDFs are nowadays known with an astonishing precision [1,2,3,4,5]. Driven by the need for accuracy, the study of unpolarised collinear fragmentation functions (FFs) has also recently been intensified, leading to accurate determinations [6,7,8,9,10,11].

Although PDFs and FFs have a broad phenomenological applicability, they encode partial information on the hadronic structure corresponding to the longitudinal-momentum distribution of partons inside hadrons. Transverse-momentum-dependent (TMD) distributions are instead also sensitive to the transverse momentum of partons, thus extending the information provided by PDFs and FFs [12, 13]. Also, due to their relevance in hot topics such as the precision determination of the mass of the W boson, TMDs are currently receiving particular attention, and much work is being invested in their determination [14,15,16,17,18,19].

Another typology of distributions relevant to the study of the hadronic structure is that of generalised parton distributions (GPDs). These distributions give us access to the energy-momentum tensor of hadrons [20, 21], providing us with a handle on important quantities like the transverse position of partons [22, 23] and their angular momentum [24]. Phenomenological determinations of GPDs do exist [25,26,27], but as of today, they are much less developed than modern analyses of PDFs/FFs and TMDs. However, the recent approval of the electron-ion colliders in China (EicC) [28] and in the USA (EIC) [29], has revived the interest in GPDs that is now a rapidly growing field. In addition, new technologies on the lattice [30, 31] have made first-principle computations of GPDs more accessible, giving additional momentum to their study.

It turns out that PDFs, TMDs and GPDs are all projections of more general quantities, usually referred to as generalised TMDs (GTMDs), that can be considered as their mother distributions [32,33,34,35,36,37,38]. As of today, very little is quantitatively known about GTMDs, and most of the existing studies are based on models [39,40,41,42,43,44,45,46,47,48,49,50,51,52]. In spite of recent proposals to access GTMDs experimentally [53,54,55], data sensitive to these distributions is presently scarce, making phenomenological studies laborious.

The goal of this paper is to exploit as much as possible our knowledge of the projections of GTMDs, namely GPDs and TMDs, to reconstruct GTMDs themselves to the best accuracy possible. For this purpose, a proper definition of GTMDs has to be devised, which enables the computation of relevant perturbative quantities. In the spirit of TMD factorisation, this was done in Ref. [56], where generalising the TMD operators to the off-forward case, a rapidity-divergence-free definition of GTMD correlators was given. We will start from this definition to compute for the first time the full set of the so-called matching functions at one-loop accuracy. These functions allow one to obtain GTMDs in terms of GPDs by accounting for the emission of partons with large transverse momentum ${\textbf{k}}_T$, whose effect is thus computable in perturbation theory. We will finally use these matching functions to reconstruct realistic GTMDs also accounting for evolution and non-perturbative effects.

The outline of the paper is as follows. In Sect. 2.1, we will give a precise operator definition of the GTMD correlators of our interest. Section 2.2 is devoted to the renormalisation of the ultraviolet (UV) divergences of these correlators and to the derivation of their evolution equations. In Sect. 2.3, we will state the GTMD matching formula that factorises large-${\textbf{k}}_T$ emissions into a set of matching functions to be convoluted with GPDs. In this section, we will exploit this matching formula to express the one-loop matching functions in terms of perturbative quantities, namely, the soft function and the unsubtracted parton-in-parton GTMD correlator, whose one-loop corrections are computed in Sects. 2.4 and 2.5, respectively. At this point, we will be in a position to obtain the explicit expression for the one-loop matching functions. In Sect. 2.6, we will show that, as expected, these expressions tend to their TMD counterpart in the forward limit. As a by-product of UV renormalisation, in Sect. 2.7, we will extract the anomalous dimensions that govern the GTMD evolution at one loop, finding that they agree with the TMD ones. A numerical implementation of the matching functions, combined with other existing perturbative and non-perturbative ingredients, puts us in a position to reconstruct realistic quark and gluon GTMDs to high accuracy. A study of the resulting distributions is presented in Sect. 3 where we will consider the behaviour of GTMDs under different points of view, commenting on some peculiar features. Finally, in Sect. 4, we will draw our conclusions.

2 Theoretical setup and results

As argued in the pioneering Ref. [56], a sound definition of GTMDs, that is, a definition free of rapidity divergences and that can thus be used for phenomenology, requires taking into proper account the soft function. In this section, we will review the reasoning of Ref. [56] that leads to a rapidity-divergence-free definition of the GTMD correlator. Working in the space of the Fourier-conjugate variable of the partonic momentum ${\textbf{k}}_T$, denoted by ${\textbf{b}}_T$, we will formulate the explicit definition of the GTMD correlator both for quarks and gluons. We will then state the matching formula that, for ${\textbf{b}}_T\simeq {\textbf{0}}$, allows us to express the GTMD correlator in terms of (collinear) GPD correlators by means of perturbatively computable matching functions. Appealing to the concept of parton-in-parton distributions (see, e.g., Refs. [13, 57]), we will finally obtain the one-loop corrections to these matching functions.

This programme requires the computation of the first perturbative correction to the soft function and the so-called unsubtracted GTMD correlators. By combining these two quantities, we will explicitly exhibit the cancellation of the rapidity divergences and obtain, for the first time, the full set of off-forward matching functions at one loop. In addition, we will show that their forward limit coincides with the one-loop TMD matching functions, as it should. Finally, as a by-product of the renormalisation of the UV divergences, we will also derive the GTMD evolution equations and extract the leading-order term of the relevant anomalous dimensions, confirming previous results.

The seminal work of Refs. [35, 38] has provided us with a thorough classification of the structures that emerge from the analysis of the GTMD correlators for spin-1/2 targets. However, if taken literally, the GTMD definitions given in those papers produce divergent results upon inclusion of radiative corrections, even after the customary renormalisation of the UV divergences. Exactly like in the TMD case [13, 58], these spurious divergences are of infrared (IR) origin and stem from the region of loop momenta k where the rapidity of the emitted partons, $y=\ln (k^+/k^-)$, becomes large [59]: hence, they are usually called rapidity divergences. This “inconvenience” can be fixed in two main steps:

• 1. removing the overlap of the GTMD correlator with the soft modes by means of the so-called zero-bin (or soft) subtraction; and

• 2. redistributing the soft function, that typically appears in a factorisation theorem along with two beam functions, between the two GTMD correlators: this usually amounts to assigning the square root of the soft function to each of them.

An implementation of these steps is guaranteed to lead to a cancellation of the rapidity divergences. The reader is referred to, for example, Refs. [60, 61] and references therein for a more detailed discussion in the TMD framework.

In order to perform an explicit perturbative calculation, rapidity divergences need to be regulated on a diagram-by-diagram basis. Several rapidity-divergence regulators have been proposed so far in the literature. However, for a specific class of regulators [61, 62], steps (1) and (2) combine in a way that the net effect is that the unsubtracted GTMD correlator has to be divided by the square root of the soft function. In this paper, we will use the principal-value regulator, originally introduced in Ref. [63], that indeed belongs to this class of regulators.

2.1 Operator definition of GTMDs

This brief introduction allows us to finally give the exact definition for the quark and gluon GTMD correlators that will be used for the computation of the one-loop matching functions. As anticipated above, it is convenient to work in ${\textbf{b}}_T$ space which simply amounts to taking the Fourier transform of the ${\textbf{k}}_T$-space GTMD correlators. In addition, for definiteness, we will limit the discussion to the specific twist-2 GTMD correlators whose forward limit gives the unpolarised quark and gluon TMDs. The corresponding operator definition for a generic hadronic target H is then

$$\begin{aligned} \hat{{\mathcal {F}}}_{i/H}(x,\xi ,{\textbf{b}}_T,t)={\hat{S}}_i^{-\frac{1}{2}}({\textbf{b}}_T) {\hat{\Phi }}_{i/H}(x,\xi ,{\textbf{b}}_T,t),\quad i=q,g,\nonumber \\ \end{aligned}$$

(1)

where ${\hat{S}}_i$ is the appropriate soft function, and ${\hat{\Phi }}_{i/H}$ is the unsubtracted GTMD correlator. The quark and gluon correlators, respectively, read as

(2)

Fig. 1

Graphical representation of the quark (left) and gluon (right) unsubtracted GTMD correlators defined in Eq. (2)

where we have used the following shorthand definitions: $\eta = yn+{\textbf{b}}_T,$¹$P_{\mathrm{in/out}}=P\pm \Delta /2$, $t=\Delta ^2$, and $\xi =2n\cdot \Delta / n\cdot P$. Here, $\psi _q$ indicates the spinor of the quark flavour q and $F_a^{\mu \nu }$ the gluon field strength, while $W_{n,i}$ is the Wilson line in the n direction. The integrals run between $-\infty$ and $+\infty$, and a summation over repeated Lorentz indices is understood. Also, the index j in ${\hat{\Phi }}_{g/H}$ is summed over and it runs over the (physical) transverse components, that is, $j=1,2$. n and ${\overline{n}}$ (the latter appearing below) are light-like four-vectors, $n^2={\overline{n}}^2=0$, and their scalar product evaluates to $n\cdot {\overline{n}}=1$.² Also notice that the parton index i, labelling both the soft function ${\hat{S}}_i$ and the Wilson line $W_{n,i}$, denotes the colour-group representation. Specifically, for $i=q(g)$, the fundamental (adjoint) representation of the SU(3) generators is to be used. Finally, the symbol $\hat{\,}$ over ${\hat{S}}_i$ and ${\hat{\Phi }}_{i/H}$ indicates that these quantities are bare, that is, they are UV divergent in four dimensions.

The definition of the GTMD correlators is still incomplete because we have not specified ${\hat{S}}_i$ and $W_i$. In fact, the precise form of these quantities depends on the choice of the gauge. Two specific gauges are usually considered in this context: the Lorentz gauge ($\partial \cdot A=0$) and the light-cone gauge ($n\cdot A=0$). When considering ${\textbf{k}}_T$-integrated quantities, which corresponds to setting ${\textbf{b}}_T={\textbf{0}}$, the light-cone gauge is often preferred. The reason is that the Wilson lines reduce to unity with a consequent reduction of the number of diagrams to be considered (see, e.g., Refs. [57, 63]). This simplification, however, comes at the price of a more complicated gluon propagator. Conversely, when ${\textbf{k}}_T$ is left unintegrated, so that ${\textbf{b}}_T\ne {\textbf{0}}$, two types of Wilson line, longitudinal and transverse, enter the game [62, 64]. It turns out that in light-cone gauge, the longitudinal Wilson lines unitarise, but the transverse ones do not [64]. The opposite happens in Lorentz gauge where the transverse Wilson lines unitarise while the longitudinal ones do not. As a consequence, using the light-cone gauge at ${\textbf{b}}_T\ne {\textbf{0}}$ brings no advantage over the Lorentz gauge in that the presence of the transverse Wilson lines invalidates the argument of a smaller number of diagrams.³ Yet, in light-cone gauge, the gluon propagator remains more complicated than in Lorentz gauge. In conclusion, in the GTMD case (as well as in the TMD one), the Lorentz gauge seems to be a preferable choice for perturbative calculations and is the one adopted here.

The choice of the gauge finally allows us to write the explicit form of the soft function:

$$\begin{aligned} {\hat{S}}_i({\textbf{b}}_T)= & {} \frac{1}{N_i}\text{ Tr}_c\langle 0| W_{{\overline{n}},i} ({\textbf{b}}_T) W_{n,i}^\dag ({\textbf{b}}_T) W_{n ,i} ({\textbf{0}}) W_{{\overline{n}},i}^\dag ({\textbf{0}})|0\rangle ,\nonumber \\{} & {} \quad i=q,g, \end{aligned}$$

(3)

with $N_q=N_c=3$ and $N_g=N_c^2-1=8$ being the dimensions of the fundamental and adjoint representations of the colour group, respectively. A trace over the colour indices is indicated by $\text{ Tr}_c$. The Wilson line, also entering the definition of the unsubtracted GTMD correlators in Eq. (2), is given by

$$\begin{aligned} W_{v,i} ({\textbf{b}}_T) = {\mathcal {P}}\exp \left[ -igt_\alpha ^{[i]} v_\mu \int _0^\infty {\textrm{d}}s\, A^\mu _\alpha ({\textbf{b}}_T+sv)\right] , \end{aligned}$$

(4)

where $t_\alpha ^{[i]}$ are the generators of SU(3) in the fundamental ($i=q$) or adjoint ($i=g$) representations. It is interesting to observe that the soft function in Eq. (3) reduces to one at ${\textbf{b}}_T={\textbf{0}}$. This can immediately be seen by plugging the Wilson line into Eq. (4) and using twice the identity $W_{v,i}^\dag ({\textbf{0}}) W_{v,i} ({\textbf{0}})={\mathbb {I}}_{N_i\times N_i}$ whose trace cancels against the factor $1/N_i$ in Eq. (3). We will explicitly verify this property up to one loop below. This explains why the soft function in ${\textbf{k}}_T$-integrated correlators plays no role.

A graphical representation of the unsubtracted GTMD correlators defined in Eq. (2) is given in Fig. 1. Here, the double lines with an arbitrary number of gluons attaching to the hadronic blobs represent the (expansion of the) Wilson line. Notice that the actual Wilson line that connects the space-time points $-\eta /2$ and $\eta /2$ runs back and forth along the light-cone direction defined by n (future- or past-pointing, according to the process) and that the two branches are connected along the transverse direction defined by ${\textbf{b}}_T$ at light-cone infinity. This latter connection is suppressed in Lorentz gauge and can thus be neglected.

A graphical representation of the soft function defined in Eq. (3) is instead given in Fig. 2. Two pairs of longitudinal Wilson lines in n and ${\overline{n}}$ directions are joined at the origin and at the space-time point with transverse displacement ${\textbf{b}}_T$. An arbitrary number of gluons is exchanged including loop corrections represented by the grey blob. Each gluon attachment to the Wilson lines comes with an SU(3) generator $t_\alpha ^{[i]}$ in the appropriate representation, as indicated by the index i associated to the Wilson lines.

Fig. 2

Graphical representation of the soft function defined in Eq. (3)

2.2 Renormalisation and evolution of GTMDs

The subtracted GTMD correlator $\hat{{\mathcal {F}}}_{i/H}$ in Eq. (1), while being free of rapidity singularities, is UV divergent. In order to renormalise it, we need to separately renormalise the soft function ${\hat{S}}_i$ and the unsubtracted GTMD correlator ${\hat{\Phi }}_{i/H}$.

The soft function is renormalised multiplicatively as follows:

$$\begin{aligned}{} & {} {S}_i({\textbf{b}}_T,\mu ,\zeta ,\delta ) \nonumber \\{} & {} \quad =\lim _{\epsilon \rightarrow 0}{Z}_{S,i}^{-1}({\textbf{b}}_T,Q, \zeta ,\mu ,\delta ,\epsilon ) {\hat{S}}_i({\textbf{b}}_T,Q,\delta ,\epsilon ). \end{aligned}$$

(5)

For completeness, in this expression, we have explicitly reported all the possible dependences that emerge from a perturbative calculation. Specifically, we have the scale Q and the regulator $\delta$ introduced by the regularisation of the rapidity divergences (to be discussed in Sect. 2.4), and the scale $\mu$ and the dimensional regulator $\epsilon$ introduced by the regularisation of the UV divergences. Also notice the presence amongst the dependences of the renormalisation constant ${Z}_{S,i}$ of the (squared) scale $\zeta$, usually dubbed rapidity scale, also discussed in Sect. 2.4.

The renormalisation of the unsubtracted GTMD correlators is also purely multiplicative, that is, it does not imply any convolution nor mixing between quark and gluon operators:

$$\begin{aligned}{} & {} \Phi _{i/H}(x,\xi ,{\textbf{b}}_T,t,\mu ,\delta )\nonumber \\{} & {} \quad =\lim _{\epsilon \rightarrow 0} Z_{\Phi ,i}^{-1}(\xi ,\mu ,\delta ,\epsilon ){\hat{\Phi }}_{i/H} (x,\xi ,{\textbf{b}}_T,t,\delta ,\epsilon ). \end{aligned}$$

(6)

This simple pattern is a direct consequence of the fact that the ${\textbf{b}}_T$ displacement prevents any contact divergence of the operators. Using Eqs. (5) and (6) in Eq. (1), one can easily deduce the renormalisation of the subtracted GTMD correlator:

$$\begin{aligned}{} & {} {\mathcal {F}}_{i/H}(x,\xi ,{\textbf{b}}_T,t,\mu ,\zeta )\nonumber \\{} & {} \quad =\lim _{\epsilon ,\delta \rightarrow 0} {Z}_{S,i}^{1/2}({\textbf{b}}_T,Q, \zeta ,\mu ,\delta ,\epsilon ) Z_{\Phi ,i}^{-1}(\xi ,\mu ,\delta ,\epsilon )\nonumber \\{} & {} \qquad \times \hat{{\mathcal {F}}}_{i/H}(x,\xi ,{\textbf{b}}_T,t,\delta ,\epsilon )\nonumber \\{} & {} \quad = \lim _{\delta \rightarrow 0}{S}_i^{-1/2}({\textbf{b}}_T,\mu ,\zeta ,\delta ) \Phi _{i/H}(x,\xi ,{\textbf{b}}_T,t,\mu ,\delta ). \end{aligned}$$

(7)

Exploiting the fact that in the limit $\epsilon ,\delta \rightarrow 0$, the bare subtracted GTMD correlator does not depend on either the renormalisation scale $\mu$ or on the rapidity scale $\zeta$, it is possible to derive the following evolution equations:

$$\begin{aligned}{} & {} \frac{{\textrm{d}}\ln {\mathcal {F}}_{i/H}(x,\xi ,{\textbf{b}}_T,t,\mu ,\zeta )}{{\textrm{d}}\ln \sqrt{\zeta }}= K_i({\textbf{b}}_T,\mu ),\nonumber \\{} & {} \frac{{\textrm{d}}\ln {\mathcal {F}}_{i/H}(x,\xi ,{\textbf{b}}_T,t,\mu ,\zeta )}{{\textrm{d}}\ln \mu }= \gamma _i(\mu ,\zeta ). \end{aligned}$$

(8)

The anomalous dimensions $K_i$ and $\gamma _i$ derive from the renormalisation constants as follows:

$$\begin{aligned}{} & {} K_i({\textbf{b}}_T,\mu )=\lim _{\epsilon ,\delta \rightarrow 0} \frac{{\textrm{d}}\ln {Z}_{S,i}({\textbf{b}}_T,Q, \zeta ,\mu ,\delta ,\epsilon )}{{\textrm{d}}\ln \zeta },\nonumber \\{} & {} \gamma _i(\mu ,\zeta )=\lim _{\epsilon ,\delta \rightarrow 0}\frac{{\textrm{d}}\ln [{Z}_{S,i}^{1/2}({\textbf{b}}_T,Q, \zeta ,\mu ,\delta ,\epsilon ) Z_{\Phi ,i}^{-1}(\xi ,\mu ,\delta ,\epsilon )]}{{\textrm{d}}\ln \mu }. \nonumber \\ \end{aligned}$$

(9)

The first equation in Eq. (8) is usually referred to as the Collins–Soper equation [66, 67], while the second is a more common renormalisation-group equation.

The anomalous dimensions in Eq. (9) can be mutually related by observing that the cross-derivative of the GTMD correlator must coincide. Indeed, it turns out that

$$\begin{aligned} \frac{{\textrm{d}} K_i({\textbf{b}}_T,\mu )}{{\textrm{d}}\ln \mu } = \frac{{\textrm{d}} \gamma _i(\mu ,\zeta )}{{\textrm{d}}\ln \sqrt{\zeta }}\equiv - \gamma _{K,i}(a_s(\mu )), \end{aligned}$$

(10)

where the anomalous dimension $\gamma _{K,i}$, usually called cusp anomalous dimension, only depends on the strong coupling $a_s=g^2/16\pi ^2=\alpha _s/4\pi$. Therefore, it is a purely perturbative quantity that, for sufficiently large scales, admits the following expansion:

$$\begin{aligned} \gamma _{K,i} (a_s(\mu )) = \sum _{n=0}^\infty a_s^{n+1}(\mu )\gamma _{K,i}^{[n]}. \end{aligned}$$

(11)

Equation (10) can be solved w.r.t. both $K_i$ and $\gamma _i$ expressing these anomalous dimensions in terms of more fundamental quantities computable in perturbation theory. To do so, we need appropriate boundary conditions. For $K_i$, the boundary condition is conveniently set at the scale $\mu =b_0/|{\textbf{b}}_T|\equiv \mu _b$, with $b_0\equiv 2{\textrm{e}}^{-\gamma _{\textrm{E}}}$ and $\gamma _{\textrm{E}}$ the Euler constant, where it admits the perturbation expansion as

$$\begin{aligned} K_i({\textbf{b}}_T,\mu _b) = \sum _{n=0}^\infty a_s^{n+1}(\mu _b)K_i^{[n]}. \end{aligned}$$

(12)

With this boundary condition, the solution to Eq. (10) reads as

$$\begin{aligned} K_i({\textbf{b}}_T,\mu ) = K_i({\textbf{b}}_T,\mu _b) - \int _{\mu _b}^{\mu } \frac{{\textrm{d}}\mu '}{\mu '} \gamma _{K,i}(a_s(\mu ')). \end{aligned}$$

(13)

We notice that, thanks to the integral in the r.h.s., this expression resums all powers in $a_s$ through its evolution, preventing the presence of potentially large logarithms. However, as we will explicitly see below, a perturbative calculation of $K_i$ at a generic fixed scale $\mu$ produces terms accompanied by powers of $\ln (\mu /\mu _b)$. Since ${\textbf{b}}_T$ is usually integrated over, these logarithms can become arbitrarily large, invalidating any fixed-order calculation. Therefore, in phenomenological applications the resummed expression of $K_i$ in Eq. (13) is to be preferred over a fixed-order calculation.

We now solve Eq. (10) for $\gamma _i$. In this case, the boundary condition is conveniently set at $\sqrt{\zeta } = \mu /\sqrt{1-\xi ^2}$ where the anomalous dimension can be expanded as

$$\begin{aligned} \gamma _i(\mu ,\mu /\sqrt{1-\xi ^2}) \equiv \gamma _{F,i}(a_s(\mu )) = \sum _{n=0}^\infty a_s^{n+1}(\mu )\gamma _{F,i}^{[n]}. \end{aligned}$$

(14)

Owing to the fact that the r.h.s. of Eq. (10) does not depend of $\zeta$, the solution to the evolution equation is particularly simple and reads

$$\begin{aligned} \gamma _i(\mu ,\zeta ) = \gamma _{F,i}(a_s(\mu )) - \gamma _{K,i}(a_s(\mu ))\ln \left( \frac{\sqrt{(1-\xi ^2)\zeta }}{\mu }\right) .\nonumber \\ \end{aligned}$$

(15)

In conclusion, the evolution equations for the GTMD correlator take the explicit form⁴

$$\begin{aligned}{} & {} \frac{{\textrm{d}}\ln {\mathcal {F}}_{i/H}(x,\xi ,{\textbf{b}}_T,t,\mu ,\zeta )}{\textrm{d}\ln \sqrt{\zeta }}\nonumber \\{} & {} \quad =K_i({\textbf{b}}_T,\mu _b) - \int _{\mu _b}^{\mu } \frac{{\textrm{d}}\mu '}{\mu '} \gamma _{K,i}(a_s(\mu ')),\nonumber \\{} & {} \frac{{\textrm{d}}\ln {\mathcal {F}}_{i/H}(x,\xi ,{\textbf{b}}_T,t,\mu ,\zeta )}{{\textrm{d}}\ln \mu }\nonumber \\{} & {} \quad =\gamma _{F,i}(a_s(\mu )) - \gamma _{K,i}(a_s(\mu ))\ln \left( \frac{\sqrt{(1-\xi ^2)\zeta }}{\mu }\right) , \end{aligned}$$

(16)

where the relevant anomalous dimensions $K_i({\textbf{b}}_T,\mu _b)$, $\gamma _{F,i}$, and $\gamma _{K,i}$ are all perturbative quantities. As a by-product of the calculation of the matching functions presented below, we will extract the leading-term coefficients $K_i^{[0]}$, $\gamma _{F,i}^{[0]}$, and $\gamma _{K,i}^{[0]}$. Unsurprisingly, they turn out to be identical to those obtained in TMD factorisation.

2.3 Matching on GPDs

As mentioned above, the transverse displacement ${\textbf{b}}_T$ is Fourier conjugated to the partonic transverse momentum ${\textbf{k}}_T$. As a consequence, when ${\textbf{b}}_T\simeq {\textbf{0}}$, which corresponds to the emission of partons with large ${\textbf{k}}_T$, it is legitimate to expect such emissions to be treatable in perturbation theory. As a matter of fact, in this regime, hard emissions can be factorised into perturbatively calculable quantities, the matching functions, that allow one to express the GTMD correlators in terms of the corresponding collinear GPD correlators. The matching formula reads

$$\begin{aligned}{} & {} {\mathcal {F}}_{i/H}(x,\xi ,{\textbf{b}}_T,t,\mu ,\zeta )\nonumber \\{} & {} \quad = \int _x^\infty \frac{{\textrm{d}}y}{y}{\mathcal {C}}_{i/k} \left( y,\frac{\xi }{x},{\textbf{b}}_T,\mu ,\zeta \right) F_{k/H} \left( \frac{x}{y},\xi ,t,\mu \right) \nonumber \\{} & {} \quad \equiv {\mathcal {C}}_{i/k}(x,\kappa ,{\textbf{b}}_T,\mu ,\zeta )\mathop {\otimes }_{x} F_{k/H}\left( x,\xi ,t,\mu \right) , \end{aligned}$$

(17)

where $F_{k/H}$, with $k=q,g$, are the renormalised GPD correlators whose explicit definition can be found, for example, in Ref. [57]. In addition, we have defined $\kappa \equiv \xi /x$.

In order to extract the matching functions ${\mathcal {C}}_{i/k}$, we make use of the concept of parton-in-parton distribution [13] (sometimes also referred to as quark-target model, see, e.g., Ref. [69]). The idea is to replace the hadronic states involved in the correlators in Eq. (2) and in their collinear analogues with partonic states. Since the action of partonic fields on partonic states is computable in perturbation theory, this enables a direct calculation. It should be pointed out that the validity of this procedure is tightly connected to factorisation and the universality of the resulting partonic distributions [13]. Since so far, to the best of my knowledge, no actual factorisation theorems involving GTMDs have been proven,⁵ here we simply assume its applicability. The parton-in-parton version of Eq. (17) is then

$$\begin{aligned}{} & {} {\mathcal {F}}_{i/j}(x,\xi ,{\textbf{b}}_T,\mu ,\zeta ) \nonumber \\{} & {} = {\mathcal {C}}_{i/k}(x,\kappa ,{\textbf{b}}_T,\mu ,\zeta )\mathop {\otimes }_{x} F_{k/j}\left( x,\xi ,\mu \right) , \end{aligned}$$

(18)

where $j=q,g$. Notice that we have dropped the dependence on t that does not participate in a partonic computation. Now all quantities involved in the matching formula admit a perturbative expansion:

$$\begin{aligned}{} & {} {\mathcal {F}}_{i/j}(x,\xi ,{\textbf{b}}_T,\mu ,\zeta ) =\sum _{n=0}^{\infty }a_s^n{\mathcal {F}}_{i/j}^{[n]}(x,\xi ,{\textbf{b}}_T,\mu ,\zeta ),\nonumber \\{} & {} {F}_{k/j}(x,\xi ,\mu ) =\sum _{n=0}^{\infty }a_s^n{F}_{k/j}^{[n]}(x,\xi ,\mu ) ,\nonumber \\{} & {} {\mathcal {C}}_{i/k}(x,\kappa ,{\textbf{b}}_T,\mu ,\zeta ) =\sum _{n=0}^{\infty }a_s^n{\mathcal {C}}_{i/k}^{[n]}(x,\kappa ,{\textbf{b}}_T,\mu ,\zeta ). \end{aligned}$$

(19)

At the lowest order in $a_s$, where no additional radiation is allowed, one finds that

$$\begin{aligned} {\mathcal {F}}_{i/j}^{[0]}(x,\xi ,{\textbf{b}}_T,\mu ,\zeta ) ={F}_{i/j}^{[0]}(x,\xi ,\mu )=D_j(\xi )\delta _{ij}\delta (1-x), \nonumber \\ \end{aligned}$$

(20)

with $D_q(\xi )=\sqrt{1-\xi ^2}$ and $D_g(\xi )=1-\xi ^2$ [57]. This immediately implies that:

$$\begin{aligned} {\mathcal {C}}_{i/k}^{[0]}(x,\kappa ,{\textbf{b}}_T,\mu ,\zeta )=\delta _{ik}\delta (1-x). \end{aligned}$$

(21)

These results allow us to determine the one-loop (${\mathcal {O}}(a_s)$) correction to ${\mathcal {C}}_{i/j}$ in terms of the one-loop GTMD and GPD parton-in-parton correlators:

$$\begin{aligned}{} & {} {\mathcal {C}}_{i/j}^{[1]}(x,\kappa ,{\textbf{b}}_T,\mu ,\zeta )\nonumber \\{} & {} \quad =D_j^{-1}(\xi )\left[ {\mathcal {F}}_{i/j}^{[1]}(x,\xi ,{\textbf{b}}_T, \mu ,\zeta )-F_{i/j}^{[1]}\left( x,\xi ,\mu \right) \right] . \end{aligned}$$

(22)

Using the perturbative expansion of the renormalised soft function $S_i$ and unsubtracted parton-in-parton GTMD correlator ${\Phi }_{i/j}$:

$$\begin{aligned}{} & {} {S}_i({\textbf{b}}_T,\mu ,\zeta ,\delta ) =\sum _{n=0}^{\infty } a_s^n {S}_i^{[n]}({\textbf{b}}_T,\mu ,\zeta ,\delta ),\nonumber \\{} & {} \Phi _{i/j}(x,\xi ,{\textbf{b}}_T,\mu ,\delta ) = \sum _{n=0}^{\infty } a_s^n \Phi _{i/j}^{[n]}(x,\xi ,{\textbf{b}}_T,\mu ,\delta ), \end{aligned}$$

(23)

with ${S}_i^{[0]}({\textbf{b}}_T,\mu ,\zeta ,\delta )=1$ (see Eq. (44)) in Eq. (7), we obtain

$$\begin{aligned}{} & {} {\mathcal {F}}_{i/j}^{[1]}(x,\xi ,{\textbf{b}}_T,\mu ,\zeta ) = \Phi _{i/j}^{[1]}(x,\xi ,{\textbf{b}}_T,\mu ,\delta ) \nonumber \\{} & {} \quad - \frac{1}{2} D_j(\xi )\delta _{ij}\delta (1-x) {S}_i^{[1]}({\textbf{b}}_T,\mu ,\zeta ,\delta ). \end{aligned}$$

(24)

Notice that, while the single terms on the r.h.s. of this equation depend on the rapidity regulator $\delta$, the l.h.s. does not. This means that this dependence must cancel in this specific combination, confirming at one loop the rapidity-divergence safety of the definition in Eq. (1). Plugging this identity into Eq. (18) finally gives

$$\begin{aligned}{} & {} {\mathcal {C}}_{i/k}^{[1]}(x,\kappa ,{\textbf{b}}_T,\mu ,\zeta )\nonumber \\{} & {} \quad =D_k^{-1}(\xi )\left[ \Phi _{i/k}^{[1]}(x,\xi ,{\textbf{b}}_T,\mu ,\delta ) -F_{i/k}^{[1]}\left( x,\xi ,\mu \right) \right] \nonumber \\{} & {} \qquad - \frac{1}{2} \delta _{ik}\delta (1-x){S}_i^{[1]}({\textbf{b}}_T,\mu ,\zeta ,\delta ). \end{aligned}$$

(25)

Therefore, the computation of the one-loop corrections to the matching functions boils down to computing parton-in-parton unsubtracted GTMD and GPD correlators, and the soft function at the same order. Moreover, the computations of $\Phi _{i/k}^{[1]}$ and of $F_{i/k}^{[1]}$ are closely related, with the latter recently presented in Ref. [57]. This similarity can be exploited to simplify the calculation. First of all, we notice that they are both made of a “real” and a “virtual” contribution:⁶

$$\begin{aligned} \Phi _{i/k}^{[1]}=\Phi _{i/k}^{[1],{\textrm{real}}}+\Phi _{i/k}^{[1],\mathrm virt},\quad F_{i/k}^{[1]} = F_{i/k}^{[1],{\textrm{real}}} + F_{i/k}^{[1],{\textrm{virt}}}. \end{aligned}$$

(26)

Since the virtual contribution to the GTMD correlator in ${\textbf{k}}_T$ space must, for kinematic reasons, be proportional to $\delta ^{(2)}({\textbf{k}}_T)$, its Fourier transform is independent from ${\textbf{b}}_T$. But this is precisely the same kinematics used to compute the virtual contribution to the GPD correlator. Therefore, one finds

$$\begin{aligned} \Phi _{i/k}^{[1],{\textrm{virt}}} = F_{i/k}^{[1],{\textrm{virt}}}, \end{aligned}$$

(27)

such that

$$\begin{aligned} \Phi _{i/k}^{[1]}-F_{i/k}^{[1]} = \Phi _{i/k}^{[1],{\textrm{real}}} -F_{i/k}^{[1], {\textrm{real}}}. \end{aligned}$$

(28)

To compute the difference in the r.h.s., we can again use the results of Ref. [57]. To do so, we observe that the only difference between $F_{i/k}^{[1], {\textrm{real}}}$, and $\Phi _{i/k}^{[1],{\textrm{real}}}$ in ${\textbf{b}}_T$ space is that the ${\textbf{k}}_T$ integral in the latter is weighted by the phase ${\textrm{e}}^{i {\textbf{b}}_T\cdot {\textbf{k}}_T}$, reflecting the displacement of the partonic fields. Specifically, one finds

$$\begin{aligned}{} & {} \Phi _{i/k}^{[1],{\textrm{real}}}(x,\xi ,{\textbf{b}}_T,\mu ,\delta ) \nonumber \\{} & {} \quad =D_k(\xi )\left[ {\mathcal {P}}_{i/k}^{[0], {\textrm{real}}}(x,\kappa ,\delta )-\epsilon {\mathcal {R}}_{i/k}^{[1]}(x,\kappa )\right] \mu ^{2\epsilon } 4\pi \nonumber \\{} & {} \qquad \times \int \frac{{\textrm{d}}^{2-2\epsilon }{\textbf{k}}_T}{(2\pi )^{2-2\epsilon }} \frac{{\textrm{e}}^{i{\textbf{b}}_T\cdot {\textbf{k}}_T}}{{\textbf{k}}_T^2}, \end{aligned}$$

(29)

where ${\mathcal {P}}_{i/k}^{[0], {\textrm{real}}}$ is the real (and rapidity divergent) contribution to the one-loop GPD splitting functions, while the “residual” functions ${\mathcal {R}}_{i/k}^{[1]}$ were so far unknown and are computed here for the first time (see Sect. 2.5). Importantly, the exponential in the integral in Eq. (29) regulates the UV divergence for ${\textbf{k}}_T\rightarrow \infty$ [70]. As a consequence, the UV divergences of the GTMD correlator can only be in the virtual (diagonal) part, hence the simple renormalisation pattern in Eq. (6).

The ${\textbf{k}}_T$ integral in Eq. (29) can be computed analytically and evaluates to

$$\begin{aligned}{} & {} 4\pi \int \frac{{\textrm{d}}^{2-2\epsilon }{\textbf{k}}_T}{(2\pi )^{2-2\epsilon }}\frac{{\textrm{e}}^{i{\textbf{b}}_T\cdot {\textbf{k}}_T}}{{\textbf{k}}_T^2} =\pi ^\epsilon b_T^{2\epsilon }\Gamma (-\epsilon )\nonumber \\{} & {} \quad = -\frac{\pi ^\epsilon b_T^{2\epsilon }(1+\gamma _\textrm{E}\epsilon )}{\epsilon _{\textrm{IR}}}+{\mathcal {O}}(\epsilon ), \end{aligned}$$

(30)

where $b_T\equiv |{\textbf{b}}_T|$. In the rightmost equality, we have highlighted the fact that the pole for $\epsilon \rightarrow 0$ is of infrared origin. This divergence is cancelled in the difference in Eq. (28) by an analogous divergence in the UV-renormalised GPD correlator. To see this, we observe that the real part of the renormalised GPD correlator can be written as

$$\begin{aligned}{} & {} F_{i/k}^{[1],{\textrm{real}}}(x,\xi ,\mu ) \nonumber \\{} & {} \quad =D_k(\xi )S_\epsilon \left[ {\mathcal {P}}_{i/k}^{[0], {\textrm{real}}}(x,\kappa ,\delta )\left( \ln \mu ^2-\frac{\mu ^{2\epsilon }}{\epsilon _{{\textrm{IR}}}}\right) \right. \nonumber \\{} & {} \qquad \left. -\epsilon {\mathcal {R}}_{i/k}^{[1]}(x,\kappa ) \left( \frac{\mu ^{2\epsilon }}{\epsilon _{{\textrm{UV}}}}-\frac{\mu ^{2 \epsilon }}{\epsilon _{\textrm{IR}}}\right) \right] , \end{aligned}$$

(31)

where

$$\begin{aligned} S_\epsilon =\frac{(4\pi )^\epsilon }{\Gamma (1-\epsilon )} = 1 + \epsilon \left( \ln 4\pi -\gamma _{\textrm{E}}\right) +{\mathcal {O}}(\epsilon ^2). \end{aligned}$$

(32)

From Eq. (31), it is apparent that the UV divergence has been removed leaving a $\ln \mu ^2$, while the infrared divergence is still present. In addition, we also retained the next term in the expansion in powers of $\epsilon$. This term, proportional to ${\mathcal {R}}_{i/k}^{[1]}$, multiplies a scaleless integral in which UV ($1/\epsilon _{\textrm{UV}}$) and IR ($1/\epsilon _{\textrm{IR}}$) poles cancel, giving a vanishing result [71] that does not contribute to $F_{i/k}^{[1],{\textrm{real}}}$. However, since the IR contribution is cancelled by the GPD correlator, the UV pole does give a finite contribution in the combination in Eq. (28). Indeed, using Eqs. (29) and (31), Eq. (28) yields

$$\begin{aligned}{} & {} \Phi _{i/k}^{[1]}(x,\xi ,{\textbf{b}}_T,\mu ,\delta ) -F_{i/k}^{[1]}(x,\xi ,\mu ) \nonumber \\{} & {} \quad =D_k(\xi )\left[ -{\mathcal {P}}_{i/k}^{[0],{\textrm{real}}}(x,\kappa ,\delta ) \ln \left( \frac{\mu ^2}{\mu _b^2}\right) +{\mathcal {R}}_{i/k}^{[1]}(x,\kappa )\right] , \nonumber \\ \end{aligned}$$

(33)

where the limit $\epsilon \rightarrow 0$ has already been taken using the equality

$$\begin{aligned} \lim _{\epsilon \rightarrow 0}\frac{S_\epsilon -\pi ^\epsilon b_T^{2\epsilon }(1+\gamma _{\textrm{E}}\epsilon )}{\epsilon } =\ln \frac{4e^{-2\gamma _{\textrm{E}}}}{b_T^{2}} =\ln \frac{b_0^2}{b_T^2}=\ln \mu _b^2.\nonumber \\ \end{aligned}$$

(34)

Remarkably, all divergences have cancelled, leaving a finite quantity. Plugging Eq. (33) into Eq. (25), we obtain

$$\begin{aligned}{} & {} {\mathcal {C}}_{i/k}^{[1]}(x,\kappa ,{\textbf{b}}_T,\mu ,\zeta )= -{\mathcal {P}}_{i/k}^{[0],\mathrm real}(x,\kappa ,\delta ) \ln \left( \frac{\mu ^2}{\mu _b^2}\right) \nonumber \\{} & {} \quad +{\mathcal {R}}_{i/k}^{[1]}(x,\kappa ) - \frac{1}{2}\delta _{ik}\delta (1-x){S}_i^{[1]}({\textbf{b}}_T,\mu ,\zeta ,\delta ). \end{aligned}$$

(35)

It is now convenient to extract the rapidity divergence from ${\mathcal {P}}_{i/k}^{[0],{\textrm{real}}}$ in a way to facilitate the cancellation against the soft function ${S}_i^{[1]}$. To this purpose, we now introduce the particular strategy that will be used to regulate the rapidity divergences. We will employ the principal-value prescription [63] that acts on eikonal propagators of the kind $1/(n\cdot k)$ and $1/({\overline{n}}\cdot k)$ by replacing them with their principal-valued version, that is

$$\begin{aligned}{} & {} \frac{1}{(n\cdot k)}\rightarrow \textrm{PV}\frac{1}{(n\cdot k)}=\frac{1}{2}\left[ \frac{1}{(n\cdot k)+i\delta (n\cdot p)}\right. \nonumber \\{} & {} \quad \left. +\frac{1}{(n\cdot k)-i\delta (n\cdot p)}\right] =\frac{(n\cdot k)}{(n\cdot k)^2+\delta ^2 (n\cdot p)^2}, \end{aligned}$$

(36)

with $\delta \rightarrow 0$, and analogously for $1/({\overline{n}}\cdot k)$ with p replaced by ${\overline{p}}$. The momenta p and ${\overline{p}}$ are to be thought as the light-like momenta ($p^2={\overline{p}}^2=0$) of two highly energetic projectiles that collide with large centre-of-mass energy $(p+{\overline{p}})^2=2p\cdot {\overline{p}}\equiv Q^2\gg \Lambda _{\textrm{QCD}}^2$. Q is precisely the scale first introduced in Eq. (5) as a consequence of the regularisation of rapidity divergences and that will appear explicitly only in the soft function. When integrating over the partonic momentum k, the longitudinal projection $(n\cdot k)$ is usually parameterised as $(n\cdot k)=(1-z) (n\cdot p)$, with the variable z running in the interval $z\in [0,1]$. One can then show that the principal-value prescription is equivalent to replacing

$$\begin{aligned} \frac{1}{1-z}\rightarrow \left( \frac{1}{1-z}\right) _+-\delta (1-z) \ln \delta , \end{aligned}$$

(37)

where the $+$-prescription is defined as

$$\begin{aligned} \int _0^1{\textrm{d}}z\left( \frac{1}{1-z}\right) _+f(z) = \int _0^1\frac{f(z)-f(1)}{1-z}, \end{aligned}$$

(38)

for a test function f well behaved at $z=1$.

It is opportune to stress that the principal-value prescription defined in Eq. (36) is not guaranteed to work beyond one-loop accuracy. As pointed out in Refs. [60, 72, 73], a better regularisation that preserves the exponentiation of the eikonal propagators and that ensures the correct implementation of the zero-bin subtraction at all orders should be adopted. To this purpose, exponentially regulated Wilson lines are introduced in the definition of the relevant objects. See also Ref. [74] for another implementation of the exponential regularisation of rapidity divergences.

Having defined the specific procedure to regularise the rapidity divergences, we can replace ${\mathcal {P}}_{i/k}^{[0],{\textrm{real}}}$ in Eq. (35) with the full splitting function ${\mathcal {P}}_{i/k}^{[0]}$ using the following equality:

$$\begin{aligned}{} & {} {\mathcal {P}}_{i/k}^{[0],{\textrm{real}}}(x, \kappa ,\delta )\nonumber \\{} & {} \quad = {\mathcal {P}}_{i/k}^{[0]}(x, \kappa ) -{\mathcal {P}}_{i/k}^{[0],{\textrm{virt}}}(x, \kappa ,\delta )\nonumber \\{} & {} \quad = {\mathcal {P}}_{i/k}^{[0]}(x,\kappa ) -\delta _{ik}\delta (1-x)2C_i \nonumber \\{} & {} \quad \left[ K_i-\ln (1-\xi ^2)-2\int _0^1\frac{{\textrm{d}}z}{1-z}\right] \nonumber \\{} & {} \quad = {\mathcal {P}}_{i/k}^{[0]}(x,\kappa ) -\delta _{ik}\delta (1-x)2C_i \nonumber \\{} & {} \quad \left[ K_i-\ln (1-\xi ^2)+2\ln \delta \right] . \end{aligned}$$

(39)

In the second line, we have used the explicit form of the virtual contribution to the splitting function extracted from Ref. [57]. The colour factors $C_i$ are defined as

$$\begin{aligned} C_g=C_A=N_c=3,\quad C_q=C_F=\frac{N_c^2-1}{2N_c}=\frac{4}{3}, \end{aligned}$$

(40)

and the values of the coefficient $K_i$ are

$$\begin{aligned} K_q = \frac{3}{2},\quad K_g = \frac{11C_A-4n_fT_R}{6C_A}, \end{aligned}$$

(41)

with $T_R=1/2$ and $n_f$ the number of active quark flavours. Finally, in the third line of Eq. (39), we have used Eq. (37) to regularise the divergent integral in z. We can now plug Eq. (39) into Eq. (35) to obtain

$$\begin{aligned}{} & {} {\mathcal {C}}_{i/k}^{[1]}(x,\kappa ,{\textbf{b}}_T,\mu ,\zeta )\nonumber \\{} & {} \quad =-{\mathcal {P}}_{i/k}^{[0]}(x,\kappa ) \ln \left( \frac{\mu ^2}{\mu _b^2}\right) +{\mathcal {R}}_{i/k}^{[1]}(x,\kappa )\nonumber \\{} & {} \qquad +\delta _{ik}\delta (1-x)\bigg [2C_i\left( K_i-\ln (1-\xi ^2) +2\ln \delta \right) \ln \left( \frac{\mu ^2}{\mu _b^2}\right) \nonumber \\{} & {} \qquad - \frac{1}{2} {S}_i^{[1]}({\textbf{b}}_T,\mu ,\zeta ,\delta )\bigg ]. \end{aligned}$$

(42)

Since the splitting functions ${\mathcal {P}}_{i/k}^{[0]}$ have been computed in Ref. [57], all we are left with to obtain the one-loop correction to the matching functions is to compute the soft function ${S}_i^{[1]}$ and the residual functions ${\mathcal {R}}_{i/k}^{[1]}$. This will be done next in Sects. 2.4 and 2.5, respectively.

2.4 The soft function

Many calculations of the soft function to one-loop order and beyond exist in the literature. For example, one-loop results can be found in Refs. [58, 61, 72, 74,75,76]. As of today, also two-loop [72, 74, 77] and three-loop [78, 79] results have been presented. As mentioned above, the soft function is affected by rapidity divergences. They are caused by the presence of eikonal propagators like $1/(n\cdot k)$ and cannot be regulated by means of dimensional regularisation. Therefore, an additional regulator needs to be introduced. We will use the principal-value prescription [63] introduced in the previous section that, to the best of my knowledge, is used here for the first time to compute the soft function. The soft function is also affected by UV divergences that we regularise by means of dimensional regularisation.

The ultimate goal of this section is to compute the first two terms, $n=0,1$, of the series in Eq. (23) for the renormalised soft function. To do so, we start from the perturbative series of the bare soft function:

$$\begin{aligned} {\hat{S}}_i({\textbf{b}}_T,Q,\delta ,\epsilon ) = \sum _{n=0}^{\infty } a_s^n {\hat{S}}_i^{[n]}({\textbf{b}}_T,Q,\delta ,\epsilon ). \end{aligned}$$

(43)

The leading order, $n=0$, is trivial. To compute it, we simply approximate all the Wilson lines in Eq. (3) with the unity operator in the appropriate colour-group representation. This immediately gives

$$\begin{aligned}{} & {} {\hat{S}}_i^{[0]}({\textbf{b}}_T,Q,\delta ,\epsilon ) \nonumber \\{} & {} \quad =\frac{1}{N_i}\text{ Tr}_c\langle 0| {\mathbb {I}}_{N_i\times N_i}{\mathbb {I}}_{N_i\times N_i}{\mathbb {I}}_{N_i\times N_i}{\mathbb {I}}_{N_i\times N_i}|0\rangle =1. \end{aligned}$$

(44)

Fig. 3

Non-vanishing one-loop diagrams contributing to the soft function. The vertical bars indicate the final-state cuts that set all partons that they cross on the mass shell

We now move on to computing the one-loop correction, that is, the term $n=1$ in the series in Eq. (43). One-loop diagrams in which the gluon attaches to Wilson lines pointing in the same light-cone direction are proportional to either $n^2=0$ or ${\overline{n}}^2=0$, and thus give no contribution. On the contrary, the diagrams displayed in Fig. 3 have a gluon that attaches to Wilson lines pointing in different directions. Therefore, they are proportional to $n\cdot {\overline{n}}=1$, giving a non-null contribution. Applying standard Feynman rules in Lorentz gauge [13] and using the principal-value prescription to regulate the eikonal propagators, diagrams (a) and (b) along with their respective Hermitian conjugates separately evaluate to

$$\begin{aligned} {\hat{S}}_i^{[a+a^\dag ]}= & {} - a_s4 C_i (4\pi \mu ^2)^{\epsilon }\frac{1+\cos (\pi \epsilon )}{2}\nonumber \\{} & {} \times \left( Q^2\delta ^2\right) ^{-\epsilon }\Gamma ^2(\epsilon ) \Gamma (1-\epsilon ),\nonumber \\ {\hat{S}}_i^{[b+b^\dag ]}= & {} - a_s4C_i (4\pi \mu ^2)^{\epsilon } \Gamma (-\epsilon )\nonumber \\{} & {} \times \Bigg [\Bigg (\frac{b_T^{2}}{4}\Bigg )^\epsilon \Bigg (\ln \frac{Q^2 \delta ^2}{\mu _b^2}-\psi (-\epsilon )-\gamma _{\textrm{E}}\Bigg )\nonumber \\{} & {} -\frac{1+\cos (\pi \epsilon )}{2}\left( Q^2\delta ^2\right) ^{-\epsilon }\Gamma (\epsilon ) \Gamma (1-\epsilon )\Bigg ]. \end{aligned}$$

(45)

The net result for the one-loop correction to the bare soft function is

$$\begin{aligned}{} & {} {\hat{S}}_i^{[1]}({\textbf{b}}_T,Q,\delta ,\epsilon ) \nonumber \\{} & {} \quad =a_s^{-1}\left[ {\hat{S}}_i^{[a+a^\dag ]}+{\hat{S}}_i^{[b+b^\dag ]}\right] \nonumber \\{} & {} \quad = - 4C_i (4\pi \mu ^2)^{\epsilon } \Gamma (-\epsilon ) \left( \frac{b_T^{2}}{4}\right) ^\epsilon \left( \ln \frac{Q^2 \delta ^2}{\mu _b^2}-\psi (-\epsilon )-\gamma _{\textrm{E}}\right) \nonumber \\{} & {} \quad =4C_i\left( -\frac{S_\epsilon ^2}{\epsilon ^2}+\frac{1}{2}\ln ^2\left( \frac{\mu ^2}{\mu _b^2}\right) -\left( \frac{S_\epsilon }{\epsilon }+\ln \left( \frac{\mu ^2}{\mu _b^2}\right) \right) \ln \left( \frac{\mu ^2}{Q^2\delta ^2}\right) \right. \nonumber \\{} & {} \qquad \left. +\frac{\pi ^2}{12}+{\mathcal {O}}(\epsilon )\right) . \end{aligned}$$

(46)

This result agrees with that of Ref. [72] despite the different (but closely related) regularisation of the rapidity divergences.

A couple of additional comments are in order. First, from the second line of Eq. (46), we immediately see that ${\hat{S}}_i^{[1]}({\textbf{b}}_T ={\textbf{0}},Q,\delta ,\epsilon )=0$. This reflects the fact that the soft function must unitarise for ${\textbf{b}}_T ={\textbf{0}}$. Considering the leading-order result in Eq. (44), this requirement is indeed fulfilled to one-loop accuracy by our calculation. It is also interesting to observe that ${\hat{S}}_i^{[b+b^\dag ]}$ is UV finite:

$$\begin{aligned} \lim _{\epsilon \rightarrow 0}{\hat{S}}_i^{[b+b^\dag ]}({\textbf{b}}_T,Q,\mu ,\delta ,\epsilon ) = 4C_i\left( \frac{1}{2}\ln ^2\left( \frac{Q^2\delta ^2 }{\mu _b^2}\right) + \frac{\pi ^2}{12}\right) .\nonumber \\ \end{aligned}$$

(47)

This was to be expected because the transverse displacement ${\textbf{b}}_T$ regulates the UV divergence in real diagrams. This is a further consistency check of the calculation.

We can now renormalise the soft function. This is done multiplicatively as in Eq. (5). Considering the first two orders of the bare soft function, Eqs. (44) and (46), the renormalisation constant ${Z}_{S,i}$ in the $\overline{\text{ MS }}$ scheme reads⁷

$$\begin{aligned}{} & {} {Z}_{S,i}({\textbf{b}}_T,Q, \zeta ,\mu ,\delta ,\epsilon )= 1 \nonumber \\{} & {} \quad - a_s4C_i\left[ \frac{S_\epsilon ^2}{\epsilon ^2}+ \frac{S_\epsilon }{\epsilon }\ln \left( \frac{\mu ^2}{Q^2\delta ^2} \right) +\ln \left( \frac{\mu ^2}{\mu _b^2}\right) \ln \left( \frac{\zeta }{Q^2}\right) \right] \nonumber \\{} & {} \quad +{\mathcal {O}}(a_s^2). \end{aligned}$$

(48)

Here, the arbitrary rapidity scale $\zeta$ is introduced as a proxy to parameterise finite ${\mathcal {O}}(a_s)$ contributions. Using this renormalisation constant, one obtains the first two coefficients of the perturbative expansion of the renormalised soft function:

$$\begin{aligned} {S}_i^{[0]}({\textbf{b}}_T,\mu ,\zeta ,\delta )= & {} 1,\nonumber \\ {S}_i^{[1]}({\textbf{b}}_T,\mu ,\zeta ,\delta )= & {} 2 C_i\left( 4\ln \left( \frac{\mu ^2}{\mu _b^2}\right) \ln \delta +\ln ^2\left( \frac{\mu ^2}{\mu _b^2}\right) \right. \nonumber \\{} & {} \left. -2\ln \left( \frac{\mu ^2}{\mu _b^2}\right) \ln \left( \frac{\mu ^2}{\zeta }\right) +\frac{\pi ^2}{6}\right) . \end{aligned}$$

(49)

Evidently, ${S}_i^{[1]}$ is still affected by a rapidity divergence for $\delta \rightarrow 0$. This divergence is precisely what is needed to cancel an analogous divergence of the unsubtracted GTMD correlator. Indeed, plugging ${S}_i^{[1]}$ into Eq. (42), one finds

$$\begin{aligned}{} & {} {\mathcal {C}}_{i/k}^{[1]}(x,\kappa ,{\textbf{b}}_T,\mu ,\zeta )\nonumber \\{} & {} \quad =-{\mathcal {P}}_{i/k}^{[0]}(x,\kappa ) \ln \left( \frac{\mu ^2}{\mu _b^2}\right) +{\mathcal {R}}_{i/k}^{[1]}(x,\kappa )\nonumber \\{} & {} \qquad -\delta _{ik}\delta (1-x)2C_i\left[ \frac{1}{2}\ln ^2\left( \frac{\mu ^2}{\mu _b^2}\right) \right. \nonumber \\{} & {} \qquad \left. -\left( K_i + \ln \left( \frac{\mu ^2}{(1-\xi ^2)\zeta }\right) \right) \ln \left( \frac{\mu ^2}{\mu _b^2}\right) +\frac{\pi ^2}{12}\right] , \end{aligned}$$

(50)

where the rapidity divergence has cancelled, leaving an explicitly finite result. We now just need to compute the residual functions ${\mathcal {R}}_{i/k}^{[1]}$ to complete the picture.

2.5 The unsubtracted GTMD correlators

In this section, we present the results for the one-loop unsubtracted parton-in-parton GTMD correlators $\Phi _{i/k}$ as defined in Eq. (2). This calculation is functional to the extraction of the residual functions ${\mathcal {R}}_{i/k}^{[1]}$ and of the renormalisation constant $Z_{\Phi , i}$ in Eq. (6).

Practically speaking, the computation is closely related to that of the splitting functions ${\mathcal {P}}_{i/k}^{[0]}$ presented in Ref. [57]. As already discussed in Sect. 2.3, the real contribution to the unsubtracted parton-in-parton GTMD $\Phi _{i/k}^{[1],{\textrm{real}}}$ is a simple generalisation of the real contribution to the respective GPD $F_{i/k}^{[1],{\textrm{real}}}$ (see Eq. (29)). In addition, as again discussed in Sect. 2.3, the virtual contribution $\Phi _{i/k}^{[1],{\textrm{virt}}}$ coincides with $F_{i/k}^{[1],{\textrm{virt}}}$. The only new element w.r.t. Ref. [57] is that we need to retain not only the pole part of the correlator but also the first finite correction in the expansion in powers of $\epsilon$. This is precisely where the residual functions ${\mathcal {R}}_{i/k}^{[1]}$ reside.

The main difference between the theoretical setup of this work and that of Ref. [57] is the choice of the gauge. While here we opted for the Lorentz gauge, in Ref. [57] the light-cone gauge was used. However, a simple analysis shows that this difference is immaterial, at least at one-loop accuracy. The reasoning runs as follows. The GPD correlators are by construction gauge-invariant, signifying that the light-cone-gauge result of Ref. [57] must coincide with a calculation in Lorenz gauge.⁸ Now, the differences between the GPD and the GTMD correlators are

• 1. a transverse displacement of the operator that causes the appearance of a phase in the real diagrams (see Sect. 2.3); and

• 2. the appearance of a transverse Wilson line at light-cone infinity.

We have already discussed in Sect. 2.3 how to take care of the phase, while, in Lorentz gauge, the transverse Wilson line gives no contribution. The conclusion is that, at least at one loop, the light-cone-gauge calculation of the GPD correlators can be entirely “recycled” to obtain the GTMD correlators in Lorentz gauge. Based on this conclusion, we refer the reader to Ref. [57] for the details of the calculation. Here, we only report the result for the bare unsubtracted GTMD correlators:

$$\begin{aligned}{} & {} {\hat{\Phi }}_{i/k}^{[1]}(x,\xi ,{\textbf{b}}_T,\mu ,\delta )\nonumber \\{} & {} \quad = D_k(\xi )\Bigg [-\frac{S_\epsilon }{\epsilon _\textrm{IR}}{\mathcal {P}}_{i/k}^{[0], \mathrm real}(x,\kappa ,\delta )-{\mathcal {P}}_{i/k}^{[0], \mathrm real}(x,\kappa ,\delta ) \nonumber \\{} & {} \qquad \times \ln \frac{\mu ^2}{\mu _b^2}+ {\mathcal {R}}_{i/k}^{[1]}(x,\kappa )\nonumber \\{} & {} \qquad + \delta _{ik}\delta (1-x) 2C_i\left( K_i-\ln (1-\xi ^2)+2\ln \delta \right) \frac{\mu ^{2\epsilon }S_{\epsilon }}{\epsilon _\textrm{UV}}\nonumber \\{} & {} \qquad +{\mathcal {O}}(\epsilon )\Bigg ], \end{aligned}$$

(51)

where we have labelled the poles as IR or UV. Using the leading-order result

$$\begin{aligned} {\hat{\Phi }}_{i/k}^{[0]}(x,\xi ,{\textbf{b}}_T,\mu ,\delta ) = D_k(\xi ) \delta _{ik}\delta (1-x), \end{aligned}$$

(52)

that descends from Eq. (20), we can extract the $\overline{\text{ MS }}$ renormalisation constant $Z_{\Phi ,i}$ as defined in Eq. (6) up to one loop:

$$\begin{aligned}{} & {} Z_{\Phi ,i}(\xi ,\mu ,\delta ,\epsilon ) \nonumber \\{} & {} \quad = 1+a_s2C_i\left( K_i-\ln (1-\xi ^2)+2\ln \delta \right) \frac{S_{\epsilon }}{\epsilon }\nonumber \\{} & {} \qquad +{\mathcal {O}}(a_s^2). \end{aligned}$$

(53)

The residual functions ${\mathcal {R}}_{i/k}^{[1]}$ in Eq. (51) are best given for non-singlet and singlet GTMD combinations, respectively defined as

$$\begin{aligned} \Phi ^- = \sum _{q}\Phi _{q/k}-\Phi _{{\overline{q}}/k},\qquad \Phi ^+= \begin{pmatrix} \sum _{q}\Phi _{q/k}+\Phi _{{\overline{q}}/k}\\ \Phi _{g/k} \end{pmatrix}.\nonumber \\ \end{aligned}$$

(54)

The sums run over the active quark flavours and the anti-quark correlators are defined as

$$\begin{aligned} \Phi _{{\overline{q}}/k}(x,\ldots ) = - \Phi _{q/k}(-x,\ldots ). \end{aligned}$$

(55)

In addition, it is convenient to introduce the following decomposition:

$$\begin{aligned}{} & {} {\mathcal {R}}^{\pm ,[1]}(y,\kappa )\nonumber \\{} & {} \quad =\theta (1-y){\mathcal {R}}_1^{\pm ,[1]}(y,\kappa )+\theta (\kappa -1){\mathcal {R}}_2^{\pm ,[1]}(y,\kappa ). \end{aligned}$$

(56)

In the non-singlet sector one finds

$$\begin{aligned} \left\{ \begin{array}{rcl} {\mathcal {R}}_1^{-,[1]}(y,\kappa ) &{}=&{} 2C_F\frac{1-y}{1-\kappa ^2y^2},\\ {\mathcal {R}}_2^{-,[1]}(y,\kappa ) &{}=&{} 2 C_F \frac{(1-\kappa )y}{1-\kappa ^2y^2}, \end{array} \right. , \end{aligned}$$

(57)

while the singlet results are

$$\begin{aligned}{} & {} \left\{ \begin{array}{l} {\mathcal {R}}_{1,qq}^{+,[1]}(y,\kappa ) = {\mathcal {R}}_1^{-,[1]}(y,\kappa ),\\ {\mathcal {R}}_{2,qq}^{+,[1]}(y,\kappa ) = 2 C_F\frac{1-\kappa }{\kappa (1-\kappa ^2y^2)}, \end{array} \right. \nonumber \\{} & {} \left\{ \begin{array}{l} {\mathcal {R}}_{1,qg}^{+,[1]}(y,\kappa ) = 4n_fT_R\frac{y(1-y)}{(1-\kappa ^2 y^2)^2},\\ {\mathcal {R}}_{2,qg}^{+,[1]}(y,\kappa ) = 4n_fT_R\frac{(1-\kappa )y^2}{(1-\kappa ^2 y^2)^2}, \end{array} \right. \nonumber \\{} & {} \left\{ \begin{array}{l} {\mathcal {R}}_{1,gq}^{+,[1]}(y,\kappa ) = 2 C_F \frac{(1-\kappa ^2)y}{1-\kappa ^2y^2},\\ {\mathcal {R}}_{2,gq}^{+,[1]}(y,\kappa ) = -2C_F \frac{1-\kappa ^2}{\kappa (1-\kappa ^2y^2)}, \end{array} \right. \nonumber \\{} & {} \left\{ \begin{array}{l} {\mathcal {R}}_{1,gg}^{+,[1]}(y,\kappa ) = 8 C_A\frac{\kappa ^2y(1-y)}{(1-\kappa ^2y^2)^2},\\ {\mathcal {R}}_{2,gg}^{+,[1]}(y,\kappa ) = C_A\frac{(1-\kappa )(1+\kappa -(1-7\kappa ) \kappa ^2 y^2)}{\kappa (1-\kappa ^2y^2)^2}. \end{array} \right. \end{aligned}$$

(58)

All the expressions for ${\mathcal {R}}_1$ and ${\mathcal {R}}_2$ above are singular at $y=1/\kappa$. It turns out that this singularity cancels out in the combination in Eq. (56) for ${\mathcal {R}}^{-,[1]}$, ${\mathcal {R}}_{qq}^{+,[1]}$, and ${\mathcal {R}}_{gq}^{+,[1]}$. Unfortunately, this does not happen for ${\mathcal {R}}_{qg}^{+,[1]}$ and ${\mathcal {R}}_{gg}^{+,[1]}$. Specifically, for $\kappa >1$ and $y<1$, one is left with undefined integrals of the following kind:

$$\begin{aligned} J=\int _x^1\frac{{\textrm{d}}y}{1-\kappa y} f(y). \end{aligned}$$

(59)

Nevertheless, appropriately interpreting these integrals as principal values, one can use a trick presented in Ref. [57] to obtain a numerically amenable representation that reads

$$\begin{aligned} J= & {} \int _x^1\frac{{\textrm{d}}y}{1-\kappa y} \left[ f(y)-f\left( \frac{1}{\kappa }\right) \left( 1+\theta \left( \kappa y-1\right) \frac{1-\kappa y}{\kappa y}\right) \right] \nonumber \\{} & {} +f\left( \frac{1}{\kappa }\right) \frac{1}{\kappa }\ln \left[ \frac{\kappa (1-\kappa x)}{\kappa -1}\right] . \end{aligned}$$

(60)

However, these integrals are clearly divergent for $\kappa \rightarrow 1^+$ because in this limit, the pole of the integrand becomes an end-point singularity. As a consequence, as we will explicitly see in Sect. 3, the GTMD correlator will manifest a divergent behaviour around $x=\xi$. We also notice that these divergences, being limited to ${\mathcal {R}}_{qg}^{+,[1]}$ and ${\mathcal {R}}_{gg}^{+,[1]}$, emerge from the convolution of the matching functions with the gluon GPD. Therefore, one may argue that a vanishing gluon GPD at $x=\xi$ would guarantee the finiteness of the GTMD correlator at $x=\xi$. However, this constraint seems to be hard to fulfil in full generality in that it should hold at all scales.

In this respect, I can foresee two possible scenarios. In the first scenario, the gluon GPD obeys a specific functional constraint such that it remains null at $x=\xi$ at all scales. Incidentally, this constraint seems to be possible to realise by means of the so-called “shadow” GPDs recently introduced in Ref. [80]. In the second scenario, the gluon GPD is different from zero at $x=\xi$, causing the GTMD correlator to diverge in this point. To the best of my knowledge, while even a discontinuity at $x=\xi$ at the level of GPDs would cause the breaking of leading-twist collinear factorisation in deeply virtual Compton scattering (DVCS) [81], a divergent GTMD correlator at $x=\xi$ is not in principle forbidden. It is interesting to observe that, beyond the leading twist, discontinuities at $x=\xi$ in DVCS partonic amplitudes do appear even in the collinear case; see for example Refs. [82,83,84]. However, in those papers, it has been shown that these discontinuities are such that they do not break factorisation within the relevant accuracy.⁹

As discussed in Ref. [85] in the context of DVCS, it can also be argued that the singularity at $x=\xi$ is a consequence of a kinematic enhancement due to the emission of soft-collinear gluons that spoil the perturbative convergence in this region. Analogously to the case of soft gluons in the threshold region $x\rightarrow 1$ in deep-inelastic scattering [86, 87], their emission can be resummed to all orders in the strong coupling $\alpha _s$ producing better-behaved results around $x=\xi$ [85].¹⁰ In any event, this point deserves further investigations.

2.6 Forward limit of the matching functions

To the best of my knowledge, this is the first calculation of GTMD matching functions to one-loop accuracy ever performed. As such, there are no previous calculations that can be used to compare these results to. However, these matching functions need to tend to their well-known forward counterpart used to match TMDs onto collinear PDFs.

To perform this check, we first observe that the structure of Eq. (50) is the same as that found in the TMD case: see for instance Eq. (4.8) of Ref. [58]. The forward limit is realised by taking the limit $\xi \rightarrow 0$ that is equivalent to $\kappa \rightarrow 0$. In this limit, we have already verified in Ref. [57] that the GPD splitting functions ${\mathcal {P}}_{i/k}^{[0]}$ tend to the Altarelli–Parisi splitting functions. We can therefore safely simplify Eq. (50) by setting $\mu =\mu _b$ in such a way that all logarithms vanish. Expressing the matching functions in the basis defined in Eq. (54), we find

$$\begin{aligned} {\mathcal {C}}^{[1],-}(y,\kappa ,{\textbf{b}}_T,\mu _b,\zeta )= & {} {\mathcal {R}}^{[1],-}(y,\kappa )-C_q\frac{\pi ^2}{6}\delta (1-y),\nonumber \\ {\mathcal {C}}_{ik}^{[1],+}(y,\kappa ,{\textbf{b}}_T,\mu _b,\zeta )= & {} {\mathcal {R}}_{ik}^{[1],+}(y,\kappa )-C_i\frac{\pi ^2}{6}\delta _{ik} \delta (1-y).\nonumber \\ \end{aligned}$$

(61)

When taking the limit $\kappa \rightarrow 0$, the term proportional to $\theta (\kappa -1)$ in the decomposition of ${\mathcal {R}}^{[1],\pm }$ in Eq. (56) drops. This leaves only the term proportional to $\theta (1-y)$ that, introduced in the convolution integral as defined in Eq. (17), turns it into a standard Mellin convolution. Therefore, all we need to do is to take the limit $\kappa \rightarrow 0$ of the functions ${\mathcal {R}}_1^{[1],\pm }$. This is easily done using Eqs. (57) and (58), and the result is

$$\begin{aligned}{} & {} \lim _{\kappa \rightarrow 0}{\mathcal {C}}^{[1],-}(y,\kappa ,{\textbf{b}}_T,\mu _b,\zeta ) = \lim _{\kappa \rightarrow 0}{\mathcal {C}}_{qq}^{[1],+}(y,\kappa ,{\textbf{b}}_T,\mu _b,\zeta )\nonumber \\{} & {} \quad = 2C_F(1-y)-C_F\frac{\pi ^2}{6}\delta (1-y),\nonumber \\{} & {} \lim _{\kappa \rightarrow 0}{\mathcal {C}}_{qg}^{[1],+}(y,\kappa ,{\textbf{b}}_T,\mu _b,\zeta ) = 4n_f T_R y(1-y),\nonumber \\{} & {} \lim _{\kappa \rightarrow 0}{\mathcal {C}}_{gq}^{[1],+}(y,\kappa ,{\textbf{b}}_T,\mu _b,\zeta ) = 2C_Fy,\nonumber \\{} & {} \lim _{\kappa \rightarrow 0}{\mathcal {C}}_{gg}^{[1],+}(y,\kappa ,{\textbf{b}}_T,\mu _b,\zeta ) = -C_A\frac{\pi ^2}{6}\delta (1-y), \end{aligned}$$

(62)

that indeed coincides with the TMD results (see, e.g., Ref. [60]).

2.7 Anomalous dimensions

As discussed in Sect. 2.2, the knowledge of the renormalisation constants of soft function and unsubtracted GTMD correlators allows us to extract the anomalous dimensions that govern the evolution of the normalised GTMD correlators (see Eq. (8)). Using Eq. (9) with the normalisation constants $Z_{S,i}$ and $Z_{\Phi ,i}$ respectively given in Eqs. (48) and (53), such that

$$\begin{aligned}{} & {} {Z}_{S,i}^{1/2}Z_{\Phi ,i}^{-1} =1 - a_s 2C_i \nonumber \\{} & {} \quad \left[ \frac{S_\epsilon ^2}{\epsilon ^2}+\frac{S_\epsilon }{\epsilon }\left( K_i+\ln \left( \frac{\mu ^2}{(1-\xi ^2)Q^2}\right) \right) \right. \nonumber \\{} & {} \quad \left. +\ln \left( \frac{\mu ^2}{\mu _b^2}\right) \ln \left( \frac{\zeta }{Q^2}\right) \right] +{\mathcal {O}}(a_s^2), \end{aligned}$$

(63)

one immediately finds

$$\begin{aligned}{} & {} K_i({\textbf{b}}_T,\mu ) = -a_s4C_i\ln \left( \frac{\mu ^2}{\mu _b^2}\right) +{\mathcal {O}}(a_s^2),\nonumber \\{} & {} \gamma _i(\mu ,\zeta ) = a_s 4C_i\left( K_i+\ln \left( \frac{\mu ^2}{(1-\xi ^2)\zeta }\right) \right) +{\mathcal {O}}(a_s^2). \nonumber \\ \end{aligned}$$

(64)

Notice that for computing the total derivative of ${Z}_{S,i}^{1/2}Z_{\Phi ,i}^{-1}$ w.r.t. $\mu$ as prescribed by the second equation in Eq. (9), we have used the identity

$$\begin{aligned}{} & {} \frac{{\textrm{d}} \ln [{Z}_{S,i}^{1/2}Z_{\Phi ,i}^{-1}]}{{\textrm{d}}\ln \mu } = \frac{\partial \ln [{Z}_{S,i}^{1/2}Z_{\Phi ,i}^{-1}]}{\partial \ln \mu }\nonumber \\{} & {} \quad +\,2(-\epsilon a_s+\beta (a_s)) \frac{\partial \ln [{Z}_{S,i}^{1/2}Z_{\Phi ,i}^{-1}]}{\partial a_s}, \end{aligned}$$

(65)

where $\beta (a_s)={\mathcal {O}}(a_s^2)$ is the four-dimensional QCD $\beta$ function that governs the running of the strong coupling:

$$\begin{aligned} \frac{{\textrm{d}}a_s}{{\textrm{d}}\ln \mu ^2}=\beta (a_s(\mu )). \end{aligned}$$

(66)

The cusp anomalous dimension is extracted using Eq. (10) and equals¹¹

$$\begin{aligned} \gamma _{K,i}(a_s) = a_s8C_i + {\mathcal {O}}(a_s^2). \end{aligned}$$

(67)

Finally, Eqs. (64) and (67) allow us to extract the leading-order coefficients:

$$\begin{aligned} \begin{array}{l} K_i^{[0]} = 0,\\ \\ \gamma _{F,i}^{[0]} = 4C_iK_i,\\ \\ \gamma _{K,i}^{[0]} = 8C_i. \end{array} \end{aligned}$$

(68)

These results coincide with those obtained in the TMD framework (see, e.g., Ref. [88] and references therein for a recent overview). It is therefore legitimate to expect that the same holds true to all orders such that we can use the known higher-order corrections to these quantities to evolve GTMDs to higher accuracies. We will do so in the next section.

3 Numerical setup

In this section, we present numerical results showing how the matching functions computed above can be used to reconstruct realistic GTMDs by combining: a sensible model for GPDs along with collinear off-forward evolution, perturbative TMD evolution, and a model for the non-perturbative transverse effects as determined by recent TMD extractions.

The limitation of this procedure is that it can only be applied to GTMDs that have, at the same time, a projection on a GPD and on a TMD.¹² This largely restricts the range of possibilities. In addition, to give an estimate of such GTMDs, it is highly desirable that the GPDs and TMDs on which they project are, to some extent, phenomenologically known.

Using the results of Refs. [35, 38], the renormalised GTMD correlator, Eq. (7), has the following tensorial decomposition:

$$\begin{aligned}{} & {} {\mathcal {F}}_{i/H} = \frac{1}{2M}{\overline{u}}(P_\textrm{out})\left[ F_{1,1}^{i}+\frac{i\sigma ^{{\textbf{k}}_Tn}}{n\cdot P}F_{1,2}^{i}+\frac{i\sigma ^{\mathbf {\Delta }_Tn}}{n\cdot P}F_{1,3}^{i}\right. \nonumber \\{} & {} \quad \left. +\frac{i\sigma ^{{\textbf{k}}_T\mathbf {\Delta }_T}}{M^2}F_{1,4}^{i}\right] u(P_{\textrm{in}}), \end{aligned}$$

(69)

with u and ${\overline{u}}$ being, respectively, the spinors of the incoming and outgoing hadron H having mass M, and where we have used the shorthand notation $a_\mu b_\nu \sigma ^{\mu \nu }\equiv \sigma ^{ab}$. The scalar coefficients $F_{1,l}$, with $l=1,\dots ,4$, are the actual GTMD distributions (GTMDs for short). Each of them can be split into a T-even and a T-odd part, both real, as follows:

$$\begin{aligned} F_{1,l}^{i}=F_{1,l}^{i,e}+i F_{1,l}^{i,o}. \end{aligned}$$

(70)

The peculiarity of $F_{1,l}^{i,e}(F_{1,l}^{i,o})$ is that it conserves (flips) the sign upon swapping of the Wilson line from past to future pointing and vice versa.

The one GTMD that fulfils all the requirements stated above is $F_{1,1}^{i,e}$. As a matter of fact, in ${\textbf{b}}_T$ space and at small ${\textbf{b}}_T$, $F_{1,1}^{i,e}$ matches on the following combination of GPDs:

$$\begin{aligned}{} & {} F_{1,1}^{i,e}(x,\xi ,b_T,t,\mu ,\zeta ) = {\mathcal {C}}_{i/j}(x,\kappa , b_T,\mu ,\zeta )\nonumber \\{} & {} \quad \mathop {\otimes }_{x}\left[ (1-\xi ^2)H_{j}(x,\xi ,t,\mu )-\xi ^2 E_{j}(x,\xi ,t,\mu )\right] , \end{aligned}$$

(71)

where we have replaced ${\textbf{b}}_T$ with $b_T$ everywhere because $F_{1,1}^{i,e}$ does not depend of the direction on ${\textbf{b}}_T$. Currently, different models for $H_j$ and $E_j$ exist. In addition, the forward limit of $F_{1,1}^{i,e}$ is the fully unpolarised TMD $f_{1,i}$, that is

$$\begin{aligned} \lim _{\xi ,t\rightarrow 0} F_{1,1}^{i,e}(x,\xi , b_T,t,\mu ,\zeta ) = f_{1,i}(x, b_T,\mu ,\zeta ). \end{aligned}$$

(72)

Accurate phenomenological determinations from data of $f_{1,i}$, with $i=q$, have recently been presented.

In order to compute $F_{1,1}^{i,e}$ numerically, we employ a standard procedure used in TMD factorisation (see, e.g., Ref. [15]). Let us first state the complete formula and then comment on it:

$$\begin{aligned}{} & {} \nonumber F_{1,1}^{i,e}(x,\xi ,b_T,t,\mu ,\zeta ) \\{} & {} \nonumber \quad = {\mathcal {C}}_{i/j}(x,\kappa , b_*,\mu _{b_*},\mu _{b_*}^2)\nonumber \\{} & {} \qquad \mathop {\otimes }_{x}\left[ (1-\xi ^2)H_{j}(x,\xi ,t,\mu _{b_*})-\xi ^2 E_{j}(x,\xi ,t,\mu _{b_*})\right] \nonumber \\{} & {} \qquad \times R_i\left[ (\mu ,\zeta )\leftarrow (\mu _{b_*},\mu _{b_*}^2)\right] \nonumber \\{} & {} \qquad \times f_{\textrm{NP}}(x,b_T,(1-\xi ^2)\zeta ). \end{aligned}$$

(73)

First, we have introduced the so-called $b_*$ prescription with the following particular functional form [14, 15, 17]:

$$\begin{aligned} b_*\equiv b_*(b_T)=\frac{b_0}{Q}\left( \frac{1-\exp \left( - \frac{b_T^4Q^4}{b_0^4}\right) }{1-\exp \left( -\frac{b_T^4}{b_0^4}\right) }\right) ^{\frac{1}{4}}, \end{aligned}$$

(74)

and the associated scale $\mu _{b_*}=b_0/b_*(b_T)$. The scale Q is usually identified with the physical hard scale of the process under consideration. In all cases, one must have $Q\simeq \mu \simeq \sqrt{\zeta }$. Since the ${\textbf{k}}_T$-space GTMD is obtained through a Fourier transform of the ${\textbf{b}}_T$-space expression, the scope of the $b_*$ prescription is to prevent the scale $\mu _{b_*}$ from becoming too low to enter the non-perturbative regime when $b_T\rightarrow \infty$. Indeed, if $\mu _{b_*}$ became of the order of $\Lambda _{\textrm{QCD}}$ or smaller, any perturbative calculation would be invalidated. To this purpose, Eq. (74) is engineered in a way that

$$\begin{aligned} \lim _{b_T\rightarrow \infty }\mu _{b_*}(b_T)= 1~\text{ GeV }. \end{aligned}$$

(75)

In addition, Eq. (74) is also such that

$$\begin{aligned} \lim _{b_T\rightarrow 0}\mu _{b_*}(b_T)= Q. \end{aligned}$$

(76)

This last property is not strictly needed, but it helps keep under better control the large-${\textbf{k}}_T$ region [89, 90].

While the $b_*$ prescription ensures the applicability of perturbation theory, it remains necessary to keep into account non-perturbative transverse effects. This is done through the function $f_{\textrm{NP}}$ in the third line of Eq. (73) that parameterises non-perturbative large-$b_T$ effects. See Refs. [15,16,17] for recent global determinations of $f_{\textrm{NP}}$ from experimental data in the TMD framework. For our numerical estimates, we will use the determination of $f_{\textrm{NP}}$ from Ref. [15], henceforth referred to as PV19 (for Pavia 2019), for both quark and gluon distributions.¹³$f_{\textrm{NP}}$ in Eq. (73) only depends on x, $b_T$, and the combination $(1-\xi ^2)\zeta$ (the latter related to the evolution pattern). However, in the GTMD case, this function should in principle also depend on $\xi$ and t. Since $f_{\textrm{NP}}$ from PV19 is obtained from a fit of TMDs where $\xi$ and t are identically zero, we do not have access to the non-perturbative transverse dependence on these variables. Nonetheless, this dependence is expected to be mild in that most of the effect is accounted for by the collinear GPDs.¹⁴

The first line of Eq. (73) corresponds to the matching onto the collinear GPDs as in Eq. (71). However, the scales at which the matching is performed are chosen to be $\mu =\sqrt{\zeta }=\mu _{b_*}$ so that the matching functions are free of potentially large logarithms. The GPDs $H_i$ and $E_i$ are computed using the Goloskokov–Kroll (GK) model [91,92,93] at the initial scale $\mu _0=2$ GeV and appropriately evolved to $\mu _{b_*}$ through collinear evolution.

The evolution to the hard scales $\mu$ and $\zeta$ is provided by the factor $R_i$ in the second line of Eq. (73), often referred to as Sudakov form factor. This factor results from the solution of the evolution equations in Eq. (16) and reads

$$\begin{aligned}{} & {} R_i\left[ (\mu ,\zeta )\leftarrow (\mu _{b_*},\mu _{b_*}^2)\right] \nonumber \\{} & {} \quad = \exp \bigg \{ K_i(b_*,\mu _{b_*}) \ln \frac{\sqrt{(1-\xi ^2)\zeta }}{\mu _{b_*}} + \int _{\mu _{b_*}}^{\mu } \nonumber \\{} & {} \qquad \frac{{\textrm{d}}\mu '}{\mu '}\left[ \gamma _{F,i} (\alpha _s(\mu ')) - \gamma _{K,i} (\alpha _s(\mu ')) \ln \frac{\sqrt{(1-\xi ^2)\zeta }}{\mu '} \right] \bigg \}.\nonumber \\ \end{aligned}$$

(77)

As discussed in Sect. 2.2, the anomalous dimensions $K_i(b_*,\mu _{b_*})$, $\gamma _{F,i}$, and $\gamma _{K,i}$ are computable in perturbation theory, and they (are expected to) coincide with their TMD counterparts.

Finally, we can obtain the $F_{1,1}^{i,e}$ in ${\textbf{k}}_T$ space by taking a 2D Fourier transform of Eq. (73) w.r.t. ${\textbf{b}}_T$ that, with abuse of notation, gives

$$\begin{aligned}{} & {} F_{1,1}^{i,e}(x,\xi ,k_T,t,\mu ,\zeta )\nonumber \\{} & {} \quad =\frac{1}{2\pi }\int _{0}^{\infty }{\textrm{d}}b_T\,b_T J_0(k_Tb_T) F_{1,1}^{i,e}(x,\xi ,b_T,t,\mu ,\zeta ), \end{aligned}$$

(78)

where $J_0$ is a Bessel function of the first kind.

The theoretical precision of the formula in Eq. (73) is usually expressed in terms of logarithmic accuracy. Table 1 of Ref. [15] tells us that the knowledge of the matching functions ${\mathcal {C}}_{i/j}$ up to ${\mathcal {O}}(\alpha _s)$ would in principle allow us to reach next-to-next-to-leading logarithmic (NNLL) accuracy. In the case of GTMDs, the only obstacle to reaching this accuracy is the evolution of collinear GPDs that, as of today, is publicly available only up to leading order [57]. However, in the following, we will use a NNLL-accurate setup, except for the GPD evolution for which the ${\mathcal {O}}(\alpha _s)$ (i.e. leading-order) splitting functions are used. This practically means that $K_i$ and $\gamma _{F,i}$ are computed at ${\mathcal {O}}(\alpha _s^2)$, $\gamma _{K,i}$ at ${\mathcal {O}}(\alpha _s^3)$, ${\mathcal {C}}_{i/j}$ at ${\mathcal {O}}(\alpha _s)$, and the evolution of the strong coupling $\alpha _s$ is computed using a $\beta$ function truncated at ${\mathcal {O}}(\alpha _s^3)$. In the following, in order to simplify the presentation, we will set $\mu =\sqrt{\zeta }=Q$.

The numerical implementation of the GTMD in Eq. (73) and its Fourier transform Eq. (78) makes use of a compound of publicly available tools. Specifically, the GK model for the GPDs $H_i$ and $E_i$ at the initial scale $\mu _0$ is provided by PARTONS [94]¹⁵; their collinear evolution and the perturbative ingredients relevant to the matching and the Sudakov form factor are provided by APFEL++ [95, 96]¹⁶; the non-perturbative function $f_{\textrm{NP}}$ and the Fourier transform are instead provided by NangaParbat [15].¹⁷ The resulting code is public¹⁸ and can be used to reproduce the results shown below.

Fig. 4

$xF_{1,1}^{u,e}$ (up quark) vs $k_T$ at $x=0.2$, $\xi =0.1$, and $t=-0.1$ GeV$^2$ for four different values of the hard scale $Q=\mu =\sqrt{\zeta }$: $Q=2$ GeV (blacked dashed curve, GPD model at the initial scale $\mu _0$), $Q=5$ GeV (solid red curve), $Q=10$ GeV (solid blue curve), and $Q=100$ GeV (solid green curve). $F_{1,1}^{u,e}$ is computed at NNLL accuracy (see text) using the GK model for the GPDs and the PV19 determination for the non-perturbative transverse component

We can finally present a selection of quantitative results. We start with Fig. 4 that shows the behaviour w.r.t. $k_T$ of the up-quark GTMD $F_{1,1}^{u,e}$ at fixed values of x, $\xi$, and t for different values of the hard scale $Q=\mu =\sqrt{\zeta }$. The curve at the GPD initial scale $Q=\mu _0$ is also shown. The typical broadening of the $k_T$ distribution caused by the Sudakov form factor as Q increases is clearly observed [56, 97]. We also note that the curve at $\mu _0$ becomes negative at large $k_T$ as a consequence of the specific functional form $f_{\textrm{NP}}$ used here. The gluon distribution behaves very similarly; therefore, we do not show the corresponding plot.

Fig. 5

$xF_{1,1}^{u,e}$ (up quark) vs $k_T$ at $x=0.2$, $t=-0.1$ GeV$^2$, and $Q=\mu =\sqrt{\zeta }=10$ GeV for five different values of the skewness: $\xi =0$ (blacked dashed curve), $\xi =0.1$ (solid red curve), $\xi =0.3$ (solid blue curve), $\xi =0.5$ (solid green curve), and $\xi =0.7$ (solid orange curve). $F_{1,1}^{u,e}$ is computed at NNLL accuracy (see text) using the GK model for the GPDs and the PV19 determination for the non-perturbative transverse component

It is also interesting to look at how the $k_T$ dependence of $F_{1,1}^{i,e}$ changes with $\xi$. To this purpose, Fig. 5 shows curves for the up-quark distribution at fixed values of x, t, and Q for different values of $\xi$. The black dashed line corresponds to $\xi =0$ that is, up to a fully non-perturbative t dependence of the underlying GPDs, the “partonic” forward limit. The behaviour in $\xi$ is clearly non-monotonic but exhibits a strong suppression as the value of $\xi$ approaches one. Again, the gluon distribution does not present any further salient features and thus is not shown.

Fig. 6

$x^2F_{1,1}^{g,e}$ (gluon) vs x at $k_T=1$ GeV, $t=-0.1$ GeV$^2$, and $Q=\mu =\sqrt{\zeta }=10$ GeV for four different values of the skewness: $\xi = 0.1$ (red curve), $\xi = 0.3$ (blue curve), $\xi = 0.5$ (green curve), and $\xi = 0.7$ (orange curve). $F_{1,1}^{g,e}$ is computed at NNLL accuracy (see text) using the GK model for the GPDs and the PV19 determination for the non-perturbative transverse component

Finally, we consider the behaviour of $F_{1,1}^{i,e}$ as a function of the longitudinal momentum fraction x. As pointed out in Sect. 2.5, ${\mathcal {R}}_{qg}^{+,[1]}$ and ${\mathcal {R}}_{gg}^{+,[1]}$ cause a divergence at $x=\xi$ in both quark and gluon GTMD correlators. To study this effect, we concentrate on the gluon GTMD $F_{1,1}^{g,e}$, and in Fig. 6, we show this distribution (weighted by $x^2$ to improve the readability of the plot) as a function of x at fixed values of $k_T$, t, and Q for different values of $\xi$. It is evident that indeed the curves exhibit a divergence at $x=\xi$ that is a manifestation of the divergence of the principal-valued integral in Eq. (60) for $\kappa \rightarrow 1^+$. The conclusion is that GTMDs obtained through perturbative matching onto GPDs beyond tree level are undefined at $x=\xi$. Whether this is an acceptable feature of GTMDs remains an open question. As already mentioned in Sect. 2.5, a possibility to remove the divergence is to require the gluon GPDs to vanish at $x=\xi$ at all scales.

4 Conclusions

The main result of this paper is the calculation of the complete set of one-loop unpolarised GTMD matching functions. These perturbative functions allow one to express GTMDs at large partonic transverse momentum ${\textbf{k}}_T$ (or equivalently at low ${\textbf{b}}_T$) in terms of GPDs. They are obtained employing a rapidity-divergence-free definition of the GTMD correlator [56] that combines the soft function with the unsubtracted GTMD correlator. In order to carry out the calculation, the principal-value regulator for rapidity divergences has been used, allowing us to obtain one-loop results for the soft function, confirming previous results, and the unsubtracted GTMD correlator, that are instead original. By renormalising the UV divergences of these objects, we have also extracted the one-loop anomalous dimensions, again confirming known results. In addition, as a consistency check, we have verified that in the limit $\xi \rightarrow 0$, the one-loop GTMD matching functions thus obtained reproduce the well-know TMD results.

In the last part of the paper, we have presented a quantitative study by implementing the matching functions and combining them with other perturbative and non-perturbative ingredients necessary to fully reconstruct the GTMDs $F_{1,1}^{i,e}$, with $i=q,g$, at NNLL accuracy. We have studied their ${\textbf{k}}_T$ dependence along with the scale evolution and the dependence on the skewness $\xi$. An interesting consequence of our calculation is that, beyond tree level, GTMD matching functions yield GTMDs that exhibit a pole at $x=\xi$. The origin of this singularity is unclear and may signal the need of additional theoretical constraints on the underlying gluon GPDs. A more extensive investigation on this point is left for a future study. Finally, we stress that the numerical results presented here are readily accessible in that all of the relevant ingredients are available in public codes. A code that consistently combines them to reproduce the plots of this paper is also released with the paper.

The calculation presented here is limited to the unpolarised twist-2 Lorentz structure (see Eq. (2)). The extension of this study to the remaining twist-2 structures is planned for the future. This will require computing the corresponding splitting functions along the lines of Ref. [57] and the one-loop corrections to the unsubtracted GTMD correlators. A different approach based on the background-field technique and working in position space at the level of operators (and not matrix elements) was proposed in Ref. [98]. Originally used to compute the operator product expansion of TMDs, this method can potentially be used also to extract the GTMD matching functions. Once this calculation will be done, it will be interesting to compare the ensuing expressions against the results obtained in this paper.

Phenomenological applications of the GTMDs presented here are also envisaged. Specifically, the computation of measurable (and possibly measured) cross sections using the GTMDs presented in this paper would allow us to gauge the accuracy of the calculation and the reliability of the underlying non-perturbative ingredients, such as the GPDs and the TMDs.