# NNLO QCD corrections to semi-inclusive DIS

Saurav Goyal<EMAIL_ADDRESS>The Institute of Mathematical Sciences, A CI
of Homi Bhabha National Institute, Taramani, 600113 Chennai, India  Sven-Olaf
Moch<EMAIL_ADDRESS>II. Institute for Theoretical Physics, Hamburg
University, D-22761 Hamburg, Germany  Vaibhav Pathak<EMAIL_ADDRESS>The
Institute of Mathematical Sciences, A CI of Homi Bhabha National Institute,
Taramani, 600113 Chennai, India  Narayan Rana<EMAIL_ADDRESS>School
of Physical Sciences, National Institute of Science Education and Research, An
OCC of Homi Bhabha National Institute, 752050 Jatni, India  V. Ravindran
<EMAIL_ADDRESS>The Institute of Mathematical Sciences, A CI of Homi
Bhabha National Institute, Taramani, 600113 Chennai, India

###### Abstract

We present the first results for the next-to-next-to leading order (NNLO)
corrections to the semi-inclusive deep-inelastic scattering process in
perturbative quantum chromodynamics. We consider the quark initiated flavor
non-singlet process and obtain the complete contributions analytically at
leading color. All relevant virtual and real emission Feynman diagrams have
been computed using integration-by-parts reduction to master integrals and two
approaches for their subsequent evaluation (parametric phase-space integration
and method of differential equations). The numerical analysis demonstrates the
significance of the NNLO corrections and their great impact on the reduction
of the residual scale dependence.

Landmark inclusive measurements of structure functions (SF) in deep-inelastic
scattering (DIS) experiments provide a wealth of information on the internal
structure of hadrons at high energies in terms of quarks and gluons and give
valuable insights into the underlying strong interaction dynamics. Collinear
factorization in quantum chromodynamics (QCD) with the systematic separation
of short- and long-distance phenomena provides the theoretical foundations for
the use of SF data in the extraction of the process-independent parton
distribution functions (PDFs), which describe the parton dynamics within the
colliding hadrons. Semi-inclusive DIS (SIDIS) measurements with an observed
hadron in the final state allows, in addition, for the study of the parton
(quark/gluon) to hadron fragmentation, encoded in universal fragmentation
functions (FFs). This makes SIDIS the most promising and valuable probe of
PDFs and FFs at the upcoming Electron-Ion collider (EIC) at the Brookhaven
National Laboratory in the USA. The EIC opens up unique opportunities to
explore the nucleon structure as well as the underlying dynamics of hadrons in
various environments. This includes a thorough understanding of the light-
quark flavor PDFs and also the spin structure of the nucleon, using polarized
beams, to measure the polarized PDFs. A quantitative assessment of the
anticipated precision of SIDIS measurements and their impact on the
determination of various hadronic observables at the EIC has been presented in
ref. Aschenauer et al. (2019).

We consider the SIDIS process $l(k_{l})+H(P)\rightarrow
l({k}^{\prime}_{l})+H^{\prime}(P_{H})+X$ with lepton momenta $k_{l}$,
${k}^{\prime}_{l}$ and space-like momentum transfer,
$q={k}_{l}-{k}^{\prime}_{l}$ with $Q^{2}=-q^{2}$. The momenta of the incoming
and outgoing hadrons are denoted as $P$, $P_{H}$. The computation of SIDIS
observables in perturbation theory can nowadays be performed with an
unprecedented accuracy, thanks to the remarkable theoretical developments in
the calculation of multi-loop and multi-leg scattering processes. In addition,
resummation techniques applicable in particular kinematical limits allow, e.g.
to obtain approximate predictions based on the summation of large threshold
logarithms to all orders, see refs. Cacciari and Catani (2001); Anderle et al.
(2013a, b); Abele et al. (2021, 2022).

The longitudinal momentum distribution of the final state hadron is sensitive
to the fragmentation functions and depends on scaling variables $x$, $y$ and
$z$,

$\displaystyle\frac{d^{2}\sigma_{e^{-}H}}{dxdydz}=\frac{2\pi
y\alpha_{e}^{2}}{Q^{4}}L^{\mu\nu}(k_{l},k^{\prime}_{l},q)W_{\mu\nu}(q,P,P_{H})\,,$
(1)

where $\alpha_{e}$ is the fine structure constant, $x=\frac{Q^{2}}{2P\cdot q}$
is the Bjorken variable, $y=\frac{P\cdot q}{P\cdot k_{l}}$ the fraction of the
initial energy transferred to the hadron and $z=\frac{P\cdot P_{H}}{P\cdot q}$
the scaling variable corresponding to the fragmenting hadron. The leptonic
tensor reads
$L^{\mu\nu}=2{k}_{l}^{\mu}{k}_{l}^{\prime\nu}+2{k}_{l}^{\prime\mu}{k}_{l}^{\nu}-Q^{2}g^{\mu\nu}$,
while the hadronic tensor $W_{\mu\nu}$ can be expressed in terms of SFs
$F_{I}(x,z,Q^{2}),I=1,2$, as

$\displaystyle
W_{\mu\nu}=F_{1}(x,z,Q^{2})T_{1,\mu\nu}+F_{2}(x,z,Q^{2})T_{2,\mu\nu}\,,$ (2)

with tensors $T_{1,\mu\nu}$ = $-g_{\mu\nu}+\frac{q_{\mu}q_{\nu}}{q\cdot q}$
and $T_{2,\mu\nu}=\frac{1}{P.q}(P_{\mu}-\frac{P\cdot q}{q\cdot
q}q_{\mu})(P_{\nu}-\frac{P\cdot q}{q\cdot q}q_{\nu})$. QCD factorization
allows us to express the SFs as

