Two-loop coefficient function for DVCS: Vector contributions

Using the approach based on conformal symmetry we calculate the two-loop coefficient function for the vector flavor-nonsinglet contribution to deeply-virtual Compton scattering (DVCS). The analytic expression for the coefficient function in momentum fraction space is presented in the $\overline{\text{MS}}$ scheme. The corresponding next-to-next-to-leading order correction to the Compton form factor $\mathcal{H}$ for a simple model of the generalized parton distribution appears to be rather large: a factor two smaller than the next-to-leading order correction, approximately $\sim 10$\% of the tree level result in the bulk of the kinematic range, for $Q^2=4$ GeV$^2$.


Introduction
With the JLAB 12 GeV upgrade completed [1] and the Electron Ion Collider (EIC) [2] proposal having received a major boost last year, there are bright perspectives for new generation of hadron physics studies in the coming decade and beyond. The foreseen very high luminosity of these new machines will allow one to study hard exclusive and semi-inclusive reactions with identified particles in the final state with unprecedented precision. Such processes are interesting as they allow one to access hadron properties on a much more detailed level as compared to totally inclusive reactions. In this way one hopes to understand the full three-dimensional proton structure by "holographic imaging" of quark and gluon distributions in distance and momentum spaces.
Deeply-virtual Compton scattering (DVCS) [3][4][5] is generally accepted to be the "gold-plated" process that would have the highest potential impact for the transverse distance imaging. The general framework for the QCD description of DVCS is based on collinear factorization in terms of generalized parton distributions (GPDs) [6,7] and is well understood at the leading-twist level. The main challenge of these studies is that the quantities of interest (GPDs) are functions of three variables (apart from the scale dependence). Their extraction requires a massive amount of data and very high precision for both experimental and theory inputs. In the ideal case one would like to reach the same level of accuracy as in inclusive reactions. The next-to-next-to leading order (NNLO) analysis of parton distributions and fragmentation functions has become the standard in this field [8], so that the NNLO precision for DVCS is necessary as well.
The NNLO accuracy implies that one needs to derive three-loop evolution equations for GPDs and also calculate the two-loop corrections to the coefficient functions (CFs) in the operator product expansion (OPE) of the DVCS amplitude. The first part of this program, the three-loop evolution equation for flavor-nonsinglet GPDs, is completed [9] and in this work we calculate the corresponding two-loop CF for vector operators. The two-loop axial vector CF involves an additional subtlety of dealing with the γ 5 matrix in off-forward reactions and will be considered in a separate publication.
From the theory point of view the main difference of DVCS from e.g. the classical case of the deep-inelastic scattering (DIS) is that the target hadron has different momenta in the initial and the final state. As a consequence, from the OPE point of view, one has to take into account additional contributions of operators containing total derivatives. Such operators contribute both to operator mixing and to the CFs. In conformal field theories the contributions of operators with total derivatives are related to the contributions of the operators without total derivatives by symmetry transformations and do not need to be calculated separately. Although conformal invariance of QCD is broken by quantum corrections, one can hope that the symmetry of the Lagrangian can be used in some way in order to simplify the calculation. The first (unsuccessful) attempt to predict the NLO evolution kernel from conformal symmetry [10] was missing an important element: the scheme-dependent difference between the dilatation and special conformal anomalies. This problem was pointed out and solved by D. Müller [11] who subsequently developed and applied (with collaborators) the conformal symmetry based technique to DVCS and GPDs. In this way the complete two-loop mixing matrix was calculated for twist-two operators in QCD [12][13][14] and the two-loop evolution kernels derived for GPDs [15][16][17]. It was also shown [14] that conformal symmetry allows one to obtain the one-loop CF in DVCS from the known CF in DIS. Thus the complete NLO calculation of the DVCS amplitudes to the leading-twist accuracy is available. An estimate of the NNLO contribution in a special renormalization scheme ("conformal scheme") was also attempted [18][19][20].
In Ref. [21] we suggested a somewhat different implementation of the same idea. Instead of studying conformal symmetry breaking in the physical theory [12][13][14] we proposed to make use of the exact conformal symmetry of large-n f QCD in d = 4 − 2 dimensions at the Wilson-Fischer fixed point at critical coupling [22]. Due to specifics of the minimal subtraction scheme (MS) the renormalization group equations (RGEs) in the physical four-dimensional theory inherit conformal symmetry of the d = 4 − 2 theory so that the evolution kernel commutes with the generators of conformal transformations. This symmetry is exact, however, the generators are modified by quantum corrections and differ from their canonical form. The consistency relations that follow from the conformal algebra can be used in order to restore the -loop off-forward kernel from the -loop anomalous dimensions and the ( − 1)-loop result for the deformation of the generators, which is equivalent to the statement in Ref. [11]. The two-loop expression for the generator of special conformal transformations for flavor-nonsinglet operators is obtained in [23] and is the main ingredient in the calculation of the three-loop evolution kernel in Ref. [9]. For the flavor-singlet case one had first to clarify the structure of the contributions of gauge non-invariant operators to conformal Ward identities, see Ref. [22]. We used this result in [24] to derive the two-loop flavorsinglet evolution equation for light-ray operators and confirmed in this way the expression for the evolution kernel obtained originally in [16].
This work is devoted to the conformal symmetry based approach to the calculation of the DVCS CFs. The presentation is organized as follows. Sect. 2 is introductory, it contains general definitions, our notation and conventions. In Sect. 3 we present the general framework for the calculation of CFs in the OPE of two electromagnetic currents using conformal symmetry of QCD at the Wilson-Fischer fixed point in non-integer dimensions. The main statement is that calculation of the -loop off-forward CF can be reduced to the -loop forward CF, known from DIS, and the ( − 1)-loop calculation of the off-forward CF in 4 − 2 dimensions, including terms O( −1 ). The one-loop example for the application of this machinery is considered in Sect. 4 where we reproduce the corresponding well-known expression [25]. The following Sect. 5 contains our main result: the derivation of the two-loop CF for vector operators. The analytic expression for the CF is presented in momentum fraction space in the MS scheme. We find that the CF in the conformal scheme satisfies the reciprocity relation [26][27][28][29] that arises in the sum of many contributions and provides one with a highly nontrivial check of the results, in particular for the two-loop conformal anomaly [23]. Numerical estimates of the size of the NNLO correction for two popular GPD models are presented in Sect. 6. We conclude in Sect. 7 with a short summary and outlook. Some useful integrals are collected in the Appendix.

