#  On evolution kernels of twist-two operators

Yao Ji<EMAIL_ADDRESS>Physik Department T31, Technische Universität München,
D-85748 Garching, Germany  Alexander Manashov<EMAIL_ADDRESS>Max-Planck-
Institut für Physik, Werner-Heisenberg-Institut, D-80805 München, Germany
Sven-Olaf Moch<EMAIL_ADDRESS>II. Institut für Theoretische Physik,
Universität Hamburg D-22761 Hamburg, Germany

###### Abstract

The evolution kernels that govern the scale dependence of the generalized
parton distributions are invariant under transformations of the
$\mathrm{SL}(2,\mathrm{R})$ collinear subgroup of the conformal group. Beyond
one loop the symmetry generators, due to quantum effects, differ from the
canonical ones. We construct the transformation which brings the full symmetry
generators back to their canonical form and show that the eigenvalues
(anomalous dimensions) of the new, canonically invariant, evolution kernel
coincide with the so-called parity respecting anomalous dimensions. We develop
an efficient method that allows one to restore an invariant kernel from the
corresponding anomalous dimensions. As an example, the explicit expressions
for NNLO invariant kernels for the twist two flavor-nonsinglet operators in
QCD and for the planar part of the universal anomalous dimension in ${\cal
N}=4$ SYM are presented.

evolution kernels, DVCS, conformal symmetry, generalized parton distribution

††preprint:  TUM–HEP–1461/23, MPP–2023–137, DESY–23–091

##  I Introduction

The study of deeply-virtual Compton scattering (DVCS) gives one access to the
generalized parton distributions [1, 2, 3] (GPDs) that encode the information
on the transverse position of quarks and gluons in the proton in dependence on
their longitudinal momentum. In order to extract the GPDs from experimental
data one has to know, among other things, their scale dependence. The latter
is governed by the renormalization group equations (RGEs) or, equivalently,
evolution equations for the corresponding twist two operators. Essentially the
same equations govern the scale dependence of the ordinary parton distribution
functions (PDFs) in the Deep Inelastic Scattering (DIS) process. In DIS one is
interested in the scale dependence of forward matrix elements of the local
twist-2 operators and therefore can neglect the operator mixing problem
between local operators from the operator product expansion (OPE). In the
nonsinglet sector, there is only one operator for a given spin/dimension. The
anomalous dimensions of such operators are known currently with the three-loop
accuracy [4, 5] and first results at four loops are becoming available [6, 7].
In contrast, the DVCS process corresponds to non-zero momentum transfer from
the initial to the final state and, as a consequence, the total derivatives of
the local twist-two operators have to be taken into consideration. All these
operators mix under renormalization and the RGE has a matrix form. The DIS
anomalous dimensions appear as the diagonal entries of the anomalous
dimensions matrix which, in general, has a triangular form for the latter.

It was shown by Dieter Müller [8, 9] that the off-diagonal part of the
anomalous dimension matrix is completely determined by a special object, the
so-called conformal anomaly. Moreover, in order to determine the off-diagonal
part of the anomalous dimension matrix with $\ell$-loop accuracy it is enough
to calculate the conformal anomaly at one loop less. This technique was used
to reconstruct all relevant evolution kernels/anomalous dimension matrices in
QCD at two loops [10, 11, 12].

A similar approach, but based on the analysis of QCD at the critical point in
non-integer dimensions, was developed in refs. [13, 14, 15]. It was shown that
the evolution kernels in $d=4$ in the $\overline{\mathrm{MS}}$-like
renormalization scheme inherit the symmetries of the critical theory in
$d=4-2\epsilon$ dimensions. As expected, the symmetry generators deviate from
their canonical form. Corrections to the generators have a rather simple form
if they are written in terms of the evolution kernel and the conformal
anomaly. It was shown in ref. [16] that by changing a renormalization scheme
one can get rid of the conformal anomaly term in the generators bringing them
into the so-called “minimal” form. Beyond computing the evolution kernels, the
conformal approach has also been employed to calculate the NNLO coefficient
(hard) functions of vector and axial-vector contributions in DVCS [17, 18],
the latter in agreement with a direct Feynman diagram calculation [19].
Moreover, the conformal technique is also applicable to computing kinematic
higher-power corrections in two-photon processes as was recently shown in
refs. [20, 21].

In this paper we construct a similarity transformation that brings the full
quantum generators back to the canonical form. Correspondingly, the
transformed evolution kernel is invariant under the canonical
$\mathrm{SL}(2,\mathrm{R})$ transformation. Moreover, we will show that the
eigenvalues of this kernel are given by the so-called parity respecting
anomalous dimension, $f(N)$ [22, 23] which is related to the PDF anomalous
dimension spectrum $\gamma(N)$ as

$\displaystyle\gamma(N)=f\left(N+{\bar{\beta}}(a)+\frac{1}{2}\gamma(N)\right),$
(1)

where $\bar{\beta}(a)=\beta(a)/a$ with $\beta(a)$ being the QCD beta function.
The strong coupling $\alpha_{s}$ is normalized as $a=\alpha_{s}/(4\pi)$. We
develop an effective approach to restore the canonically invariant kernel from
its eigenvalues $\gamma(N)$. As an example, we present explicit expressions
for three-loop invariant kernels in QCD and ${\cal N}=4$ supersymmetric Yang-
Mills (SYM) theory. The answers are given by linear combinations of harmonic
polylogarithms [24], up to weight four in QCD and up to weight three in ${\cal
N}=4$ SYM. We also compare our exact result with the approximate expression
for the three-loop kernels in QCD given in ref. [16].

The paper is organized as follows: in section II we describe the general
structure of the evolution kernels of twist-two operators. In section III we
explain how to effectively recover the evolution kernel from the known
anomalous dimensions and present our results for the invariant kernels in QCD
and ${\cal N}=4$ SYM. Sect. IV contains the concluding remarks. Some technical
details are given in the Appendices.

##  II Kernels & symmetries

We are interested in the scale dependence of the twist-two light-ray flavor
nonsinglet operator [25]

$\displaystyle\mathcal{O}(z_{1},z_{2})$
$\displaystyle=[\bar{q}(z_{1}n)\gamma_{+}[z_{1}n,z_{2}n]q(z_{2}n)]_{\overline{\text{MS}}},$
(2)

where $n^{\mu}$ is an auxiliary light-like vector, $n^{2}=0$, $z_{1,2}$ are
real numbers, $\gamma_{+}=n^{\mu}\gamma_{\mu}$ and $[z_{1}n,z_{2}n]$ stands
for the Wilson line ensuring gauge invariance, and the subscript
$\overline{\rm MS}$ denotes the renormalization scheme. This operator can be
viewed as the generating function for local operators,
${\mathcal{O}}^{\mu_{1}\dots\mu_{N}}$ that are symmetric and traceless in all
Lorentz indices $\mu_{1}\dots\mu_{N}$.

The renormalized light-ray operator (2) satisfies the RGE

$\displaystyle\Big{(}\mu\partial_{\mu}+\beta(a)\partial_{a}+\mathbb{H}(a)\Big{)}{\mathcal{O}}(z_{1},z_{2})=0,$
(3)

where $\beta(a)$ is $d$-dimensional beta function

$\displaystyle\beta(a)=-2a\big{(}\epsilon+\beta_{0}a+\beta_{1}a^{2}+O(a^{3})\big{)},$
(4)

$\beta_{0}=11/3N_{c}-2/3n_{f}$, etc., and
$\mathbb{H}(a)=a\mathbb{H}_{1}+a^{2}\mathbb{H}_{2}+\ldots$ is an integral
operators in $z_{1},z_{2}$.

It follows from the invariance of the classical QCD Lagrangian under conformal
transformations that the one-loop kernel $\mathbb{H}_{1}$ commutes with the
canonical generators of the collinear conformal subgroup, $S_{0},S_{\pm}$,

$\displaystyle S_{-}$ $\displaystyle=-\partial_{z_{1}}-\partial_{z_{2}}\,,$
$\displaystyle S_{0}$
$\displaystyle=z_{1}\partial_{z_{1}}+z_{2}\partial_{z_{2}}+2\,,$
$\displaystyle S_{+}$
$\displaystyle=z_{1}^{2}\partial_{z_{1}}+z_{2}^{2}\partial_{z_{2}}+2z_{1}+2z_{2}\,.$
(5)

This symmetry is preserved beyond one loop albeit two of the generators,
$S_{0},S_{+}$ receive quantum corrections,
$S_{\alpha}\mapsto\widetilde{S}_{\alpha}(a)=S_{\alpha}+\Delta S_{\alpha}(a)$.
The explicit form of these corrections can be found in ref. [15].

It is quite useful to bring the generators to the following form using the
similarity transformation [16],

$\displaystyle\mathbb{H}(a)=e^{-X(a)}\mathrm{H}(a)e^{X(a)}\,,$
$\displaystyle\widetilde{S}_{\alpha}(a)=e^{-X(a)}\mathrm{S}_{\alpha}(a)e^{X(a)}\,,$
(6)

where $X(a)=aX_{1}+a^{2}X_{2}+\ldots$ is an integral operator known up to
terms of $O(a^{3})$ [11, 16]. This transformation can be thought of as a
change in a renormalization scheme.

The shift operator $\mathrm{S}_{-}$ is not modified and hence identical to
$S_{-}$ in Eq. (II), and the quantum corrections to $\rm S_{0}$ and $\rm
S_{+}$ come only through the evolution kernel

$\displaystyle\mathrm{S}_{0}(a)$
$\displaystyle=S_{0}+\bar{\beta}(a)+\frac{1}{2}\mathrm{H}(a)\,,$ (7a)
$\displaystyle\mathrm{S}_{+}(a)$
$\displaystyle=S_{+}+(z_{1}+z_{2})\left(\bar{\beta}(a)+\frac{1}{2}\mathrm{H}(a)\right)\,,$
(7b)

where $\bar{\beta}(a)=\beta_{0}a+\beta_{1}a^{2}+\cdots$ is the beta function
in four dimensions, cf. Eq. (1). The form of the generator $\mathrm{S}_{0}(a)$
is completely fixed by the scale invariance of the theory, while Eq. (7b) is
the “minimal” ansatz consistent with the commutation relation
$[\mathrm{S}_{+},\mathrm{S}_{-}]=2\mathrm{S}_{0}$. Since the operator
$\mathrm{H}(a)$ commutes with the generators,
$[\mathrm{H}(a),\mathrm{S}_{\alpha}(a)]=0$ its form is completely determined
by its spectrum (anomalous dimensions). However, since the generators do not
have the simple form as in Eq. (II), it is yet necessary to find a way to
recover the operator from its spectrum.

To this end we construct a transformation which brings the generators
$\mathrm{S}_{\alpha}(a)$ to the canonical form $S_{\alpha}$, Eq. (II). Let us
define an operator $\mathrm{T}(\mathrm{H})$:

$\displaystyle\mathrm{T}(\mathrm{H})=\sum_{n=0}^{\infty}\frac{1}{n!}\mathrm{L}^{n}\left(\bar{\beta}(a)+\frac{1}{2}\mathrm{H}(a)\right)^{n}\,,$
(8)

where $\mathrm{L}=\ln z_{12}$, $z_{12}\equiv z_{1}-z_{2}$. Recall that
$z_{1},z_{2}$ are real variables, so for $z_{12}<0$ it is necessary to choose
a specific branch of the logarithm function. Although this choice is
irrelevant for further analysis we chose the $+i0$ recipe for concreteness,
i.e., $\mathrm{L}=\ln(z_{12}+i0)$. It can be shown that the operator
$\mathrm{T}(\mathrm{H})$ intertwines the symmetry generators
$\mathrm{S}_{\alpha}(a)$ and the canonical generators, $S_{\alpha}$. Namely,

$\displaystyle\mathrm{T}(\mathrm{H})\,\mathrm{S}_{\alpha}(a)=S_{\alpha}\,\mathrm{T}(\mathrm{H}),$
(9)

see Appendix A for details. Let us also define a new kernel
$\widehat{\mathrm{H}}$ as

$\displaystyle\mathrm{T}(\mathrm{H})\,\mathrm{H}(a)=\widehat{\mathrm{H}}(a)\,\mathrm{T}(\mathrm{H}).$
(10)

It follows from Eqs. (9), (10) that the operator $\widehat{\mathrm{H}}$
commutes with the canonical generators in Eq. (II)

$\displaystyle[S_{\alpha},\widehat{\mathrm{H}}(a)]=0.$ (11)

The problem of restoring a canonically invariant operator $\widehat{\rm H}(a)$
from its spectrum is much easier than that for the operator $\mathrm{H}(a)$
and will be discussed in the next section. It can be shown that the inverse of
$\mathrm{T}(\mathrm{H})$ takes the form

$\displaystyle\mathrm{T}^{-1}(\mathrm{H})=\sum_{n=0}^{\infty}\frac{(-1)^{n}}{n!}\mathrm{L}^{n}\left(\bar{\beta}(a)+\frac{1}{2}\widehat{\mathrm{H}}(a)\right)^{n}\,,$
(12)

see Appendix A. Further, it follows from Eq. (10) that