$\displaystyle F_{I}$
$\displaystyle=x^{I-1}\sum_{a,b}\int_{x}^{1}\frac{dx_{1}}{x_{1}}f_{a}(x_{1},\mu_{F}^{2})\int_{z}^{1}\frac{dz_{1}}{z_{1}}D_{b}(z_{1},\mu_{F}^{2})$
$\displaystyle\times\mathcal{F}_{I,ab}\left(\frac{x}{x_{1}},\frac{z}{z_{1}},\mu_{F}^{2},Q^{2}\right)\,,$
(3)

where $f_{a}$ and $D_{b}$ denote the PDFs and the FFs, subject to summation
over initial state partons ‘$a$’ from the incoming hadron and final state
partons ‘$b$’ that fragments into the observed hadron. The coefficient
functions (CFs) $\mathcal{F}_{I,ab}$ can be computed in a perturbative
expansion in powers of the strong coupling
$a_{s}(\mu_{R}^{2})=\frac{\alpha_{s}(\mu_{R}^{2})}{4\pi}$,

$\displaystyle{\cal F}_{I,ab}$ $\displaystyle=$
$\displaystyle\sum_{i=0}^{\infty}\,a_{s}^{i}(\mu_{R}^{2})\,{\cal
F}_{I,ab}^{(i)}(\mu_{R}^{2})\,,$ (4)

with $\mu_{F}$, $\mu_{R}$ the factorisation and renormalisation scales.

To date, exact results for CFs are available only at the next-to-leading order
(NLO) in QCD Altarelli et al. (1979); Furmanski and Petronzio (1982), covering
all parton combinations $ab=qq,\bar{q}\bar{q},qg,\bar{q}g,gq,g\bar{q}$. In
this letter, we present the first results at next-to-next-to leading order
(NNLO) for the non-singlet channels $ab=qq,\bar{q}\bar{q}$, i.e., for SIDIS
processes with an initial quark or an anti-quark, which fragments into the
observed final state hadron. The results are complete to leading color $N_{c}$
for QCD as an $SU(N_{c})$ gauge theory. Previously, refs. Abele et al. (2021),
Abele et al. (2022) have derived the threshold enhanced logarithms using
resummed results for these partonic channels ($ab=qq,\bar{q}\bar{q}$) at NNLO
and even up to next-to-next-to-next-to leading order (N3LO). We confirm that
the resummation of ref. Abele et al. (2021) correctly predicts the dominant
threshold logarithms.

The computation of the CFs, $\mathcal{F}_{I,ab}$ in perturbative QCD in powers
of $a_{s}$, see eq. (4), starts from the parton level cross sections denoted
by $d\hat{\sigma}_{I,ab}$, appropriately set up with the corresponding
projectors $\mathcal{P}^{\mu\nu}_{I}$,

$\displaystyle
d\hat{\sigma}_{I,ab}=\frac{\mathcal{P}_{I}^{\mu\nu}}{4\pi}\int\text{dPS}_{X+b}\,\overline{\Sigma}|{M}_{ab}|^{2}_{\mu\nu}\,\delta\Big{(}\frac{z}{z_{1}}-\frac{p_{a}\cdot
p_{b}}{p_{a}\cdot q}\Big{)}\,,$ (5)

where $|M_{ab}|^{2}$ is the squared amplitude for the process
$a(p_{a})+\gamma^{*}(q)\rightarrow b(p_{b})+X$, with the parton ‘$b$’ tagged
to fragment into hadron $H^{\prime}$. $dPS_{X+b}$ stands for the phase space
of the final state consisting of $X$ and the fragmenting parton $b$ and
$\overline{\Sigma}$ denotes the average over initial and summation over final
state spin/polarization and color quantum numbers. In $D$ dimensions, the
projectors $\mathcal{P}^{\mu\nu}_{I}$ are given by
$\mathcal{P}^{\mu\nu}_{1}$=$\frac{1}{(D-2)}(T^{\mu\nu}_{1}+2xT^{\mu\nu}_{2})$
and
$\mathcal{P}^{\mu\nu}_{2}$=$\frac{2x}{(D-2)x_{1}}(T^{\mu\nu}_{1}+2x(D-1)T^{\mu\nu}_{2})$.
The partonic scaling variables are $x_{1}=\frac{p_{a}}{P}$ for the momentum
fraction carried by the initial parton ‘$a$’ of incident hadron $H$ and
$z_{1}=\frac{P_{H}}{p_{b}}$ for the corresponding fraction of hadron
$H^{\prime}$ with respect to the final state parton ‘$b$’.

Beyond leading order in perturbation theory, both ultraviolet (UV) and
infrared divergences resulting from soft and collinear partons appear in the
computation of $d\hat{\sigma}_{I,ab}$, and we work in $D=4+\varepsilon$
dimension to regulate them. The ultraviolet divergences are removed by
renormalisation of the strong coupling at the scale $\mu_{R}$, while all
infrared divergences cancel among virtual and real emission processes, except
for the collinear divergences related to the ‘$ab$’ partons in the initial
state and the fragmentation into hadrons. The structure of those divergences
is described by mass factorisation in QCD. The partonic cross sections in eq.
(5) factorise into the Altarelli-Parisi (AP) kernels $\Gamma_{c\leftarrow a}$
of PDFs and $\tilde{\Gamma}_{b\leftarrow d}$ of FFs, which capture the
collinear divergences in $1/\varepsilon$, and the CFs (${\cal F}_{I,cd}$),
which are finite as $\varepsilon\rightarrow 0$. Mass factorisation at the
scale $\mu_{F}$ reads

$\displaystyle\frac{d\hat{\sigma}_{I,ab}(\varepsilon)}{x^{\prime
I-1}}=\Gamma_{c\leftarrow a}(\mu_{F}^{2},\varepsilon)\otimes{\cal
F}_{I,cd}(\mu_{F}^{2},\varepsilon)\tilde{\otimes}\tilde{\Gamma}_{b\leftarrow
d}(\mu_{F}^{2},\varepsilon)\,,\,$ (6)