DVCS kinematics, notation and conventions
The amplitude of the DVCS process is given by a matrix element of the time-ordered product of two electromagnetic currents Here q, q are the momenta of the virtual (incoming) and real (outgoing) photons and p, p are the target (nucleon) momenta in initial and final states. We use the photon momenta to define the longitudinal plane spanned by two light-like vectors [30], In the leading-twist approximation the DVCS amplitude can be written as a sum of vector and axial-vector contributions 1 The ellipses stand for the higher twist contributions of helicity-flip amplitudes and terms ∼ q µ , which do not contribute to physics observables thanks to electromagnetic Ward identities. In this work we will consider flavor-nonsinglet contributions to the vector amplitude V in the leading-twist approximation. To this accuracy it can be written as a convolution of the CF with the corresponding GPD Here and in what follows q is a quark flavor, q = u, d, s, . . ., µ stands for the factorization scale in the MS scheme and we use established conventions for the kinematical variables The GPD F q is defined by the appropriate matrix element of the light-ray quark-antiquark operator, where [z 1 n, z 2 n] is a Wilson line. The matrix element can, in general, be parameterized as The expression for F q depends on the spin of the target (e.g. nucleon vs. pion or 4 He nucleus) whereas the CF is the same in all cases. For the most important case of DVCS off the nucleon, F q can further be decomposed in contributions of the standard GPDs, H, E, as follows [6] F q (x, ξ) Note that only the charge conjugation even C = +1 part of the GPDs F can contribute to the vector amplitude, which is reflected in antisymmetry of the CF (2.10) The CF can be calculated in perturbation theory where a s = α s (µ) 4π , (2.12) and the first two terms in this series are well known (see e.g. [7]) In this work we calculate the two-loop expression C (2) (x/ξ, Q 2 /µ 2 ). It is very important that the CFs only depend on the ratio x/ξ and are real functions at |x/ξ| < 1 (ERBL region), and can be continued analytically to the DGLAP region |x/ξ| ≥ 1 using the ξ → ξ − i prescription. This property holds to all orders and allows one to simplify the calculation assuming ξ = 1 in which case the CF in DVCS coincides (after a redefinition of kinematic variables) with the CF in the transition form factors γ * γ → meson with appropriate quantum numbers. If the latter is known, the corresponding CF in DVCS is obtained with the substitution x → x/(ξ − i ).