$\displaystyle\mathrm{H}(a)$
$\displaystyle=\mathrm{T}^{-1}(\mathrm{H})\,\widehat{\mathrm{H}}(a)\,\mathrm{T}(\mathrm{H})$
$\displaystyle=\widehat{\mathrm{H}}(a)+\sum_{n=1}^{\infty}\frac{1}{n!}\mathrm{T}_{n}(a)\left(\bar{\beta}(a)+\frac{1}{2}{\mathrm{H}}(a)\right)^{n}\,.$
(13)

The operators $\mathrm{T}^{(n)}$ are defined by recursion

$\displaystyle\mathrm{T}_{n}(a)=[\mathrm{T}_{n-1}(a),\mathrm{L}]$ (14)

with the boundary condition $\mathrm{T}_{0}(a)=\widehat{\mathrm{H}}(a)$. The
$n$-th term in the sum in Eq. (II) is of order $\mathcal{O}(a^{n+1})$ so that
one can easily work out an approximation for $\mathrm{H}(a)$ with arbitrary
precision, e.g.,

$\displaystyle\mathrm{H}(a)$
$\displaystyle=\widehat{\mathrm{H}}(a)+\mathrm{T}_{1}(a)\left(1+\frac{1}{2}\mathrm{T}_{1}(a)\right)\left(\bar{\beta}(a)+\frac{1}{2}\widehat{\mathrm{H}}(a)\right)$
$\displaystyle+\frac{1}{2}\mathrm{T}_{2}(a)\left(\bar{\beta}(a)+\frac{1}{2}\widehat{\mathrm{H}}(a)\right)^{2}+\mathcal{O}(a^{4})\,.$
(15)

It can be checked that this expression coincides with that obtained in ref.
[16, Eq. (3.9)] ***The notations adopted here and in ref. [16] differ
slightly. To facilitate a comparison we note that the operators
$\mathrm{T}_{n}$ defined here satisfy the equation
$[S_{+},\mathrm{T}_{n}]=n[\mathrm{T}_{n-1},z_{1}+z_{2}]$..

The evolution kernel $\widehat{\mathrm{H}}(a)$ can be realized as an integral
operator. It acts on a function of two real variables as follows

$\displaystyle\widehat{\mathrm{H}}(a)\,f(z_{1},z_{2})$
$\displaystyle=\\!Af(z_{1},z_{2})\\!+\\!\\!\int_{+}h(\tau)f(z_{12}^{\alpha},z_{21}^{\beta}),$
(16)

where $A$ is a constant, $z_{12}^{\alpha}\equiv
z_{1}\bar{\alpha}+z_{2}\alpha$, $\bar{\alpha}\equiv 1-\alpha$, and

$\displaystyle\int_{+}\equiv\int_{0}^{1}d\alpha\int_{0}^{\bar{\alpha}}d\beta.$
(17)

$\tau=\alpha\beta/\bar{\alpha}\bar{\beta}$ is called conformal ratio. The
weight function $h(\tau)$ in Eq. (16) only depends on this particular
combination of the variables $\alpha,\beta$ as a consequence of invariance
properties of $\widehat{\mathrm{H}}$, Eq. (11).

It is easy to find that the operators $\mathrm{T}_{n}$ take the form

$\displaystyle\mathrm{T}_{n}(a)f(z_{1},z_{2})$
$\displaystyle=\\!\int_{+}\ln^{n}(1\\!-\\!\alpha\\!-\\!\beta)h(\tau)f(z_{12}^{\alpha},z_{21}^{\beta})\,,$
(18)

that again agrees with the results of ref. [16]. Note, that this expression
does not depend on the choice of the branch of the logarithm defining the
function $\mathrm{L}=\ln z_{12}$ in Eq. (8), see Appendix A for more
discussion.

##  III Anomalous dimensions vs kernels

First of all let us establish a connection between the eigenvalues of the
operators $\mathrm{H}$ and $\widehat{\mathrm{H}}$. Since both of them are
integral operators of the functional form in Eqs. (16), (18), both operators
are diagonalized by functions of the form
$\psi_{N}(z_{1},z_{2})=(z_{1}-z_{2})^{N-1}$, where $N$ is an arbitrary complex
number. One may worry that the continuation of the function $\psi_{N}$ for
negative $z_{12}$ is not unique and requires special care. But it does not
matter for our analysis. Indeed,
$z_{12}^{\alpha}-z_{21}^{\beta}=(1-\alpha-\beta)z_{12}$ with $\alpha+\beta<1$,
therefore the operators do not mix the regions $z_{12}\gtrless 0$. For
definiteness let us suppose that

$\displaystyle\psi_{N}(z_{1},z_{2})=\theta(z_{12})z_{12}^{N-1}.$ (19)

Let $\gamma(N),\widehat{\gamma}(N)$ be eigenvalues (anomalous dimensions) of
the operators $\mathrm{H}$, $\widehat{\mathrm{H}}$ corresponding to the
function $\psi_{N}$, respectively,

$\displaystyle\mathrm{H}(a)\psi_{N}$ $\displaystyle=\gamma(N)\psi_{N},$ (20)
$\displaystyle\widehat{\mathrm{H}}(a)\psi_{N}$
$\displaystyle=\widehat{\gamma}(N)\psi_{N}.$ (21)

The anomalous dimensions $\gamma(N),\widehat{\gamma}(N)$ are analytic
functions of $N$ in the right complex half-plane, $\mathrm{Re}(N)>0$. For
integer even (odd) $N$, $\gamma(N)$ gives the anomalous dimensions of the
local (axial)vector operators †††As usual one has to consider the operators of
certain parity,
$\mathcal{O}_{\pm}(z_{1},z_{2})=\mathcal{O}(z_{1},z_{2})\mp\mathcal{O}(z_{2},z_{1})$,
then the functions $\gamma_{\pm}(N)$ give the anomalous dimensions of local
operators, for even and odd $N$ respectively..

Now let us note that the operator $\mathrm{T}(\mathrm{H})$ acts on $\psi_{N}$
as follows

$\displaystyle\mathrm{T}(\mathrm{H})\psi_{N}(z_{1},z_{2})$
$\displaystyle=\sum_{n=0}^{\infty}\frac{\mathrm{L}^{n}}{n!}\left(\bar{\beta}(a)+\frac{1}{2}\gamma(N)\right)^{n}\\!\psi_{N}(z_{1},z_{2})$
$\displaystyle=z_{12}^{\bar{\beta}(a)+\frac{1}{2}\gamma(N)}\psi_{N}(z_{1},z_{2})$
$\displaystyle=\psi_{N+\bar{\beta}+\frac{1}{2}\gamma(N)}(z_{1},z_{2}).$ (22)

Thus, it follows from Eq. (II) that the anomalous dimensions $\gamma(N)$ and
$\widehat{\gamma}(N)$ satisfy the relation (cf. also Eq. (1))

$\displaystyle\gamma(N)=\widehat{\gamma}\left(N+\bar{\beta}(a)+\frac{1}{2}\gamma(N)\right).$
(23)

This relation appeared first in refs. [23, 22] as an generalization of the
Gribov-Lipatov reciprocity relation [26, 27]. It was shown that the asymptotic
expansion of the function $\widehat{\gamma}(N)$ for large $N$ is invariant
under the reflection $N\to-N-1$, see e.g., refs. [22, 28, 29, 30]. This
property strongly restricts harmonics sums which can appear in the
perturbative expansion of the anomalous dimension $\widehat{\gamma}(N)$ [29].
Explicit expressions for $\widehat{\gamma}(N)$ are known at four loops in QCD
[6] and at seven loops in the ${\cal N}=4$ SYM, see refs. [31, 29, 32, 33,
34].

###  III.1 Kernels from anomalous dimensions

For large $N$ the anomalous dimension $\widehat{\gamma}(N)$ grows as $\ln N$.
This term enters with a coefficient $2\Gamma_{\text{cusp}}(a)$ where
$\Gamma_{\text{cusp}}(a)$ is the so-called cusp anomalous dimension [35, 36]
known to the four-loop order in QCD [37, 38] and in ${\cal N}=4$ SYM [37].
Thus, we write $\widehat{\gamma}(N)$ in the following form

$\displaystyle\widehat{\gamma}(N)=2\Gamma_{\text{cusp}}(a)S_{1}(N)+A(a)+\Delta\widehat{\gamma}(N)\,,$
(24)

where $S_{1}(N)=\psi(N+1)-\psi(1)$ is the harmonic sum responsible for the
$\ln N$ behavior at large $N$, and $A(a)$ is a constant term. The remaining
term, $\Delta\widehat{\gamma}(N)$, vanishes at least as $O(1/N(N+1))$ at large
$N$. The constant $A(a)$ is exactly the same which appears in Eq. (16). The
first term in Eq. (24) comes from a special $\mathrm{SL}(2,\mathbb{R})$
invariant kernel

$\displaystyle\widehat{\mathcal{H}}f$
$\displaystyle=\int_{0}^{1}\frac{d\alpha}{\alpha}\Big{\\{}2f(z_{1},z_{2})-\bar{\alpha}\big{(}f(z_{12}^{\alpha},z_{2})+f(z_{1},z_{21}^{\alpha})\big{)}\Big{\\}},$
(25)

which in momentum space gives rise to the so-called plus-distribution. The
eigenvalues of this kernel are $2S_{1}(N)$
($\widehat{\mathcal{H}}z_{12}^{N-1}=2S_{1}(N)z_{12}^{N-1}$). It corresponds to
a singular contribution of the form $-\delta_{+}(\tau)$ to the invariant
kernel $h(\tau)$, see ref. [16, Eq. (2.19)] for detail. Thus the evolution
kernel can be generally written as

$\displaystyle\widehat{\mathrm{H}}$
$\displaystyle=\Gamma_{\text{cusp}}(a)\widehat{\mathcal{H}}+A(a)+\Delta\widehat{\mathrm{H}}\,.$
(26)

Here $\Delta\widehat{\mathrm{H}}$ is an integral operator,

$\displaystyle\Delta\widehat{\mathrm{H}}f(z_{1},z_{2})$
$\displaystyle=\int_{+}h(\tau)f(z_{12}^{\alpha},z_{21}^{\beta})\,,$ (27)

where the weight function $h(\tau)$ is a regular function of $\tau\in(0,1)$.
The eigenvalues of $\Delta\widehat{\mathrm{H}}$ are equal to
$\Delta\widehat{\gamma}(N)$ and are given by the following integral

$\displaystyle\Delta\widehat{\gamma}(N)=\int_{+}h(\tau)(1-\alpha-\beta)^{N-1}\,.$
(28)

The inverse transformation takes the form [14]

$\displaystyle h(\tau)$ $\displaystyle=\int_{C}\frac{dN}{2\pi
i}(2N+1)\Delta\widehat{\gamma}(N)P_{N}\left(\frac{1+\tau}{1-\tau}\right),$
(29)

where $P_{N}$ are the Legendre polynomials. The integration path $C$ goes
along the line parallel to the imaginary axis, $\text{Re}(N)>0$, such that all
poles of $\Delta\widehat{\gamma}(N)$ lie to the left of this line. Some
details of the derivation can be found in Appendix B.

One can hardly hope to evaluate the integral (29) in a closed form for an
arbitrary function $\Delta\widehat{\gamma}(N)$. However, as was mentioned
before, the anomalous dimensions $\Delta\widehat{\gamma}(N)$ in quantum field
theory are rather special functions. Most of the terms in the perturbative
expansion of $\Delta\widehat{\gamma}(N)$ have the following form

$\displaystyle\eta^{k}(N)\,\Omega_{\vec{m}}(N),$
$\displaystyle\eta^{k}(N)\,\Omega_{1}^{p}(N)$ (30)

where $\eta(N)=1/(N(N+1))$, and the functions
$\Omega_{\vec{m}}(N)=\Omega_{{m_{1},\ldots,m_{p}}}(N)$ are the parity
respecting harmonic sums [29],
($\Omega_{\vec{m}}(N)\sim\Omega_{\vec{m}}(-N-1)$ for $N\to\infty$). We will
assume that the sums $\Omega_{\vec{m}}(N)$ are “subtracted”, i.e.
$\Omega_{\vec{m}}(N)\to 0$ at $N\to\infty$. The second structure occurs only
for $k>0$, since $\Omega_{1}(N)=S_{1}(N)$ grows as $\ln N$ for large $N$.

Since all $\mathrm{SL}(2,\mathbb{R})$ invariant operators share the same
eigenfunctions, the product of two invariant operators $H_{1}$ and $H_{2}$,
$H_{1}H_{2}(=H_{2}H_{1})$ with eigenvalues $H_{1}(N)$ and $H_{2}(N)$
respectively, has eigenvalues $H_{1}(N)H_{2}(N)$. One can use this property to
reconstruct an operator with the eigenvalue (30).

First, we remark that the operator with the eigenvalues $\eta(N)$, (we denote
it as $\mathcal{H}_{+}$), has (as follows from Eq. (28)) a very simple weight
function, $h_{+}(\tau)=1$. This can also be derived from Eq. (29). Since
$P_{N}(x)=P_{-N-1}(x)$ the integral in Eq. (29) vanishes for the integration
path $\mathrm{Re}(N)=-1/2$ due to antisymmetry of the integrand. Therefore,
the integral (29) can be evaluated by the residue theorem ‡‡‡This trick allows
one to calculate the integral (29) for any function
$\Delta\widehat{\gamma}(N)$ with exact symmetry under $N\to-1-N$ reflection.