where $x^{\prime}=x/x_{1}$, summation over $c,d$ is implied and $\otimes$ and
$\tilde{\otimes}$ denote convolutions over scaling variables corresponding to
the PDFs and FFs, respectively, cf. eq. (NNLO QCD corrections to semi-
inclusive DIS). The AP kernels are determined by the evolution equations for
PDFs and FFs in terms of space- and time-like splitting functions, fully known
to third order in $a_{s}$ Moch et al. (2004); Vogt et al. (2004); Almasy et
al. (2012); Chen et al. (2021). Due to the collinear poles in $1/\varepsilon$
of the AP kernels, the extraction of the CFs (${\cal F}_{I,ab}$) to NNLO
accuracy from the convolutions in eq. (6) requires the computation of NLO
$d\hat{\sigma}_{I,ab}$ up to positive powers of $\varepsilon$.

At NLO, the partonic cross sections in eq. (5) get contributions from the one-
loop corrections to the Born process $\gamma^{*}+q(\bar{q})\rightarrow
q(\bar{q})$ and the real emission $\gamma^{*}+q(\bar{q})\rightarrow
q(\bar{q})+g$. In addition, there is the gluon-initiated channel
$\gamma^{*}+g\rightarrow q+\bar{q}$. In the new NNLO computation of
$d\hat{\sigma}_{i,qq}$, we restrict ourselves to the quark flavor non-singlet
case, i.e., we consider the following three classes, namely, two-loop
corrections to $\gamma^{*}+q(\bar{q})\rightarrow q(\bar{q})$, one-loop
contributions to the single gluon real emission
$\gamma^{*}+q(\bar{q})\rightarrow q(\bar{q})+g$, and double real emissions
$\gamma^{*}+q(\bar{q})\rightarrow q(\bar{q})+g+g$.

One- and two-loop virtual corrections to the Born process can be obtained
using the vector form factor which is known to fourth order in $a_{s}$ Lee et
al. (2022). For the real-virtual and double real emission processes we follow
the standard diagrammatic approach. We generate Feynman diagrams using QGRAF
Nogueira (1993) and use a set of in-house routines in FORM Kuipers et al.
(2013); Ruijl et al. (2017), which convert the output of QGRAF into a suitable
format that allows one to use the Feynman rules and to perform Lorentz
contractions, color as well as Dirac algebra. The real-virtual processes
contain one-loop and two-body phase-space integrals, while double real
emission processes contain three-body phase-space integrals. The phase-space
integrals are performed with the constraint $\frac{z}{z_{1}}=\frac{p_{a}\cdot
p_{b}}{p_{a}\cdot q}$, because of which, the computation of these phase-space
integrals is technically challenging compared to the fully inclusive ones.
Using the method of reverse unitarity Anastasiou et al. (2004, 2012), we
convert all the phase-space integrals into loop integrals and apply
integration-by-parts identities (IBP) Chetyrkin and Tkachov (1981); Laporta
(2000) to reduce them to a smaller number of the master integrals (MI). In the
reverse unitary method, we replace the on-shell Dirac delta functions by the
corresponding propagators. The constraint ${z\over z_{1}}=\frac{p_{a}\cdot
p_{b}}{p_{a}\cdot q}$ is imposed by an additional delta function
$\delta\left(z^{\prime}-\frac{p_{a}\cdot p_{b}}{p_{a}\cdot q}\right)$, where
$z^{\prime}=\frac{z}{z_{1}}$, which is then replaced by a propagator-like term
$-\frac{1}{\pi}\text{Im}(1/(z^{\prime}-\frac{p_{a}\cdot p_{b}}{p_{a}\cdot
q}+i\epsilon))$ with $p_{2}=p_{1}+q-k_{1}$ or $p_{2}=p_{1}+q-k_{1}-k_{2}$ for
two- or three-body final states. We perform the IBP reduction with the
Mathematica package LiteRed Lee (2014) and obtain at the NNLO level at leading
$N_{c}$, 7 MIs for both the real-virtual and the double real emission sub-
processes. We note that the computation of the double real emission process
$\gamma^{*}+q(\bar{q})\rightarrow q(\bar{q})+q+\bar{q}$ as a part of the
flavor pure-singlet SFs (not considered in this paper) requires an additional
13 MIs for the color non-planar contributions, i.e. those suppressed as
$1/N_{c}$.

The calculation of the MIs poses one of the primary challenges in this
computation. We have followed two different approaches, applying either
parametric integration (PI) Matsuura et al. (1989); Zijlstra and van Neerven
(1992); Rijken and van Neerven (1997); Ravindran et al. (2003) or differential
equations (DE) Kotikov (1991); Argeri and Mastrolia (2007); Remiddi (1997);
Henn (2013); Ablinger et al. (2016). In the PI method, a convenient choice of
the Lorentz frame helps to solve the integrals. For example, in the case of a
two-body phase-space, we use the center-of-mass (COM) frame of the photon and
the incoming parton, so that only the one-loop integrals need to be done,
which reduce to hypergeometric functions of the scaling variables
$x^{\prime}$, $z^{\prime}$ and $\varepsilon$. The three-body phase-space is
conveniently integrated in the COM frame of two outgoing partons that do not
fragment and we need to perform integrals over one of the parametric variables
and two angles. Upon expansion in $\varepsilon$ the three-body phase-space
integrals then reduce to multiple polylogarithms (MPLs) and Nielsen
polylogarithms.

In the second approach, we have used the method of differential equations to
evaluate the MIs. The corresponding system of differential equations has been
generated for each topology by taking derivatives with respect to the
independent variables $x^{\prime}$ and $z^{\prime}$, using LiteRed for the
differentiation and further IBP reduction. We have obtained a system without
any coupled differential equations, which can easily be arranged in an upper-
triangular form, and allows for the bottom-up approach to solve the MIs one by
one. While solving for each MI, we also perform a Taylor series expansion in
$\varepsilon$ and obtain the results in terms of generalized harmonic
polylogarithms (GPLs). The boundary conditions for the MIs are either derived
from regularity conditions or by calculating explicitly the threshold limit.
We have used HarmonicSums Ablinger (2009); Ablinger et al. (2011); Ablinger
(2014) and PolyLogTools Duhr and Dulat (2019) in various intermediate steps of
the computation including the conversion of MPLs and Nielsen polylogarithms
into GPLs. We have cross-checked the results obtained from both approaches and
found perfect agreement.