General framework
The OPE for the product of currents has a generic form, schematically where O N (0) are local operators of increasing dimension and C N k are the corresponding CFs. For forward matrix elements, like in DIS, contributions of the operators with total derivatives vanish identically and can be omitted. Thus only the sum over N remains and the necessary CFs C N ≡ C N,k=0 are known to three-loop accuracy [31]. In off-forward reactions, like DVCS, operators with total derivatives have to be included and one needs to calculate their coefficients, C N,k with k = 0, as well. In conformal field theories a direct calculation is not needed since the CFs of operators with total derivatives are related to the CFs of the operators without total derivatives by the symmetry transformations (for the special choice of the operator basis O N ). We will show how to restore C N,k from C N,0 using conformal algebra in what follows. QCD in four space-time dimensions, however, is not a conformal theory. The idea is to consider QCD in non-integer d = 4 − 2 dimensions at the intermediate step, for the specially chosen (critical) value of the coupling α * s such that the β(α * s ) = 0 (Wilson-Fisher fixed point). This theory is conformally invariant [22] and all renormalization constants/anomalous dimensions for composite operators in this theory in minimal subtraction schemes coincide with the renormalization constants/anomalous dimensions of the corresponding operators for the "real" QCD in d = 4 [21].
One can consider, formally, the DVCS process in a generic 4 − 2 -dimensional theory. All definitions in Sec. 2 can be taken over without modifications except for that the CFs acquire an -dependence so that and their perturbative expansion involves -dependent coefficients: Note that the tree-level CF C (0) does not depend on . We are interested in the CF in four dimensions (2.11) as a function of the coupling, whereas methods of conformal field theories allow one to calculate C * = C(α * s , ) on the line in the ( , α s ) plane where β(α * s ) = 0 so that α * s = α * s ( ) or, equivalently, with N c and n f the numbers of colors and light flavors, respectively. Trading the -dependence for the a * s dependence one can write the CF at the critical point as an expansion in the coupling alone, where, obviously, (3.5) The coefficients C (k) * can be related to the known CFs for DIS (not without effort) thanks to conformal invariance. Thus In other words, the one-loop CF in d = 4 coincides with the one-loop CF in conformal QCD in 4−2 dimensions and in order to find the two-loop CF in d = 4 one needs to know the corresponding result in conformal QCD and, in addition, terms of order in the one-loop CF. Since all necessary one-loop integrals can be calculated in terms of Γ-functions for arbitrary space-time dimensions, the latter calculation is rather straightforward and in what follows we will present the final result only. The calculation of C (2) * presents the main challenge and will be discussed in detail. The expansion in Eq. (3.5) or (3.6) can obviously be continued to higher orders. The general statement is that the -loop off-forward CF in QCD in d = 4 in the MS scheme can be obtained from the corresponding result in conformal theory (alias from the CF in the forward limit), adding terms proportional to the QCD beta-function. Such extra terms require the calculation of the corresponding − 1-loop off-forward CF in d = 4 − 2 dimensions to the ∼ O( −1 ) accuracy.

Conformal OPE
Retaining the contributions of vector operators only, the most general expression for the OPE of the product of two electromagnetic currents to the twist-two accuracy has the form where and where O µ1...µ N N (y) are the leading-twist conformal operators that transform in the proper way under conformal transformations Here and below N is the operator spin, ∆ N is its scaling dimension, is the conformal spin. We have separated in Eq. (3.7) the scale factor µ γ N to make the invariant functions A N (u), . . . , D N (u) dimensionless. Note that only vector operators with even spin N contribute to the expansion.
The conditions of conformal invariance and current conservation ∂ µ j µ = 0 lead to constraints on the functional form and also certain relations between the invariant functions A N (u), . . . , D N (u) in Eq. (3.7). One obtains The remaining functions C N (u) and D N (u) are given by the following expressions: The coefficients c N and d N are not independent and are given in terms of a N and b N by linear relations Thus the form of the OPE of the product of two conserved spin-one currents in a conformal theory is fixed up to two constants, a N (a s ) and b N (a s ), for each (even) spin N . In QCD the expansion of a N (a s ) starts at order O(a s ), so that also c N − b N = O(a s ) and d N = O(a s ).