$\displaystyle
h_{+}(\tau)=\frac{2N+1}{N+1}P_{N}\left(\frac{1+\tau}{1-\tau}\right)\Big{|}_{N=0}=1.$
(31)

Let us consider the product $H_{2}=H_{+}\,H_{1}(=H_{1}\,H_{+})$, where $H_{1}$
is an integral operator with the weight function $h_{1}(\tau)$. Then the
weight function $h_{2}(\tau)$ of the operator $H_{2}$ is given by the
following integral

$\displaystyle
h_{2}(\tau)=\int_{0}^{\tau}\frac{ds}{\bar{s}^{2}}\ln(\tau/s)h_{1}(s),$ (32)

see Appendix B for details. Thus the contribution to the anomalous dimension
of type (30) can be evaluated with the help of this formula if the weight
function corresponding to the harmonic sums $\Omega_{\vec{m}}$ is known.

We also give an expression for another product of the operators:
$H_{2}=\widehat{\mathcal{H}}H_{1}$,

$\displaystyle h_{2}(\tau)$
$\displaystyle=-\ln\tau\,h_{1}(\tau)+2\bar{\tau}\int_{0}^{\tau}\frac{ds}{\bar{s}}\frac{h_{1}(\tau)-h_{1}(s)}{(\tau-s)}\,,$
(33)

which appears to be useful in the calculations as well.

###  III.2 Recurrence procedure

Let us consider the integral (29) with
$\Delta\widehat{\gamma}=\Omega_{\vec{m}}$,

$\displaystyle h_{\vec{m}}(\tau)$ $\displaystyle=\int_{C}\frac{dN}{2\pi
i}(2N+1)\Omega_{\vec{m}}(N)P_{N}\left(z\right),$ (34)

where $z=(1+\tau)/(1-\tau)$. Using a recurrence relation for the Legendre
functions

$\displaystyle(2N+1)P_{N}(z)$
$\displaystyle=\frac{d}{dz}\Big{(}P_{N+1}(z)-P_{N-1}(z)\Big{)}$ (35)

we obtain

$\displaystyle h_{\vec{m}}(\tau)$
$\displaystyle=-\frac{d}{dz}\int_{C}\frac{dN}{2\pi i}P_{N}(z)F_{\vec{m}}(N),$
(36)

where

$\displaystyle F_{\vec{m}}(N)$
$\displaystyle=\Big{(}\Omega_{\vec{m}}(N+1)-\Omega_{\vec{m}}(N-1)\Big{)}.$
(37)

It is easy to see that the function $F_{\vec{m}}(N)$ has the negative parity
under $N\to-N-1$ transformation and can be represented in the form

$\displaystyle F_{m_{1},\ldots,m_{p}}(N)$
$\displaystyle=\sum_{k=2}^{p}r_{k}(N)\Omega_{m_{k},\ldots,m_{p}}(N)+r(N),$
(38)

where $r_{k}(N)$ are rational functions of $N$. The harmonic sums
$\Omega_{m_{k},\ldots,m_{p}}(N)$ in Eq. (38) can be either of positive or
negative parity. Therefore the coefficient $r_{k}(N)$ accompanying the
positive parity function $\Omega_{m_{k},\ldots,m_{p}}(N)$ has the form
$r_{k}(N)=(2N+1)P_{k}(\eta)$, where $P_{k}$ is some polynomial, while
$r_{k}=P_{k}(\eta)$ for the harmonic sums of negative parity. The free term
has the form $r(N)=(2N+1)P(\eta)$. Together, they make
$F_{m_{1},\ldots,m_{p}}(N)$ with negative parity. For example, for the
harmonic sum $\Omega_{1,3}$ (see appendix C for a definition), one gets

$\displaystyle F_{1,3}(N)$
$\displaystyle=(2N+1)\underline{\eta\left(\Omega_{3}+\zeta_{3}-\eta^{2}-\frac{1}{2}\eta^{3}\right)},$
(39)

while for the harmonic sum $\Omega_{2,2}$

$\displaystyle F_{2,2}(N)$
$\displaystyle=(2N+1)\frac{1}{2}\eta^{3}(3+\eta)+\underline{\eta(2+\eta)\Omega_{2}}.$
(40)

Note the reappearance of the common factor $(2N+1)$ in the first case, (39).
This implies that, up to the derivative $d/dz$, the integral (36) has the form
(29). Hence, if the kernel corresponding to the underlined terms in Eq. (39)
is known, the kernel corresponding to $\Omega_{1,3}$ can be easily obtained.
Thus the problem of finding the invariant kernel with the eigenvalues
$\Omega_{1,3}(N)$ is reduced to the problem of finding the kernel with the
eigenvalues $\Omega_{3}(N)$ $(\Omega_{1,3}\mapsto\Omega_{3})$.

However, as it seen from our second example, not all parity preserving
harmonic sums share this property. Indeed, the underlined term on the right
hand side (rhs) of Eq. (40) does not have the factor $(2N+1)$. Hence, all
these transformations do not help to solve the problem for $\Omega_{2,2}$.

It is easy to see that the above recurrence procedure works only if all the
harmonic sums $\Omega_{m_{k},\ldots,m_{p}}$ appearing in Eq. (38) are of
positive parity. It was proven in ref. [29, Theorem 2] that any harmonic sum,
$\Omega_{\vec{m}}$, with all indices $\vec{m}$ positive odd or negative even
has positive parity (see Appendix C for explicit examples of the harmonic sums
satisfying these conditions). Therefore, the rhs of Eq. (38) only contains
harmonic sums of the same type. Thus the invariant kernels corresponding to
the harmonic sums of positive parity can always be calculated recursively,
using Eqs. (36), (38) and (32), (33). Crucially, only such harmonic sums
appear in the anomalous dimensions $\widehat{\gamma}(N)$ in QCD and ${\cal
N}=4$ SYM. All convolution integrals (32) and (33) can in turn be
systematically calculated with the packages HyperInt [39] or PolyLogTools
[40].

The explicit expressions for the kernels corresponding to the lowest harmonic
sums are given in Appendix C for references.

###  III.3 Invariant kernels: QCD

Below we give an explicit expression for the invariant kernel of twist-two
flavor nonsinglet operator in QCD. We will not split the operator
$\mathcal{O}(z_{1},z_{2})$ into positive (negative) parity operators. The
evolution operator still takes the form (26), with
$\Delta\widehat{\mathrm{H}}$ given by the following integral

$\displaystyle\Delta\widehat{\mathrm{H}}f(z_{1},z_{2})$
$\displaystyle=\int_{+}(h(\tau)+\bar{h}(\tau)P_{12})f(z_{12}^{\alpha},z_{21}^{\beta})\,,$
(41)

where $P_{12}$ is a permutation operator,
$P_{12}f(z_{1},z_{2})=f(z_{2},z_{1})$§§§ In order to avoid possible
misunderstandings we write down it explicitly,
$P_{12}f(z_{12}^{\alpha},z_{21}^{\beta})=f(z_{21}^{\alpha},z_{12}^{\beta})$..
For (anti)symmetric functions $f(z_{1},z_{2})$ the operator (41) takes a
simpler form (27) with the kernel $h\pm\bar{h}$.

Our expression for the constant term $A(a)$ agrees with the constant term
$\chi$ given in ref. [16, Eq. (5.5)], $A=\chi-2\Gamma_{\text{cusp}}$. For
completeness, we provide explicit expressions for the constant
$A=aA_{1}+a^{2}A_{2}+a^{3}A_{3}+\cdots$,

$\displaystyle A_{1}$ $\displaystyle=-6C_{F}\,,$ $\displaystyle A_{2}$
$\displaystyle=C_{F}\left[n_{f}\left(\frac{16}{3}\zeta_{2}+\frac{2}{3}\right)-N_{c}\left(\frac{52}{3}\zeta_{2}+\frac{43}{6}\right)+\frac{1}{N_{c}}\left(24\zeta_{3}-12\zeta_{2}+\frac{3}{2}\right)\right]\,,$
$\displaystyle A_{3}$
$\displaystyle=C_{F}\bigg{[}n_{f}^{2}\left(\frac{32}{9}\zeta_{3}-\frac{160}{27}\zeta_{2}+\frac{34}{9}\right)+n_{f}N_{c}\left(-\frac{256}{15}\zeta_{2}^{2}+\frac{8}{9}\zeta_{3}+\frac{2492}{27}\zeta_{2}-17\right)+\frac{n_{f}}{N_{c}}\left(\frac{232}{15}\zeta_{2}^{2}-\frac{136}{3}\zeta_{3}+\frac{20}{3}\zeta_{2}-23\right)$
$\displaystyle\qquad+N_{c}^{2}\left(-80\zeta_{5}+\frac{616}{15}\zeta_{2}^{2}+\frac{266}{9}\zeta_{3}-\frac{5545}{27}\zeta_{2}+\frac{847}{18}\right)+\left(-120\zeta_{5}-16\zeta_{2}\zeta_{3}-\frac{124}{15}\zeta_{2}^{2}+\frac{1048}{3}\zeta_{3}-\frac{356}{3}\zeta_{2}+\frac{209}{4}\right)$
$\displaystyle\qquad+\frac{1}{N_{c}^{2}}\left(120\zeta_{5}+16\zeta_{2}\zeta_{3}-\frac{144}{5}\zeta_{2}^{2}-34\zeta_{3}-9\zeta_{2}-\frac{29}{4}\right)\bigg{]}\,,$
(42)

where $C_{F}=(N_{c}^{2}-1)/(2N_{c})$ is the quadratic Casimir in the
fundamental representation of $SU(N_{c})$ and we take $T_{F}=1/2$. Note that
we are adopting a different color basis compared to ref. [16].

The explicit expressions for the cusp anomalous dimensions $\Gamma_{\rm
cusp}(a)=a\Gamma_{\rm cusp}^{(1)}+a^{2}\Gamma_{\rm
cusp}^{(2)}+a^{3}\Gamma_{\rm cusp}^{(3)}$ up to three loops are provided in
Eq. (D). Finally we give answers for the kernels
$h(\bar{h})(a)=\sum_{k}a^{k}h_{k}(\bar{h}_{k})$. Explicit one- and two-loop
expressions are known [14, 16] but for completeness we give them here

$\displaystyle h_{1}=-4C_{F}\,,$ $\displaystyle\bar{h}_{1}=0\,,$ (43)

and

$\displaystyle h_{2}$
$\displaystyle=C_{F}\biggl{\\{}n_{f}\frac{88}{9}+N_{c}\left(-2\mathrm{H}_{1}+8\zeta_{2}-\frac{604}{9}\right)$
$\displaystyle\quad+\frac{1}{N_{c}}\left(-8\Big{(}\mathrm{H}_{11}+\mathrm{H}_{2}\Big{)}+2\left(1-\frac{4}{\tau}\right)\mathrm{H}_{1}\right)\biggr{\\}}\,,$
$\displaystyle\bar{h}_{2}$
$\displaystyle=-\frac{8C_{F}}{N_{c}}\left(\mathrm{H}_{11}+\tau\,\mathrm{H}_{1}\right)\,,$
(44)

where $\mathrm{H}_{\vec{m}}=\mathrm{H}_{\vec{m}}(\tau)$ are the harmonic
polylogarithms (HPLs) [24]. The three-loop expression¶¶¶A file with our main
results can be obtained from the preprint server http://arXiv.org by
downloading the source. Furthermore, they are available from the authors upon
request. is more involved