The phase-space integrals in $D=4+\varepsilon$ generally result into functions
of the form

$(1-x^{\prime})^{-1+a\frac{\varepsilon}{2}}(1-z^{\prime})^{-1+b\frac{\varepsilon}{2}}f_{1}(x^{\prime},z^{\prime},\varepsilon)\,,$

which contain simple poles at the thresholds i.e. $x^{\prime}\rightarrow 1$
and $z^{\prime}\rightarrow 1$ for integers $a$, $b$. In addition, we find
singularities at $x^{\prime}=z^{\prime}$ and at $x^{\prime}+z^{\prime}=1$ in
the form

$\frac{(1-x^{\prime})^{a\varepsilon}(1-z^{\prime})^{b\varepsilon}(z^{\prime}-x^{\prime})^{c\varepsilon}(1-z^{\prime}-x^{\prime})^{d\varepsilon}}{(1-x^{\prime})(1-z^{\prime})(z^{\prime}-x^{\prime})(1-z^{\prime}-x^{\prime})}f_{2}(x^{\prime},z^{\prime},\varepsilon)\,.$

The $f_{i}(x^{\prime},z^{\prime},\varepsilon)$ are regular functions in the
threshold limits $x^{\prime}\rightarrow 1$ and/or $z^{\prime}\rightarrow 1$.
The sign of the imaginary part of the term
$(z^{\prime}-x^{\prime})^{a\varepsilon}$ depends on whether
$x^{\prime}>z^{\prime}$ or $x^{\prime}<z^{\prime}$, while for the term
$(1-z^{\prime}-x^{\prime})^{b\varepsilon}$, it depends on whether
$x^{\prime}+z^{\prime}>1$ or $x^{\prime}+z^{\prime}<1$. We introduce the
identities $\theta(x^{\prime}-z^{\prime})+\theta(z^{\prime}-x^{\prime})=1$ and
$\theta(1-x^{\prime}-z^{\prime})+\theta(z^{\prime}+x^{\prime}-1)=1$ to
separate the different sectors. Note that these scaling variables should be
understood as $x^{\prime}-i\epsilon$ and $z^{\prime}-i\epsilon$ with Feynman’s
$i\epsilon$ prescription. The divergences resulting from the threshold region
can be isolated by using the following identity

$\displaystyle(1-w)^{-1+n\varepsilon}=\frac{1}{n\varepsilon}\delta(1-w)+\sum_{k=0}^{\infty}\frac{(n\varepsilon)^{k}}{k!}\mathcal{D}_{w,k}\,,$
(7)

where $w=x^{\prime},z^{\prime}$ and
$\mathcal{D}_{w,k}=\left[\frac{\log^{k}(1-w)}{(1-w)}\right]_{+}$ denote the
usual ‘plus’-distributions, as defined by

$\int_{0}^{1}dw~{}g_{+}(w)~{}f(w)=\int_{0}^{1}dw~{}g(w)~{}(f(w)-f(1))\,.$ (8)

Due to the presence of these singularities, it is henceforth necessary to
compute the phase-space integrals in closed form in $\varepsilon$, or at least
partially. At NLO level, the leading double poles ($1/\varepsilon^{2}$) from
the virtual contributions and the real emission diagrams cancel each other and
the remaining collinear divergence ($1/\varepsilon$) is removed by mass
factorisation with the AP kernels $\Gamma$ and $\tilde{\Gamma}$, see eq. (6).
At NNLO level, the leading $1/\varepsilon^{4}$ and $1/\varepsilon^{3}$ terms
cancel among the pure virtual, real-virtual and double real emission
contributions. All remaining $1/\varepsilon^{2}$ and $1/\varepsilon$ terms
cancel against those of the AP kernels during mass factorization.

After mass factorization given by eq. (6) we obtain the finite CFs ${\cal
F}_{I}^{(i)}$ ($\equiv{\cal F}_{I,ab}^{(i)}$ as a short-hand now) for
$i=0,1,2$. The Born contribution ${\cal F}_{I}^{(0)}$ is proportional to
$\delta(1-x^{\prime})\delta(1-z^{\prime})$ and our NLO results for ${\cal
F}_{I}^{(1)}$ are in complete agreement with refs. Anderle et al. (2013a);
Altarelli et al. (1979). It is convenient to separate the dependence on
distributions
$\delta(1-x^{\prime}),\delta(1-z^{\prime}),\mathcal{D}_{x^{\prime},k},\mathcal{D}_{z^{\prime},k}$
and on regular functions in $x^{\prime}$, $z^{\prime}$ such as MPLs as
follows,

${\cal F}_{I}^{(i)}={\cal F}_{I,2}^{(i)}+{\cal F}_{I,1}^{(i)}+{\cal
F}_{I,r}^{(i)}\,.$ (9)

At NNLO, the first term ${\cal F}_{I,2}^{(2)}$ contains only distributions in
$x^{\prime}$ and $z^{\prime}$, often called soft-plus-virtual (SV) terms and
is in agreement with ref. Abele et al. (2021), which provides the first
verification of those results from an explicit computation. The remaining
contributions in eq. (9), consisting of single distributions denoted by ${\cal
F}_{I,1}^{(2)}$ and pure regular terms ${\cal F}_{I,r}^{(2)}$, are new. As the
latter are too lengthy, we present here only the single distributions ${\cal
F}_{1,1}^{(2)}$ and introduce the following abbreviations for compact
presentation:
$\bar{x}=(1-x^{\prime})\,,\bar{z}=(1-z^{\prime})\,,\tilde{x}=(1+x^{\prime})\,,\tilde{z}=(1+z^{\prime})\,,\delta_{\bar{x}}=\delta(1-x^{\prime})\,,\delta_{\bar{z}}=\delta(1-z^{\prime})\,,l_{x}=\log(x^{\prime})\,,l_{z}=\log(z^{\prime})\,,l_{\bar{x}}=\log(1-x^{\prime})\,,l_{\bar{z}}=\log(1-z^{\prime})\,$.
The full results for ${\cal F}_{I}^{(i)}$ in eq. (9) are attached to this
publication as a supplementary file. We also note that, although our NNLO
results consist of the leading color contributions, we retain their full
Casimir structure in the presentation, since all sub-leading color
contributions (to be presented in future work) are proportional to
$C_{F}(C_{A}-2C_{F})$. Thus, the NNLO result for the single distributions of
${\cal F}_{1}^{(2)}$ reads:

$\displaystyle{\cal F}_{1,1}^{(2)}$
$\displaystyle=C_{F}^{2}\Big{[}\delta_{\bar{x}}\Big{\\{}2l_{z}^{2}(1-4\bar{z})+4(1-8\bar{z})-8\,\text{Li}_{3}(\bar{z})\tilde{z}+\frac{25}{3}l_{z}^{3}\tilde{z}-4l_{z}l_{\bar{z}}^{2}\tilde{z}-4l_{\bar{z}}^{3}\tilde{z}+52\,\text{S}_{12}(\bar{z})\tilde{z}+\,\text{Li}_{2}(\bar{z})(4(1-6\bar{z})$
$\displaystyle~{}~{}+40l_{z}\tilde{z})+\frac{1}{\bar{z}}\big{(}8\,\text{Li}_{3}(\bar{z})-64\,\text{Li}_{2}(\bar{z})l_{z}-\frac{40}{3}l_{z}^{3}+12l_{z}l_{\bar{z}}^{2}-88\,\text{S}_{12}(\bar{z})+l_{\bar{z}}\big{(}-8\,\text{Li}_{2}(\bar{z})-12l_{z}^{2}\big{)}+l_{z}\big{(}-64+24\zeta_{2}\big{)}\big{)}$
$\displaystyle~{}~{}+l_{\bar{z}}\big{(}14+24\tilde{z}+4l_{z}(1-2\bar{z})+8\,\text{Li}_{2}(\bar{z})\tilde{z}+10l_{z}^{2}\tilde{z}+16\tilde{z}\zeta_{2}\big{)}+l_{z}\big{(}-2+38\tilde{z}-16\tilde{z}\zeta_{2}\big{)}+8\bar{z}\zeta_{2}-16\tilde{z}\zeta_{3}\Big{\\}}$
$\displaystyle~{}~{}+{\mathcal{D}_{x,0}}\big{\\{}12+24\tilde{z}+4l_{z}(1-3\bar{z})+12\,\text{Li}_{2}(\bar{z})\tilde{z}+16l_{z}^{2}\tilde{z}-4l_{z}l_{\bar{z}}\tilde{z}-12l_{\bar{z}}^{2}\tilde{z}-\frac{1}{\bar{z}}\big{(}16\,\text{Li}_{2}(\bar{z})+24l_{z}^{2}-16l_{z}l_{\bar{z}}\big{)}+16\tilde{z}\zeta_{2}\big{\\}}$
$\displaystyle~{}~{}+{\mathcal{D}_{x,1}}\big{\\{}(4l_{z}\tilde{z}-24l_{\bar{z}}\tilde{z})\big{\\}}-12\tilde{z}~{}{\mathcal{D}_{x,2}}+\delta_{\bar{z}}\Big{\\{}-4-48\bar{x}-2l_{x}^{2}+\frac{11}{3}l_{x}^{3}\tilde{x}+16l_{x}l_{\bar{x}}^{2}\tilde{x}-4l_{\bar{x}}^{3}\tilde{x}-24\,\text{S}_{12}(\bar{x})\tilde{x}$
$\displaystyle~{}~{}+\,\text{Li}_{2}(\bar{x})(4+8\bar{x}-12l_{x}\tilde{x})+\frac{1}{\bar{x}}\big{(}-8\,\text{Li}_{3}(\bar{x})+16\,\text{Li}_{2}(\bar{x})l_{x}-4l_{x}^{3}-28l_{x}l_{\bar{x}}^{2}+48\,\text{S}_{12}(\bar{x})+l_{\bar{x}}\big{(}8\,\text{Li}_{2}(\bar{x})+32l_{x}^{2}\big{)}$
$\displaystyle~{}~{}+l_{x}\big{(}64+32\zeta_{2}\big{)}\big{)}+l_{\bar{x}}\big{(}14+26\tilde{x}+4l_{x}-20l_{x}^{2}\tilde{x}+16\tilde{x}\zeta_{2}\big{)}+l_{x}\big{(}-8-34\tilde{x}-20\tilde{x}\zeta_{2}\big{)}+8\bar{x}\zeta_{2}-16\tilde{x}\zeta_{3}\Big{\\}}$
$\displaystyle~{}~{}+{\mathcal{D}_{z,0}}\big{\\{}12+28\tilde{x}+4l_{x}(1+\bar{x})-4\,\text{Li}_{2}(\bar{x})\tilde{x}-12l_{x}^{2}\tilde{x}+28l_{x}l_{\bar{x}}\tilde{x}-12l_{\bar{x}}^{2}\tilde{x}+\frac{1}{\bar{x}}\big{(}16\,\text{Li}_{2}(\bar{x})+16l_{x}^{2}-48l_{x}l_{\bar{x}}\big{)}+16\tilde{x}\zeta_{2}\big{\\}}$
$\displaystyle~{}~{}+{\mathcal{D}_{z,1}}\Big{\\{}-\frac{32}{\bar{x}}l_{x}+20l_{x}\tilde{x}-24l_{\bar{x}}\tilde{x}\Big{\\}}-12\tilde{x}{\mathcal{D}_{z,2}}\Big{]}+C_{A}C_{F}\Big{[}\delta_{\bar{x}}\Big{\\{}4\,\text{Li}_{2}(\bar{z})(1-2\tilde{z})+\frac{2}{3}l_{z}^{2}(3-11\tilde{z})+\frac{1}{27}(396+179\tilde{z})$
$\displaystyle~{}~{}+\frac{1}{9}l_{z}(1+70\bar{z})-4\,\text{Li}_{3}(\bar{z})\tilde{z}-\frac{5}{3}l_{z}^{3}\tilde{z}+\frac{11}{3}l_{\bar{z}}^{2}\tilde{z}+6\,\text{S}_{12}(\bar{z})\tilde{z}+\frac{1}{\bar{z}}\big{(}6\,\text{Li}_{2}(\bar{z})+8\,\text{Li}_{3}(\bar{z})+\frac{62}{9}l_{z}+\frac{49}{6}l_{z}^{2}+\frac{10}{3}l_{z}^{3}+6l_{z}l_{\bar{z}}$
$\displaystyle~{}~{}-12\,\text{S}_{12}(\bar{z})\big{)}+l_{\bar{z}}\Big{(}\frac{4}{9}(15-41\tilde{z})-6l_{z}\tilde{z}+4\tilde{z}\zeta_{2}\Big{)}-\frac{4}{3}(3+4\tilde{z})\zeta_{2}-14\tilde{z}\zeta_{3}\Big{\\}}+{\mathcal{D}_{x,0}}\Big{\\{}\frac{2}{9}(39-82\tilde{z})-4\,\text{Li}_{2}(\bar{z})\tilde{z}-6l_{z}\tilde{z}$
$\displaystyle~{}~{}-2l_{z}^{2}\tilde{z}+\frac{22}{3}l_{\bar{z}}\tilde{z}+\frac{1}{\bar{z}}\big{(}8\,\text{Li}_{2}(\bar{z})+6l_{z}+4l_{z}^{2}\big{)}+4\tilde{z}\zeta_{2}\Big{\\}}+\frac{22}{3}\tilde{z}{\mathcal{D}_{x,1}}+\delta_{\bar{z}}\Big{\\{}\frac{46}{3}+\frac{197}{27}\tilde{x}+8\,\text{Li}_{3}(\bar{x})\tilde{x}+\frac{55}{6}l_{x}^{2}\tilde{x}+\frac{11}{3}l_{\bar{x}}^{2}\tilde{x}$
$\displaystyle~{}~{}+2\,\text{S}_{12}(\bar{x})\tilde{x}+\,\text{Li}_{2}(\bar{x})\Big{(}-4l_{x}\tilde{x}-\frac{4}{3}(3+4\tilde{x})\Big{)}+l_{x}\Big{(}\frac{1}{3}(-13+77\tilde{x})-4\tilde{x}\zeta_{2}\Big{)}+\frac{1}{\bar{x}}\Big{(}-16\,\text{Li}_{3}(\bar{x})-\frac{83}{6}l_{x}^{2}$
$\displaystyle~{}~{}+\frac{2}{3}\,\text{Li}_{2}(\bar{x})(13+12l_{x})+\frac{70}{3}l_{x}l_{\bar{x}}-4\,\text{S}_{12}(\bar{x})+l_{x}\Big{(}-\frac{116}{3}+8\zeta_{2}\Big{)}\Big{)}+l_{\bar{x}}\Big{(}-\frac{44}{3}l_{x}\tilde{x}-\frac{4}{9}(-15+41\tilde{x})+4\tilde{x}\zeta_{2}\Big{)}$
$\displaystyle~{}~{}-\frac{4}{3}(3+4\tilde{x})\zeta_{2}-14\tilde{x}\zeta_{3}\Big{\\}}+{\mathcal{D}_{z,0}}\Big{\\{}4\,\text{Li}_{2}(\bar{x})\tilde{x}-\frac{44}{3}l_{x}\tilde{x}+\frac{22}{3}l_{\bar{x}}\tilde{x}+\frac{26}{9}(3-7\tilde{x})+\frac{1}{\bar{x}}\big{(}-8\,\text{Li}_{2}(\bar{x})+\frac{70}{3}l_{x}\big{)}+4\tilde{x}\zeta_{2}\Big{\\}}$
$\displaystyle~{}~{}+\frac{22}{3}\tilde{x}{\mathcal{D}_{z,1}}\Big{]}+\frac{1}{3}C_{F}n_{F}\Big{[}\delta_{\bar{x}}\Big{\\{}4-\frac{2}{3}\frac{1}{\bar{z}}l_{z}(10+3l_{z})-\frac{74}{9}\tilde{z}+l_{z}^{2}\tilde{z}-2l_{\bar{z}}^{2}\tilde{z}+\frac{8}{3}l_{\bar{z}}(-3+4\tilde{z})+\frac{2}{3}l_{z}(-12+11\tilde{z})+4\tilde{z}\zeta_{2}\Big{\\}}$
$\displaystyle~{}~{}+{\mathcal{D}_{x,0}}\Big{\\{}\frac{32}{3}\tilde{z}-8-4l_{\bar{z}}\tilde{z}\Big{\\}}-4\tilde{z}{\mathcal{D}_{x,1}}+\delta_{\bar{z}}\Big{\\{}4\,\text{Li}_{2}(\bar{x})\tilde{x}-4-\frac{38}{9}\tilde{x}-5l_{x}^{2}\tilde{x}-2l_{\bar{x}}^{2}\tilde{x}+2l_{x}(2-7\tilde{x})+\frac{1}{\bar{x}}\big{(}20l_{x}+10l_{x}^{2}$
$\displaystyle~{}~{}-8\,\text{Li}_{2}(\bar{x})-16l_{x}l_{\bar{x}}\big{)}+l_{\bar{x}}\Big{(}\frac{32}{3}\tilde{x}-8+8l_{x}\tilde{x}\Big{)}+4\tilde{x}\zeta_{2}\Big{\\}}+{\mathcal{D}_{z,0}}\Big{\\{}\frac{32}{3}\tilde{x}-8-\frac{16}{\bar{x}}l_{x}+8l_{x}\tilde{x}-4l_{\bar{x}}\tilde{x}\Big{\\}}-4\tilde{x}{\mathcal{D}_{z,1}}\Big{]}\,.$
(10)