Relating DIS and DVCS
It is convenient to fix the normalization of the leading-twist conformal operators such that where D µ = ∂ µ + igA µ and {. . .} denotes the symmetrization of all enclosed Lorentz indices and the subtraction of traces. In this way the forward matrix elements of these operators are related to moments of quark parton distributions (PDFs) Using this parametrization and taking the Fourier transform of the forward matrix element of Eq. (3.7) one obtains the OPE for the forward Compton tensor in a generic conformal theory is the Bjorken scaling variable and Here and below B(j N , j N ) is the Euler Beta function. The same expansion in QCD is usually written as [31] T µν (p, q) = N,even (we drop electromagnetic charges and the sum over flavors), so that we can identify The coefficient functions C 2 and C L contribute to the OPE for the structure functions F 2 and F L , respectively, and are known to third order in the coupling. Next, let us consider the DVCS process (2.1). In comparison to forward scattering there are two modifications. First, the position of the operator O N (ux) in Eq. (3.7) (we assume here 20) and effectively results in a shift of the momentum in the Fourier integral q → q − u∆. Second, the matrix element becomes more complicated and can be parameterized as Hence one needs a more general Fourier integral where we neglected all power-suppressed corrections ∆ 2 /Q 2 and also used that to this accuracy Up to the obvious replacement p → P there are two differences to the forward case (DIS): an extra factorū −N − 1 2 γ N , and the matrix element f N → f N (ξ). Note that the leading-twist DVCS vector amplitude (2.3) corresponds to the Lorentz structure g ⊥ µν in Eq. (2.4) in the Compton tensor (2.1), and can be traced by contributions ∼ g µν (in momentum space). Starting from the position space expression in Eq. (3.7), such terms can only originate from structures ∼ g µν and ∼ x µ x ν which involve the invariant functions A N (u) and B N (u) with the same u-dependence (3.11). The extra factorū −N − 1 2 γ N , therefore, results in both cases in the following modification: where (see above) j N = N + 1 − * + 1 2 γ N . Since this modification affects the contribution of the structures A N (u) and B N (u) in the same way, we actually do not need to consider them separately. Thus the OPE for the DVCS amplitude V in conformal QCD in non-integer dimensions can be written as and is completely determined by the forward-scattering coefficients C 1 (N ) (in d = 4 − 2 * ). For d = 4 the above expression agrees with [14,Eq. (22)]. In the next section we show how the DVCS coefficient function in momentum fraction space (2.5) can be obtained starting from this representation.

Coefficient function in momentum fraction space
GPDs are defined as matrix elements (2.8) of nonlocal light-ray operators (2.7) so that in order to relate the CFs in position or momentum fraction space to the CFs in the OPE one needs an expansion of the type We will tacitly assume that conformal operators are normalized as in Eq. (3.14). The light-ray operator [O(z 1 , z 2 )], in turn, satisfies the RGE of the form where H (evolution kernel) is an integral operator acting on the coordinates z 1 , z 2 . Conformal symmetry ensures that H commutes with the generators of SL(2, R) subgroup of the conformal group. Translation-invariant polynomials z N 12 ≡ (z 1 − z 2 ) N are eigenfunctions of the evolution kernel and the corresponding eigenvalues define the anomalous dimensions (3.28) The expansion coefficients Ψ N k in Eq. (3.25) are homogeneous polynomials in z 1 , z 2 of degree N + k − 1 and are given by a repeated application of the generator of special conformal transformations S + to the coefficient of the conformal operator, The problem is that S + in the interacting theory in the MS scheme contains a rather complicated conformal anomaly term z 12 ∆ + [23] so that finding explicit expression for S k + z N 12 is difficult. The way out is to do a rotation to the "conformal scheme" at the intermediate step using a similarity transformation defined in [9] [O(z 1 , Note that H and H obviously have the same eigenvalues (anomalous dimensions). Going over to the "boldface" operators can be thought of as a change of the renormalization scheme. The GPD in the conformal scheme, F q , is related to the GPD F q in the MS scheme by the U-"rotation" 2 and, similarly, for the CF in DVCS (2.5) The "rotated" light-ray operator O(z 1 , z 2 ) in d = 4 satisfies the RGE Looking for the operator U in the form we require that the "boldface" generators do not include conformal anomaly terms, are the canonical generators. Explicit expressions for X (1) and X (2) are given in [9]. With this choice, the generators S α on the subspace of the eigenfunctions of the operator H with a given anomalous dimension γ N take the canonical form with shifted conformal spin 36) and the eigenfunctions of the rotated kernel, HΨ N k = γ N Ψ N k , can be constructed explicitly [32]: For the forward matrix element of the light-ray operator one obtains, in our normalization, and for the rotated light-ray operator 39) where σ N are the eigenvalues of U: The generalization of this expansion to include off-forward matrix elements is completely fixed by conformal algebra and effectively amounts to the operator relation where This expression can be derived applying ∂ + to Eq. (3.39). On the one hand, taking a derivative amounts to a shift k → k − 1. On the other hand, it corresponds to an application of (−S − ) and using the commutation relation one gets a recurrence relation a N,k−1 = k(2j N + k − 1)a N,k . The overall normalization (function of N ) is fixed by the condition a N,k=0 = 1.
Taking a matrix element of Eq. (3.41) between states with fixed momenta and using that The sum over k can be evaluated with the help of Eq. (B.10) in Ref. [33] (for the special case n = 2): where P (λ N )