$\displaystyle h_{3}$
$\displaystyle=C_{F}\biggl{\\{}-\frac{64}{9}n_{f}^{2}+n_{f}N_{c}\frac{8}{3}\biggl{[}\mathrm{H}_{3}-\mathrm{H}_{110}-\mathrm{H}_{20}+\mathrm{H}_{12}+\frac{1}{\tau}\mathrm{H}_{2}-\frac{1}{\tau}\mathrm{H}_{10}-\frac{19}{12}\mathrm{H}_{1}+8\zeta_{3}-\frac{32}{3}\zeta_{2}+\frac{5695}{72}\biggr{]}$
$\displaystyle\quad+\frac{n_{f}}{N_{c}}\frac{16}{3}\biggr{[}3\zeta_{3}-\frac{75}{16}+\mathrm{H}_{3}+\mathrm{H}_{21}+\mathrm{H}_{12}+\mathrm{H}_{111}+\left(\frac{16}{3}+\frac{1}{\tau}\right)\Big{(}\mathrm{H}_{2}+\mathrm{H}_{11}\Big{)}+\left(\frac{31}{24}+\frac{10}{3\tau}\right)\mathrm{H}_{1}\biggr{]}$
$\displaystyle\quad+N_{c}^{2}4\biggl{[}\mathrm{H}_{13}+\mathrm{H}_{112}-\mathrm{H}_{120}-\mathrm{H}_{1110}+2\mathrm{H}_{4}-2\mathrm{H}_{30}-2\mathrm{H}_{210}+2\mathrm{H}_{22}+\left(\frac{8}{3}-\frac{2}{\tau}\right)\Big{(}\mathrm{H}_{20}-\mathrm{H}_{3}+\mathrm{H}_{110}-\mathrm{H}_{12}\Big{)}$
$\displaystyle\quad-\frac{5}{4}\Big{(}\mathrm{H}_{10}+\mathrm{H}_{11}\Big{)}+\frac{2}{3\tau}\Big{(}\mathrm{H}_{10}-\mathrm{H}_{2}\Big{)}-\frac{5}{2}\mathrm{H}_{0}+\left(\frac{115}{72}+\zeta_{2}+\frac{1}{\tau}\right)\mathrm{H}_{1}-\frac{44}{5}\zeta_{2}^{2}-\frac{22}{3}\zeta_{3}+\frac{436}{9}\zeta_{2}-\frac{4783}{27}\biggr{]}$
$\displaystyle\quad+16\,\biggl{[}\mathrm{H}_{4}-\mathrm{H}_{30}+\mathrm{H}_{13}+\mathrm{H}_{121}-\frac{3}{2}\mathrm{H}_{120}+\frac{3}{2}\mathrm{H}_{22}+\frac{3}{2}\mathrm{H}_{112}+2\mathrm{H}_{31}+2\mathrm{H}_{1111}+3\mathrm{H}_{211}-\frac{1}{2}\mathrm{H}_{1110}-\left(\frac{1}{\tau}+1\right)\mathrm{H}_{20}$
$\displaystyle\quad-\left(\frac{11}{6}-\frac{1}{\tau}\right)\mathrm{H}_{3}-\mathrm{H}_{110}+\left(-\frac{37}{12}+\frac{3}{2\tau}\right)\mathrm{H}_{12}-\left(\frac{7}{3}-\frac{2}{\tau}\right)\mathrm{H}_{21}+\left(-\frac{43}{12}+\frac{3}{\tau}\right)\mathrm{H}_{111}+\Big{(}\frac{13}{8}+\frac{1}{2}\zeta_{2}\Big{)}\mathrm{H}_{10}$
$\displaystyle\quad-\left(\frac{1}{2}\zeta_{2}+\frac{127}{9}+\frac{11}{6\tau}\right)\mathrm{H}_{2}-\left(\frac{899}{72}+\frac{1}{3\tau}\right)\mathrm{H}_{11}+\Big{(}\zeta_{2}-1\Big{)}\mathrm{H}_{0}+\left(\frac{7}{4}\zeta_{2}-\frac{143}{36}-\frac{1}{\tau}\left(\frac{1}{2}\zeta_{2}+\frac{67}{9}\right)\right)\mathrm{H}_{1}+\frac{5}{2}\zeta_{2}-\frac{47}{24}\biggr{]}$
$\displaystyle\quad+\frac{8}{N_{c}^{2}}\biggl{[}\mathrm{H}_{4}-\mathrm{H}_{30}-\mathrm{H}_{210}+\mathrm{H}_{112}-\mathrm{H}_{1111}-2\,\mathrm{H}_{120}+2\,\mathrm{H}_{13}+2\mathrm{H}_{31}-2\mathrm{H}_{1110}-2\mathrm{H}_{211}+3\mathrm{H}_{121}$
$\displaystyle\quad-\left(\frac{1}{2}+\frac{1}{\tau}\right)\Big{(}\mathrm{H}_{20}-\mathrm{H}_{3}+\mathrm{H}_{110}\Big{)}+2\left(1+\frac{1}{\tau}\right)\,\mathrm{H}_{21}+\left(\frac{3}{2}-\frac{2}{\tau}\right)\,\mathrm{H}_{111}+\left(\frac{7}{8}+\frac{3}{2\tau}\right)\,\mathrm{H}_{10}-\left(\zeta_{2}-\frac{1}{2}+\frac{3}{2\tau}\right)\,\mathrm{H}_{2}$
$\displaystyle\quad+\Big{(}\frac{11}{8}-\zeta_{2}\Big{)}\,\mathrm{H}_{11}-\frac{11}{4}\,\mathrm{H}_{0}+\left(\zeta_{2}-\frac{107}{16}-\frac{\zeta_{2}}{\tau}-\frac{1}{2\tau}\right)\,\mathrm{H}_{1}+\frac{7}{2}\biggr{]}\biggr{\\}}$
(45) and $\displaystyle\bar{h}_{3}$
$\displaystyle=-8C_{F}\Biggl{\\{}-\frac{2n_{f}}{3N_{c}}\left[\mathrm{H}_{111}+\mathrm{H}_{110}+\tau\,\mathrm{H}_{10}+\left(\frac{16}{3}+\tau\right)\mathrm{H}_{11}+\left(\frac{1}{2}+\frac{10}{3}\tau\right)\mathrm{H}_{1}+\frac{1}{2}\right]$
$\displaystyle\quad+\mathrm{H}_{120}+\mathrm{H}_{22}-\mathrm{H}_{1110}-\mathrm{H}_{112}-2\mathrm{H}_{121}+2\mathrm{H}_{211}-4\mathrm{H}_{1111}+\tau\mathrm{H}_{20}+\left(\frac{13}{6}-\tau\right)\mathrm{H}_{110}+\Big{(}\frac{1}{2}-2\tau\Big{)}\mathrm{H}_{12}$
$\displaystyle\quad+\Big{(}\frac{5}{2}-2\tau\Big{)}\mathrm{H}_{21}+\left(\frac{43}{6}-6\tau\right)\mathrm{H}_{111}-\left(\zeta_{2}-\frac{13}{6}\tau\right)\mathrm{H}_{10}-\Big{(}3+\zeta_{2}+\frac{3}{2}\tau\Big{)}\mathrm{H}_{2}+\left(\frac{236}{9}+\frac{2}{3}\tau\right)\mathrm{H}_{11}-\zeta_{2}\tau\mathrm{H}_{0}$
$\displaystyle\quad+\left(\frac{53}{6}+\zeta_{2}+3\zeta_{3}+\frac{134}{9}\tau+\zeta_{2}\tau\right)\mathrm{H}_{1}+\frac{11}{6}+3\left(\zeta_{2}+\zeta_{3}-\frac{1}{2}\right)\tau$
$\displaystyle\quad+\frac{1}{N_{c}^{2}}\biggl{[}\mathrm{H}_{1111}-\mathrm{H}_{22}-\mathrm{H}_{211}+3\mathrm{H}_{120}+3\mathrm{H}_{112}-3\mathrm{H}_{1110}+4\mathrm{H}_{121}+3\,\tau\,\mathrm{H}_{20}+3\,\left(\frac{1}{2}-\tau\right)\,\mathrm{H}_{110}-\left(\frac{7}{2}-4\tau\right)\,\mathrm{H}_{12}$
$\displaystyle\quad+\left(\frac{1}{2}+4\tau\right)\,\mathrm{H}_{21}+\left(-\frac{3}{2}+2\tau\right)\,\mathrm{H}_{111}-3\,\left(\zeta_{2}-\frac{1}{2}\tau\right)\,\mathrm{H}_{10}+\left(\zeta_{2}-2-\frac{3}{2}\tau\right)\,\mathrm{H}_{2}+2\Big{(}\zeta_{2}-1\Big{)}\,\mathrm{H}_{11}-3\zeta_{2}\,\tau\,\mathrm{H}_{0}$
$\displaystyle\quad+\Big{(}5+2\zeta_{2}+3\zeta_{3}+\zeta_{2}\tau\Big{)}\,\mathrm{H}_{1}+3\tau\left(\zeta_{3}-\frac{1}{2}\right)\biggr{]}\Biggr{\\}}.$
(46)

The kernels are smooth functions of $\tau$ except for the endpoints $\tau=0$
and $\tau=1$. For $\tau\to 1$ the three-loop kernel functions behave as
$\sum_{0\leq k\leq 4}\sum_{m>0}r_{km}\bar{\tau}^{m}\ln^{k}\bar{\tau}$. For
small $\tau$ – which determines the large $N$ asymptotic of the anomalous
dimensions – the kernels (for each color structure) have the form $\sum_{k\geq
0}(a_{k}+b_{k}\ln\tau)\tau^{k}$. We note here that the reciprocity property of
the anomalous dimension is equivalent to the statement that the small $\tau$
expansion of the kernels does not involve non-integer powers of $\tau$, namely
$h(\tau)\sim\sum_{m,k\geq 0}a_{mk}\tau^{m}\ln^{k}\tau$.

Below we compare our exact three-loop results with the approximate expressions
constructed in ref. [16]. The approximate expressions reproduce the asymptotic
behaviors of the exact kernels at both $\tau\to 0,1$. We therefore subtract
the logarithmically divergent pieces (see Eqs. (D) and (D) for explicit
expressions) from both the exact and the approximated expressions to highlight
their (small) deviations as shown in Figs. 1 and 2. For illustrative purposes,
we plot the planar contribution ($C_{F}N_{c}^{2}$ and $C_{F}$ in $h_{3}$ and
$\bar{h}_{3}$ respectively) and the subsubplanar contribution
($C_{F}/N_{c}^{2}$). The former is numerically dominant and generates the
leading contribution in the large-$N_{c}$ limit whereas the latter shows the
worst-case scenario for the previous approximation using a simple HPL function
ansatz. The error of other color structures all fall between the planar and
subsubplanar cases, hence are numerically small.

Figure 1: Comparison of two distinct color contributions ($C_{F}N_{c}^{2}$ and
$C_{F}/N_{c}^{2}$) in the exact (black solid) and approximated (red dashed)
three-loop kernel $h_{3}$ (see ref. [16] for explicit expressions of the
latter). The inset curves show the relative percentage errors ($(1-h_{3,{\rm
appr}}^{(c)}/h_{3,{\rm exact}}^{(c)})\times 100\%$) of the approximation.
Figure 2: Same as Fig. 1 for the color structures $C_{F}$ and
$C_{F}/N_{c}^{2}$ in $\bar{h}_{3}$.

###  III.4 Invariant kernels: ${\cal N}=4$ SYM

In this section we present the invariant kernels for the universal anomalous
dimensions of the planar ${\cal N}=4$ SYM, see e.g., refs. [41, 31] for
expressions up to NNLO. They are rather short so that we quote them here. We
use the parametrization (24), where $\Gamma_{\text{cusp}}(a)$ can be found in
ref. [37] and the constant term $A(a)$ is

$\displaystyle
A(a)=-24a^{2}\zeta_{3}+32a^{3}\big{(}\zeta_{2}\zeta_{3}+5\zeta_{5}\big{)}+O(a^{4}),$
(47)

where $a=\frac{N_{c}g^{2}_{\text{SYM}}}{16\pi^{2}}$, and

$\displaystyle\Delta\widehat{\gamma}(N)$
$\displaystyle=-a^{2}16\Big{(}\Omega_{3}-2\,\Omega_{-2,1}+2\,\Omega_{1}\,\Omega_{-2}\Big{)}$
$\displaystyle\quad+a^{3}{64}\Big{(}\Omega_{5}+2\,\Omega_{3,-2}-8\,\Omega_{1,1,-2,1}+2\Omega_{1,-4}$
$\displaystyle\quad+\Omega_{1}(\Omega_{-4}+\Omega_{-2}^{2}+\zeta_{2}\Omega_{-2})\\!-2\zeta_{2}\Omega_{-2,1}\Big{)}.$
(48)

For the kernels we find $h_{1}=\bar{h}_{1}=0$,

$\displaystyle h_{2}$
$\displaystyle=8\frac{\bar{\tau}}{\tau}\mathrm{H}_{1}\,,$
$\displaystyle\bar{h}_{2}$ $\displaystyle=-8{\bar{\tau}}\mathrm{H}_{1}$ (49)

and

$\displaystyle h_{3}$
$\displaystyle=-16\frac{\bar{\tau}}{\tau}\Big{(}4\,\mathrm{H}_{111}+\mathrm{H}_{21}+\mathrm{H}_{12}+\mathrm{H}_{110}\Big{)}\,,$
(50) $\displaystyle\bar{h}_{3}$
$\displaystyle=16\bar{\tau}\Big{(}4\mathrm{H}_{111}+3\big{(}\mathrm{H}_{21}+\mathrm{H}_{12}\big{)}-\mathrm{H}_{110}+\mathrm{H}_{20}-\zeta_{2}\mathrm{H}_{0}\Big{)}.$

These expressions are extremely simple in comparison with the expressions in
QCD of the same order. Let us notice that the two-loop kernels contain only
HPLs of weight one with the three-loop kernels involving HPLs of weight three,
while in QCD the corresponding kernels require HPLs of weight two and four,
respectively. Note also that the kernel $h$ is proportional to the factor
$\bar{\tau}/\tau$ and the kernel $\bar{h}$ to the factor $\bar{\tau}$. It
would be interesting to see if these properties persist in higher loops.

##  IV Summary

We have constructed a transformation that brings the evolution kernels of
twist-two operators to the canonically conformal invariant form. The
eigenvalues of these kernels are given by the parity respecting anomalous
dimensions. We have developed a recurrence procedure that allows one to
restore the weight functions of the corresponding kernels. It is applicable to
a subset of the harmonic sums (with positive odd and negative even indices).
It is interesting to note that exactly only such harmonic sums appear in the
expressions for the reciprocity respecting anomalous dimensions.

We have calculated the three-loop invariant kernels in QCD and in ${\cal N}=4$
SYM (in the planar limit). In QCD it was the last missing piece to obtain the
three-loop evolution kernels for the flavor-nonsinglet twist-two operators in
a fully analytic form, see ref. [16].

