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

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 GTMDs are presented, highlighting some salient features.


Introduction
At present, the study of the hadronic structure is a particularly lively field of research. If on the one hand this is a very interesting topic in itself, 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, also the study of unpolarised collinear fragmentation functions (FFs) has 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 energymomentum 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 US (EIC) [29], has revived the interest in GPDs that is now a rapidly growing field. In addition, new technologies on lattice [30,31] have made firstprinciple computations of GPDs more accessible, giving additional momentum to their study.
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. To 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 rapiditydivergence-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 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. Sect. 2.2 is devoted to the renormalisation of the 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-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 behavior of GTMDs under different points of view, commenting on some peculiar features. Finally, in Sect. 4 we will draw our conclusions.

Theoretical setup and results
As argued in the pioneering Ref. [56], a sound definition of GTMDs, i.e. 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 Fourierconjugate variable of the partonic momentum k T , denoted by 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 b T 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. Ref. [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 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: zero-bin (or soft) subtraction, 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, e.g., 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.

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 b T space which simply amounts to taking the Fourier transform of the 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: whereŜ i is the appropriate soft function andΦ i/H is the unsubtracted GTMD correlator. The quark and gluon correlators respectively read: where we have used the following shorthand definitions: η = yn + b T 1 , P in/out = P ± ∆/2, t = ∆ 2 , ξ = 2n · ∆/n · P . Here, ψ q indicates the spinor of the quark flavour q and F µν a the gluon field strength, while W n,i is the Wilson line in the n direction. The integrals run between −∞ and +∞ and a summation over repeated Lorentz indices is understood. Also the index j in Φ g/H is summed over and it runs over the (physical) transverse components, i.e. j = 1, 2. n and n (the latter appearing below) are light-like four-vectors, n 2 = n 2 = 0, and their scalar product evaluates to n · n = 1. 2 Also notice that the parton index i, labeling both the soft functionŜ 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 overŜ i andΦ i/H indicates that these quantities are bare, i.e. they are UV divergent in four dimensions.
The definition of the GTMD correlators is still incomplete because we have not specifiedŜ 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 (∂ · A = 0) and the light-cone gauge (n · A = 0). When considering k T -integrated quantities, which corresponds to setting b T = 0, the light-cone gauge is often preferred. The reason is that the Wilson lines 1 Here, the transverse vector bT is to be interpreted as a four-dimenesional embedding with null time and longitudinal components and transverse component coinciding with bT . This notation is repeatedly used below.
2 In light-cone coordinates defined as a µ = (a + , a − , aT ) with a ± = (at ± az)/ √ 2, the components of n and n are usually chosen to be n µ = (0, 1, 0) and n µ = (1, 0, 0). Figure 1: Graphical representation of the quark (left) and gluon (right) unsubtracted GTMD correlators defined in Eq. (2). 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 k T is left unintegrated, so that b T = 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 b T = 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. 3 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: with N q = N c = 3 and N g = N 2 c − 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 Tr c . The Wilson line, also entering the definition of the unsubtracted GTMD correlators in Eq. (2), is given by: where t [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 b T = 0. This can immediately be seen by plugging the Wilson line into Eq. (4) and using twice the identity W † v,i (0)W v,i (0) = I N i ×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 in k T -integrated correlators the soft function 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 −η/2 and η/2 runs back and forth along the light-cone direction defined by n (future-or past-pointing, according to the process) and the two branches are connected along the transverse direction defined by 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.

Renormalisation and evolution of GTMDs
The subtracted GTMD correlatorF 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Ŝ i and the unsubtracted GTMD correlatorΦ i/H . The soft function is renormalised multiplicatively as follows: 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 δ introduced by the regularisation of the rapidity divergences (to be discussed in Sect. 2.4), the scale µ and the dimensional regulator 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 ζ, usually dubbed rapidity scale, also discussed in Sect. 2.4. The renormalisation of the unsubtracted GTMD correlators is also purely multiplicative, i.e. it does not imply any convolution nor mixing between quark and gluon operators: This simple pattern is a direct consequence of the fact that the 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: Exploiting the fact that, in the limit , δ → 0, the bare subtracted GTMD correlator does not depend on either the renormalisation scale µ or on the rapidity scale ζ, it is possible to derive the following evolution equations: The anomalous dimensions K i and γ i derive from the renormalisation constants as follows: The first equation in Eq. (8) is usually referred to as 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: where the anomalous dimension γ K,i , usually called cusp anomalous dimension, only depends on the strong coupling a s = g 2 /16π 2 = α s /4π. Therefore, it is a purely perturbative quantity that, for sufficiently large scales, admits the expansion: Eq. (10) can be solved w.r.t. both K i and γ 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 where it admits the perturbation expansion: With this boundary condition, the solution to Eq. (10) reads: 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 µ produces terms accompanied by powers of ln(µ/µ b ). Since 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 γ i . In this case, the boundary condition is conveniently set at √ ζ = µ/ 1 − ξ 2 where the anomalous dimension can be expanded as: Owing to the fact that the r.h.s. of Eq. (10) does not depend of ζ, the solution to the evolution equation is particularly simple and reads: In conclusion, the evolution equations for the GTMD correlator take the explicit form: where the relevant anomalous dimensions K i (b T , µ b ), γ F,i , and γ 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 K,i . Unsurprisingly, they turn out to be identical to those obtained in TMD factorisation.

Matching on GPDs
As mentioned above, the transverse displacement b T is Fourier conjugated to the partonic transverse momentum k T . As a consequence, when b T 0, which corresponds to the emission of partons with large 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: where F k/H , with k = q, g, are the renormalised GPD correlators whose explicit definition can be found, e.g., in Ref. [57]. In addition, we have defined κ ≡ ξ/x. In order to extract the matching functions C i/k , we make use of the concept of parton-inparton distribution [13] (sometimes also referred to as quark-target model, see e.g. Ref. [68]). 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, 4 here we simply assume its applicability. The parton-in-parton version of Eq. (17) is then: 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: k/j (x, ξ, µ) , 4 However, factorisation is often assumed. See, for example, Refs. [52][53][54][55].
At the lowest order in a s , where no additional radiation is allowed, one finds that: with D q (ξ) = 1 − ξ 2 and D g (ξ) = 1 − ξ 2 [57]. This immediately implies that: These results allow us to determine the one-loop (O(a s )) correction to C i/j in terms of the one-loop GTMD and GPD parton-in-parton correlators: Using the perturbative expansion of the renormalised soft function S i and unsubtracted partonin-parton GTMD correlator Φ i/j : (44)) in Eq. (7), we obtain: Notice that, while the single terms on r.h.s. of this equation depend on the rapidity regulator δ, 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. (22), finally gives: 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 Φ [1] i/k and of F [1] i/k 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: 5 Since the virtual contribution to the GTMD correlator in k T space must, for kinematic reasons, be proportional to δ (2) (k T ), its Fourier transform is independent from b T . But this is precisely the same kinematics used to compute the virtual contribution to the GPD correlator. Therefore, one finds: Φ such that: Φ [1] i/k − F [1] i/k = Φ 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 [1],real i/k and Φ [1],real in the latter is weighted by the phase e ib T ·k T , reflecting the displacement of the partonic fields. Specifically, one finds: where P [0],real i/k is the real (and rapidity divergent) contribution to the one-loop GPD splitting functions, while the "residual" functions R [1] i/k 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 k T → ∞ [69]. 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 k T integral in Eq. (29) can be computed analytically and evaluates to: where b T ≡ |b T |. In the rightmost equality we have highlighted the fact that the pole for → 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: where: From Eq. (31), it is apparent that the UV divergence has been removed leaving a ln µ 2 , while the infrared divergence is still present. In addition, we also retained the next term in the expansion in powers of . This term, proportional to R [1] i/k , multiplies a scaleless integral in which UV (1/ UV ) and IR (1/ IR ) poles cancel giving a vanishing result [70] that does not contribute to F [1],real i/k . 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: where the limit → 0 has already been taken using the equality: Remarkably, all divergences have cancelled leaving a finite quantity. Plugging Eq. (33) into Eq. (25), we obtain: It is now convenient to extract the rapidity divergence from P i . To this purpose, we now introduce the particular strategy that will be used to regulate the rapidity divergences. We will employ the principalvalue prescription [63] that acts on eikonal propagators of the kind 1/(n · k) and 1/(n · k) by replacing them with their principal-valued version, that is: with δ → 0, and analogously for 1/(n · k) with p replaced by p. The momenta p and p are to be thought as the light-like momenta (p 2 = p 2 = 0) of two highly energetic projectiles that collide with large centre-of-mass energy (p + p) 2 = 2p · p ≡ Q 2 Λ 2 QCD . 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 · k) is usually parameterised as (n · k) = (1 − z)(n · p), with the variable z running in the interval z ∈ [0, 1]. One can then show that the principal-value prescription is equivalent to replacing: where the +-prescription is defined as: for a test function f well-behaved at z = 1.
With this at hand, we can replace P i/k using the following equality: In the second line we have used the explicit form of the virtual contribution to the splitting function extracted from Ref. [57]. The color factors C i are defined as: and the values of the coefficient K i are: 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: Since the splitting functions P [0] i/k 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 [1] i and the residual functions R [1] i/k . This will be done next in Sects. 2.4 and 2.5, respectively.

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,[71][72][73][74]. As of today, also twoloop [72,73,75] and three-loop [76,77] 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 · 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:Ŝ 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: We now move on to computing the one-loop correction, i.e. 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 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 · n = 1, giving a non-null contribution. Applying standard Feynman rules in Lorentz gauge [13] and using the principalvalue prescription to regulate the eikonal propagators, diagrams (a) and (b) along with their respective hermitian conjugates separately evaluate to: The net result for the one-loop correction to the bare soft function is: 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Ŝ [1] i (b T = 0, Q, δ, ) = 0. This reflects the fact that the soft function must unitarise for b T = 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Ŝ This was to be expected because the transverse displacement 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 MS scheme reads: 6 Here, the arbitrary rapidity scale ζ is introduced as a proxy to parameterise finite O(a s ) contributions. Using this renormalisation constant, one obtains the first two coefficients of the perturbative expansion of the renormalised soft function: Evidently, S [1] i is still affected by a rapidity divergence for δ → 0. This divergence is precisely what is needed to cancel an analogous divergence of the unsubtracted GTMD correlator. Indeed, plugging S [1] i into Eq. (42), one finds: where the rapidity divergence has cancelled, leaving an explicitly finite result. We now just need to compute the residual functions R [1] i/k to complete the picture. 6 Notice that in the MS scheme poles in S / , rather than in 1/ , are subtracted. The expansion in the second line of Eq. (46) is purposely written in terms of S / in view of renormalisation in the MS scheme.