45)
C λ N are Gegenbauer polynomials, and . (3.46) Using this representation and the parametrization of the matrix element in Eq. (3.21) we obtain This expression should be matched to the definition of the GPD in the rotated scheme Changing variables x → x/ξ one can bring the exponential factor in Eq. (3.47) to the same form as in Eq. (3.48) and then try to interchange the order of summation and integration to obtain the answer for the GPD as a series in contributions of local conformal operators. Attempting this one would find, however, that F(x, ξ, t) vanishes outside the ERBL region |x| ≤ ξ, which is certainly wrong. This problem is well known and is caused by non-uniform convergence of a sum representation for GPDs in the DGLAP region ξ < |x|. It can be avoided, however, because the CF in DVCS only depends on the ratio x/ξ so that for our purposes we can set ξ = 1 and eliminate the DGLAP region completely. In this way we obtain and the DVCS (vector) amplitude is then given by On the other hand, from the conformal OPE (3.24) .

(3.51)
Comparing the coefficients in front of f N (ξ = 1) for these two representations we obtain , (3.52) where d = 4 − 2 * and σ N (3.40) are the eigenvalues of U on z N −1 12 . It remains to solve this equation to obtain an explicit expression for the CF in momentum fraction space and, as the last step, to go over from the "rotated" to the MS scheme. In the remaining part of this section we outline the general procedure for this calculation. To leading order (tree level) everything is simple. To this accuracy γ N = 0, λ N = 3/2, j N = N + 1 and C (0) 1 (N ) = 1 for even N and zero otherwise. Eq. (3.52) reduces to N −1 (x) = 1 , N = 2, 4, . . . , (3.53) and is solved by in agreement with Eq. (2.13). The problem is that beyond the leading order λ N depends on N in a nontrivial way and the functions P N −1 (x) are not orthogonal for different N with some simple weight function. Note, however, that they are eigenfunctions of the (exact) "rotated " evolution kernel and also eigenfunctions of the "rotated" SL(2) Casimir operator. This property suggests the following ansatz for the CF: (3.57) Using the above ansatz one obtains Comparing this expression with Eq. (3.52) we obtain i.e., the spectrum of K is given directly in terms of moments of the DIS CF and the spectrum of eigenvalues of the rotation operator U. An SL(2)-invariant operator, i.e., an operator that commutes with the generators S ±,0 (3.34) of SL(2, R) transformations, is fixed uniquely by its spectrum. Therefore, Eq. (3.59) unambiguously defines the operator K and by virtue Eqs. (3.56), (3.31) also the coefficient function C(x). In the next two sections we describe this calculation for the one-loop and the two-loop CFs, respectively.