In the case of ${\cal N}=4$ SYM the lowest order expressions for the kernels
are rather simple and exhibit some regularities, $h\sim\bar{\tau}/\tau$,
$\bar{h}\sim\bar{\tau}$. It would be interesting to check if these properties
survive at higher loops. We expect that at $\ell$-loops the kernels
$h^{(\ell)}(\tau)$ will be given by linear combinations (up to common
prefactors) of HPLs of weight $2\ell-3$ with positive indices. Therefore going
over to the invariant kernel can lead to a more compact representation of the
anomalous dimensions than representing the anomalous dimension spectrum
$\gamma(N)$ in terms of harmonic sums. The much smaller function basis in
terms of HPLs ($\tau/\bar{\tau}H_{\vec{m}}$ and $\bar{\tau}H_{\vec{m}}$) opens
the possibility of extracting the analytical expressions of the higher-order
evolution kernels from minimal numerical input through the PSLQ algorithm.

## Acknowledgments

We are grateful to Vladimir M. Braun and Gregory P. Korchemsky for
illuminating discussions and comments on the manuscript. This work is
supported by Deutsche Forschungsgemeinschaft (DFG) through the Collaborative
Research Center TRR110/2, grant 409651613 (Y.J.), and the Research Unit FOR
2926, project number 40824754 (S.M.).

## Appendices

## Appendix A

In this appendix, we describe in detail the derivations of some of the
equations presented in section II. Let us start with Eq. (9). For the
generator $\mathrm{S}_{-}(a)=S_{-}$ the statement is trivial. Next, making use
of Eq. (8) for the operator $\mathrm{T}(\mathrm{H})$ and, taking into account
that $\mathrm{H}(a)$ commutes with the generators $\mathrm{S}_{\alpha}(a)$,
one can write the left hand side (lhs) of Eq. (9) in the form

$\displaystyle\sum_{n=0}^{\infty}\frac{1}{n!}\mathrm{L}^{n}\mathrm{S}_{\alpha}(a)\mathrm{X}^{n}\,,$
(A.51)

where $\mathrm{X}=\bar{\beta}(a)+\frac{1}{2}\mathrm{H}(a)$. Using the
representation (7) for the generators and taking into account that
$[S_{0},\mathrm{L}]=1$ and $[S_{+},\mathrm{L}]=z_{1}+z_{2}$ (we recall that
$\mathrm{L}=\ln z_{12}$) one obtains

$\displaystyle\mathrm{L}^{n}\mathrm{S}_{0}$
$\displaystyle=S_{0}\mathrm{L}^{n}-n\mathrm{L}^{n-1}+\mathrm{L}^{n}X,$
$\displaystyle\mathrm{L}^{n}\mathrm{S}_{+}$
$\displaystyle=S_{+}\mathrm{L}^{n}+(z_{1}+z_{2})\left(-n\mathrm{L}^{n-1}+\mathrm{L}^{n}X\right).$
(A.52)

Substituting these expressions back into Eq. (A.51) one finds that the
contributions of the last two terms on the rhs of Eq. (A) cancel each other.
Hence Eq. (A.51) takes the form

$\displaystyle
S_{\alpha}\sum_{n=0}^{\infty}\frac{1}{n!}\mathrm{L}^{n}\mathrm{X}^{n}=S_{\alpha}\mathrm{T}(\mathrm{H}),$
(A.53)

that finally results in Eq. (9).

Let us now show that the inverse to $\mathrm{T}(\mathrm{H})$ has the form
(12). The product
$\mathcal{I}=\mathrm{T}^{-1}(\mathrm{H})\mathrm{T}(\mathrm{H})$ can be written
as

$\displaystyle\mathcal{I}$
$\displaystyle=\sum_{n=0}^{\infty}\frac{(-1)^{n}}{n!}\mathrm{L}^{n}\left(\bar{\beta}(a)+\frac{1}{2}\widehat{\mathrm{H}}(a)\right)^{n}\mathrm{T}(\mathrm{H}).$
(A.54)

Moving $\mathrm{T}(\mathrm{H})$ to the left with help of the relation (10) and
then using Eq. (8) for $\mathrm{T}(\mathrm{H})$ one gets
($\mathrm{X}=\bar{\beta}(a)+\frac{1}{2}\mathrm{H}(a)$)

$\displaystyle\mathcal{I}=\sum_{n=0}^{\infty}\frac{(-1)^{n}}{n!}\mathrm{L}^{n}\mathrm{T}(\mathrm{H})\mathrm{X}^{n}=\sum_{n,k=0}^{\infty}\frac{(-1)^{n}}{n!k!}\mathrm{L}^{n+k}\mathrm{X}^{n+k}=1\,.$

Finally, we consider the product of operators $\mathrm{T}$ with a differently
defined function $\mathrm{L}$. Namely, let us take
$\mathrm{T}_{\pm}(\mathrm{H})\equiv\mathrm{T}(\mathrm{L}_{\pm},\mathrm{H})$,
where $L_{\pm}=\ln(z_{12}\pm i0)$ so that
$\mathrm{L}_{+}-\mathrm{L}_{-}=2\pi\theta(z_{2}-z_{1})$. In order to calculate
the product $\mathrm{U}=\mathrm{T}_{+}(H)\mathrm{T}_{-}(H)$ one proceeds as
before: use expansion (8) for $\mathrm{T}_{+}(H)$, move $\mathrm{T}_{-}(H)$ to
the left and then expand it into a power series. It yields

$\displaystyle\mathrm{U}$
$\displaystyle=\sum_{n,k=0}^{\infty}\frac{(-1)^{n}}{n!k!}\mathrm{L}_{+}^{n}\mathrm{L}_{-}^{k}\widehat{\mathrm{X}}^{n+k}\,,$
(A.55)

where
$\widehat{\mathrm{X}}=\bar{\beta}(a)+\frac{1}{2}\widehat{\mathrm{H}}(a)$. Let
$L_{+}=L_{-}+2\pi i\theta(z_{2}-z_{1})$ one can get for the sum in Eq. (A.55)

$\displaystyle\mathrm{U}$ $\displaystyle=\sum_{m=0}^{\infty}\frac{(2\pi
i\theta)^{m}}{m!}\widehat{\mathrm{X}}^{m}(a)=1-\theta\left(1-e^{2\pi
i\left(\bar{\beta}+\frac{1}{2}\widehat{\mathrm{H}}\right)}\right),$

where $\theta\equiv\theta(z_{2}-z_{1})$. Since $S_{0,+}\theta(z_{21})\sim
z_{21}\delta(z_{21})=0$ one concludes that $\mathrm{U}$ commutes with the
canonical generators $S_{\alpha}$ and hence
$\mathrm{U}\widehat{\mathrm{H}}=\widehat{\mathrm{H}}\mathrm{U}$.

## Appendix B

Let us check that the kernel $h(\tau)$ given by Eq. (29) has the eigenvalues
$\Delta\widehat{\gamma}(N)$. First, after some algebra, the integral in Eq.
(28) can be brought to the following form

$\displaystyle\Delta\widehat{\gamma}(N)=\int_{1}^{\infty}dt\,h\left(\frac{t-1}{t+1}\right)\,Q_{N}(t)\,,$
(B.56)

where $Q_{N}(t)$ is the Legendre function of the second kind [42]. Inserting
$h$ in the form of Eq. (29) into Eq. (B.56) one gets

$\displaystyle\int_{C}\frac{dN^{\prime}}{2\pi
i}(2N^{\prime}+1)\Delta\gamma(N^{\prime})\int_{1}^{\infty}dt\,P_{N^{\prime}}(t)\,Q_{N}(t)\,.$
(B.57)

The $t$-integral of the product of the two Legendre functions gives [42]

$\displaystyle\left((N-N^{\prime})(N+N^{\prime}+1)\right)^{-1}.$ (B.58)

Then closing the integration contour in the right half-plane one evaluates the
$N^{\prime}$ integral with the residue theorem at $N^{\prime}=N$ yielding the
desired lhs of Eq. (B.56).

Finally, in order to verify Eq. (32) one can check that the integral (B.56)
with the kernel $h_{2}$, $\Delta\widehat{\gamma}_{2}(N)$, is equal to
$\Delta\widehat{\gamma}_{1}(N)/N/(N+1)$. The simplest way to do it is to
substitute the Legendre function in the form

$\displaystyle Q_{N}(t)=-\partial_{t}(1-t^{2})\partial_{t}Q_{N}(t)/N/(N+1),$
(B.59)

and perform integration by parts.

## Appendix C

In this appendix, we collect the harmonic sums and the corresponding kernels
which we have used. We split them into two parts: the first one includes the
harmonic sums $\Omega_{m_{1},\ldots,m_{k}}$ such that
$\prod_{i}^{k}\mathrm{sign}(m_{i})=1$.

$\displaystyle\Omega_{3}=S_{3}-\zeta_{3},$
$\displaystyle\Omega_{3,1}=S_{3,1}-\frac{1}{2}S_{4}-\frac{3}{10}\zeta_{2}^{2}$
$\displaystyle\Omega_{-2,-2}=S_{-2,-2}-\frac{1}{2}S_{4}+\frac{1}{2}\zeta_{2}S_{-2}+\frac{1}{8}\zeta_{2}^{2},$
$\displaystyle\Omega_{1,3,1}=S_{1,3,1}-\frac{1}{2}S_{1,4}-\frac{1}{2}S_{4,1}+\frac{1}{4}S_{5}-\frac{3}{10}\zeta_{2}^{2}S_{1}+\frac{3}{4}\zeta_{5},$
$\displaystyle\Omega_{-2,-2,1}=S_{-2,-2,1}-\frac{1}{2}S_{4,1}-\frac{1}{2}S_{-2,-3}+\frac{1}{4}\zeta_{3}S_{-2}+\frac{5}{16}\zeta_{5},$
$\displaystyle\Omega_{5}=S_{5}-\zeta_{5}.$ (C.60)

Here $S_{\vec{m}}$ are the harmonic sums with argument $N$. We define the sums
of negative signature, $\prod_{i}^{k}\mathrm{sign}(m_{i})=-1$, with an
additional sign factor:

$\displaystyle\Omega_{-2}=(-1)^{N}\left[S_{-2}+\frac{\zeta_{2}}{2}\right],$
$\displaystyle 0$
$\displaystyle\Omega_{-2,1}=(-1)^{N}\left[S_{-2,1}-\frac{1}{2}S_{-3}+\frac{1}{4}\zeta_{3}\right],$
$\displaystyle\Omega_{1,-2,1}=(-1)^{N}\left[S_{1,-2,1}-\frac{1}{2}S_{1,-3}-\frac{1}{2}S_{-3,1}+\frac{1}{4}S_{-4}\right.$
$\displaystyle\quad\left.+\frac{1}{4}\zeta_{3}S_{1}-\frac{1}{80}\zeta_{2}^{2}\right],$
$\displaystyle\Omega_{-4,1}=(-1)^{N}\left[S_{-4,1}-\frac{1}{2}S_{-5}+\frac{11}{8}\zeta_{5}-\frac{1}{2}\zeta_{2}\zeta_{3}\right],$
$\displaystyle\Omega_{3,-2}=(-1)^{N}\left[S_{3,-2}-\frac{1}{2}S_{-5}+\frac{1}{2}\zeta_{2}S_{3}+\frac{9}{8}\zeta_{5}-\frac{3}{4}\zeta_{2}\zeta_{3}\right],$
$\displaystyle\Omega_{1,1,-2,1}=(-1)^{N}\biggl{[}S_{1,1,-2,1}-\frac{1}{2}S_{1,1,-3}-\frac{1}{2}S_{1,-3,1}$
$\displaystyle\qquad\quad-\frac{1}{2}S_{2,-2,1}+\frac{1}{4}S_{2,-3}+\frac{1}{4}S_{-4,1}+\frac{1}{4}S_{1,-4}-\frac{1}{8}S_{-5}$
$\displaystyle\qquad\quad+\frac{1}{4}\zeta_{3}S_{1,1}-\frac{1}{80}\zeta_{2}^{2}S_{1}-\frac{1}{8}\zeta_{3}S_{2}+\frac{1}{8}\zeta_{5}-\frac{1}{16}\zeta_{2}\zeta_{3}\biggr{]},$
$\displaystyle\Omega_{1,-4}=(-1)^{N}\biggl{[}S_{1,-4}-\frac{1}{2}S_{-5}+\frac{7}{20}\zeta_{2}^{2}S_{1}-\frac{11}{8}\zeta_{5}+\frac{1}{2}\zeta_{2}\zeta_{3}\biggr{]}.$
(C.61)

These combinations of harmonic sums are generated by the following kernels,