The unsubtracted GTMD correlators
In this section, we present the results for the one-loop unsubtracted parton-in-parton GTMD correlators Φ i/k as defined in Eq. (2). This calculation is functional to the extraction of the residual functions R [1] i/k as well as of the renormalisation constant Z Φ,i in Eq. (6). Practically speaking, the computation is closely related to that of the splitting functions P [0] i/k presented in Ref. [57]. As already discussed in Sect. 2.3, the real contribution to the unsubtracted parton-in-parton GTMD Φ [1],real i/k is a simple generalisation of the real contribution to the respective GPD F [1],real i/k (see Eq. (29)). In addition, as again discussed in Sect. 2.3, the virtual contribution Φ [1],virt i/k coincides with F [1],virt i/k . 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 . This is precisely where the residual functions R [1] i/k 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 oneloop 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. 7 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),

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: where we have labelled the poles as IR or UV. Using the leading-order result: that descend from Eq. (20), we can extract the MS renormalisation constant Z Φ,i as defined in Eq. (6) up to one loop: The residual functions R [1] i/k in Eq. (51) are best given for non-singlet and singlet GTMD combinations, respectively defined as: The sums run over the active quark flavours and the anti-quark correlators are defined as: In addition, it is convenient to introduce the following decomposition: In the non-singlet sector one finds: while the singlet results are: All the expressions for R 1 and R 2 above are singular at y = 1/κ. It turns out that this singularity cancels out in the combination in Eq. (56) for R −, [1] , R +, [1] qq , and R +, [1] gq . Unfortunately this does not happen for R +, [1] qg and R +, [1] gg . Specifically, for κ > 1 and y < 1, one is left with undefined integrals of this kind: 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: However, these integrals are clearly divergent for κ → 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 = ξ. We also notice that these divergences, being limited to R +, [1] qg and R +, [1] gg , emerge from the convolution of the matching functions with the gluon GPDs. Therefore, one may argue that vanishing gluon GPDs at x = ξ would guarantee the finiteness of the GTMD correlator at x = ξ. However, this constraint on the gluon GPDs 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 GPDs obey a specific functional constraint such that they remain null at x = ξ at all scales. Incidentally, this constraint seems to be possible to realise by means of the so-called "shadow" GPDs recently introduced in Ref. [78]. In the second scenario, the gluon GPDs are different from zero at x = ξ causing the GTMD correlator to diverge in this point. To the best of my knowledge, while even a discontinuity at x = ξ at the level of GPDs would cause the breaking of leading-twist collinear factorisation in deeply-virtual Compton scattering (DVCS) [79], a divergent GTMD correlator at x = ξ is not in principle forbidden. It is interesting to observe that, beyond leading twist, discontinuities at x = ξ in DVCS partonic amplitudes do appear even in the collinear case, see e.g. Ref. [80][81][82]. However, in those papers it has been shown that these discontinuities are such that they do not break factorisation within the relevant accuracy. 8 As discussed in Ref. [83] in the context of DVCS, it can also be argued that the singularity at x = ξ 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 → 1 in deep-inelastic scattering [84,85], their emission can be resummed to all orders in the strong coupling α s producing better-behaved results around x = ξ [83]. 9 In any event, this point deserves further investigations.