where $S_{n,p}(z)$= ${(-1)^{n+p-1}\over(n-1)!p!}\int_{0}^{1}{dy\over
y}\log^{n-1}(y)\log^{p}(1-zy)$, $Li_{n}(z)$ = $S_{n-1,1}(z)$, and $\zeta_{n}$
denote values of the Riemann zeta-function.

A few points are in order. As a check of our results, we have reproduced the
inclusive SF results of refs. Zijlstra and van Neerven (1992); van Neerven and
Zijlstra (1991) for the channels considered by integrating over the scaling
variable $z^{\prime}$, including all the scale-dependant terms. The terms that
are proportional to the Dirac delta functions and ‘plus’ distributions are in
complete agreement with ref. Abele et al. (2021), which uses the framework of
threshold resummation. Our results from an explicit Feynman diagrammatic
approach confirm those predictions including the next-to-SV (NSV) terms quoted
in ref. Abele et al. (2021). However, for the scale dependent term
proportional to $\delta_{\bar{x}}\delta_{\bar{z}}\log({\mu_{F}^{2}\over
Q^{2}})$, we find $\big{(}{17\over 24}C_{F}n_{F}-{193\over
48}C_{A}C_{F}\big{)}$ instead of $\big{(}{1\over 24}C_{F}n_{F}-{17\over
48}C_{A}C_{F}\big{)}$ as in ref. Abele et al. (2021).