$\displaystyle{\mathcal{H}}_{3}=-\frac{1}{2}\frac{\bar{\tau}}{\tau}\mathrm{H}_{1},$
$\displaystyle{\mathcal{H}}_{3,1}=\frac{1}{4}\frac{\bar{\tau}}{\tau}\left(\mathrm{H}_{11}+\mathrm{H}_{10}\right)$
$\displaystyle{\mathcal{H}}_{-2,-2}=\frac{1}{4}\frac{\bar{\tau}}{\tau}\mathrm{H}_{11},$
$\displaystyle{\mathcal{H}}_{1,3,1}=-\frac{1}{8}\frac{\bar{\tau}}{\tau}\left(\mathrm{H}_{20}+\mathrm{H}_{110}+\mathrm{H}_{21}+\mathrm{H}_{111}\right),$
$\displaystyle{\mathcal{H}}_{-2,-2,1}=\frac{1}{8}\frac{\bar{\tau}}{\tau}\left(\mathrm{H}_{12}-\mathrm{H}_{110}\right),$
$\displaystyle\mathcal{H}_{5}=-\frac{1}{2}\frac{\bar{\tau}}{\tau}\left(\mathrm{H}_{111}+\mathrm{H}_{12}\right)$
(C.62)

and

$\displaystyle{\mathcal{H}}_{-2}=\frac{1}{2}\bar{\tau},$
$\displaystyle{\mathcal{H}}_{-2,1}=-\frac{1}{4}\bar{\tau}(\mathrm{H}_{1}+\mathrm{H}_{0}),$
$\displaystyle{\mathcal{H}}_{1,-2,1}=\frac{1}{8}\bar{\tau}\left(\mathrm{H}_{10}+\mathrm{H}_{11}\right),$
$\displaystyle{\mathcal{H}}_{-4,1}=-\frac{1}{4}\bar{\tau}\left(\mathrm{H}_{21}+\mathrm{H}_{20}+\mathrm{H}_{111}+\mathrm{H}_{110}\right),$
$\displaystyle{\mathcal{H}}_{3,-2}=-\frac{1}{4}\bar{\tau}\left(\mathrm{H}_{21}+\mathrm{H}_{111}\right),$
$\displaystyle{\mathcal{H}}_{1,1,-2,1}=-\frac{1}{16}\bar{\tau}\,\left(\mathrm{H}_{111}+\mathrm{H}_{110}\right),$
$\displaystyle{\mathcal{H}}_{1,-4}=-\frac{1}{4}\bar{\tau}\left(\mathrm{H}_{12}+\mathrm{H}_{111}\right),$
(C.63)

where all HPLs have argument $\tau$. These functions serve as a basis and more
complicated structures can be generated as products of $\Omega_{\vec{m}}$.

## Appendix D

Here we give the small ($\tau\to 0$) and large ($\tau\to 1)$ expansions of the
invariant kernels $h_{3},\bar{h}_{3}$. By $h_{3}^{(A)}$ ($\bar{h}_{3}^{(A)}$)
we denote the function which appears in the expression for $h_{3}$
($\bar{h}_{3}^{(A)}$) with the color factor $C_{F}\times A$. We will keep the
logarithmically enhanced and constant terms in both limits. The former is
subtracted from both the exact and approximated three-loop kernel to obtain
the two figures in Eqs. 1 and 2. At $\tau\to 0$ one gets

$\displaystyle h_{3}^{(n_{f}N_{c})}$
$\displaystyle=\frac{5839}{27}-\frac{256}{9}\zeta_{2}+\frac{64}{3}\zeta_{3}-\frac{8}{3}\ln\tau\,,$
$\displaystyle h_{3}^{(n_{f}/N_{c})}$
$\displaystyle=-\frac{17}{9}+16\zeta_{3}\,,$
$\displaystyle\bar{h}_{3}^{(n_{f}/N_{c})}$ $\displaystyle=\frac{8}{3}\,,$
$\displaystyle h_{3}^{(N_{c}^{2})}$
$\displaystyle=-\frac{18520}{27}-\frac{88}{3}\zeta_{3}-\frac{176}{5}\zeta_{2}^{2}+\frac{1744}{9}\zeta_{2}-\frac{46}{3}\ln\tau$
$\displaystyle h_{3}^{(N_{c}^{0})}$
$\displaystyle=-\frac{1186}{9}+32\zeta_{2}+(-32+16\zeta_{2})\ln\tau\,,$
$\displaystyle h_{3}^{(N_{c}^{-2})}$
$\displaystyle=24-8\zeta_{2}-18\ln\tau\,,$
$\displaystyle\bar{h}_{3}^{(N_{c}^{0})}$ $\displaystyle=-\frac{44}{3}\,,$
$\displaystyle\bar{h}_{3}^{(N_{c}^{-2})}$
$\displaystyle=-48\tau\left(\zeta_{2}+\zeta_{3}+\frac{1}{4}-\zeta_{2}\ln\tau\right)\,,$
(D.64)

and for $\tau\to 1$ one obtains

$\displaystyle h_{3}^{(n_{f}N_{c})}$
$\displaystyle=\frac{5695}{27}-\frac{208}{9}\zeta_{2}+\frac{64}{3}\zeta_{3}+\left(-\frac{16}{3}\zeta_{2}+\frac{38}{9}\right)\ln\bar{\tau}\,,$
$\displaystyle h_{3}^{(n_{f}/N_{c})}$
$\displaystyle=\frac{304}{9}\zeta_{2}+16\zeta_{3}-25-\left(\frac{16}{3}\zeta_{2}+\frac{74}{3}\right)\ln\bar{\tau}$
$\displaystyle\quad+\frac{152}{9}\ln^{2}\bar{\tau}-\frac{8}{9}\ln^{3}\bar{\tau}\,,$
$\displaystyle\bar{h}_{3}^{(n_{f}/N_{c})}$
$\displaystyle=\frac{16}{3}\left(\frac{1}{2}-\zeta_{2}+\zeta_{3}\right)+\left(\frac{16}{3}\zeta_{2}-\frac{184}{9}\right)\ln\bar{\tau}$
$\displaystyle\quad+\frac{152}{9}\ln^{2}\bar{\tau}-\frac{8}{9}\ln^{3}\bar{\tau}\,,$
$\displaystyle h_{3}^{(N_{c}^{2})}$
$\displaystyle=-\frac{72}{5}\zeta_{2}^{2}+\frac{1741}{9}\zeta_{2}-\frac{88}{3}\zeta_{3}-\frac{19132}{27}$
$\displaystyle\quad+\left(\frac{4}{3}\zeta_{2}-\frac{187}{18}\right)\ln\bar{\tau}+\left(-\frac{5}{2}+4\zeta_{2}\right)\ln^{2}\bar{\tau}\,,$
$\displaystyle h_{3}^{(N_{c}^{0})}$
$\displaystyle=\frac{136}{5}\zeta_{2}^{2}-\frac{2170}{9}\zeta_{2}+80\zeta_{3}-\frac{94}{3}$
$\displaystyle\quad+\left(-\frac{32}{3}\zeta_{2}-24\zeta_{3}+\frac{548}{3}\right)\ln\bar{\tau}$
$\displaystyle\quad+\left(16\zeta_{2}-\frac{923}{9}\right)\ln^{2}\bar{\tau}+\frac{14}{9}\ln^{3}\bar{\tau}+\frac{4}{3}\ln^{4}\bar{\tau}\,,$
$\displaystyle h_{3}^{(N_{c}^{-2})}$
$\displaystyle=-28\zeta_{2}^{2}-27\zeta_{2}+56\zeta_{3}+28$
$\displaystyle\quad+\left(\frac{115}{2}-12\zeta_{2}-40\zeta_{3}\right)\ln\bar{\tau}$
$\displaystyle\quad+\left(8\zeta_{2}+\frac{11}{2}\right)\ln^{2}\bar{\tau}+\frac{2}{3}\ln^{3}\bar{\tau}-\frac{1}{3}\ln^{4}\bar{\tau}\,,$
$\displaystyle\bar{h}_{3}^{(N_{c}^{0})}$
$\displaystyle=-\frac{136}{5}\zeta_{2}^{2}+\frac{88}{3}\zeta_{2}-\frac{136}{3}\zeta_{3}-\frac{8}{3}$
$\displaystyle\quad+\left(-\frac{16}{3}\zeta_{2}+\frac{1708}{9}\right)\ln\bar{\tau}$
$\displaystyle\quad-\frac{968}{9}\ln^{2}\bar{\tau}+\frac{14}{9}\ln^{3}\bar{\tau}+\frac{4}{3}\ln^{4}\bar{\tau}\,,$
$\displaystyle\bar{h}_{3}^{(N_{c}^{-2})}$
$\displaystyle=-\frac{216}{5}\zeta_{2}^{2}+40\zeta_{2}+8\zeta_{3}+12$
$\displaystyle\quad+(40\zeta_{2}-64\zeta_{3}+40)\ln\bar{\tau}$
$\displaystyle\quad-8(4\zeta_{2}-1)\ln^{2}\bar{\tau}+\frac{2}{3}\ln^{3}\bar{\tau}-\frac{1}{3}\ln^{4}\bar{\tau}\,.$
(D.65)

Here we quote the cusp anomalous dimensions up to three loops for reference
[35, 36, 4],

$\displaystyle\Gamma_{\rm cusp}^{(1)}$ $\displaystyle=4C_{F}\,,$
$\displaystyle\Gamma_{\rm cusp}^{(2)}$
$\displaystyle=C_{F}\left[N_{c}\left(\frac{536}{9}-16\zeta_{2}\right)-\frac{40}{9}n_{f}\right]\,,$
$\displaystyle\Gamma_{\rm cusp}^{(3)}$
$\displaystyle=C_{F}\bigg{[}N_{c}^{2}\left(\frac{176}{5}\zeta_{2}^{2}+\frac{88}{3}\zeta_{3}-\frac{1072}{9}\zeta_{2}+\frac{490}{3}\right)$
$\displaystyle\quad+N_{c}n_{f}\left(-\frac{64}{3}\zeta_{3}+\frac{160}{9}\zeta_{2}-\frac{1331}{27}\right)$
$\displaystyle\quad+\frac{n_{f}}{N_{c}}\left(-16\zeta_{3}+\frac{55}{3}\right)-\frac{16}{27}n_{f}^{2}\bigg{]}\,.$
(D.66)