One-loop example
In this section we take µ = Q as logarithmic terms ln k (µ/Q) in the CF can easily be restored from the evolution equation. To one-loop accuracy one needs to expand Eq. (3.59) to order O(a s ) taking into account that * = −β 0 a s + . . .. Since the tree-level CF does not depend on the space-time dimension, the -dependence in C 1 N, Q 2 µ 2 , a s , * starts at order O(a s ) and can be neglected here. Thus we only need the one-loop result for C 1 (N ) in physical d = 4 dimensions, which can be taken from [31,34]: where S k1...kn (N ) are harmonic sums [35]. We also need the one-loop flavor-nonsinglet anomalous dimension and the (one-loop) eigenvalues σ N = 1 + a s σ (1) N + . . . of the rotation matrix U = 1l + a s X (1) + . . . in Eq. (3.33). This is the only new element that requires a calculation. From Ref. [9] [X (1) In order to calculate σ N we take f (z 1 , z 2 ) = z N −1

12
and get Collecting everything we obtain K(N ) = 1 + a s K (1) (N ) + . . . , whereγ (1) N in the last line is the one-loop anomalous dimension (4.2) stripped of the color factor, γ (1) N . Note that the asymptotic expansion of the anomalous dimensionγ (1) N and therefore also K (1) (N ) at large j N is symmetric under the substitution j N → 1 − j N , alias N → −N − 1 (reciprocity relation). We will find that this relation holds to two-loops as well, in agreement with the general argumentation in [27][28][29], see the next Section.
As already mentioned, a SL(2)-invariant operator is completely determined by its spectrum. It is easy to do in the case under consideration, because an invariant operator with eigenvaluesγ (1) N is, obviously, the one-loop evolution kernel, and 1/(N (N + 1)) is nothing else but the inverse Casimir operator. The corresponding explicit expressions in position space are well known: where These kernels commute with the canonical generators S (0) ±,0 . Thus the operator K (1) can be written as The same expression holds in momentum fraction space, apart from that the kernels have to be taken in the corresponding representation. For the one-loop example considered here the transformation from position to momentum fraction space is not difficult to do and the results are available from Ref. [36]: where ω, ω are rescaled momentum fractions, ω = (1 − x)/2, and the plus distribution is defined as It remains to calculate the convolution of K (1) (x, x ) with the leading-order CF (3.56) and "rotate" the result to the MS scheme: The one-loop rotation kernel in momentum fraction space X (1) (x, x ) can also be found explicitly, Here, as above, ω = 1 2 (1 − x), ω = 1 2 (1 − x ) are rescaled momentum fractions. Collecting all terms one reproduces after some algebra the well-known expression for the one-loop CF in Eq. (2.13).
Beyond one loop, the last part of this strategy -restoration of momentum fraction kernels from the known position space results and taking the remaining convolution integrals -becomes impractical because of very complicated expressions. It can be avoided, however, using the following algorithm.
Let f (x) be a function of the momentum fraction x so that its position space analogue is The convolution of f (x) with the leading order CF C (0) (x) can be rewritten as a position space integral (4.14) Assume that the invariant operator K in position space can be written in the following form where k(α, β) is a certain weight function. Then The (momentum fraction space) convolution of the leading order CF and K can therefore be written directly in terms of the weight function k(α, β) If K is given by a product of several kernels of the type (4.15), then the right hand side of Eq. (4.17) can be written as a manyfold integral of the same type, e.g., for K = K 1 K 2 one gets Integrals of this kind can be evaluated with the help of the Maple HyperInt package by E. Panzer [37] in terms of harmonic polylogarithms, see e.g. [38]. In this way a very time consuming transformation of beyond-one-loop kernels to the momentum fraction representation can be avoided. For instance, instead of using the momentum fraction expression for X (1) in Eq. (4.12), its convolution with C (0) (x) can be calculated using Eq. (4.17) directly from the position space representation in Eq. (4.3). This leads to simple integrals where ω = (1 − x)/2. Beyond one loop, this simplification proves to be crucial.

Two-loop coefficient function
The spectrum of the invariant operator K(N ) to two-loop accuracy is obtained by expanding Eq. (3.59) to second order in the coupling. Since * = O(a s ), to this end one needs the two-loop CF in DIS in d = 4, and also terms O( ) in the one-loop DIS CF as inputs. The corresponding expressions are available from Ref. [31,39]. In addition we need to calculate the spectrum of eigenvalues σ N = 1 + a s σ N + . . . of the rotation operator (3.33) to the two-loop accuracy. Explicit expressions for the corresponding kernels X are collected in Appendix B in Ref. [9] 3 . The necessary integrals can be done analytically in terms of harmonic sums up to fourth order using computer algebra packages [35,[40][41][42][43]. The resulting expressions are rather cumbersome and we do not present them here. The final expressions for the CFs turn out to be considerably shorter thanks to many cancellations. The next step is to restore the invariant kernel K from its spectrum. This is not as simple as at one loop, because the invariant kernel has to commute with deformed SU (3) generators (3.34) (including O(a s ) terms) rather than the canonical generators (3.35). In other words, we are looking now for the integral operator (with the given spectrum) with eigenfunctions P (λ N ) N ) + . . . rather than λ N = 3/2. All expressions can of course be truncated at order a 2 s so that for the second-order contributions to the spectrum, K(N ) = . . . + a 2 s K(N ) (2) , it is sufficient to require canonical conformal invariance. However, we need to modify the first-order kernel (4.9) K (1) → K (1) = K (1) + δK (1) in such a way that K (1) has eigenfunctions P (λ N ) it commutes with deformed generators S α in Eq. (3.34) (up to terms O(a 2 s )). This can be achieved by replacing . The spectrum of eigenvalues of K (1) will, of course, differ from the spectrum of K (1) , K (1) (N ) = K (1) (N ) + O(a s ) and this difference will have to be compensated by the corresponding change in K (2) → K (2) , which is, however, straightforward. The two-loop evolution kernel can be written as [9], so that they have the same eigenfunctions up to O(a 2 s ). In other words, by throwing away the canonically invariant part H (2,inv) of the two-loop evolution kernel the eigenfunctions remain the same up to terms O(a 2 s ) that are not relevant to our accuracy. Thus, we can replace the full evolution kernel in the expression for K 1 in Eq. (5.1) by its (canonically) non-invariant part where T (1) ≡ 4C F T (1) and In addition, we need to find the inverse of the Casimir operator where One can show that to the required accuracy where H + is defined in Eq. (4.8) and , The term ∼ H 2 + is a canonically invariant operator and can be dropped for the same reasons as H (2,inv) in the evolution kernel. Thus, we choose where The terms O(a s ) in Eqs. (5.5) and (5.11) modify the spectrum of eigenvalues of K (1) as compared to K (1) and have to be subtracted from K (2) (N ). Note, that the two-loop kernel has three color structures A (N ) . (5.12) and only the ∼ C 2 F term is affected by this subtraction. We obtain where S m ≡ S m (N ). It can be checked that these expressions satisfy the reciprocity relation [27][28][29]: their asymptotic expansion at N → ∞ is symmetric under the substitution N → −N − 1. 4 It remains to find an invariant operator K (2) with the spectrum K (2) (N ). This is not very hard to do since we only need canonical SL(2)-invariance, i.e. the operators in question have to be diagonal in the basis of P With the invariant operators at hand, the DVCS coefficient function in the rotated scheme is obtained by the convolution with the leading-order CF 14) and the CF in the MS scheme (so far still in conformal QCD at the critical point) recovered as In both cases one can follow the procedure described in the previous section and avoid a Fourier transformation of the kernels to the momentum fraction space. All necessary integrals can be computed using the HyperInt package [37] in terms of the harmonic polylogarithms [38]. The two-loop CF contains contributions of three color structures We obtain Here σ = Q 2 /µ 2 and C (0) , C (1) ( ), C (2) ( ) are the CFs in d dimensions (3.2) at µ 2 = Q 2 , alias σ = 1. Note that the contribution in the second line vanishes at the critical point, β(a s , * ) = 0. For the physical case d = 4 one obtains The functions C DVCS CF in the ERBL region (|x| < ξ) is given directly by the above expressions with an obvious substitution x → x/ξ and is obtained by the analytic continuation ξ → ξ − i in the DGLAP region |x| > ξ. We have used the Mathematica package HPL-2.0 by D. Maitre [46,47] for a numerical evaluation of the harmonic polylogarithms at complex arguments.
The results are shown in Figs. 1 and 2, respectively. In the first figure, we also show on the right panel the ratios of NLO and NNLO to the leading order (LO) contribution, NLO/LO and NNLO/LO, and the ratio NNLO/NLO. It is seen that the NNLO (two-loop) and NLO (one-loop) contributions to the CF have the same sign and are negative with respect to the LO (tree-level) result in the bulk of the kinematic region apart from the end points |x| → |ξ| where the loop corrections are positive and dominated by the contributions of threshold double-logarithms (5.19). We observe that the NNLO contribution is significant. In the ERBL region, it is generally about 10% of the LO result (factor two below NLO). In the DGLAP region it is less important and in fact negligible for the real part at x/ξ > 2, and for the imaginary part at x > 4ξ.
The relative size of the contributions of different origin to the NNLO CF is illustrated in Fig. 3. We repeat the definitions for convenience: Here C (2) * (x) defines the CF in conformal QCD at the critical point and the β 0 C (1,1) (x) term describes the shift to integer d = 4 dimension. We show in Fig. 3 the ratio β 0 C (1,1) /C (2) by the black solid, β 0 C F C (2β) * /C (2) by the dashed, and C 2 F C (2P ) * /C (2) by the dash-dotted curves, respectively. The C (2) * /C (2) ratio is shown by the red dashed curve; it includes also the contribution of nonplanar diagrams ∼ C (2A) * which is very small numerically. We observe that C (2) (x) is dominated in almost the entire range of x by the simplest contribution, β 0 C (1,1) (x), that comes from the O( ) correction to the one-loop diagrams. The CF in the conformal theory at = * is small as the result of a strong cancellation between the contributions ∼ C 2 F of planar diagrams and the term proportional to the QCD β function, β 0 C F C  The dominance of the β 0 C (1,1) (x) term does not hold for very large x/ξ 0.95. In this region C (2) changes sign so that the representation in Fig. 3 is not very informative. Asymptotically, for x/ξ → 1, the NNLO CF C (2) is dominated by the Sudakov-type double-logarithmic term ∼ a 2 s C 2 F ln 4 (1 − x/ξ) in Eq. (5.19) that is part of the C and allows for a simple analytic representation: (An overall normalization is irrelevant for our purposes so we omit it). We use the value of the parameter n = 1/2 which corresponds to a valence-like PDF q(x) ∼ x −1/2 (1 − x) 3 in the forward limit. The C-even part of the GPD (6.3), H(x, ξ) − H(−x, ξ), is shown in Fig. 4 for several values of ξ.
For a numerical evaluation of the contribution to the integral in Eq. (6.2) from the DGLAP region it proves to be convenient to shift the integration contour in the complex plane. We have checked that the results do not depend on the shape of the integration contour, which is a good test of numerical accuracy. The results are presented in Fig. 5. Following [20] we show the ratios for the absolute value and the phase of the Compton form factor, H(ξ) = R(ξ) e iΦ(ξ) , (6.4) calculated to NNLO and NLO accuracy and normalized to the LO. One sees that the NNLO correction to the absolute value of the Compton form factor H is quite large: it is only about factor two smaller than the NLO correction and decreases the Compton form factor by about 10% in the whole kinematic range. The NNLO correction for the phase proves to be much smaller.

Summary
Using an approach based on conformal symmetry [21] we have calculated the two-loop CF in DVCS in MS scheme for the flavor-nonsinglet vector contributions. Analytic expressions for the CF in momentum fraction space at µ = Q are presented in Eqs. (5.17), (5.18) and in Sect. 5.1 for an arbitrary scale. Numerical estimates in Sect. 6 suggest that the two-loop contribution gives rise to a ∼ 10% correction to the Compton form factor, which is significantly above the projected accuracy at the JLAB 12 GeV facility and the Electron Ion Collider. We find an interesting hierarchy of different contributions to the two-loop CF, suggesting that the perturbative series in conformal QCD at the critical coupling is converging much faster than in the physical case. It would be interesting to check whether a similar hierarchy holds in other QCD examples, where the first few terms in perturbative expansion are known. The corresponding study goes beyond the tasks of this work. and M n τ Li 2 (τ ) + 1 2 ln 2τ = 2S −4 (n) + 7 10 ζ 2 2 , M n τ 4τ Li 2 (τ ) + 1 2 ln 2τ = S 1,3 (n) − 1 2 S 4 (n) − ζ 3 S 1 (n) + 3 10