text:
#1ζ_ #1^
ζ_2^ 2
ζ_3^ 2

n_c^ 2
n_c^ 3

C^ 2_A
C^ 3_A
C^ 4_A

C^_ F
C^  2_ F
C^  3_ F
C^  4_ F

n^_ f
n^ 0_ f
n^ 1_ f
n^ 2_ f
n^ 3_ f
n^ 4_ f

T^ 2_F
T^ 3_F
T^ 4_F

d^ abcd_abc/n_c

μ_r^ 2
μ_f^ 2

(#1)α_s^ #1
(#1)a_s^ #1

(#1)ln^ #1(1-x)

Q^  2

(#1,#2) 𝐒_#1,#2
(#1,#2,#3) 𝐒_#1,#2,#3
(#1,#2,#3,#4) 𝐒_#1,#2,#3,#4
(#1,#2,#3,#4,#5) 𝐒_#1,#2,#3,#4,#5
(#1,#2,#3,#4,#5,#6) 𝐒_#1,#2,#3,#4,#5,#6

(#1) D_#1
(#1,#2)  D_#1^ #2

(#1,#2)D_#1^ #2

Nikhef 2022-021 November 2022

LTH 1324

Four-loop large-n$_{\:\!\bf f}$ contributions to the

non-singlet structure functions F$_{\bf 2}$ and F$_{\bf L}$

A. Basdew-Sharma, A. Pelloni

Nikhef Theory Group

Science Park 105, 1098 XG Amsterdam, The Netherlands

F. Herzog

Higgs Centre for Theoretical Physics, School of Physics and

University of Edinburgh, Edinburgh EH9 3FD, Scotland, UK

A. Vogt

Department of Mathematical Sciences, University of Liverpool

Liverpool L69 3BX, United Kingdom

We have calculated the $\nfs$ and $\nft$ contributions to the flavour
non-singlet structure functions $F_2$ and $F_L$ in inclusive deep-inelastic
scattering at the fourth order in the strong coupling $\als\:\!$.
The coefficient functions have been obtained by computing a very large
number of Mellin-$N$ moments using the method of differential equations,
and then  determining the analytic forms in $N$ and Bjorken-$x$ from these.

Our new $\nfs$ terms are numerically much larger than the $\nft$ leading
large-$\nf$ parts which were already known; they agree with predictions of
the threshold and high-energy resummations. Furthermore our calculation
confirms the earlier determination of the four-loop $\nfs$ part of the
corresponding anomalous dimension.

Via the no-$\pi^2$ conjecture$/$theorem for Euclidean physical quantities,
we predict the $\z4\:\!\nft$ part of the fifth-order anomalous
dimension for the evolution of non-singlet quark distributions.

§ INTRODUCTION

Inclusive deep-inelastic lepton nucleon scattering (DIS) via the exchange
of an electro-weak gauge boson is an experimental and theoretical benchmark
process of perturbative QCD.

Data on its main structure functions provide a rather direct determination
of (linear combinations of) the quark momentum distributions of the nucleon.
The structure function $F_2(x,\Qs)$, in particular, has been determined
in the past decades over a wide range of the scale $\Qs$ ($\, = - q^2$,
where $q$ is the momentum of the exchanged boson) and the Bjorken variable
$x$ ($\,= \Qs/(2\,p\cdot q)$, where $p$ is the momentum of the nucleon)
in fixed-target experiments and at the electron-proton collider HERA,
see ref. [1] and references therein.
Further measurements of inclusive DIS are planned for future
facilities, in particular the Electron Ion Collider (EIC) at Brookhaven
National Lab [2, 3]
and the Large Hadron Electron Collider (LHeC)
[4, 5].

Precise determinations of the quark momentum distributions
$q_{\:\! \sf i}^{}(\xi,\Qs)$
(with $\xi = x$ at the leading-order of perturbative QCD) as well as, less
directly, of the gluon distribution $g(\xi,\Qs)$ and the strong coupling
$\als$ from structure-function data require higher-order calculations of
the corresponding coefficient functions (partonic structure functions).
These coefficient functions are of relevance also beyond the cross sections
for inclusive DIS, see, e.g., refs. [6, 7]
on Higgs production in vector-boson fusion and ref. [8]
on jet production in DIS.

For the quantities under consideration in this article, the flavour
non-singlet contributions to $F_2$ and the longitudinal structure function
$F_L$, the second-order corrections have been calculated and verified long
ago <cit.>. The corresponding three-loop expressions were obtained in
ref. [13] and recently re-calculated in
ref. [14].

At the fourth order, only the lowest five Mellin-$N$ moments have been
computed so far [15, 16, 17] using the
Forcer program [18], in addition to the leading
terms in the limit of a large number of flavours $\nf$
[19, 20].

In the present article, we take the next step towards the determination of
the fourth-order non-singlet coefficient functions $c_{2,\rm ns}^{\,(4)}(x)$
and $c_{L,\rm ns}^{\,(4)}(x)$ and compute their doubly fermionic $\nfs$
These results are obtained by a new method which allows the determination
of their moments up to very high (even) values of $N$, beyond $N = 1000$,
from which the analytic dependence on $N$, and hence on $x$, can be
re-constructed in terms of harmonic sums [21]
and harmonic polylogarithms (HPLs) [22], respectively.
As a by-product, we have checked the $\nfs$-contributions to the four-loop
non-singlet splitting function $P_{\rm ns}^{\,(3)+}(x)$ of
ref. [23].

The remainder of this article is organized as follows:
In section 2 we briefly recall the theoretical framework for the
coefficient functions in inclusive DIS and their determination to the
fourth order in $\als$.
In section 3 we describe our method of the calculation based on
iteratively solving a system of recurrence relations which is derived
via the method of differential equations
[24, 25, 26, 27].
The analytic results for the coefficient functions in $N$-space are
presented in section 4, which also includes a resulting partial
prediction for the five-loop non-singlet splitting function
$P_{\rm ns}^{\,(4)}(N)$.
The corresponding $x$-space coefficient functions and their threshold
and high-energy limits are written down and discussed in section 5.
We summarize our method and results and give a brief outlook in section 6.

§ THEORETICAL FRAMEWORK AND NOTATIONS

The subject of our computations is unpolarized inclusive lepton-nucleon DIS

(k)  +  (p)

(k^ ')  +   X

at the lowest order of QED (i.e., via the exchange of one photon with
momentum $q = k - k^{\,\prime}$). $X$ stands for all hadronic states
allowed by quantum number conservation.

The double-differential cross section
in $\Qs = - q^2$ and $x = \Qs/(2\,p\cdot q)$
for this process can be expressed as the product of a calculable and
well-known leptonic tensor and the hadronic tensor

W_ μν(p,q)   =
( q_μ q_ν/q^2 - g_μν )
- ( q_μ + 2x p_μ ) ( q_ν + 2x p_ν )
1/2xq^2  F_2(x,Q^2)    .

Neglecting contributions that are suppressed at large scales by powers of
$1/\Qs$, the structure functions
$\,{\cal F}_2= 1/x\, F_2$ and $\,{\cal F}_{\!L}= 1/x\, (F_2 - 2x\:\! F_1)$
can be expressed in terms of the universal but perturbatively incalculable
quark and gluon parton distribution functions (PDFs),
$q_{\:\! \sf i}^{}(\xi,\Qs)$ and $g(\xi,\Qs)$,
and the perturbative coefficient functions
${\cal C}_{a,\rm p}(x,\Qs)$. In the present article, we are specifically
interested in non-singlet (combinations of) structure functions
$F_{a,\rm ns\,}$, such as $F_a^{\,\rm proton} - F_a^{\,\rm neutron}$, which
decouple from the gluon distribution, viz

F_a,ns(x,)    =
[   C_a,ns ⊗q_ ns ] (x,)

where $\otimes$ abbreviates the Mellin convolution. The non-singlet
combinations $q_{\:\!\rm ns}$ of quark distributions are normalized such
that the expansion of the coefficient functions in powers of
$\ars \equiv \als(\Qs) /(4\:\!\pi)$ is given by

C_a,ns(x,)   =   (1-δ_aL)   δ(1-x)
+  ∑_n=1 (n ) c_a,ns^ (n)(x)    .

Here and below we identify the  renormalization and mass-factorization
scales $\muRs$ and $\muFs$, at which the strong coupling and the PDFs
are evaluated, with $\Qs$. The dependence on $\muRs$ and $\muFs$ can be
readily reconstructed a posteriori, see, e.g., eqs. (2.17) and (2.18)
of ref. [28].

In terms of the general framework, our determination of the coefficient
functions uses the method set out (and applied to the lowest
moments at the third order) in
refs. [29, 30],
see also refs. [12, 13, 14]:

The cross section for inclusive DIS for quark external states, projected
onto the structure functions ${\cal F}_{a}$, is related by the optical
theorem to the imaginary parts of the corresponding amplitudes for
photon-quark forward scattering. Via a dispersion relation the
coefficients of $[(2\,p\cdot q)/\Qs]^N = 1/x^N$ lead to the even-integer
(see also ref. [31]) Mellin-$N$ moments

F_a,ns(N,)    =
∫_0^1  dx  x^ N-1   F_a,ns(x,)

of the bare partonic structure functions. These are computed from Feynman
diagrams in dimensional regularization with $D = 4 -2\:\!\ep$ dimensions.

After the  renormalization of the coupling constant to the fourth order,

()    =    ( 1
-  β_0^/   +  ( β_0^2/^2
- β_1^/2)  (2)
-  ( β_0^3/^3
- 7β_1^β_0^/6^2
+ β_2^/3 )  (3) +  …)

with $\,\beta_0 = 11 - 2/3\,\nf\,$ etc in QCD
[32, 33],
the left-hand-side of eq. (<ref>) can be written as

F_a,ns(N,)    =
C_a,ns(N,)   Z_ns(N,)

The $D$-dimensional coefficient function $\widetilde{\cal C}_{a,\rm ns}$
includes additional terms with positive powers of $\ep$ on top of
the Mellin transform of eq. (<ref>), i.e.,

c̃_a,ns^ (n)(N)    =
c_a,ns^ (n)(N)
+    a_a,ns^(n)(N)
+  ^2   b_a,ns^ (n)(N)
+  ^3   d_a,ns^ (n)(N)
+  …   .

The quantity $Z_{\rm ns}$ which renormalizes the non-singlet quark
distributions is given by

Z_ns     =
+   1/ γ_ns^(0)
+ (2) [
1/2 ^2 { (γ_ns^(0)
- β_0 ) γ_ns^(0) }
+ 1/2   γ_ns^(1) ]

+ (3)  [   1/6 ^3
{ ( γ_ns^(0) - 2 β_0 )
( γ_ns^(0) - β_0 )
γ_ns^(0) }

.      + 1/6 ^2
{ ( 3 γ_ns^(0) - 2 β_0 )
- 2 β_1 γ_ns^(0) }
+ 1/3   γ_ns^(2) ]

+ (4)  [
1/24 ^4
{ ( γ_ns^(0) - 3 β_0 )
( γ_ns^(0) - 2 β_0 )
( γ_ns^(0) - β_0 )
γ_ns^(0) }

.      + 1/12 ^3
{ ( ( 3 γ_ns^(0) - 7 β_0)
γ_ns^(0) + 3 β_0^2 ) γ_ns^(1)
- 2 ( 2 γ_ns^(0) - 3 β_0 )
β_1 γ_ns^(0) }

.      + 1/24 ^2
{ 2 ( 4 γ_ns^(0) - 3 β_0)
+ 3 ( γ_ns^(1) - 2 β_1) γ_ns^(1)
- 6 β_2 γ_ns^(0) }
+ 1/4   γ_ns^(3)

to the fourth order. Here $\gamma_{\rm ns}^{(n)}(N)$ – the arguments $N$
have been suppressed in eq. (<ref>) for brevity – are the N$_{}^{n}$LO
non-singlet anomalous dimensions related by

γ_ns^(n)(N)    = - ∫_0^1  dx   x^N-1

to the expansion coefficients of the -scheme splitting function for
the evolution of flavour differences of the sums (hence `+') of quark and
antiquark PDFs,

P_ns^+(x,)    =
∑_n=0 (n+1) P_ns^+(n)(x)    .

Inserting the expansions (<ref>), (<ref>) and (<ref>)
into eq. (<ref>), the anomalous dimension and ($D$-dimensional)
coefficient functions can be extracted order by order from the results
of the diagram calculations. In order to obtain the fourth-order
coefficient functions $c_{a,\rm ns}^{\,(4)}$, the lower-order calculations
need to include terms up to $\ep^{\:\!4-n}$ at order $\as(n)$. In particular,
the determination of the $\nfs$ contributions to $c_{a,\rm ns}^{\,(4)}$
requires the $\nf$ parts of $a_{a,\rm ns}^{(3)}$ which were beyond the
scope of ref. [13] – at the time only the integrals
required for one simpler Lorentz projection of $W_{\:\!\!\mu\nu}$ were
extended to this accuracy for ref. [34].

§ METHOD AND  COMPUTATIONS

In terms of the diagram sets and the treatment to the point at which
the Feynman integrals are evaluated, our computation is closely related
to that of third-order fermionic ($\nf$ and $\nfs$) contributions
in ref. [35] and the non-singlet part of
ref. [23].
Our evaluation of the Feynman integrals is entirely different, though,
from both.
In ref. [35] the analytic $N$-dependence was determined
by setting up and solving, in a far from fully automated manner,
complicated systems of difference equations.
In ref. [23] the even moments were computed to
$N \!=\! 22$ for the $\cf \ca\, \nfs$ terms and to $N \!=\! 42$ for the
$\cfs\, \nfs$ terms using Forcer [18]. From these
it was possible, just, to reconstruct the analytic $N$-dependence of
the four-loop anomalous dimensions using  all available physics
constraints and systems of Diophantine equations.

Below we give details on the methods used in the present computation.
We believe that some of the techniques we have used here have been
employed for the first time in a multi-loop calculation, and that these
should be useful not only for tackling DIS at four loops but also for
other multi-loop calculations.
Using these techniques we have been able to generate a very large number
of Mellin moments, up to $N = 1500$, by evaluating the series
expansion of the forward scattering amplitude around
$1/x \equiv \omega = 0$, recall the discussion above eq. (<ref>).
With that many moments, it is possible to reconstruct the analytic
$N$-dependence of the fourth-order coefficient functions by a direct
(and over-constrained) Gaussian elimination for a sufficiently
general ansatz in terms of harmonic sums.

Our basic setup relies on Qgraf [36] and
Form [37, 38, 39],
employing the program Minos [40] as a diagram database
Many of the Form libraries we use have been employed in a
substantial number of earlier calculations, e.g., in
refs. <cit.> and have been highly optimized for multi-loop perturbative
QCD calculations.
In a particular, as in refs. [15, 23], the
database combines diagrams whose underlying graph topology is equivalent
and whose colour factors are the same.  Such sets of diagrams lead to
faster evaluation times as they allow to realize algebraic cancellations
between the individual diagrams.

As the $\nf$ contributions at three loops, the present $\nfs$ contributions
at four loops are special since they do not yet involve the more difficult
topologies in their respective orders. In fact, the most difficult
four-loop cases derive from the hardest three-loop diagrams shown in
fig. <ref> by simply inserting an additional quark loop into one of
the gluon propagators, i.e., no `genuine' (non-insertion) four-loop
self-energy topologies are required [15].

The hardest diagrams contributing to the
$\nf$ contributions to the third-order coefficient functions for
$F_{2,\rm ns}$ and $F_{L,\rm ns}$. In the notation of mincer
[41, 42], the two on the left are of the BE
topology and the two on the right are O4 cases. All four diagrams have
the colour factor $C_{\:\!\!F}\, C_{\!A}\, \nf$.

The best route to determine the all-$N$ expressions for the
$\nfs$ four-loop contributions is via the large-$\nc$ limits and the
$\cfs \nfs$ terms. The former do not involve alternating harmonic sums,
as even and odd $N$-values must lead to the same $x$-space function.
For the latter one expects a somewhat simpler form than for the
$\cf\ca \nfs$ terms, since only two-loop diagrams with two one-loop or
one two-loop self-energy insertion(s) contribute.
In practice, the non-$\zeta$ part of $c_{2,\rm ns}^{\,(4)}$ was the most
difficult case, which we solved in Form using an ansatz with a
little less than 600 even-$N$ coefficients.
All other cases required fewer than 400 coefficients.

§.§ Topologies for the non-singlet $\nfs$ structure functions

As discussed at the end of section 2, we also need the lower-order
corrections to a sufficient power in $\ep$; in particular we have to
compute the $\ep^{\,1}$ terms for the $\nf$ parts at three loops.

The diagrams contributing to the final result are organized into
topologies, and described in terms of a linear independent set of
propagators and scalar products ready to be reduced to a smaller set
of master integrals [43].
In order to obtain a more efficient reduction, especially at four loops,
we build our topologies by keeping the number of propagators within a
topology as low as possible.
In practice, a four-loop DIS topology has 18 linearly independent scalar
products; instead of describing them in terms of as many linearly
independent propagators, we opt for at most 12 propagators and 6 scalar
products that will appear only in the numerator. At three-loop, this results
in a topology described by 10 propagators and 2 scalar products.

We further simplify the expressions of the diagrams by rewriting
multi-loop self-energy corrections in terms of chained bubbles which has
the additional benefit of keeping all the propagator powers as whole
numbers more suited for the Laporta reductions,

q SelfEnergy2q   .
In fig. <ref> we show the range of topologies that are
required for the  non-singlet $\nfs$ correction at the fourth order,
with diagrams ranging from 12 to 9 distinct scalar propagators, one
propagator short from the most complicated case one can encounter
at four loops.

in 1,...,6[width=0.1]Topologies3L

in 7,16,18,22[width=0.1]Topologies3L
in 1,...,2[width=0.1]TopologiesNf2

in 3,...,8[width=0.1]TopologiesNf2

in 9,...,13[width=0.1]TopologiesNf2

Three- and four-loop propagator topologies used for the reduction to master
integrals and the expansion in $\omega$ using differential equations.
The external double line represent the off-shell photon while the simple
lines the parton.
Each diagram allows for two independent flows of the external momenta
resulting in the same topology up to a sign inversion $\omega \leftrightarrow

The number of scalar propagators in the definition of each topology is not
the only factor to take into account when considering the problem of reducing
our integrals to master integrals (very important is also, for example,
the dimensionality of the numerator).
However, it remains a key aspect since it drastically reduces the degrees of
freedom one has to consider during the reductions by effectively vetoing some
of the sectors. To perform the reductions to master integrals we have used
both the fire [44] and kira
[45] programs.

§.§ Series expansion and differential equations

Once all the diagrams appearing in the process are reduced to master
integrals, we build a system of differential equations by acting on them
with the differential operator $\partial / \partial\:\! \omega$
and further reduce the resulting integrals.
This step may also involve extending the initial set of master integrals to
achieve a closed system.
The resulting differential equations can be written as

ωM⃗ (ω,)  =
A(ω,)·M⃗(ω,) ,

where $A\in\mathbb{Q}^{\:\!n\times n}(\omega,\ep)$ is a square matrix whose
elements are rational functions  in $\omega$ and the dimensional regulator
$\ep$ over the field $\mathbb{Q}$, with $n$ being the number of master
integrals. The general ansatz for this system of differential equations can
be written as a linear combination of meromorphic functions,

M⃗(ω,)  =   ω^  α⃗
∑_s ⊂ ℤ
∑_k = 0^∞ m⃗^(s)_k(ω,)
ω^ k+s·  ,
α⃗∈ℤ^n  ,

where $\vec{\alpha}$ defines the leading power of each master and the sum over
$s$ goes over a finite set of integers which multiply $\epsilon$ in the power
of $\omega$, and represent independent (not necessarily physical) solutions to
the differential equations. The space of solutions is reduced by noticing that
the boundary conditions for the master integrals in the limit $\omega=0$ are
finite ($\alpha \geq 0$) and support only the regular sector $s=0$.
This may in fact not be true for individual integrals, but it is well known
that the forward amplitudes for DIS are regular at $p=0$ and hence $w=0$.
We can thus safely ignore solutions with $s \neq 0$.
The general problem of obtaining series expansions from differential
equations is of course not new in the context of Feynman integrals and
indeed has been employed successfully in, e.g., refs.
and public implementations exist
[54, 55, 56].

Here we present a new method which takes advantage of the specific analytic
properties of the problem and which is well suited for obtaining many

coefficients in the expansion, while at the same time
being simple enough to be implemented purely in Form. We describe
this procedure in the following.

We make use of the information coming from the boundary conditions to simplify
the ansatz for the master integral expansion, and at the same time we consider
the differential matrix to contain at most simple poles at $\omega=0$, viz

M_i  =  ω^ α_i ∑_k=0^∞m_ik  ω^k
,   A  =  A_-1/w + ∑_k=0^∞A_k  ω^k  .

In general, the matrix A can contain higher order poles around
$\omega \!=\! 0$, depending on the choice of master integrals.
While it is known in the Mathematics literature since the 1950s that a basis
transformation to remove higher poles exists [59], finding the
transformation to bring the differential equations into a so-called canonical
form is non-trivial in practice.
By now there exist nevertheless several strategies/algorithms to solve this
problem [57, 58], and some implementations are now
publicly available [60, 61, 62].
The feasibility to run these algorithms on systems containing hundreds master
integrals is however questionable.
For our purpose, we do not actually require the complete reduction of the
system in to the canonical form – as it is sufficient to simply remove higher
order poles at $\omega=0$.

We achieve this by applying a rescaling transformation $T$ to the master

M  →  T·M   ,  A  →  Tω T^-1 + T ·A ·T^ -1,

with the transformation matrix taking the form

T  =  (ω^β )  ,    β∈ℕ_0^n
⇒  A_ij  →  δ_ij β_i/ω
+ ω^ β_i-β_j A_ij   .

In general we found that the form of $T$ required to remove all poles is not
We explored several algorithms allowing to construct a suitable $T$;
the probably simplest was based on constructing $T$ by iteratively removing
a pole of order $k$ in the $j^{\:\!\rm th}$ row by setting
While this particular pole would be removed, new poles could be created in
other rows in the transformed matrix $A$.
In all cases we found that iterating this procedure allowed to eventually
remove all non-simple poles from $A$.
We actually have no prove for this simple procedure to terminate and it
remains to be seen whether it will also work for even more complicated cases.
The advantage of the procedure in comparison to the much more involved
algorithms to bring $A$ into a Fuchsian form is that it is computationally
far simpler and can be applied with

ease also to comparably large systems of sizes around the 100s or even 1000s.
It is not clear to us that the same holds for the algorithms mentioned above.

By using the definitions in eq. (<ref>) it is possible
to rewrite eq. (<ref>) into a recursive expression that allows for an
efficient extraction of the expansion coefficients of the master integrals,

((k+1)𝕀- A_-1)_:= B_k· m⃗_k+1^
=  ∑_j=0^k A_j m⃗_k-j^

where $\id$ is the identity matrix.
Note that we can have at most $n$ cases where $\det(B_k)=0$, corresponding
to positive integer values of the eigenvalue of $A_{-1}$.
In these cases the system cannot be inverted and is solved by means of
Gaussian elimination.
For high enough $k$ we have $\det(B_k)\neq 0$ and we can write

m⃗_k+1^  =  B_k^ -1·(
∑_j=0^k A_j  m⃗_k-j^)   .

The inversion of the matrix $B$ is performed for generic values of $k$,
avoiding the expensive procedure of performing a new inversion for every
step of the expansion, especially for high values of $k$.

We have implemented the Gaussian elimination for the arbitrary steps of
eq. (<ref>) in C while the case of non-singular
$B_k$ matrices in eq. (<ref>) has been implemented in
The required boundary conditions have been computed using Forcer
in the limit of $\omega=0$ for all the master integrals,
where the external parton is taken to be soft.
The performances of eq. (<ref>) can be further improved
by truncating the expansion in the dimensional regulator to the required order
to obtain an amplitude known up to $\ep^{4-n}$ at the $n$-th order.

One can benefit from the combined expansion only when the two limits
$\ep,\omega=0$ of the differential matrix $A$ commute.
For example, the series expansion in $\omega$ of a denominator of the form
$(\omega-\ep)$ has a convergence radius that goes to zero as $\ep$ vanishes.
These kind of poles are unphysical and will lead to arbitrary high poles in
$\ep$ when expanding the differential matrix and will cancel in the final
expression of the expansion coefficients $\vec{m}_{k+1}$.

To circumvent this problem one can transform to a basis of master integrals
in which the reduction coefficients – and therefore also the coefficients
of the differential equation matrix $A$ – have the property that their
dependence on $\epsilon$ factorizes from their dependence on $\omega$.
The existence of such a basis was proven in ref. [63],
and two independent implementations which construct the corresponding basis
of master integrals have been published [64, 65]
as complementary codes to both the fire and kira reduction
programs respectively. We used the former one in our
reductions with fire.

With this implementation, and after obtaining a factorized form for the
differential equations, we were able to push the recursive algorithm in
eq. (<ref>) up to $\order{\omega^{\:\!1500}}$ within
no more than a few days for all of the masters.
Such a large number of coefficients was necessary to be certain to fully
constrain the expression of the coefficient functions in $N$-space from an
ansatz in terms of harmonic sums and rational coefficients.

§.§ Four-loop rescaling example

Here we expand on the method used to bring the differential system into
the form with at most simple poles in $\omega$.
We will illustrate this method with one of the four-loop topologies we
encountered in our computation.

\begin{equation}
\begin{gathered}
\includegraphics[width=.1\textwidth]{./img/TopologiesNf2_1.eps}
\end{gathered}:\qquad
\Itop{\nu_1^{},\nu_2^{},\nu_3^{},\dots,\nu_{18}^{}} \,:=\,
\int\left(\prod_{j=1}^4\frac{\dd^d k_j}{(2\pi)^d}\right)
\frac{1}{\prod_{n=1}^{\:\!18}D_n^{\:\!\nu_n}}
\end{equation}

\begin{align}
\nonumber
&&D_{13}=2 k_1\cdot k_4,\quad
\\&
&&D_{14}=2 k_2\cdot k_4\,,\quad
\nonumber
\\&
&&D_{15}=2 k_1\cdot p\,,\quad
\nonumber
\\&
&&D_{16}=2 k_4\cdot p\,,\quad
\nonumber
\\&
&&D_{17}=2 k_1\cdot q\,,\quad
\nonumber
\\&
&&D_{12}=\left(k_3+k_4+p-q \right){}^2,\quad
&&D_{18}=2 k_4\cdot q\,.\quad
\end{align}

The reduction of this topology produced a set of $55$ distinct
master integrals:

\begin{align*}
\{\;
\Itop{0,0,1,1,0,1,1,0,1,1,0,1,0,0,0,0,0,0},\;
\Itop{0,0,1,1,1,1,1,0,1,1,0,1,0,0,0,0,0,0},\;
\\&
\Itop{0,1,0,1,0,1,1,0,0,1,0,1,0,0,0,0,0,0},\;
\Itop{0,1,0,1,0,1,1,0,1,1,0,1,0,0,0,0,0,0},\;
\Itop{0,1,0,1,1,1,1,0,0,1,0,1,0,0,0,0,0,0},\;
\\&
\Itop{0,1,0,1,1,1,1,0,0,2,0,1,0,0,0,0,0,0},\;
\Itop{0,1,1,1,0,1,1,0,1,1,0,1,0,0,0,0,0,0},\;
\Itop{1,0,0,1,0,1,1,0,0,0,0,1,0,0,0,0,0,0},\;
\\&
\Itop{1,0,0,1,0,1,1,0,1,0,0,1,0,0,0,0,0,0},\;
\Itop{1,0,0,1,1,1,1,0,0,0,0,1,0,0,0,0,0,0},\;
\Itop{1,0,0,1,1,1,1,0,0,0,1,1,0,0,0,0,0,0},\;
\\&
\Itop{1,0,0,1,1,1,1,0,0,1,0,1,0,0,0,0,0,0},\;
\Itop{1,0,0,1,1,1,1,0,1,1,0,1,0,0,0,0,0,0},\;
\Itop{1,0,1,1,0,0,1,0,0,1,1,1,0,0,0,0,0,0},\;
\\&
\Itop{1,0,1,1,0,0,1,0,1,1,0,1,0,0,0,0,0,0},\;
\Itop{1,0,1,1,0,1,0,0,1,0,0,1,0,0,0,0,0,0},\;
\Itop{1,0,1,1,0,1,0,0,1,0,0,2,0,0,0,0,0,0},\;
\\&
\Itop{1,0,1,1,0,1,0,0,1,1,0,1,0,0,0,0,0,0},\;
\Itop{1,0,1,1,0,1,1,0,0,0,1,1,0,0,0,0,0,0},\;
\Itop{1,0,1,1,0,1,1,0,0,1,0,1,0,0,0,0,0,0},\;
\\&
\Itop{1,0,1,1,0,1,1,0,1,0,0,1,0,0,0,0,0,0},\;
\Itop{1,0,1,1,0,1,1,0,1,1,0,1,0,0,0,0,0,0},\;
\Itop{1,0,1,1,0,1,1,0,1,1,0,2,0,0,0,0,0,0},\;
\\&
\Itop{1,0,1,1,0,1,1,0,1,2,0,1,0,0,0,0,0,0},\;
\Itop{1,0,1,1,0,1,2,0,1,0,0,1,0,0,0,0,0,0},\;
\Itop{1,0,1,1,0,1,2,0,1,1,0,1,0,0,0,0,0,0},\;
\\&
\Itop{1,0,1,1,1,0,1,0,0,1,1,1,0,0,0,0,0,0},\;
\Itop{1,0,1,1,1,0,1,0,1,1,0,1,0,0,0,0,0,0},\;
\Itop{1,0,1,1,1,1,0,0,1,1,0,1,0,0,0,0,0,0},\;
\\&
\Itop{1,0,1,1,1,1,1,0,0,0,1,1,0,0,0,0,0,0},\;
\Itop{1,0,1,1,1,1,1,0,0,0,1,2,0,0,0,0,0,0},\;
\Itop{1,0,1,1,1,1,1,0,0,1,0,1,0,0,0,0,0,0},\;
\\&
\Itop{1,0,1,1,1,1,1,0,0,1,0,2,0,0,0,0,0,0},\;
\Itop{1,0,1,1,1,1,1,0,0,1,1,1,0,0,0,0,0,0},\;
\Itop{1,0,1,1,1,1,1,0,0,1,1,2,0,0,0,0,0,0},\;
\\&
\Itop{1,0,1,1,1,1,1,0,1,0,0,1,0,0,0,0,0,0},\;
\Itop{1,0,1,1,1,1,1,0,1,1,0,1,0,0,0,0,0,0},\;
\Itop{1,0,1,1,1,1,1,0,1,1,0,2,0,0,0,0,0,0},\;
\\&
\Itop{1,0,1,1,1,1,1,0,1,2,0,1,0,0,0,0,0,0},\;
\Itop{1,0,1,1,1,1,1,0,2,1,0,1,0,0,0,0,0,0},\;
\Itop{1,1,0,1,0,0,1,0,0,1,0,1,0,0,0,0,0,0},\;
\\&
\Itop{1,1,0,1,0,0,1,0,1,1,0,1,0,0,0,0,0,0},\;
\Itop{1,1,0,1,1,0,1,0,0,1,0,1,0,0,0,0,0,0},\;
\Itop{1,1,0,1,1,0,1,0,0,1,1,1,0,0,0,0,0,0},\;
\\&
\Itop{1,1,0,1,1,1,1,0,0,1,0,1,0,0,0,0,0,0},\;
\Itop{1,1,0,1,1,1,1,0,0,1,0,2,0,0,0,0,0,0},\;
\Itop{1,1,1,1,0,0,0,0,0,1,1,1,0,0,0,0,0,0},\;
\\&
\Itop{1,1,1,1,0,0,0,0,1,1,0,1,0,0,0,0,0,0},\;
\Itop{1,1,1,1,0,0,0,0,1,1,0,2,0,0,0,0,0,0},\;
\Itop{1,1,1,1,0,0,1,0,1,1,0,1,0,0,0,0,0,0},\;
\\&
\Itop{1,1,1,1,0,1,0,0,1,1,0,1,0,0,0,0,0,0},\;
\Itop{1,1,1,1,0,1,0,0,1,1,0,2,0,0,0,0,0,0},\;
\Itop{1,1,1,1,1,0,0,0,0,1,1,1,0,0,0,0,0,0},\;
\\&
\Itop{1,1,1,1,1,0,1,0,0,1,1,1,0,0,0,0,0,0}\;\}.
\end{align*}

For the purpose of the rescaling, the only information we are concerned about
is the leading behaviour in $\omega$ of each of the coefficients of the
matrix $A$ appearing in the differential euqations.
We represent the leading behaviour for the topology of eq. (<ref>) using
the following colour coding

\begin{equation*}
\begin{gathered}\tikz{\draw[fill=yellow!0!red ,draw=none]  rectangle ++(0.5em,0.5em);}\end{gathered},
\begin{gathered}\tikz{\draw[fill=yellow!25!red,draw=none]  rectangle ++(0.5em,0.5em);}\end{gathered},
\begin{gathered}\tikz{\draw[fill=yellow!50!red,draw=none]  rectangle ++(0.5em,0.5em);}\end{gathered},
\begin{gathered}\tikz{\draw[fill=yellow!75!red,draw=none]  rectangle ++(0.5em,0.5em);}\end{gathered},
\begin{gathered}\tikz{\draw[fill=green!90!blue,draw=none]  rectangle ++(0.5em,0.5em);}\end{gathered},
\begin{gathered}\tikz{\draw[fill=green!60!blue,draw=none]  rectangle ++(0.5em,0.5em);}\end{gathered},
\begin{gathered}\tikz{\draw[fill=green!30!blue,draw=none]  rectangle ++(0.5em,0.5em);}\end{gathered}
\} \; .
\end{equation*}

We also use a circle, instead of a square, to highlight the position of
the deepest poles at each step of the transformation.
The elimination is performed row-wise, where we rescale each master whose
row contains at least one of the deepest poles by a factor of $\omega$.
In the current example the iterative rescaling procedure leads to the following sequence of transformations:

\begin{align*}
&\phantom{ \xrightarrow{\;T_1\;\;\;} }
\left(\begin{gathered}
\includegraphics[width=\matrixplotwidth]{./img/RescaleDE_1.eps}
\end{gathered}\right)
\\
\left(\begin{gathered}
\includegraphics[width=\matrixplotwidth]{./img/RescaleDE_2.eps}
\end{gathered}\right)
\\
\left(\begin{gathered}
\includegraphics[width=\matrixplotwidth]{./img/RescaleDE_3.eps}
\end{gathered}\right)
\\
\left(\begin{gathered}
\includegraphics[width=\matrixplotwidth]{./img/RescaleDE_4.eps}
\end{gathered}\right)
\\
& \xrightarrow{T_4\;\;}
\left(\begin{gathered}
\includegraphics[width=\matrixplotwidth]{./img/RescaleDE_5.eps}
\end{gathered}\right)
\end{align*}

where the transformation map is defined as

\begin{equation}
A \;\xrightarrow{\;T_a\;}\; T_a\cdot\left(A \cdot T_a^{\:\!-1}
- \pdv{\:\! T_a^{-1}}{\:\!\omega}\right)
\end{equation}

\begin{align}
\left(T_{a}\right)_{ij} \;=\;\delta_{ij}\times
\begin{cases}
\;\omega,\quad a=1,\;i\in\{39,40,41\}\\
\;\omega,\quad a=2,\;i\in\{25,39,40,41\}\\
\;\omega,\quad a=3,\;i\in\{23,24,25,27,32,34,\dots,41,47,52,53\}\\
\;\omega,\quad a=4,\;i\in\{38\}\\
\;1 \; .
\end{cases}
\end{align}

Here we use $\dots$ as a short-hand for sequence of consecutive integers.
The complete transformation for the matrix $A$ can then be build by combining
the intermediate steps

\begin{align}
\nonumber
T= T_4\cdot T_3\cdot T_2\cdot T_1= \text{diag}\Big(&1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, \omega, \omega,
\\&
\nonumber
\omega^2, 1, \omega, 1, 1, 1, 1, \omega, 1, \omega, \omega, \omega, \omega,
\omega^2, \omega^3, \omega^3, \omega^3, 1, 1, 1,
\\&
1, 1, \omega, 1, 1, 1, 1, \omega, \omega, 1, 1\Big) \; .
\end{align}

Alternatively, one could decide remove the deepest pole by rescaling
column-wise, where each master, whose column contains one of the deepest
poles, is rescaled by a factor of $\omega^{\:\!-1}$ or, equivalently,
all the others are rescaled by $\omega$.
We use the latter because it does not affect the ability to express each
master integral as a Taylor expansion.

If we were to start with the same matrix but cancel all the poles column-wise,
the resulting transformation matrix would have been

\begin{align}
\left(T_{a}\right)_{ij} \;=\; \delta_{ij}\times
\begin{cases}
\;\omega,\quad a=1,\;i\notin\{9,10\}\\
\;\omega,\quad a=2,\;i\notin\{39,40,41\}\\
\;\omega,\quad a=3,\;i\notin\{1,\dots,6,9,\dots,19,21,22,23,24,26,\dots,
\;1
\end{cases}
\end{align}

and $T = T_3\cdot T_2 \cdot T_1$ with

\begin{align}
\nonumber
T = \text{diag}\Big(&\omega^2, \omega^2, \omega^2, \omega^2, \omega^2,
\omega^2, \omega^3, \omega^3, 1, 1, \omega^2, \omega^2, \omega^2, \omega^2,
\omega^2, \omega^2, \omega, \omega, \omega^2,
\\&
\nonumber
\omega^3, \omega^2, \omega^2, \omega^2, \omega^2, \omega^3, \omega^2,
\omega^2, \omega^2, \omega^2, \omega^2, \omega^3, \omega^3, \omega^3,
\omega^3, \omega^3, \omega^3, \omega^2,
\\&
\omega^3, \omega^3, \omega^3, \omega^3, \omega^2, \omega^3, \omega^2,
\omega^3, \omega^3, \omega^3, \omega^3, \omega^3, \omega^3, \omega^3,
\omega^3, \omega^3, \omega^3, \omega^3\Big) \; .
\end{align}

In general, one could combine the two approaches by alternating between
column- and row-wise rescaling at any step of the procedure.
Such a combination will generate different rescaling matrices $T$.
We have observed that a column-wise approach is generally preferred for our
computation because it requires the rescaling of fewer integrals, resulting in
the ability to truncate the expansion for a larger set of master integrals.

§ RESULTS IN N-SPACE

We are now ready to present our new analytic even-$N$ expressions for
the $\nfs$ contributions to the fourth-order coefficient functions
$c_{2,\rm ns}^{\,(4)}$ and $c_{L,\rm ns}^{\,(4)}$ in eq. (<ref>).
The corresponding $\nft$ results were derived long ago
[19, 20]. We confirm those results and
include them below for completeness; for the case of ${\cal C}_2^{}$
in a more transparent form than given in ref. [20].

The quantities under consideration can be expressed in terms of
harmonic sums. Following the notation of ref. [21],
these sums are recursively defined by

𝐒_±m(M)  =
∑_i = 1^M  (±1)^i/i^  m

𝐒_±m_1^, m_2^,…, m_k^(M)  =
∑_i = 1^M  (±1)^i/i^  m_1^
S_m_2^,…, m_k^(i)    .

The sum of the absolute values of the indices $m_k$ defines the weight $w$
of the harmonic sum. In the $n$-loop coefficient functions one encounters
sums up to weights $2n$ for ${\cal C}_2$ and $2n-1$ for ${\cal C}_L$.
The present non-singlet $\nfs$ contributions only include sums op to
$w = 6$ for ${\cal C}_2$ and $w = 5$ for ${\cal C}_L^{}$, for the $\nft$
terms the corresponding maximal weights are lower by 1.
Below all harmonic sums have the argument $N$, which is omitted in the
formulae for brevity, and we use the short-hand

D_a  =   (N+a)^-1    .

We first present the result for the coefficient function for $F_L$
which we write as

c_L,ns^ (4)(N)     =

+      16/9
c_L,ns^ (4)L(N)
+   (-12  )     16/9
c_L,ns^ (4)N(N)

+    16/27
c_L,ns^ (4)F(N)

where $\ca = \nc$ and $\cf = (\nc - n_c^{-1})/2$ in SU($\nc$), with
$\nc = 3$ colours in QCD.
For compactness, we have decomposed the $\nfs$ part into leading (L)
and non-leading (N) contributions in the large-$\nc$ limit. The factors
16/9 and 16/27 have been put in order to shorten, on average, the lengths
of the fractions in the expressions below.
The $\nfs$ contributions are given by

The $\cfs$ contribution (<ref>) includes, as the corresponding
lower-order quantities, only the denominator $D_1 = 1/(N\!+\!1)$ at the
maximal overall weight $w = 5$ of the sums and Riemann $\zeta$-values.
The $\cf\ca$ part, and hence the large-$\nc$ coefficient (<ref>),
does not have this expected property, but instead involves the linear

2 D_{-2} - D_{-1} + D_1 - 2 D_2 + 3 D_3
\quad \mbox{ and } \quad
3 D_1 - 12 D_2 + 10 D_3
\:\: .

The presence of terms with $D_{-2} = 1/(N\!-2)$ and $D_{-1} = 1/(N\!-1)$
does not imply poles at $N=2$ or $N=1$, as the corresponding numerators also
vanish at this point. This feature already occurred in the second-order
coefficient functions of
refs. <cit.>. At $N\!=\!2$ these functions are given by

c_L,ns^ (4)L(N=2)     =
- 6713/135  3
- 32 5
+ 24    ,

c_L,ns^ (4)N(N=2)     =
+ 247/270 3
+ 30 5

Note the $\zts$ term which does not occur at higher $N$.
The $\nft$ contribution to eq. (<ref>) reads

The corresponding result for ${\cal C}_2$ can be decomposed as

c_2,ns^ (4)(N)     =

+      16/9
c_2,ns^ (4)L(N)
+   (-12  )     16/9
c_2,ns^ (4)N(N)

+   (-)     1/3   4
c_2,ns^ (4)Z(N)
+    16/27
c_2,ns^ (4)F(N)

In addition to the structures present in eq. (<ref>), the $\nfs$
coefficient function for $F_2$ includes a $\z4$ contribution,
which is proportional to $(\cf\! - \ca)$ and hence vanishes vanish for
$ \ca = \cf$  which is part of the choice of the colour factors that leads
to a ${\cal N}\! =\! 1$ supersymmetric theory.

The three $\nfs$ coefficients in eq. (<ref>) read

Unlike the above expressions for ${\cal C}_L$, eqs. (<ref>) and
(<ref>) include harmonic sums (up to $w=6$) without prefactors
$D_a$; these are the terms that do not vanish in the soft-gluon limit
$N \!\ra \infty$. As in eq. (<ref>), only non-alternating sums
occur in the large-$\nc$ result (<ref>), and additional denominators
(here in addition to $D_0$ and $D_1$) occur with $w=5$ sums in the form

D_{-2} -  6D_2 + 9 D_3 \:\: .

Again as for ${\cal C}_L$, terms with $D_a$, $a \neq 0,\,1$, are not
present with $w=5$ sums in the $\cfs\, \nfs$ coefficient (<ref>).

The situation at $N=2$ is the same as that for ${\cal C}_L$ above, with

c_2,ns^ (4)L(N=2)     =
- 4613/540  3
- 290/3  5
+ 6    ,

c_2,ns^ (4)N(N=2)     =
+ 247/270  3
+ 30 5

Finally the $\nft$ coefficient in eq. (<ref>) is given by

The numerical size of the above fourth-order results is shown in
figs. <ref> and <ref> for QCD, i.e., $\ca = 3$ and
$\cf = 4/3$, with $\nf = 4$ light flavours, together with the
corresponding third-order contributions for the physically very
wide range $2 \leq N \leq 50$.
The coefficient functions $c_{L,\rm ns}^{\,(n)}$ vanish for
$N \!\ra \infty$, yet only slowly: the size of the $\nfs$ parts of
$c_{L,\rm ns}^{\,(4)}$ in the right part of fig. <ref>
decreases only by a factor of 2 from $N \!=\! 10$ to $N \!=\! 50$.
A further reduction by another factor of 2 and 4 is only reached
at $N \!=\! 175$ and $N = 540$, respectively.

The shape of the leading large-$\nf$ contributions
($ \sim n_{\! f}^{\,n-1}$ at order $\as(n)$) in fig. <ref> is
similar to that of the subleading large-$\nf$ contributions
($\sim n_{\! f}^{\,n-2}$ at order $\as(n)$) in the $N$-range of the
figure; its relative size is decreasing for $\nf = 4$ from about
1/10 at $n=3$ to about 1/15 at $n=4$.
This pattern of a reduced numerical significance of the leading
large-$\nf$ term at the fourth order is also seen for
$c_{2,\rm ns}^{\,(n)}$ in fig. <ref>. Here the subleading
large-$\nf$ contributions are larger by factors between 10 and 15
at $4 \leq N \leq 50$ for $n=3$; the corresponding range for $n=4$
is 19 to 24.

The leading (dashed) and sub-leading (solid) large-$\nf$ contributions
to the three-loop (left panel) and four-loop (right panel) coefficient
functions for the structure function $F_{L,\rm ns}$ for QCD with
$\nf = 4$ flavours. The results in eqs. (<ref>) – (<ref>)
have been converted to an expansion in $\als$, and the leading
large-$\nf$ curves have been scaled up by a factor of 10 for better

As fig. <ref>, but for the structure function $F_{2,\rm ns}$.
In addition the $\ln^\ell \!\wN$ and $\ln^\ell \! \wN \!+\!
\mbox{const.}$ threshold contributions are shown here, respectively,
by the upper and lower dotted curves.

Also shown in fig. <ref> are the dominant contributions to
${\cal C}_{2,\rm ns}$ in the large-$N$ threshold limit, where our
new contributions (<ref>) – (<ref>) to
$c_{2,\rm ns}^{\,(4)}$ lead to the numerical expansion

c_2,ns^ (4)(N) |_    =
0.7233196  ln^  6 + 12.339095  ln^  5

+ 87.721224  ln^  4 + 293.04552  ln^  3

+ 233.48456  ln^  2 + 65.035706  ln- 12175.412

+ N^ -1 (
2.1069959  ln^  5 + 67.193416  ln^  4 .
+ 691.16782  ln^  3 + 3429.9787  ln^  2 .
+ 8755.5832  ln+ 12282.167 )
+ O ( N^-2 )

with $\ln \wN = \ln N + \gamma_e$, where $\gamma_e$ is the
Euler-Mascheroni constant.
Keeping only the $\ln^{\:\!\ell} \!\wN$ terms in the first three lines,
and the corresponding contributions at the third order, one arrives
at the upper dotted curves in the figure. Adding to these the constant-$N$
contributions yields the lower dotted curves. Neither of these results
can quantitatively replace the exact expressions at physically interesting
moderate values of $N$.

Analytic $x$-space expression corresponding to eq. (<ref>) and
a further discussion of the threshold and high-energy limits can be found
in section 4 below.

A five-loop prediction

It is possible to predict a small part, the $\z4$ coefficients of $\nft$
and $\nff$, of the $n\!=\!4$ five-loop non-singlet anomalous dimension
in eq. (<ref>), from our result (<ref>) and eq. (<ref>)
for $c_{2,\rm ns}^{\,(4)}(N)$. This possibility arises from the
no$\:\!$-$\:\!\pi^2$ conjecture$/$theorem [66, 67]
for Euclidean (space-like) physical quantities in a suitable renormalization
scheme, which was investigated in the context of inclusive DIS in
ref. [68].
For this one considers the physical evolution kernels which were
expressed in terms of the splitting functions and coefficient functions
for the non-singlet case to the fifth order in eqs. (2.7) – (2.9) of
ref. [69].

Keeping only the terms with $\z4$ simplifies the fifth-order  kernel
for $F_{2,\rm ns}$ to

_2 ns^ (4)(N)  =
- _ ns^ (4)(N)
- 3  β_1  _a^ (3)(N)
- 4  β_0 ( _2 ns^ (4)(N)
- c_2 ns^ (1)(N)  _2 ns^ (3)(N) )

where the tilde indicates the contribution with $\z4\!=\! 1/90\:\pi^{\,4}$.
At this order, a scheme transformation removing the $\z4$ term of the
five-loop beta function [70, 71, 72]
needs to be performed, and the prediction of the no$\:\!$-$\:\!\pi^2$
theorem becomes [68]

β_0  _a^ (4)
+ 13  β_ 4   K_a^ (0)  =  0   .

Since $\beta_0$ includes $\nf$ and $K_{2,\rm ns}^{\,(0)} \sim \cf$,
only the $\zeta_4 \,\nff$ terms in $\beta_4$ can contribute, but
there is no such term. Hence already the quantity (<ref>)
has to vanish for the $\nft$ and $\nff$ terms.

The three-loop coefficient function $c_{2,\rm ns}^{\,(3)}(N)$ does not
include a $\z4$ term with $\nfs$. Since the prefactors $\beta_1$ and
$\beta_0 \, c_{2,\rm ns}^{\,(1)}$ include no more than one power in
$\nf$, the two terms with $c_{2,\rm ns}^{\,(3)}$ do not contribute to
the $\nft$ and $\nff$ parts of eq. (<ref>).

So we end up with a simple relation for the $\z4\,\nft$ and
$\z4\,\nff$ contributions,

_ns^ (4)(N)
=  -  _ns^ (4)(N)
=  4  β_0  _2 ns^ (4)(N)

which leads to

At $N \!=\!2 $ and $N \!=\!3 $ this result agrees with those obtained
by diagram calculations in ref. [73] using the program
of ref. [74] for the $R^{\:\!\ast}$ operation.
Its last line was derived long ago as part of the complete leading
large-$\nf$ result [75].
Due to eqs. (<ref>) and (<ref>), the whole of (<ref>)
is proportional to $\beta_0$ and the lowest-order anomalous dimension
$\gamma_{\rm qq}^{\,(0)}(N)$ for $\cf = \ca$.

§ THE X-SPACE COEFFICIENT FUNCTIONS

The coefficient functions $c_{a,\rm ns}^{\,(n)}(x)$ are obtained from
the $N$-space results of the previous section by an inverse Mellin
transformation, which expresses these functions in terms of harmonic
polylogarithms (HPLs) $H_{m_1,...,m_w}(x)$, $m_j = 0,\pm 1$. Following
ref. [22], to which the reader is referred for a
detailed discussion, the lowest-weight ($w = 1$) functions $H_m(x)$ are
given by

H_0(x)         =   lnx    ,     H_±1(x)   =   ∓  ln(1 ∓x)    ,

and the higher-weight ($w \geq 2$) functions are recursively defined as

H_m_1,...,m_w(x)   =
1/w! ln^w x    ,
if     m^_1,...,m^_w = 0,…,0

∫_0^x  dz  f_m_1(z)   H_m_2,...,m_w(z)
,        else

f_0(x)         =   1/x    ,     f_±1(x)   =   1/1 ∓x    .

The inverse Mellin transformation exploits an isomorphism between
the set of harmonic sums for even or odd $N$ and the set of HPLs.
Hence it can be performed by a completely algebraic procedure
[12, 22], based on the fact that harmonic
sums occur as coefficients of the Taylor expansion of harmonic
A Fortran program for the HPLs up to weight $w=4$ has been
provided in ref. [76], later this was extended
to $w=5$ and $w=6$.

In our results below, the argument $x$ of the HPLs is suppressed
for brevity, and we use the abbreviations

=   1-x        =   1+x    .

Analogous to eq. (<ref>) above, the $x$-space coefficient
function for $F_{L,\rm ns}$ is decomposed as

c_L,ns^ (4)(x)     =

+      16/9
c_L,ns^ (4)L(x)
+   (-12  )     16/9
c_L,ns^ (4)N(x)

+    16/27
c_L,ns^ (4)F(x)
The two $\nfs$ contributions in the second line are given by

The $\nft$ coefficient in the last line of eq. (<ref>) reads

The corresponding results for $F_{2,\rm ns}$ are written in a very
similar manner as

c_2,ns^ (4)(x)     =

+      49
c_2,ns^ (4)L(x)
+   (-12  )     49
c_2,ns^ (4)N(x)

+    427
c_2,ns^ (4)F(x)

The structure of these expressions reflects the $N$-space results
presented above. The large-$\nc$ coefficients (<ref>),
(<ref>), (<ref>) and (<ref>) include only HPLs
with non-negative indices corresponding to the non-alternating harmonic
sums in eqs. (<ref>), (<ref>), (<ref>) and
At the highest weight, $w=5$, the $\nfs$ coefficients
(<ref>) and (<ref>) do not only appear with factors $x$
for ${\cal C}_L$ and $\xm^{-1} = 1/(1-x)$, $1$ and $x$ for ${\cal C}_2$,
but also include additional fixed combinations of $x^a$ with
$a = -2,\,-1,\,2,\,3$.
These do not occur with $w=5$ HPLs in the $\cfs \nfs$ coefficients
(<ref>) and (<ref>) corresponding to eqs. (<ref>)
and (<ref>), which however include HPLs with negative indices
and terms with $\xp^{-1} = 1/(1+x)$.

Up to terms that vanish in the corresponding limit,
the large-$x$ and small-$x$ behaviour of the coefficient functions is
given by plus-distributions,

(n)  =  [ (n)/1-x ]_+

the $\delta$-function $\delta(1\!-\!x)$ and powers of the logarithms

L_1   =   ln(1-x)
L_0   =   lnx

The $\nfs$ and $\nft$ coefficients of $\Dplus(n)$ and $\delta(1\!-\!x)$ for
${\cal C}_2$ are

The coefficients of $\Dplus(n)$ have been obtained from the soft-gluon
exponentiation in eqs. (5.6) - (5.9) of ref. [77] for
$n = 2,\,\ldots,\,5$, and in eq. (A.4) of ref. [78] for
$n=1$ and $n=0$; our results in eqs. (<ref>) - (<ref>)
agree with these predictions.
On the other hand, the coefficient of $\delta(1\!-\!x)$ in
eq. (<ref>) is a new result of the present article.

The corresponding terms with powers of $\LntO$, which are subleading for
${\cal C}_2$ but leading for ${\cal C}_L$ in the threshold limit,
are given by

The coefficients (<ref>) and (<ref>) have been
predicted in ref. [79] up to one number that was
determined as $\xi_{\rm DIS_4}^{} = 100/3$ a little later
[80, 81]; eq. (<ref>)
provides the first check of this result by a four-loop diagram
Also the coefficient (<ref>) has been predicted before
[82, 79]; all other results in
eqs. (<ref>) – (<ref>) are new.

The small-$x$ limit of the non-singlet splitting functions and
coefficient function is given by the terms with
$x^{\:\!0} \ln^{\:\!\ell} x = x^{\:\!0} L_0^{\:\!\ell}$.
Their $\nfs$ and $\nft$ contributions read

Eqs. (<ref>) and (<ref>) agree with the resummation
predictions (4.5) and (4.6) of ref. [83] after using
that $\ln^{\:\!\ell\!} x$ transforms to $(-1)^\ell\, \ell\, !\,
N^{\,-\ell-1}$.  All other small-$x$ coefficients are new.

In the case of QCD, eqs. (<ref>) – (<ref>) lead
to the numerical small-$x$ expansions

c_2,ns^ (4)(x) |_    =
- 2.1410151  ln^  5  x
- 57.187471  ln^  4  x
- 358.88192  ln^  3  x

- 945.87933  ln^  2  x
- 1327.8897  lnx
- 688.88341
+ O ( x )

c_L,ns^ (4)(x) |_    =
- 20.192044  ln^  3  x
- 347.47874  ln^  2  x
- 1539.0328  lnx

- 2177.6994
+ O ( x )

which yield the successive approximations shown in fig. <ref>.
The pattern seen in this figure is the same seen for the complete
third-order non-singlet coefficient functions in the right parts
of fig. 2 and fig. 7 in ref. [13] and for the
four- and five-loop resummation predictions in figs. 1 – 3 of
ref. [83]: the leading, next-to-leading and
next-to-next-to-leading logarithmic small-$x$ approximations wildly
oscillate and thus, in general, cannot be used to obtain any reliable
prediction for even the rough shape of high-order coefficient
functions at any physically relevant small values of $x$.

Successive small-$x$ approximations for the $\nfs$ contributions to
the fourth-order coefficient functions for $F_{2,\rm ns}$ (left panel)
and $F_{L,\rm ns}$ (right panel). The coefficients in eqs. (<ref>)
and (<ref>) have been converted to an expansion in $\als$.

§ SUMMARY AND OUTLOOK

As a step towards the determination of the fourth-order QCD coefficient
functions $c_{a}^{\,(4)}(x)$ in inclusive deep-inelastic scattering
(DIS), we have computed the double-fermionic ($\nfs$) non-singlet
contributions to the structure functions $F_2$ and $F_L$. Our results
are applicable to electromagnetic DIS (where we have ignored the very
small [30, 13] contributions in which the
photon couples to different quark lines) and to charged-current
$(W^+ \!+\! W^-)$ exchange.
The analytic dependence of these coefficient functions on Mellin-$N$,
and hence Bjorken-$x$, has been reconstructed from a very large number
of even values of $N$, up to $N \simeq 1200$.

Our calculations verify the corresponding contributions to the
next-to-next-to-leading order splitting function $P_{\rm ns}^{\,+(3)}(x)$
obtained in ref. [23], and the much simpler (and much
smaller) $\cf \nft$ leading large-$\nf$ contributions to $c_{L}^{\,(4)}(x)$
and $c_{\:\!2}^{\,(4)}(x)$ [19, 20].
The coefficients of all large-$x$ plus-distributions in ${\cal C}_2$
at this order,
$[ (1\!-\!x)^{-1} \ln^{\:\!\ell} (1\!-\!x)]_+$ with $0 \leq \ell \leq 5$,
have been predicted by the soft-gluon exponentiation
[77, 78]; we agree with these predictions.
Our results for both structure functions make contact, just, with the
resummations of $\ln (1\!-\!x)$ and $\ln x$ double logarithms in refs.
[79, 80, 81, 82, 83],
as the leading terms of the $\nfs$ parts contribute to the overall
next-to-next-to-leading logarithms; also here we find full agreement.
Finally we have predicted the analytic form of the $\z4 \:\! \nft$
contributions to the Mellin-$N$ moments of the five-loop splitting function
$P_{\rm ns}^{\,+(4)}(x)$, and verified the corresponding $\z4 \:\! \nff$
parts of ref. [75].

To compute the large number of moments required for the reconstruction of the
all-$N$ coefficient functions, we have developed a new approach based on the
method of integration by parts, differential equations and the optical theorem.
By employing the knowledge that the Mellin moments of the forward scattering
amplitude correspond to the coefficients of the Taylor series around
$\omega=1/x=0$, we derive a system of recurrence relations which holds for any
linear system of differential equations, whose matrix contains no higher-order
poles in $\omega$.
While such higher poles are generically present, we find a simple algorithm
which allows to change the basis of master integrals into a form where only
simple poles appear in the differential matrix. The thereby obtained
recurrence relations for the master integrals allow to obtain the Mellin
moments at, in principle, arbitrary high $N$ from the knowledge of the
boundary conditions, which we compute using the Forcer program

In practice, we find that the algorithm starts to become computationally
demanding for the problem under consideration at values of $N$ around about
One element which allowed us to speed up the calculation considerably was to
convert the basis of master integrals into one where the $\ep$-dependence of
the reduction coefficients factorizes from the dependence on $x$.
An important feature of the method is the simplicity with which the
recurrences can be solved - a procedure which we have implemented in
Form [37, 38, 39].

In spirit the method is not dissimilar from the `method of arbitrary high
moments' [47], which was used recently to recompute the
3-loop DIS coefficient functions [14].
Their method also derives recurrence relations from differential equations,
but differs from our method in several ways.
It relies on more advanced combinatorial algorithms to construct and
recursively solve the differential equations for each master integral (while
our method is based on solving a recurrence relation for the entire system).
Furthermore it requires that the differential equations are first-order
factorizable, a property not needed in our approach.
While it is possible that their approach leads to a faster algorithm at
very large $N$, it is difficult at this stage to comment on timings without
explicit performance benchmarks.

We plan to present the $\nfs$ fourth-order contributions to the structure
function $F_3$ in $(W^+ \!+\! W^-)$ charged-current DIS in a forthcoming
publication. The computation of its `standard' $\cf \ca\:\! \nfs$ and
$\cfs\:\!\nfs$ parts is not much more difficult than that presented in
this article, but some additional contributions with the group invariant
$d^{\:\!abc}d_{abc\:\!}$, for the third-order results see
refs. [84, 85, 14], prove to be more
We have also started  to work on a first set of $\nfo$ contributions, the
$\cft\, \nf$ terms, where we hope to be able to provide the first exact
results of $\nfo$ parts of $P_{\rm ns}^{\,+(3)}(x)$ beyond the large-$\nc$
limit of ref. [86].

Form files with our results can be obtained from the preprint server
https://arXiv.org by downloading the source of this article.
They are also available from the authors upon request.
These files include also the analytic $N$-dependence of the third-order
quantities $a_{2,\rm ns}^{(n)}(N)$ and $a_{L,\rm ns}^{(n)}(N)$ in
eq. (<ref>) which we have not included in section 4 above.

§.§ Acknowledgements

We would like to thank Jos Vermaseren, Johann Usovitsch and Sven Moch
for useful discussions.
This work has been supported by the Vidi grant 680-47-551 of the
Dutch Research Council (NWO), the UK Research and Innovation (UKRI)
FLF Mr/S03479x/1 and the Consolidated Grant ST/T000988/1 of the UK
Science and Technology Facilities Council (STFC).
A.B.S. is grateful to Nikhef for the additional support provided
to him during the final stages of this project.

[1]
R.L. Workman et al., Particle Data Group,
Review of Particle Physics,
Prog. Theor. Exp. Phys. 2022, 083C01

[2]
A. Accardi et al.,
Electron Ion Collider: The Next QCD Frontier: Understanding the glue
that binds us all,
Eur. Phys. J. A52 (2016) 268, arXiv:1212.1701

[3]
R. Abdul Khalek et al.,
Science Requirements and Detector Concepts for the Electron-Ion
Collider: EIC Yellow Report,
Nucl. Phys. A1026 (2022) 122447, arXiv:2103.05419

[4]
J.L. Abelleira Fernandez et al., LHeC Study Group,
A Large Hadron Electron Collider at CERN: Report on the Physics and
Design Concepts for Machine and Detector.
J. Phys. G39 (2012) 075001, arXiv:1206.2913

[5]
P. Agostini et al., LHeC and FCC-he Study Group,
The Large Hadron–Electron Collider at the HL-LHC,
J. Phys. G48 (2021) 110501, arXiv:2007.14491

[6]
P. Bolzoni, F. Maltoni, S.O. Moch and M. Zaro,
Higgs production via vector-boson fusion at NNLO in QCD,
Phys. Rev. Lett. 105 (2010) 011801, arXiv:1003.4451

[7]
F.A. Dreyer and A. Karlberg,
Vector-Boson Fusion Higgs Production at Three Loops in QCD,
Phys. Rev. Lett. 117 (2016) 072001, arXiv:1606.00840

[8]
J. Currie, T. Gehrmann, E.W.N. Glover, A. Huss, J. Niehues and A. Vogt,
N$^{\,3}$LO corrections to jet production in deep-inelastic scattering
using the Projection-to-Born method,
JHEP 05 (2018) 209, arXiv:1803.09973

[9]
J. Sanchez Guillen, J. Miramontes, M. Miramontes, G. Parente and O.A. Sampayo,
Next-to-leading order analysis of the deep inelastic
$R = \sigma_{\!L}^{} / \sigma_{\rm total}^{}$,
Nucl. Phys. B353 (1991) 337

[10]
W.L. van Neerven and E.B. Zijlstra,
Order $\as(2)$ contributions to the deep inelastic Wilson coefficient,
Phys. Lett. B272 (1991) 127

[11]
E.B. Zijlstra and W.L. van Neerven,
Order $\as(2)$ QCD corrections to the deep inelastic proton structure
functions $F_2$ and $F_L$,
Nucl. Phys. B383 (1992) 525

[12]
S. Moch and J.A.M. Vermaseren,
Deep inelastic structure functions at two loops,
Nucl. Phys. B573 (2000) 853, hep-ph/9912355

[13]
J.A.M. Vermaseren, A. Vogt and S. Moch,
The third-order QCD corrections to deep-inelastic scattering by
photon exchange,
Nucl. Phys. B724 (2005) 3, hep-ph/0504242

[14]
J. Blümlein, P. Marquard, C. Schneider and K. Schönwald,
The massless three-loop Wilson coefficients for the deep-inelastic
structure functions $F_2$, $F_L,$ $xF_3$ and $g_1^{}$,

[15]
B. Ruijl, T. Ueda, J.A.M. Vermaseren, J. Davies and A. Vogt,
First Forcer results on deep-inelastic scattering
and related quantities,
PoS LL 2016 (2016) 071, arXiv:1605.08408

[16]
S. Moch, B. Ruijl, T. Ueda, J.A.M. Vermaseren and A. Vogt,
DIS coefficient functions at four loops in QCD and beyond,
PoS LL2022 (2022) 047, arXiv:2208.11067

[17]
S. Moch, B. Ruijl, T. Ueda, J. A. M. Vermaseren and A. Vogt,
to appear

[18]
B. Ruijl, T. Ueda and J.A.M. Vermaseren,
Forcer, a FORM program for the parametric reduction of
four-loop massless propagator diagrams,
Comput. Phys. Commun. 253 (2020) 107198, arXiv:1704.06650

[19]
J.A. Gracey,
Large $\nf$ methods for computing the perturbative structure of
deep inelastic scattering,
talk at the 4th International Workshop on Software Engineering
and Artificial Intelligence for High-energy and Nuclear Physics,
Pisa (Italy), April 1995, hep-ph/9509276

[20]
L. Mankiewicz, M. Maul and E. Stein,
Perturbative part of the non-nsinglet structure function $F_2$ in the
large $\nf$ limit,
Phys. Lett. B404 (1997) 345, hep-ph/9703356

[21]
J.A.M. Vermaseren,
Harmonic sums, Mellin transforms and integrals,
Int. J. Mod. Phys. A14 (1999) 2037, hep-ph/9806280

[22]
E. Remiddi and J.A.M. Vermaseren,
Harmonic polylogarithms,
Int. J. Mod. Phys. A15 (2000) 725, hep-ph/9905237

[23]
J. Davies, B. Ruijl, T. Ueda, J.A.M. Vermaseren and A. Vogt,
Large-$\nf$ contributions to the four-loop splitting functions in QCD,
Nucl. Phys. B915 (2017) 335, arXiv:1610.07477

[24]
A.V. Kotikov,
Differential equations method: New technique for massive Feynman
diagrams calculation,
Phys. Lett. B254 (1991) 158

[25]
A.V. Kotikov,
Differential equation method: The Calculation of N-point Feynman
Phys. Lett. B267 (1991) 123
[erratum: Phys. Lett. B295 (1992) 409]

[26]
E. Remiddi,
Differential equations for Feynman graph amplitudes,
Nuovo Cim. A110 (1997) 1435

[27]
T. Gehrmann and E. Remiddi,
Differential equations for two loop four point functions,
Nucl. Phys. B580 (2000) 485, hep-ph/9912329

[28]
W.L. van Neerven and A. Vogt,
NNLO evolution of deep-inelastic structure functions:
The Singlet case,
Nucl. Phys. B588 (2000) 345, hep-ph/0006154

[29]
S.A. Larin, T. van Ritbergen and J. Vermaseren,
The next next-to-leading QCD approximation for non-singlet moments
of deep inelastic structure functions,
Nucl. Phys. B427 (1994) 41

[30]
S.A. Larin, P. Nogueira, T. van Ritbergen and J.A.M. Vermaseren,
The three loop QCD calculation of the moments of deep inelastic
structure functions,
Nucl. Phys. B492 (1997) 338, hep-ph/9605317

[31]
S. Moch and M. Rogal,
Charged current deep-inelastic scattering at three loops,
Nucl. Phys. B782 (2007) 51, arXiv:0704.1740

[32]
O.V. Tarasov, A.A. Vladimirov and A.Y. Zharkov,
The Gell-Mann-Low Function of QCD in the Three Loop Approximation,
Phys. Lett. 93B (1980) 429

[33]
S.A. Larin and J.A.M. Vermaseren,
The Three loop QCD beta function and anomalous dimensions,
Phys. Lett. B303 (1993) 334, hep-ph/9302208

[34]
S. Moch, J.A.M. Vermaseren and A. Vogt,
Three-loop results for quark and gluon form-factors,
Phys. Lett. B625 (2005) 245, hep-ph/0508055

[35]
S. Moch, J.A.M. Vermaseren and A. Vogt,
Non-singlet structure functions at three loops:
Fermionic contributions,
Nucl. Phys. B646 (2002) 181, hep-ph/0209100

[36]
P. Nogueira,
Automatic Feynman graph generation,
J. Comput. Phys. 105 (1993) 279

[37]
J.A.M. Vermaseren, New features of FORM, math-ph/0010025

[38]
J. Kuipers, T. Ueda, J.A.M. Vermaseren and J. Vollinga,
FORM version 4.0,
Comput. Phys. Commun. 184 (2013) 1453, arXiv:1203.6543

[39]
B. Ruijl, T. Ueda and J.A.M. Vermaseren,
FORM version 4.2,

[40]
J.A.M. Vermaseren,
The minos database facility,

<https://www.nikhef.nl/ form/maindir/others/minos/minos.html>.

[41]
S.G. Gorishnii, S.A. Larin, L.R. Surguladze and F.V. Tkachov,
Mincer: Program for Multiloop Calculations in Quantum Field Theory
for the Schoonschip System,
Comput. Phys. Commun. 55 (1989) 381

[42]
S.A. Larin, F.V. Tkachov and J.A.M. Vermaseren,
The FORM version of MINCER, NIKHEF-H-91-18

[43]
S. Laporta,
High precision calculation of multiloop Feynman integrals by
difference equations,
Int. J. Mod. Phys. A15 (2000) 5087, hep-ph/0102033

[44]
A.V. Smirnov and F.S. Chuharev,
FIRE6: Feynman Integral REduction with Modular Arithmetic,
Comput. Phys. Commun. 247 (2020) 106877

[45]
P. Maierhöfer, J. Usovitsch and P. Uwer,
Kira—A Feynman integral reduction program,
Comput. Phys. Commun. 230 (2018) 99

[46]
R. Boughezal, M. Czakon and T. Schutzmeier,
NNLO fermionic corrections to the charm quark mass dependent
matrix elements in $\bar B \to X_s \gamma$,
JHEP 09 (2007) 072, arXiv:0707.3090

[47]
J. Blümlein and C. Schneider,
The Method of Arbitrarily Large Moments to Calculate Single Scale
Processes in Quantum Field Theory,
Phys. Lett. B771 (2017) 31, arXiv:1701.04614

[48]
R.N. Lee, A.V. Smirnov and V.A. Smirnov,
Solving differential equations for Feynman integrals by expansions
near singular points,
JHEP 03 (2018) 008, arXiv:1709.07525

[49]
X. Liu, Y.Q. Ma and C.Y. Wang,
A Systematic and Efficient Method to Compute Multi-loop Master
Phys. Lett. B779 (2018) 353, arXiv:1711.09572

[50]
B. Mistlberger,
Higgs boson production at hadron colliders at N$^{3}$LO in QCD,
JHEP 05 (2018) 028, arXiv:1802.00833

[51]
F. Moriello,
Generalised power series expansions for the elliptic planar families
of Higgs$\,+\,$jet production at two loops,
JHEP 01 (2020) 150, arXiv:1907.13234

[52]
I. Dubovyk, A. Freitas, J. Gluza, K. Grzanka, M. Hidding and J. Usovitsch,
Evaluation of multi-loop multi-scale Feynman integrals for precision

[53]
M. Fael, F. Lange, K. Schönwald and M. Steinhauser,
Singlet and non-singlet three-loop massive form factors,
Phys. Rev. D106 (2022) 034029, arXiv:2207.00027

[54]
M. Hidding,
DiffExp, a Mathematica package for computing Feynman
integrals in terms of one-dimensional series expansions,
Comput. Phys. Commun. 269 (2021) 108125, arXiv:2006.05510

[55]
X. Liu and Y. Q. Ma,
AMFlow: A Mathematica package for Feynman integrals computation via
auxiliary mass flow,
Comput. Phys. Commun. 283 (2023) 108565, arXiv:2201.11669

[56]
T. Armadillo, R. Bonciani, S. Devoto, N. Rana and A. Vicini,
Evaluation of Feynman integrals with arbitrary complex masses via
series expansions,
Comput. Phys. Commun. 282 (2023) 108545, arXiv:2205.03345

[57]
J.M. Henn,
Multiloop integrals in dimensional regularization made simple,
Phys.  Rev. Lett. 110 (2013) 251601, arXiv:1304.1806

[58]
R.N. Lee,
Reducing differential equations for multiloop master integrals,
JHEP 04 (2015) 108, arXiv:1411.0911

[59]
J. Moser,
The order of a singularity in Fuchs' theory,
Mathematische Zeitschrift 72 (1959) 379

[60]
R.N. Lee,
Libra: A package for transformation of differential systems
for multiloop integrals,
Comput. Phys. Commun. 267 (2021) 108058, arXiv:2012.00279

[61]
O. Gituliar and V. Magerya,
Fuchsia: a tool for reducing differential equations for Feynman
master integrals to epsilon form,
Comput. Phys. Commun. 219 (2017) 329, arXiv:1701.04269

[62]
M. Prausa,
epsilon: A tool to find a canonical basis of master integrals,
Comput. Phys. Commun. 219 (2017) 361, arXiv:1701.00725

[63]
C. Sabbah,
Lieu des pôles d'un système holonome d'équations aux
différences finies,
Bulletin de la Société Mathématique de France 120 (1992) 371

[64]
A.V. Smirnov and V.A. Smirnov,
How to choose master integrals,
Nucl. Phys. B960 (2020) 115213, arXiv:2002.08042

[65]
J. Usovitsch,
Factorization of denominators in integration-by-parts reductions,

[66]
M. Jamin and R. Miravitllas,
Absence of even-integer $\zeta$-function values in Euclidean
physical quantities in QCD,
Phys. Lett. B779 (2018) 452, arXiv:1711.00787

[67]
P.A. Baikov and K.G. Chetyrkin,
The structure of generic anomalous dimensions and no-$\pi$ theorem
for massless propagators,
JHEP 06 (2018) 141, arXiv:1804.10088

[68]
J. Davies and A. Vogt,
Absence of $\pi^2$ terms in physical anomalous dimensions in DIS:
Verification and resulting predictions,
Phys. Lett. B776 (2018) 1894, arXiv:1711.05267

[69]
W.L. van Neerven and A. Vogt,
Non-singlet structure functions beyond the next-to-next-to-leading
Nucl. Phys. B603 (2001) 42, arXiv:hep-ph/0103123

[70]
P.A. Baikov, K.G. Chetyrkin and J.H. Kühn,
Five-Loop Running of the QCD coupling constant,
Phys. Rev. Lett. 118 (2017) 082002, arXiv:1606.08659

[71]
F. Herzog, B. Ruijl, T. Ueda, J.A.M. Vermaseren and A. Vogt,
The five-loop beta function of Yang-Mills theory with fermions,
JHEP 02 (2017) 090, arXiv:1701.01404

[72]
T. Luthe, A. Maier, P. Marquard and Y. Schröder,
The five-loop beta function for a general gauge group and
anomalous dimensions beyond Feynman gauge,
JHEP 10 (2017) 166, arXiv:1709.07718

[73]
F. Herzog, S. Moch, B. Ruijl, T. Ueda, J. Vermaseren, A. Vogt,
Five-loop contributions to low-N non-singlet anomalous dimensions
in QCD,
Phys. Lett. B790 (2019) 436, arXiv:1812.11818

[74]
F. Herzog and B. Ruijl,
The R$^{\:\!\ast}$-operation for Feynman graphs with generic
numerators, JHEP 05 (2017) 037, arXiv:1703.03776

[75]
J.A. Gracey,
Anomalous dimension of non-singlet Wilson operators at $O (1/N_f)$
in deep inelastic scattering,
Phys. Lett. B322 (1994) 141, hep-ph/9401214

[76]
T. Gehrmann and E. Remiddi,
Numerical evaluation of harmonic polylogarithms,
Comput. Phys. Commun. 141 (2001) 296, hep-ph/0107173

[77]
S. Moch, J.A.M. Vermaseren and A. Vogt,
Higher-order corrections in threshold resummation,

Nucl. Phys. B726 (2005) 317, hep-ph/0506288v2

[78]
G. Das, S. Moch, A. Vogt,
Soft corrections to inclusive deep-inelastic scattering at four loops
and beyond,
JHEP 03 (2020) 116, arXiv:1912.12920

[79]
S. Moch and A. Vogt,
On non-singlet physical evolution kernels and large-x coefficient
in perturbative QCD,
JHEP 11 (2009) 099, arXiv:0909.2124

[80]
G. Grunberg,
Large-x structure of physical evolution kernels in Deep Inelastic
Phys. Lett. B687 (2010) 405, arXiv:0911.4471v5

[81]
A.A. Almasy, G. Soar and A. Vogt,
Generalized double-logarithmic large-x resummation in inclusive
deep-inelastic scattering,
JHEP 1103 (2011) 030, arXiv:1012.3352

[82]
S. Moch and A. Vogt,
Threshold Resummation of the Structure Function $F_{\!L\,}$,
JHEP 04 (2009) 081, arXiv:0902.2342

[83]
J. Davies, C.H. Kom, S. Moch and A. Vogt,
Resummation of small-$x$ double logarithms in QCD:
inclusive deep-inelastic scattering,
JHEP 08 (2022) 135, arXiv:2202.10362

[84]
A. Retey and J.A.M. Vermaseren,
Some higher moments of deep inelastic structure functions at
next-to-next-to-leading order of perturbative QCD,
Nucl. Phys. B604 (2001) 281, hep-ph/0007294

[85]
S. Moch, J.A.M. Vermaseren and A. Vogt,
Third-order QCD corrections to the charged-current structure function
Nucl. Phys. B813 (2009) 220, arXiv:0812.4168

[86]
S. Moch, B. Ruijl, T. Ueda, J.A.M. Vermaseren and A. Vogt,
Four-Loop Non-Singlet Splitting Functions in the Planar Limit and
JHEP 10 (2017) 041, arXiv:1707.08315

authors:
A. Basdew-Sharma, A. Pelloni (Nikhef, Amsterdam), Franz Herzog (U.
  Edinburgh, Higgs Ctr. Theor. Phys.) and A. Vogt (Liverpool U., Dept. Math.)

submitter: Andreas Vogt