## References

  * Müller  _et al._ [1994] D. Müller, D. Robaschik, B. Geyer, F. M. Dittes, and J. Hořejši, Wave functions, evolution equations and evolution kernels from light ray operators of QCD, Fortsch. Phys. 42, 101 (1994), arXiv:hep-ph/9812448 .
  * Radyushkin [1996] A. V. Radyushkin, Scaling limit of deeply virtual Compton scattering, Phys. Lett. B380, 417 (1996), arXiv:hep-ph/9604317 [hep-ph] .
  * Ji [1997] X.-D. Ji, Deeply virtual Compton scattering, Phys. Rev. D55, 7114 (1997), arXiv:hep-ph/9609381 [hep-ph] .
  * Moch  _et al._ [2004] S. Moch, J. A. M. Vermaseren, and A. Vogt, The Three loop splitting functions in QCD: The Nonsinglet case, Nucl. Phys. B 688, 101 (2004), arXiv:hep-ph/0403192 .
  * Vogt  _et al._ [2004] A. Vogt, S. Moch, and J. A. M. Vermaseren, The Three-loop splitting functions in QCD: The Singlet case, Nucl. Phys. B 691, 129 (2004), arXiv:hep-ph/0404111 .
  * Moch  _et al._ [2017] S. Moch, B. Ruijl, T. Ueda, J. A. M. Vermaseren, and A. Vogt, Four-Loop Non-Singlet Splitting Functions in the Planar Limit and Beyond, JHEP 10, 041, arXiv:1707.08315 [hep-ph] .
  * Falcioni  _et al._ [2023] G. Falcioni, F. Herzog, S. Moch, and A. Vogt, Four-loop splitting functions in QCD – The quark-quark case –, Phys. Lett. B 842, 137944 (2023), arXiv:2302.07593 [hep-ph] .
  * Müller [1991] D. Müller, Constraints for anomalous dimensions of local light cone operators in $\phi^{3}$ in six-dimensions theory, Z. Phys. C 49, 293 (1991).
  * Müller [1998] D. Müller, Restricted conformal invariance in QCD and its predictive power for virtual two photon processes, Phys. Rev. D 58, 054005 (1998), arXiv:hep-ph/9704406 .
  * Belitsky and Müller [1998] A. V. Belitsky and D. Müller, Predictions from conformal algebra for the deeply virtual Compton scattering, Phys. Lett. B417, 129 (1998), arXiv:hep-ph/9709379 [hep-ph] .
  * Belitsky and Müller [1999] A. V. Belitsky and D. Müller, Broken conformal invariance and spectrum of anomalous dimensions in QCD, Nucl. Phys. B 537, 397 (1999), arXiv:hep-ph/9804379 .
  * Belitsky  _et al._ [2000] A. V. Belitsky, A. Freund, and D. Müller, Evolution kernels of skewed parton distributions: Method and two loop results, Nucl. Phys. B 574, 347 (2000), arXiv:hep-ph/9912379 .
  * Braun and Manashov [2013] V. M. Braun and A. N. Manashov, Evolution equations beyond one loop from conformal symmetry, Eur. Phys. J. C 73, 2544 (2013), arXiv:1306.5644 [hep-th] .
  * Braun and Manashov [2014] V. M. Braun and A. N. Manashov, Two-loop evolution equations for light-ray operators, Phys. Lett. B 734, 137 (2014), arXiv:1404.0863 [hep-ph] .
  * Braun  _et al._ [2016] V. M. Braun, A. N. Manashov, S. Moch, and M. Strohmaier, Two-loop conformal generators for leading-twist operators in QCD, JHEP 03, 142, arXiv:1601.05937 [hep-ph] .
  * Braun  _et al._ [2017] V. M. Braun, A. N. Manashov, S. Moch, and M. Strohmaier, Three-loop evolution equation for flavor-nonsinglet operators in off-forward kinematics, JHEP 06, 037, arXiv:1703.09532 [hep-ph] .
  * Braun  _et al._ [2020] V. M. Braun, A. N. Manashov, S. Moch, and J. Schönleber, Two-loop coefficient function for DVCS: vector contributions, JHEP 09, 117, arXiv:2007.06348 [hep-ph] .
  * Braun  _et al._ [2021a] V. M. Braun, A. N. Manashov, S. Moch, and J. Schönleber, Axial-vector contributions in two-photon reactions: Pion transition form factor and deeply-virtual Compton scattering at NNLO in QCD, Phys. Rev. D 104, 094007 (2021a), arXiv:2106.01437 [hep-ph] .
  * Gao  _et al._ [2022] J. Gao, T. Huber, Y. Ji, and Y.-M. Wang, Next-to-Next-to-Leading-Order QCD Prediction for the Photon-Pion Form Factor, Phys. Rev. Lett. 128, 062003 (2022), arXiv:2106.01390 [hep-ph] .
  * Braun  _et al._ [2021b] V. M. Braun, Y. Ji, and A. N. Manashov, Two-photon processes in conformal QCD: resummation of the descendants of leading-twist operators, JHEP 03, 051, arXiv:2011.04533 [hep-ph] .
  * Braun  _et al._ [2023] V. M. Braun, Y. Ji, and A. N. Manashov, Next-to-leading-power kinematic corrections to DVCS: a scalar target, JHEP 01, 078, arXiv:2211.04902 [hep-ph] .
  * Basso and Korchemsky [2007] B. Basso and G. P. Korchemsky, Anomalous dimensions of high-spin operators beyond the leading order, Nucl. Phys. B 775, 1 (2007), arXiv:hep-th/0612247 .
  * Dokshitzer  _et al._ [2006] Y. L. Dokshitzer, G. Marchesini, and G. P. Salam, Revisiting parton evolution and the large-x limit, Phys. Lett. B 634, 504 (2006), arXiv:hep-ph/0511302 .
  * Remiddi and Vermaseren [2000] E. Remiddi and J. A. M. Vermaseren, Harmonic polylogarithms, Int. J. Mod. Phys. A 15, 725 (2000), arXiv:hep-ph/9905237 .
  * Balitsky and Braun [1989] I. I. Balitsky and V. M. Braun, Evolution Equations for QCD String Operators, Nucl. Phys. B 311, 541 (1989).
  * Gribov and Lipatov [1972a] V. N. Gribov and L. N. Lipatov, Deep inelastic e p scattering in perturbation theory, Sov. J. Nucl. Phys. 15, 438 (1972a).
  * Gribov and Lipatov [1972b] V. N. Gribov and L. N. Lipatov, $e^{+}e^{-}$ pair annihilation and deep inelastic $ep$ scattering in perturbation theory, Sov. J. Nucl. Phys. 15, 675 (1972b).
  * Dokshitzer and Marchesini [2007] Y. L. Dokshitzer and G. Marchesini, $\mathcal{N}=4$ SUSY Yang-Mills: three loops made simple(r), Phys. Lett. B 646, 189 (2007), arXiv:hep-th/0612248 .
  * Beccaria and Forini [2009] M. Beccaria and V. Forini, Four loop reciprocity of twist two operators in $\mathcal{N}=4$ SYM, JHEP 03, 111, arXiv:0901.1256 [hep-th] .
  * Alday  _et al._ [2015] L. F. Alday, A. Bissi, and T. Lukowski, Large spin systematics in CFT, JHEP 11, 101, arXiv:1502.07707 [hep-th] .
  * Kotikov  _et al._ [2009] A. V. Kotikov, A. Rej, and S. Zieme, Analytic three-loop Solutions for $\mathcal{N}=4$ SYM Twist Operators, Nucl. Phys. B 813, 460 (2009), arXiv:0810.0691 [hep-th] .
  * Beccaria  _et al._ [2009] M. Beccaria, V. Forini, T. Lukowski, and S. Zieme, Twist-three at five loops, Bethe Ansatz and wrapping, JHEP 03, 129, arXiv:0901.4864 [hep-th] .
  * Marboe  _et al._ [2015] C. Marboe, V. Velizhanin, and D. Volin, Six-loop anomalous dimension of twist-two operators in planar $\mathcal{N}=4$ SYM theory, JHEP 07, 084, arXiv:1412.4762 [hep-th] .
  * Marboe and Velizhanin [2016] C. Marboe and V. Velizhanin, Twist-2 at seven loops in planar $\mathcal{N}=4$ SYM theory: full result and analytic properties, JHEP 11, 013, arXiv:1607.06047 [hep-th] .
  * Polyakov [1980] A. M. Polyakov, Gauge Fields as Rings of Glue, Nucl. Phys. B 164, 171 (1980).
  * Korchemsky and Radyushkin [1987] G. P. Korchemsky and A. V. Radyushkin, Renormalization of the Wilson Loops Beyond the Leading Order, Nucl. Phys. B 283, 342 (1987).
  * Henn  _et al._ [2020] J. M. Henn, G. P. Korchemsky, and B. Mistlberger, The full four-loop cusp anomalous dimension in $\mathcal{N}=4$ super Yang-Mills and QCD, JHEP 04, 018, arXiv:1911.10174 [hep-th] .
  * von Manteuffel  _et al._ [2020] A. von Manteuffel, E. Panzer, and R. M. Schabinger, Cusp and collinear anomalous dimensions in four-loop QCD from form factors, Phys. Rev. Lett. 124, 162001 (2020), arXiv:2002.04617 [hep-ph] .
  * Panzer [2015] E. Panzer, Algorithms for the symbolic integration of hyperlogarithms with applications to Feynman integrals, Comput. Phys. Commun. 188, 148 (2015), arXiv:1403.3385 [hep-th] .
  * Duhr and Dulat [2019] C. Duhr and F. Dulat, PolyLogTools — polylogs for the masses, JHEP 08, 135, arXiv:1904.07279 [hep-th] .
  * Kotikov  _et al._ [2003] A. V. Kotikov, L. N. Lipatov, and V. N. Velizhanin, Anomalous dimensions of Wilson operators in $\mathcal{N}=4$ SYM theory, Phys. Lett. B 557, 114 (2003), arXiv:hep-ph/0301021 .
  * Gradshteyn and Ryzhik [1994] I. S. Gradshteyn and I. M. Ryzhik,  _Table of integrals, series, and products_ , russian ed. (Academic Press, Inc., Boston, MA, 1994) pp. xlviii+1204, translation edited and with a preface by Alan Jeffrey.


# Measurement Based Evaluation and Mitigation of Flood Attacks on a LAN Test-
Bed††thanks: This research has been supported by the European Commission H2020
Program through the IoTAC Research and Innovation Action, under Grant
Agreement No. 952684.

Mohammed Nasereddin and Mert Nakıp  Institute of Theoretical & Applied
Informatics
Polish Acad. Sci., 44-100 Gliwice, Poland
<EMAIL_ADDRESS>Erol Gelenbe<EMAIL_ADDRESS>Institute of
Theoretical & Applied Informatics
Polish Acad. Sci., 44-100 Gliwice, Poland
& Yaşar University, Bornova, Turkey,
& Lab. I3S Université Côte d’Azur, Nice, France

###### Abstract

The IoT is vulnerable to network attacks, and Intrusion Detection Systems
(IDS) can provide high attack detection accuracy and are easily installed in
IoT Servers. However, IDS are seldom evaluated in operational conditions which
are seriously impaired by attack overload. Thus a Local Area Network test-bed
is used to evaluate the impact of UDP Flood Attacks on an IoT Server, whose
first line of defence is an accurate IDS. We show that attacks overload the
multi-core Server and paralyze its IDS. Thus a mitigation scheme that detects
attacks rapidly, and drops packets within milli-seconds after the attack
begins, is proposed and experimentally evaluated.

###### Index Terms:

Internet of Things, Local Area Networks, Cybersecurity, Random Neural
Networks, G-Networks, UDP Flood Attacks, Intrusion Detection and Mitigation

##  I Introduction

Denial of service (DoS) disables systems or networks by flooding them with
huge streams of requests, causing reputational damage, with financial and
productivity losses [1]. In the last year, a $150\%$ increase in such attacks
occurred worldwide [2], targeting the IoT, industrial control systems, power
grids and transportation systems [3, 4, 5, 6], with DoS and Botnet attacks
spreading via their victims, who then become attackers [7, 8, 9]. Flood
attacks [10] overwhelm networks with large numbers of forged-source address
packets, causing delayed or lost data, inaccurate or incomplete readings and
overload [11]. Thus these threats require effective Intrusion Detection
Systems (IDS).

While much of the literature on IDS evaluates them under ideal conditions
where attack traffic is treated as data, this paper compares “ideal” results
about attack detection (AD) algorithms, with system performance measurements
in a LAN environment. Section I-A reviews related work, Section II describes
the experimental test-bed that we use and Section II-A briefly discusses the
IDS and its performance under ideal and real-world conditions. Section III
presents measurements during different UDP Flood attacks, and summarizes
improvements resulting from simple attack mitigation. Finally, Section IV
concludes the paper and outlines future work directions.

###  I-A Related Work

Using test-beds to evaluate IDSs was recommended in early work [12], and
several test-beds for cyber-physical systems, industrial control and IoT
environments have been described [13], such as real-time test-bed for cyber-
physical systems [14], power systems [15], and wind farms [16]. In [17], a
semi-physical test-bed for ICSs was proposed, while in [18], a low-cost Smart
Grid test-bed for SCIDS systems was evaluated for TCP flood attacks. Test-beds
for SCADA systems are discussed in [19, 20, 21]. In [22], a test-bed using six
NetFlow tools for collecting, analyzing, and displaying data was proposed for
HTTP-GET flood attacks in a WAN. In [23], the impact of attack datasets on IoT
systems is discussed and real-time data collection for DNS amplification is
investigated, while DoS attacks on software defined networks are discussed in
[24]. In [25], DoS attacks on an autonomous vehicle test-bed are described.
Attack datasets are reviewed in [22, 23], while the present paper uses the
MHDDoS repository [26] of real-world DoS attacks.

##  II Experimental Setup

Figure 1: Testing Environment using Ethernet for communications, with
Raspberry Pi machines acting as forwarders of normal and attack traffic, and
an Intel 8-Core Processer used as a Server to process incoming packet traffic
and run the IDS algorithm.

In this paper, a physical test-bed environment is constructed to evaluate IDS
algorithms in more realistic conditions, with an arbitrary number of linked
devices, multiple sources of normal and attack traffic, and a Server that
supports the UDP Protocol for incoming traffic, runs the IDS algorithm, and
processes the incoming packets’ contents. The devices that generate benign or
attack traffic, are embodied by Raspberry Pi $4$ Model B Rev $1.2$ machines
(RPi$1$ and RPi$2$), each with a $1.5$GHz ARM Cortex-A$72$ quad-core processor
and $2$GB LPDDR$4-3200$ SDRAM, running Raspbian GNU/Linux $11$ (bullseye), a
Debian-based operating system for the Raspberry Pi hardware. A Server with
eight Intel Core i$7-8705$G processors acts as the receiver of the packet
traffic and is responsible for detecting the attack and storing the arriving
packets. It has $16$GB of RAM, a $500$GB hard drive, and runs Linux
$5.15.0-60-$ generic $66-$Ubuntu SMP, an Ubuntu-based operating system. Its
cores run at $3.10$GHz. UDP traffic is carried over Ethernet connections
between all devices via the Hub in Figure 1, without ACKs or error recovery
[27].

###  II-A The IDS and its Ideal Performance

The IDS used in this paper, based on the Deep Random Neural Network (DRNN)
[28], is shown in Figure 2. It learns from the  first $500$ benign packets
received by the Server, with metrics $x^{i}=[x_{i}^{1},x_{i}^{2},x_{i}^{3}]$
related to successive sets of packets. The IDS predicts the expected metrics:
$\hat{x}_{i}=[\hat{x}_{i}^{1},\hat{x}_{i}^{2},\hat{x}_{i}^{3}]$, and the
difference between its input and the prediction yields the decision variables
$y_{i}$ (attack or non-attack), which are also used to update the algorithm’s
weights. The IDS uses the DRNN [29], a Random Neural Network [30] with soma-
to-soma triggering between neurons. This IDS provides accurate detection with
different datasets [31, 32, 33, 34], with the excellent statistical
performance shown in Figure 3, which reports the Accuracy, TPR, and TNR, for a
$10$ second attack. This IDS attains high accuracy both when a fixed threshold
$\gamma=0.3$ is used, and for the best threshold $\gamma=0.3787$, exhibiting
$99.7\%$ Accuracy and TPR, while TNR is $98.48\%$. These results are not
significantly different when the attack lasts for $60$ seconds. Note that
Figure 4 shows that, due to the decision delay, the IDS may raise an alarm
just after an attack ceases.