In the following we illustrate the numerical impact of the new results for the
${\cal F}_{1,qq}$ CF for a future EIC with a COM energy $\sqrt{s}$ = 140 GeV.
We have convoluted the CFs with a set of order independent, but sufficiently
realistic model distributions for both, PDFs and FFs,

$\displaystyle xq(x,\mu_{F}^{2})$ $\displaystyle=$ $\displaystyle
0.6x^{-0.3}(1-x)^{3.5}(1+5.0x^{0.8})\,,$ $\displaystyle xg(x,\mu_{F}^{2})$
$\displaystyle=$ $\displaystyle 1.6x^{-0.3}(1-x)^{4.5}(1-0.6x^{0.3})\,.$ (11)

In Fig. 1 we present the $K$-factor, defined by the ratio $K$=NiLO/LO for
$i=1,2$ as function of $z$ for EIC after integrating $x$ between $0.1$ to
$0.8$ and $y$ between $0.1$ to $0.9$. We also show the variation of the
renormalisation scale $\mu_{R}^{2}$ in the range
$\mu_{R}^{2}\in[Q^{2}/2,2Q^{2}]$, keeping $\mu_{F}^{2}$ = $Q^{2}$ fixed. We
use $n_{F}=5$ for the number of active quarks, the value of $\alpha_{e}$=1/128
and for the strong coupling constant $\alpha_{s}(M_{Z})=0.120$ at NLO and
$\alpha_{s}(M_{Z})=1.118$ at NNLO.

Figure 1: The upper panel contains the $K$-factor as a function of $z$ for the
NLO and NNLO results, using kinematics of the EIC at $\sqrt{s}=140$ GeV and
different approximations: the SV terms at NLO (blue dashed), full NLO (blue
solid), SV terms at NNLO (red dashed) and full (non-singlet, leading color)
NNLO (red solid). The lower panel contains the corresponding uncertainty due
to the renormalisation scale variation in the range
$\mu_{R}^{2}\in[Q^{2}/2,2Q^{2}]$.

We find, for example at $z=0.5$, that the new NNLO contributions increase the
$K$-factor from $1.26$ at NLO to $1.32$ at NNLO level, while reducing the
renormalisation scale dependence from variations of $\mu_{R}^{2}$ by a factor
of 2 around $Q^{2}$ from $\\{1.51\%,-1.31\%\\}$ at NLO to
$\\{-0.32\%,0.06\%\\}$ at NNLO.

In this letter, we report the first results for the complete NNLO QCD
corrections to the quark-initiated SIDIS process, based on calculating the
complete set of Feynman diagrams including all relevant interference terms.
However, as mentioned earlier, we have selectively excluded the flavor pure-
singlet parts of the double real emission process
$\gamma^{*}+q(\bar{q})\rightarrow q(\bar{q})+q+\bar{q}$. The contributions
from these diagrams along with the detailed description of the computational
method will be presented in a future publication.

The primary challenge resides in the computation of the Feynman integrals,
especially the phase-space integrals. To obtain them analytically, we have
used state-of-the-art loop computation techniques, namely the reverse
unitarity method, the IBP reduction technique, and the method of differential
equations. The parametric integration delivers the MIs in terms of MPLs and
Nielsen polylogarithms, although their arguments are not simple, due to the
presence of two variables, $x^{\prime}$ and $z^{\prime}$. Simplification arise
from a conversion to GPLs. Our results have been exposed to a number of cross-
checks, in particular the comparison of the SV and NSV limits with ref. Abele
et al. (2021), finding complete agreement.

The new NNLO results display a moderate increase of the $K$-factor over a
large range of kinematics, except for the threshold limits, where additional
resummations need to be performed. At the same time, the residual scale
uncertainties are significantly reduced. As such, they mark a milestone in the
precision study of the SIDIS process, and soon will play a very important role
in shedding light on the physics of the hadron structure and the mechanism of
fragmentation, facilitating new determinations of PDFs and FFs.

A Mathematica notebook with all results for the CFs ${\cal F}_{I}^{(i)}$ has
been deposited at the preprint server https://arXiv.org with the sources of
this letter. They are also available from the authors upon request.

###### Acknowledgements.