Forward limit of the matching functions
To the best of my knowledge, this is the first calculation of GTMD matching functions to oneloop 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 ξ → 0 that is equivalent to κ → 0. In this limit, we have already verified in Ref. [57] that the GPD splitting functions P i/k tend to the Altarelli-Parisi splitting functions. We can therefore safely simplify Eq. (50) by setting µ = µ b in such a way that all logarithms vanish. Expressing the matching functions in the basis defined in Eq. (54), we find: When taking the limit κ → 0, the term proportional to θ(κ − 1) in the decomposition of R [1],± in Eq. (56) drops. This leaves only the term proportional to θ(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 κ → 0 of the functions R . This is easily done using Eqs. (57) and (58) and the result is: that indeed coincide with the TMD results (see, e.g., Ref. [60]).

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 Φ,i respectively given in Eqs. (48) and (53), such that: one immediately finds: Notice that, for computing the total derivative of Z µ as prescribed by the second equation in Eq. (9), we have used the identity: where β(a s ) = O(a 2 s ) is the four-dimensional QCD β -function that governs the running of the strong coupling: da s d ln µ 2 = β(a s (µ)) .
The cusp anomalous dimension is extracted using Eq. (10) and equals: 10 Finally, Eqs. (64) and (67) allow us to extract the leading-order coefficients: These results coincide with those obtained in the TMD framework (see, e.g., Ref. [86] 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.
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: with u and u being respectively the spinors of the incoming and outgoing hadron H having mass M , and where we have used the shorthand notation a µ b ν σ µν ≡ σ ab . The scalar coefficients F 1,l , with l = 1, . . . , 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: The peculiarity of F i,e 1,l (F i,o 1,l ) is that it conserves (flips) the sign upon swapping of the Wilson line from past to future pointing and viceversa.
The one GTMD that fulfils all the requirements stated above is F i,e 1,1 . As a matter of fact, in b T space and at small b T , F i,e 1,1 matches on the following combination of GPDs: where we have replaced b T with b T everywhere because F i,e 1,1 does not depend of the direction on b T . Currently, different models for H j and E j exist. In addition, the forward limit of F i,e 1,1 is the fully unpolarised TMD f 1,i , that is: Accurate phenomenological determinations from data of f 1,i , with i = q, have recently been presented.
In order to compute F i,e 1,1 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: First, we have introduced the so-called b * -prescription with this particular functional form [14,15,17]: and the associated scale The scale Q is usually identified with the physical hard scale of the process under consideration. In all cases, one must have Q µ √ ζ. Since the k T -space GTMD is obtained through a Fourier transform of the b T -space expression, the scope of the b * -prescription is to prevent the scale µ b * from becoming too low such to enter the non-perturbative regime when b T → ∞. Indeed, if µ b * became of the order of Λ QCD or smaller, any perturbative calculation would be invalidated. To this purpose, Eq. (74) is engineered in a way that: lim In addition, Eq. (74) is also such that: This last property is not strictly needed but it helps keep under better control the large-k T region [87,88]. 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 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 NP from experimental data in the TMD framework. For our numerical estimates, we will use the determination of f NP from Ref. [15], henceforth referred to as PV19 (for Pavia 2019), for both quark and gluon distributions. 12 f NP in Eq. (73) only depends on x, b T , and the combination (1−ξ 2 )ζ (the latter related to the evolution pattern). However, in the GTMD case this function should in principle also depend on ξ and t. Since f NP from PV19 is obtained from a fit of TMDs where ξ 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. 13 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 µ = √ ζ = µ 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 [89][90][91] at the initial scale µ 0 = 2 GeV and appropriately evolved to µ b * through collinear evolution.
The evolution to the hard scales µ and ζ 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: As discussed in Sect. 2.2, the anomalous dimension K i (b * , µ b * ), γ F,i , and γ K,i are computable in perturbation theory and they (are expected to) coincide with their TMD counterparts. Finally, we can obtain the F i,e 1,1 in k T space by taking a 2D Fourier transform of Eq. (73) w.r.t. b T that, with abuse of notation, gives: 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. Tab. 1 of Ref. [15] tells us that the knowledge of the matching functions C i/j up to O(α 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 O(α s ) (i.e. leading-order) splitting functions are used. This practically means that K i and γ , and the evolution of the strong coupling α s is computed using a β function truncated at O(α 3 s ). In the following, in order to simplify the presentation, we will set µ = √ ζ = 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 µ 0 is provided by PARTONS [92] 14 ; their collinear evolution as well as the perturbative ingredients relevant to the matching and the Sudakov form factor are provided by APFEL++ [93,94] 15 ; the non-perturbative function f NP and the Fourier transform are instead provided by NangaParbat [15] 16 . The resulting code is public 17 and can be used to reproduce the results shown below.
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 u,e 1,1 at fixed values of x, ξ, and t for different values of the hard scale Q = µ = √ ζ. The curve at the GPD initial scale Q = µ 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,95]. We also note that the curve at µ 0 becomes negative at large k T as a consequence of the specific functional form f NP used here. The gluon distribution behaves very similarly, therefore we do not show the corresponding plot.
It is also interesting to look at how the k T dependence of F i,e 1,1 changes with ξ. To this purpose, Fig. 5 shows curves for the up-quark distribution at fixed values of x, t, and Q for different values of ξ. The black dashed line corresponds to ξ = 0 that is, up to a fully nonperturbative t dependence of the underlying GPDs, the "partonic" forward limit. The behaviour in ξ is clearly non monotonic but exhibits a strong suppression as the value of ξ approaches one. Again, the gluon distribution does not present any further salient features and thus is not shown.
Finally, we consider the behaviour of F i,e 1,1 as a function of the longitudinal momentum fraction x. As pointed out in Sect. 2.5, R +, [1]  is computed at NNLL accuracy (see text) using the GK model for the GPDs and the PV19 determination for the non-perturbative transverse component. 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 ξ. It is evident that indeed the curves exhibit a divergence at x = ξ that is a manifestation of the divergence of the principal-valued integral in Eq. (60) for κ → 1 + . The conclusion is that GTMDs obtained through perturbative matching onto GPDs beyond tree level are undefined at x = ξ. 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 = ξ at all scales.

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 k T (or equivalently at low 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 soft function, confirming previous results, and 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 ξ → 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 i,e 1,1 , with i = q, g, at NNLL accuracy. We have studied their k T dependence along with the scale evolution and the dependence on the skewness ξ. An interesting consequence of our calculation is that, beyond tree level, GTMD matching functions yield GTMDs that exhibit a pole at x = ξ. 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] as well as the one-loop corrections to the unsubtracted GTMD correlators. Phenomenological applications of the GTMDs presented here are also envisaged. Specifically, the computation of measurable (and possible measured) cross sections using the GTMDs presented in this paper would allow us to gauge the accuracy of the calculation as well as the reliability of the underlying non-perturbative ingredients, such as the GPDs and the TMDs.