Figure 2: The structure of the IDS system that computes the decision variable
$y_{i}$ from the network traffic metrics $[x_{i}^{1},x_{i}^{2},x_{i}^{3}]$
with the DRNN based Auto-Associative Random Neural Network (AAD RNN) and the
postprocessing module. Figure 3: The performance of the IDS with $\gamma=0.3$,
and compared with the best value of $\gamma=0.3787$, is evaluated for
Accuracy, TPR, and TNR, in an experiment where RPi2 starts a UDP Flood attack
lasting $10$ seconds. Figure 4: The IDS’s binary decisions are shown for
$\gamma=0.3$, when the RPi2 starts a UDP Flood attack lasting $10$ seconds.

##  III Experiments with Normal and Attack Traffic

In Figure 5, the Server receives packets from linked devices on port $5555$,
which are then passed to the buffer manager by the network protocol, and
queued for analysis at the IDS, whose decisions are based on the average of a
batch of $10$ successive packets, that are being classified as normal or
attack traffic. Packets classified as “normal” are forwarded to the packet
content processor for the rest of the Server’s operations.

When there is no attack, each RPi device generates normal IP packet traffic
containing the device’s own CPU temperature and transmits it to the Server
every ($1$) second using UDP. RPi$2$ is programmed to send both normal and
attack traffic via random sampling, and generates attacks from the public
repository MHDDoS [26]: each $1$ second, it initiates a UDP Flood attack with
a probability of $0.10$, or sends normal traffic packet with a probability of
$0.90$. RPi$1$ only sends normal traffic.

Figure 5: Schematic organization of the Server that supports the IDS.
Mitigation is based on triggering “packet drop” decisions for all packets in
the IDS Input Buffer, when it detects a majority of attack packets among the
most recent $20$ packets. The IDS then resumes testing the incoming packets,
and the decision and mitigatiion process is repeated.

Figure 6, displays an intense flow of $1032$ byte attack packets during $10$
seconds, while “normal” traffic typically consists of two small packets sent
each second, for each IoT device. Figure 7 shows that the packet queue length
at the Server rises sharply infront of the IDS, with a gradual decrease after
the attack. Figure 8 shows (Above) the effect of a $60$ second attack on the
Server’s intermittently paralyzed packet processing rate ($y$-axis in
packets/sec), and (Below) the resulting huge input queue length. Figures 7 and
8 show that a $10$ second attack floods the packet queue, and the IDS
completes the analysis of the accumulated packets over a long $15$ minute
period, while for a $60$ second attack, the IDS is intermittently paralyzed
and its analysis can last $5.85$ hours, since the Server’s cores are busy
handling the incoming traffic. Thus the effect of an attack typically lasts
far longer than the attacker’s activity.

Figure 6: The difference between the normal and attack traffic on the Server
that is targeted by a UDP Flood attack.

Figure 7: The top figure shows the queue length infront of the IDS in a $10$
second attack, and the vertical red dashed lines show the active duration of
the attack originating in the compromised device RPi$2$. The bottom figure
shows the delay before the packet is processed by the IDS.

Figure 8: At the top, the effect of a $60$ second UDP Flood attack on the IDS
traffic processing rate in packets per second, is shown when the attack
duration is $60$ seconds. The corresponding packet queue length infront of the
IDS is shown at the bottom.

We now report system measurements when a novel mitigation action occurs: if
the IDS decides that the majority of the $20$ most recent packets are attack
packets, then the input buffer contents and all incoming packets within the
next $30$ seconds are dropped. This action is repeated at the end of the $30$
second window. Figure 9 displays the queue length in the input buffer when the
attack mitigation is performed during a UDP Flood attack which lasts $10$
seconds: the queue length increases until the IDS processes $20$ packets and
decides to empty the buffer, and we observe that the attack is mitigated
successfully.

Figure 9: During the $10$ second attack, the decision to drop packets results
in a very short packet queue length, avoiding Server and IDS paralysis.

Figure 10 shows the queue length when an attack that lasts $60$ seconds is
mitigated: the buffer length increases to $22$ packets, which is small
compared to the value without mitigation shown in Figure 8: the mitigation
decision was taken twice, the second time between $162$ to $192$ seconds after
the start of the experiment, and the Server could then operate normally
without being paralyzed.

Figure 10: The figure shows that during the attack’s $60$ seconds, the
mitigation decision occurs twice, with the second mitigation occurring after a
detection that takes place between $162$ and $192$ seconds.

##  IV Conclusions and Future Work

Despite the high accuracy of an IDS installed on a Server that receives
traffic from devices in a LAN network test-bed subjected to UDP Flood attacks,
we observe that while short attacks are accurately detected, longer attacks
may be detected with greater delay due to Server overload. Thus fast
mitigation is proposed to discard attacking traffic with rapidly taken
decisions. Future work will study mitigation to optimize the traffic that is
discarded, the frequency with which the IDS analyzes incoming traffic, as well
as the minimization of benign traffic loss, and the system’s energy
consumption [35].

Acknowledgment Support from the European Commission’s H2020 IoTAC Project GA
No. 952684, and H2020 DOSS Project GA No. 101120270, is gratefully
acknowledged.

## References

  * [1] E. Gelenbe, P. Campegiani, T. Czachórski, S. K. Katsikas, I. Komnios, L. Romano, and D. Tzovaras, Eds., _Security in Computer and Information Sciences, First International ISCIS Security Workshop, Euro-CYBERSEC 2018, London, February 26-27, 2018_. Springer, Cham, 2018. [Online]. Available: http://library.oapen.org/handle/20.500.12657/23295
  * [2] S. Staff, “Organizations fought an average of 29.3 attacks daily in late 2022,” Feb 2023. [Online]. Available: https://www.securitymagazine.com/articles/98958-organizations-fought-an-average-of-293-attacks-daily-in-late-2022
  * [3] G. Öke, et al., “Detecting denial of service attacks with bayesian classifiers and the random neural network,” in _IEEE International Fuzzy Systems Conf._ IEEE, 2007, pp. 1–6.
  * [4] S. Evmorfos, et al, “Neural network architectures for the detection of SYN flood attacks in IoT systems,” in _Proceedings of the 13th ACM International Conference on PErvasive Technologies Related to Assistive Environments_ , 2020, pp. 1–4.
  * [5] L. Rajesh and P. Satyanarayana, “Detecting flooding attacks in communication protocol of industrial control systems,” _International Journal of Advanced Computer Science and Applications_ , vol. 11, no. 1, pp. 396–401, 2020\.
  * [6] Y. Al-Hadhrami and F. K. Hussain, “Ddos attacks in iot networks: a comprehensive systematic literature review,” _World Wide Web_ , vol. 24, no. 3, pp. 971–1001, 2021.
  * [7] N. Statt, “How an army of vulnerable gadgets took down the web today,” October 2016. [Online]. Available: https://www.theverge.com/2016/10/21/13362354/dyn-dns-ddos-attack-cause-outage-status-explained
  * [8] B. Tushir, H. Sehgal, R. Nair, B. Dezfouli, and Y. Liu, “The impact of DoS attacks on resource-constrained IoT Devices: A study on the Mirai attack,” _arXiv preprint arXiv:2104.09041_ , 2021.
  * [9] H. Sinanović and S. Mrdovic, “Analysis of Mirai malicious software,” in _2017 25th International Conference on Software, Telecommunications and Computer Networks (SoftCOM)_. IEEE, 2017, pp. 1–5.
  * [10] Cloudflare. [Online]. Available: https://www.cloudflare.com/learning/ddos/udp-flood-ddos-attack/
  * [11] J. Mirkovic, et al., “Accurately measuring denial of service in simulation and testbed experiments,” _IEEE Transactions on Dependable and Secure Computing_ , vol. 6, no. 2, pp. 81–95, 2008.
  * [12] J. Mirkovic, S. Fahmy, P. Reiher, and R. K. Thomas, “How to test dos defenses,” in _2009 Cybersecurity App. & Technology Conf._ IEEE, 2009, pp. 103–117.
  * [13] O. A. Waraga, M. Bettayeb, Q. Nasir, and M. A. Talib, “Design and implementation of automated iot security testbed,” _Computers & Security_, vol. 88, p. 101648, 2020.
  * [14] C. B. Vellaithurai, S. S. Biswas, and A. K. Srivastava, “Development and application of a real-time test bed for cyber-physical system,” _IEEE Systems Journal_ , vol. 11, no. 4, pp. 2192–2203, 2015.
  * [15] A. Ashok, et al., “Testbed-based performance evaluation of attack resilient control for agc,” in _2016 Resilience Week (RWS)_. IEEE, 2016, pp. 125–129.
  * [16] V. K. Singh, R. Sharma, and M. Govindarasu, “Testbed-based performance evaluation of attack resilient control for wind farm scada system,” in _2020 IEEE Power & Energy Society General Meeting (PESGM)_. IEEE, 2020, pp. 1–5.
  * [17] M. Kaouk, F.-X. Morgand, and J.-M. Flaus, “A testbed for cybersecurity assessment of industrial and iot-based control systems,” in _Lambda Mu 2018-21è Congrès de Maîtrise des Risques et Sûreté de Fonctionnement_ , 2018.
  * [18] M. Annor-Asante and B. Pranggono, “Development of smart grid testbed with low-cost hardware and software for cybersecurity research and education,” _Wireless Personal Communications_ , vol. 101, pp. 1357–1377, 2018.
  * [19] A. Ghaleb, S. Zhioua, and A. Almulhem, “Scada-sst: a scada security testbed,” in _2016 World Congress on Industrial Control Systems Security (WCICSS)_. IEEE, 2016, pp. 1–6.
  * [20] A. Tesfahun and D. L. Bhaskari, “A scada testbed for investigating cyber security vulnerabilities in critical infrastructures,” _Automatic Control and Computer Sciences_ , vol. 50, pp. 54–62, 2016.
  * [21] B. Reutimann and I. Ray, “Simulating measurement attacks in a scada system testbed,” in _Critical Infrastructure Protection XV: 15th IFIP WG 11.10 International Conference, ICCIP 2021, Virtual Event, March 15–16, 2021, Revised Selected Papers 15_. Springer, 2022, pp. 135–153.
  * [22] S.-U. Park and S.-M. Hwang, “Test bed construction for apt attack detection,” _International Journal of Control and Automation_ , vol. 11, no. 4, pp. 175–186, 2018.
  * [23] R. Arthi and S. Krishnaveni, “Design and development of iot testbed with ddos attack for cyber security research,” in _2021 3rd International Conference on Signal Processing and Communication (ICPSC)_. IEEE, 2021, pp. 586–590.
  * [24] A. P. Wright and N. Ghani, “A testbed for the evaluation of denial of service attacks in software-defined networks,” in _2019 SoutheastCon_. IEEE, 2019, pp. 1–6.
  * [25] P. V. Sontakke and N. B. Chopade, “Impact and analysis of denial-of-service attack on an autonomous vehicle test bed setup,” in _Proceedings of Third International Conference on Intelligent Computing, Information and Control Systems: ICICCS 2021_. Springer, 2022, pp. 221–236.
  * [26] “MHDDoS - DDoS Attack Script With 56 Methods,” Online, May 2022, accessed: 2023-02-22. [Online]. Available: https://github.com/MatrixTM/MHDDoS
  * [27] S. Kumar and S. Rai, “Survey on transport layer protocols: Tcp & udp,” _International Journal of Computer App._ , vol. 46, no. 7, pp. 20–25, 2012\.
  * [28] O. Brun, et al., “Iot attack detection with deep learning,” in _Security in Computer and Information Sciences, First International ISCIS Security Workshop, Euro-CYBERSEC 2018, London, February 26-27, 2018_. Springer, Cham, 2018. [Online]. Available: http://library.oapen.org/handle/20.500.12657/23295
  * [29] E. Gelenbe and Y. Yin, “Deep learning with random neural networks,” in _2016 International Joint Conference on Neural Networks (IJCNN)_. IEEE, 2016, pp. 1633–1638.
  * [30] E. Gelenbe, “Random neural networks with negative and positive signals and product form solution,” _Neural computation_ , vol. 1, no. 4, pp. 502–510, 1989.
  * [31] M. Nakip and E. Gelenbe, “MIRAI Botnet attack detection with Auto-Associative Dense Random Neural Network,” in _IEEE Global Communications Conference (Globecom)_ , 2021, pp. 1–6.
  * [32] M. Nakıp and E. Gelenbe, “Botnet attack detection with incremental online learning,” in _Security in Computer and Information Sciences: Second International Symposium, EuroCybersec 2021, Nice, France, October 25–26, 2021, Revised Selected Papers_ , E. Gelenbe, et al., Ed. Springer, 2022, pp. 51–60.
  * [33] E. Gelenbe and M. Nakıp, “Traffic based sequential learning during botnet attacks to identify compromised iot devices,” _IEEE Access_ , vol. 10, pp. 126 536–126 549, 2022.
  * [34] E. Gelenbe and M. Nakıp, “G-networks can detect different types of cyberattacks,” in _MASCOTS’22: 30th International Symposium on the Modeling, Analysis, and Simulation of Computer and Telecommunication Systems, IEEE Computer Society_ , pp. 1–6.
  * [35] B. Pernici, et al., “What IS can do for environmental sustainability: a report from caise’11 panel on green and sustainable IS,” _Comm. Association for Information Systems_ , vol. 30, no. 1, p. 18, 2012.