_Acknowledgements:_ This work has been supported through a joint Indo-German
research grant by the Department of Science and Technology
(DST/INT/DFG/P-03/2021/dtd.12.11.21) and Deutsche Forschungsgemeinschaft
(project number 443850114). S.M. acknowledges the ERC Advanced Grant 101095857
Conformal-EIC.

## References

  * Aschenauer et al. (2019) E. C. Aschenauer, I. Borsa, R. Sassot, and C. Van Hulse, Phys. Rev. D 99, 094004 (2019), eprint 1902.10663.
  * Cacciari and Catani (2001) M. Cacciari and S. Catani, Nucl. Phys. B 617, 253 (2001), eprint hep-ph/0107138.
  * Anderle et al. (2013a) D. P. Anderle, F. Ringer, and W. Vogelsang, Phys. Rev. D 87, 034014 (2013a), eprint 1212.2099.
  * Anderle et al. (2013b) D. P. Anderle, F. Ringer, and W. Vogelsang, Phys. Rev. D 87, 094021 (2013b), eprint 1304.1373.
  * Abele et al. (2021) M. Abele, D. de Florian, and W. Vogelsang, Phys. Rev. D 104, 094046 (2021), eprint 2109.00847.
  * Abele et al. (2022) M. Abele, D. de Florian, and W. Vogelsang, Phys. Rev. D 106, 014015 (2022), eprint 2203.07928.
  * Altarelli et al. (1979) G. Altarelli, R. K. Ellis, G. Martinelli, and S.-Y. Pi, Nucl. Phys. B 160, 301 (1979).
  * Furmanski and Petronzio (1982) W. Furmanski and R. Petronzio, Z. Phys. C 11, 293 (1982).
  * Moch et al. (2004) S. Moch, J. A. M. Vermaseren, and A. Vogt, Nucl. Phys. B 688, 101 (2004), eprint hep-ph/0403192.
  * Vogt et al. (2004) A. Vogt, S. Moch, and J. A. M. Vermaseren, Nucl. Phys. B 691, 129 (2004), eprint hep-ph/0404111.
  * Almasy et al. (2012) A. A. Almasy, S. Moch, and A. Vogt, Nucl. Phys. B 854, 133 (2012), eprint 1107.2263.
  * Chen et al. (2021) H. Chen, T.-Z. Yang, H. X. Zhu, and Y. J. Zhu, Chin. Phys. C 45, 043101 (2021), eprint 2006.10534.
  * Lee et al. (2022) R. N. Lee, A. von Manteuffel, R. M. Schabinger, A. V. Smirnov, V. A. Smirnov, and M. Steinhauser, Phys. Rev. Lett. 128, 212002 (2022), eprint 2202.04660.
  * Nogueira (1993) P. Nogueira, J. Comput. Phys. 105, 279 (1993).
  * Kuipers et al. (2013) J. Kuipers, T. Ueda, J. A. M. Vermaseren, and J. Vollinga, Comput. Phys. Commun. 184, 1453 (2013), eprint 1203.6543.
  * Ruijl et al. (2017) B. Ruijl, T. Ueda, and J. Vermaseren (2017), eprint 1707.06453.
  * Anastasiou et al. (2004) C. Anastasiou, K. Melnikov, and F. Petriello, Phys. Rev. D 69, 076010 (2004), eprint hep-ph/0311311.
  * Anastasiou et al. (2012) C. Anastasiou, S. Buehler, C. Duhr, and F. Herzog, JHEP 11, 062 (2012), eprint 1208.3130.
  * Chetyrkin and Tkachov (1981) K. G. Chetyrkin and F. V. Tkachov, Nucl. Phys. B 192, 159 (1981).
  * Laporta (2000) S. Laporta, Int. J. Mod. Phys. A 15, 5087 (2000), eprint hep-ph/0102033.
  * Lee (2014) R. N. Lee, J. Phys. Conf. Ser. 523, 012059 (2014), eprint 1310.1145.
  * Matsuura et al. (1989) T. Matsuura, S. C. van der Marck, and W. L. van Neerven, Nucl. Phys. B 319, 570 (1989).
  * Zijlstra and van Neerven (1992) E. B. Zijlstra and W. L. van Neerven, Nucl. Phys. B 383, 525 (1992).
  * Rijken and van Neerven (1997) P. J. Rijken and W. L. van Neerven, Nucl. Phys. B 487, 233 (1997), eprint hep-ph/9609377.
  * Ravindran et al. (2003) V. Ravindran, J. Smith, and W. L. van Neerven, Nucl. Phys. B 665, 325 (2003), eprint hep-ph/0302135.
  * Kotikov (1991) A. V. Kotikov, Phys. Lett. B 254, 158 (1991).
  * Argeri and Mastrolia (2007) M. Argeri and P. Mastrolia, Int. J. Mod. Phys. A 22, 4375 (2007), eprint 0707.4037.
  * Remiddi (1997) E. Remiddi, Nuovo Cim. A 110, 1435 (1997), eprint hep-th/9711188.
  * Henn (2013) J. M. Henn, Phys. Rev. Lett. 110, 251601 (2013), eprint 1304.1806.
  * Ablinger et al. (2016) J. Ablinger, A. Behring, J. Blümlein, A. De Freitas, A. von Manteuffel, and C. Schneider, Comput. Phys. Commun. 202, 33 (2016), eprint 1509.08324.
  * Ablinger (2009) J. Ablinger, Master’s thesis, Linz U. (2009), eprint 1011.1176.
  * Ablinger et al. (2011) J. Ablinger, J. Blumlein, and C. Schneider, J. Math. Phys. 52, 102301 (2011), eprint 1105.6063.
  * Ablinger (2014) J. Ablinger, PoS LL2014, 019 (2014), eprint 1407.6180.
  * Duhr and Dulat (2019) C. Duhr and F. Dulat, JHEP 08, 135 (2019), eprint 1904.07279.
  * van Neerven and Zijlstra (1991) W. L. van Neerven and E. B. Zijlstra, Phys. Lett. B 272, 127 (1991).