#  DESY 16-210, DO-TH 16/29, TTK-16-46
Heavy flavour corrections to polarised and unpolarised deep-inelastic
scattering at 3-loop order††thanks: This work was supported in part by the
Austrian Science Fund (FWF) grant SFB F50 (F5009-N15), the European Commission
through contract PITN-GA-2012-316704 (HIGGSTOOLS) and by FP7 ERC Starting
Grant 257638 PAGAP.

J. Ablingera, b,c, J. Blümleinb, A. De Freitasb, A. Hasselhuhnd, A. von
Manteuffele,f, M. Rounda, C. Schneidera, F. Wißbrockg
aResearch Institute for Symbolic Computation (RISC), Johannes Kepler
University,
Altenbergerstraße 69, A–4040, Linz, Austria
bDeutsches Elektronen-Synchrotron, DESY, Platanenallee 6, D–15738 Zeuthen,
Germany
cInstitut für Theoretische Teilchenphysik und Kosmologie, RWTH Aachen
University,
D–52056 Aachen, Germany
dInstitut für Theoretische Teilchenphysik, Karlsruher Institut für Technologie
(KIT),
D–76128 Karlsruhe, Germany
ePRISMA Cluster of Excellence, Institute of Physics, J. Gutenberg University,
D–55099 Mainz,
Germany
fDepartment of Physics and Astronomy, Michigan State University, East Lansing,
MI 48824, USA
gIHES, 35 Route de Chartres, F–91440 Bures-sur-Yvette, France
E-mail

###### Abstract:

We report on progress in the calculation of 3-loop corrections to the deep-
inelastic structure functions from massive quarks in the asymptotic region of
large momentum transfer $Q^{2}$. Recently completed results allow us to obtain
the $O(a_{s}^{3})$ contributions to several heavy flavour Wilson coefficients
which enter both polarised and unpolarised structure functions for lepton-
nucleon scattering. In particular, we obtain the non-singlet contributions to
the unpolarised structure functions $F_{2}(x,Q^{2})$ and $xF_{3}(x,Q^{2})$ and
the polarised structure function $g_{1}(x,Q^{2})$. From these results we also
obtain the heavy flavour contributions to the Gross-Llewellyn-Smith and the
Bjorken sum rules.

##  1 Introduction

Deep-inelastic scattering provides a valuable way to both test the theory of
quantum chromodynamics (QCD) and to extract theory parameters from
experiments. Among these are in particular the strong coupling constant
$\alpha_{s}$ [1, *Moch:2014tta, *Alekhin:2016evh], the parton distribution
functions (PDFs) [4, 5, *Dulat:2015mca, *Harland-Lang:2014zoa, *Ball:2014uwa,
*Abramowicz:2015mha, 10] and the masses of the charm and bottom quarks [11,
12]. To harness the full potential of the experimental data, it is necessary
to have predictions at hand for which the theoretical uncertainties keep up
with the experimental accuracy. Currently, the corrections from massive quarks
are still missing for a complete next-to-next-to-leading order (NNLO) analysis
of the deep-inelastic scattering World data. They can be calculated
analytically at NNLO in the kinematic limit $Q^{2}\gg m^{2}$ [13], where
$Q^{2}$ is the virtuality of the electro-weak gauge boson and $m$ is the mass
of the heavy quark. In this paper, we report on progress in the calculation of
these heavy flavour corrections.

In Section 2 we describe the framework for the heavy flavour corrections to
deep-inelastic scattering in the limit $Q^{2}\gg m^{2}$. In this limit, the
massive operator matrix elements of the light-cone operators are the key
quantities which have not been completely computed yet. Thus, we sketch the
steps involved in their calculation. Section 3 contains several applications
of the results which have been obtained so far. In particular, we illustrate
the impact of the heavy flavour Wilson coefficients on the structure function
$F_{2}(x,Q^{2})$, $g_{1}(x,Q^{2})$ and $xF_{3}(x,Q^{2})$, as well as their
consequences for the polarised Bjorken sum rule and the Gross-Llewellyn-Smith
sum rule. Finally, we comment on our results for the massive operator matrix
element of the non-singlet operator for transversity in Section 4 and conclude
in Section 5

##  2 Framework of calculation

The structure functions of deep-inelastic scattering can in general be written
as convolutions of PDFs and Wilson coefficients (cf., e.g., [14, *Reya:1979zk,
*Blumlein:2012bf]). The Wilson coefficients carry the process-dependent
information about the particular scattering process at hand and can be
calculated in perturbation theory. They receive contributions from massless
quarks and gluons as well as from massive quarks. In the following, we will
refer to the contributions of a single massive quark species (i.e. charm or
bottom quarks) as the heavy flavour contributions. Starting at 3-loop order,
there are also contributions from diagrams with two different massive quarks.
Their calculation poses quite different challenges and is discussed elsewhere
[17, *Ablinger:2012qj, *DESY-14-019]. The heavy flavour contributions to the
structure function $F_{2}(x,Q^{2})$, for example, can be written as [13, 20,
21]

$\displaystyle{F_{2}^{\text{h}}(x,N_{F}+1,Q^{2},m^{2})=}$ $\displaystyle
x\Biggl{\\{}\sum_{k=1}^{N_{F}}e_{k}^{2}\Biggl{[}L_{q,2}^{\text{NS}}\left(x,N_{F}+1,\frac{Q^{2}}{\mu^{2}},\frac{m^{2}}{\mu^{2}}\right)\otimes\Bigl{[}f_{k}(x,\mu^{2},N_{F})+\bar{f}_{k}(x,\mu^{2},N_{F})\Bigr{]}$
$\displaystyle+\frac{1}{N_{F}}L_{q,2}^{\text{PS}}\left(x,N_{F}+1,\frac{Q^{2}}{\mu^{2}},\frac{m^{2}}{\mu^{2}}\right)\otimes\Sigma(x,\mu^{2},N_{F})$
$\displaystyle+\frac{1}{N_{F}}L_{g,2}^{\text{S}}\left(x,N_{F}+1,\frac{Q^{2}}{\mu^{2}},\frac{m^{2}}{\mu^{2}}\right)\otimes
G(x,\mu^{2},N_{F})\Biggr{]}$
$\displaystyle+e_{Q}^{2}\Biggl{[}H_{q,2}^{\text{PS}}\left(x,N_{F}+1,\frac{Q^{2}}{\mu^{2}},\frac{m^{2}}{\mu^{2}}\right)\otimes\Sigma(x,\mu^{2},N_{F})$
$\displaystyle+H_{g,2}^{\text{S}}\left(x,N_{F}+1,\frac{Q^{2}}{\mu^{2}},\frac{m^{2}}{\mu^{2}}\right)\otimes
G(x,\mu^{2},N_{F})\Biggr{]}\Biggr{\\}}\,.$ (1)

Here, $f_{k}$, $\bar{f}_{k}$ and $G$ refer to the PDFs for quarks and anti-
quarks of flavour $k$ and the gluon PDF, respectively. The singlet PDF
combination is defined as $\Sigma=\sum_{k=1}^{N_{F}}(f_{k}+\bar{f}_{k})$ and
$\mu^{2}$ is the factorisation scale. The number of light quark flavours is
denoted by $N_{F}$ and $e_{k}$ and $e_{Q}$ represent the charges of the light
and heavy quarks. The heavy flavour Wilson coefficients are denoted by
$L_{i,a}$ and $H_{i,a}$, where the $L$ and $H$ distinguish the cases where the
electro-weak gauge boson couples to a light or a heavy quark, respectively,
and the subscripts $i$ and $a$ label the initial state parton ($q,g$) and the
structure function under consideration. The symbol $\otimes$ represents the
Mellin convolution, which turns into a simple product of moments if we apply
to it a Mellin transformation

$\displaystyle M[f](N)$
$\displaystyle=\int_{0}^{1}\mathrm{d}x\,x^{N-1}f(x)\,.$ (2)

This introduces the Mellin variable $N$, to which we will refer in several
places in the following.

The massless Wilson coefficients have been calculated up to 3-loop order [22],
while the massive ones are available only semi-numerically up to 2-loop order
[23, *Laenen:1992xs, *Riemersma:1994hv].111For a precise implementation in
Mellin space see [26]. However, it was observed in [13] that the heavy-flavour
Wilson coefficients factorise into the massless Wilson coefficients and
massive operator matrix elements (OME) in the kinematic limit where $Q^{2}\gg
m^{2}$. In this limit, the heavy flavour Wilson coefficients
$\mathcal{C}^{\text{asymp}}_{i,a}$ can be written schematically as [13]

$\displaystyle\mathcal{C}^{\text{asymp}}_{i,a}\left(x,N_{F}+1,\frac{Q^{2}}{\mu^{2}},\frac{m^{2}}{\mu^{2}}\right)=\sum_{j}C_{j,a}\left(x,N_{F}+1,\frac{Q^{2}}{\mu^{2}}\right)\otimes
A_{ij,Q}\left(x,N_{F}+1,\frac{m^{2}}{\mu^{2}}\right)+\mathcal{O}\left(\frac{m^{2}}{Q^{2}}\right)\,.$
(3)

The massless Wilson coefficients $C_{i,a}$ are evaluated for $N_{F}+1$
massless flavours and the massive OMEs $A_{ij,Q}$ are calculated with $N_{F}$
massless and one massive quark. More details on the formalism can be found in
[13, 27, 20, 21]. Below, we will frequently refer to the expansion
coefficients of the OMEs in powers of $a_{s}=\tfrac{\alpha_{s}}{4\pi}$, where
the coefficient of $a_{s}^{k}$ is denoted by $A_{ij,Q}^{(k)}$.

The massive OMEs are known analytically up to 2-loop order [13, 28, 27, 29,
*Bierenbaum:2008yu, *Bierenbaum:2009zt] including linear terms in the
dimensional regulator $\varepsilon=D-4$, where $D$ is the dimensionality of
space-time in dimensional regularisation [32, *Ashmore:1972uj, *Cicuta:1972jf,
*Bollini:1972ui]. The extension of these results to 3-loop order is the topic
of our project and here we report on progress in this regard.

The massive OMEs can be extracted from calculating two-point functions with
external on-shell partons and operators, which introduce additional Feynman
rules beyond those of QCD. Our calculation follows a diagrammatic approach,
where all relevant diagrams are generated using QGRAF [36]. After inserting
the Feynman rules, we simplify the colour and Dirac algebra using FORM [37,
*Tentyukov:2007mu, *Kuipers:2012rf] and color.h [40]. In this way, we express
the diagrams in terms of roughly $10^{5}$ scalar loop integrals. In order to
compute those, we first reduce them to a smaller number of master integrals
using integration-by-parts relations [41, *Gauss:1813, *Green:1828,
*Ostrogradski:1831, 45, *Chetyrkin:1981qh, *Tkachov:1981wb]. For this task we
use the program Reduze 2 [48, *vonManteuffel:2012np].222Reduze 2 uses the
packages GiNaC [50] and Fermat [51]. A major task is then to actually
calculate the master integrals. Over the years a number of techniques have
proven to be very useful for this task:

  * •

_Higher hypergeometric functions_ [52, *Kummer:1836, *Bailey:1935,
*Slater:1966, *Gradstein:1981, *Appell:1926, *Appell:1925,
*KampeDeFeriet:1937, *Exton:1976yx, *Exton:1978, *Srivastava:1985]: After
deriving a Feynman parameter representation for the loop integrals, they can
sometimes be brought in a form which can be integrated in terms of generalised
hypergeometric (${}_{p}F_{q}$) or Appell functions that have convergent series
representations. If this is the case, we can expand them in $\varepsilon$ and
the resulting sum representations can be simplified using the summation
algorithms [63, 64, *Schneider:2004den, *Schneider:2005rec,
*Schneider:2005pol, *Schneider:2007sim, *Schneider:2010cla,
*Schneider:2010alg, *Schneider:2015lnc, *Schneider:2008jsc,
*Schneider:2016jsc, *Schneider:2016two] implemented in Sigma [63, 75, 76],
EvaluateMultiSums and SumProduction [77, *Blumlein:2012hg, *Schneider:2013,
*Schneider:2013zna] with support from HarmonicSums [81, 82, 83, 84, 85, 86]
for dealing with the nested sums that arise.

  * •

_Mellin-Barnes integrals_ [87, *Barnes:1908, *Barnes:1910, 90, *Tausk:1999vh]:
Even if a direct evaluation of the Feynman parameter integrals in terms of
hypergeometric functions is not possible, it may still be feasible to derive a
sum representation that can be simplified using Sigma and related packages: We
split up suitably formed sums of Feynman parameters at the cost of complex
contour integrals [90, *Tausk:1999vh, 92, *Czakon:2005rk, *Smirnov:2009up].
Afterwards, the Feynman parameter integrals can usually be done in terms of
Beta and Gamma functions, while the contour integrals give rise to infinite
sums via the residue theorem.

  * •

_Almkvist-Zeilberger algorithm_ [95, *Apagodu:2006, 82]: If it is possible to
confine the Mellin variable $N$ in the integrand of the Feynman parameter
integrals to only one or a small number of places, it can be advantageous to
employ the multi-variable Almkvist-Zeilberger algorithm, implemented in the
package MultiIntegrate [82]. It allows to derive recurrence relations for the
integrals which can subsequently be solved using the algorithms implemented in
Sigma.

  * •

_Differential equations and difference equations_ [97, *Caffo:1998yd,
*Caffo:1998du, *Gehrmann:1999as, 101]: Based on the integration-by-parts
relations it is possible to derive coupled systems of differential equations
for the master integrals. We translate these into coupled systems of
difference equations and uncouple them using Zürcher’s algorithm [102], which
is implemented in OreSys [103]. Once suitable initial conditions are available
(e.g. from direct calculations using other methods), the difference equations
can be solved using the package SolveCoupledSystem [104, *Ablinger:2016yjz].

More details on the application of these techniques to the calculation of
massive OMEs can be found in [101, 106, 107, 108, 109, 110, 111]. The results
obtained so far for the OMEs can be expressed in terms of nested sums. In
particular harmonic sums [112, *Blumlein:1998if, *Blumlein:2000hw,
*Blumlein:2003gb, *Blumlein:2009ta, *Blumlein:2009fz], generalised harmonic
sums [118, 84], cyclotomic sums [83] and binomially weighted sums [119,
*Davydychev:2003mv, *Weinzierl:2004bn, 85] appear both in the intermediate
steps and in the results. These structures are related to corresponding
iterated integrals [122, 118, 84, 83, 85] via an inverse Mellin
transformation. Finally, the results for the master integrals can be inserted
into the diagrams, yielding the unrenormalised expressions for the operator
matrix elements.

The renormalisation procedure for the massive OMEs at $\mathcal{O}(a_{s}^{3})$
was worked out in [20]. Since we calculate matrix elements of the local light-
cone operators, it comes at no surprise that their renormalisation involves
the anomalous dimensions of the operators. Using known results for the beta
function, mass anomalous dimensions and lower order OMEs, we can use the pole
terms of our 3-loop results to calculate the $N_{F}$-dependent part of the
anomalous dimensions.

##  3 Results for structure functions

Over the course of the recent years, a number of analytic results for the
operator matrix elements have been completed. In particular, the OMEs
$A_{qg,Q}^{(3)}$ [106], $A_{gq,Q}^{(3)}$ [109], $A_{qq,Q}^{\mathrm{PS},(3)}$
[106], $A_{Qq}^{\mathrm{PS},(3)}$ [123], $A_{qq,Q}^{\mathrm{NS},(3)}$ [124]
have been calculated. Moreover, the gluonic OME $A_{gg,Q}^{(3)}$ is known for
even values of the Mellin variable $N$ [125]. These results allow for a
numerical illustration of their impact on the heavy flavour Wilson
coefficients of different structure functions. For the structure function
$F_{2}(x,Q^{2})$ there are five different Wilson coefficients (see Eq. (1)),
each of which requires the knowledge of the 3-loop term of one OME. Since
$A_{Qg}^{(3)}$ is not yet fully known, the Wilson coefficient $H_{g,2}$ cannot
be given yet at 3-loop order at this point. Nevertheless, we can illustrate
the impact of the remaining Wilson coefficients, supplemented by the
contribution of $H_{g,2}$ up to 2-loop order for comparison.

Figure 1: Illustration of the contributions to the structure function
$F_{2}(x,Q^{2})$ from the individual Wilson coefficients. The plot shows the
values for $Q^{2}=100\,\mathrm{GeV}^{2}$ and the scale choice $\mu^{2}=Q^{2}$
with PDFs from [4]. The charm quark mass is
$m_{c}=$1.59\text{\,}\mathrm{G}\mathrm{eV}^{2}$$ in the on-shell scheme [11].
The inset enlarges the region of $0.05\leq x\leq 1$.

Figure 1 shows the heavy flavour contributions to $F_{2}(x,Q^{2})$ for a fixed
value of $Q^{2}=$100\text{\,}\mathrm{G}\mathrm{eV}^{2}$$. The biggest
contribution comes from $H_{g,2}$. Contributions to this Wilson coefficient
start at $\mathcal{O}(a_{s})$ and are the only contribution at that order. Due
to this and the fact that it involves the gluon PDF, it is quite large in the
small-$x$ region and we have scaled down the curve by a factor $20$. The
second largest contribution in the small-$x$ region is the pure-singlet Wilson
coefficient $H_{q,2}^{\text{PS}}$. It starts at $\mathcal{O}(a_{s}^{2})$ and
is negative, except for very large values of $x$ (not visible in the plot).
The large-$x$ region is dominated by the non-singlet Wilson coefficient, which
also starts at $\mathcal{O}(a_{s}^{2})$ and is negative throughout the whole
$x$-range. Here, the even moments of the non-singlet OME
$A_{qq,Q}^{\mathrm{NS},(3)}$ enter. Somewhat smaller contributions arise from
the gluon- and quark-initiated singlet Wilson coefficients
$L_{g,2}^{\text{S}}$ and $L_{q,2}^{\text{PS}}$, which start at
$\mathcal{O}(a_{s}^{2})$ and $\mathcal{O}(a_{s}^{3})$, respectively.

Figure 2: Ratio of the heavy over the light flavour contributions to the
polarised structure functions $g_{1}(x,Q^{2})$ (left panel) and
$g_{2}(x,Q^{2})$ (right panel) for different values of $Q^{2}$. Here, the PDFs
from [126] were used; plots from [127].

The non-singlet OME $A_{qq,Q}^{\mathrm{NS},(3)}$ was calculated also for odd
moments in [124]. This allows for applications to other structure functions
besides $F_{2}(x,Q^{2})$: In particular the non-singlet contribution to the
polarised structure function $g_{1}(x,Q^{2})$ was explored in [127]. Figure 2
illustrates the impact of the charm quarks on this structure function.
Compared to the contributions from massless quarks and gluons, the charm quark
constitutes around +1% to 2% to -5% of the non-singlet structure function.
This is below the current experimental accuracy but may be of interest at
future high-luminosity colliders [128, *Accardi:2012qut,
*AbelleiraFernandez:2012cc].

At the level of twist 2, the structure function $g_{2}(x,Q^{2})$ is also
related to $g_{1}(x,Q^{2})$ by the Wandzura-Wilczek relation [131],

$\displaystyle g_{2}(x,Q^{2})$
$\displaystyle=-g_{1}(x,Q^{2})+\int_{x}^{1}\frac{\mathrm{d}y}{y}g_{1}(y,Q^{2})\,.$
(4)

This allows us to show also the impact of the charm quarks there. As can be
seen in the right panel of Figure 2, the heavy quark contributions are about
$1\%$ to $4\%$ the size of the massless contribution. The pole in the plot is
due to a change of sign of $g_{2}^{\text{light}}$.

Another application of the non-singlet OME is the polarised Bjorken sum rule
[132]. It is defined as the difference of the first moments of
$g_{1}(x,Q^{2})$ in electron-proton and electron-neutron scattering,

$\displaystyle\int_{0}^{1}\mathrm{d}x\left[g_{1}^{ep}(x,Q^{2})-g_{1}^{en}(x,Q^{2})\right]$
$\displaystyle=\frac{1}{6}\left|\frac{g_{A}}{g_{V}}\right|C_{\text{pBj}}(a_{s})\,,$
(5)

where $g_{A}$ and $g_{V}$ denote the axial-vector and vector decay constants.
The perturbative coefficient $C_{\text{pBj}}$ arises from first moment of the
Wilson coefficients. For the massless contributions it has been calculated up
to $\mathcal{O}(a_{s}^{4})$ [133, 134, 135]. Our result for the non-singlet
heavy flavour Wilson coefficient leads us to the conclusion that in the limit
$Q^{2}\gg m^{2}$ the polarised Bjorken sum rule for $N_{F}$ massless quarks
and one massive quark is given completely by the massless contributions for
$N_{F}+1$ quarks: The massive non-singlet OME, which could modify the sum rule
compared to the completely massless case, has a vanishing first moment due to
fermion number conservation. However, away from the limit $Q^{2}\gg m^{2}$,
the behaviour described above no longer holds and genuine heavy flavour
contributions to the polarised Bjorken sum rule arise, cf. [136].

Figure 3: _Left panel_ : Illustration of the structure function combination
$xF_{3}^{W^{+}-W^{-}}(x,Q^{2})$ including the contributions from a massive
charm quark with a mass of $m_{c}=$1.59\text{\,}\mathrm{G}\mathrm{eV}^{2}$$ in
the on-shell scheme [11]. _Right panel_ : Ratio of the heavy quark
contributions over the massless part of the same structure function; plots
from [137].

Moreover, the same non-singlet OME enters also the charged-current structure
function combination

$\displaystyle xF_{3}^{W^{+}-W^{-}}(x,Q^{2})$
$\displaystyle=xF_{3}^{W^{+}}(x,Q^{2})+xF_{3}^{W^{-}}(x,Q^{2})\,.$ (6)

For an explanation of the notation we refer to [138, *Blumlein:2014fqa, 137].
This structure function receives contributions from two non-singlet Wilson
coefficients: On the one hand, there is $L_{q,3}^{\text{NS}}$, which describes
reactions in which the $W$ boson couples to a light quark and mediates a
flavour transition between light quark species. This is analogous to the case
of photon exchange, except for the flavour change. On the other hand, there is
$H_{q,3}^{\text{NS}}$, which describes flavour excitation reactions (e.g.
$s\to c$). Here, the $W$ boson couples to a heavy quark. This part has no
analogy in the photon-mediated case. The impact of both Wilson coefficients on
$xF_{3}^{W^{+}-W^{-}}(x,Q^{2})$ is illustrated in Figure 3. The size of the
heavy quark contribution is again of the order of about $3\%$. Also in the
charged current sector, there is a sum rule arising from the first moment of
the structure functions: The Gross-Llewellyn-Smith sum rule [140] is given by

$\displaystyle\int_{0}^{1}\mathrm{d}x\left[F_{3}^{\bar{\nu}p}(x,Q^{2})+F_{3}^{\nu
p}(x,Q^{2})\right]$ $\displaystyle=6\,C_{\text{GLS}}(a_{s})\,.$ (7)

The massless QCD corrections to $C_{\text{GLS}}$ are again known to
$\mathcal{O}(a_{s}^{4})$ [133, 141, 142] and the situation in the heavy quark
sector is similar to the polarised case: Due to the vanishing first moment of
$A_{qq,Q}^{\text{NS}}$, the influence of the heavy quark in the asymptotic
region $Q^{2}\gg m^{2}$ reduces to incrementing the number of massless
flavours ($N_{F}\to N_{F}+1$). For the power corrections in $m^{2}/Q^{2}$ up
to 2-loop order, we refer to [136].

##  4 Operator matrix element for transversity

The non-singlet results discussed so far involve the flavour non-singlet
vector or axial-vector operators. A very similar calculation can be performed
for the non-singlet tensor operator

$\displaystyle O_{q,r}^{\text{TR,NS},\mu\mu_{1}\dots\mu_{N}}(z)$
$\displaystyle=\mathrm{i}^{N-1}S\left[\bar{\psi}(z)\sigma^{\mu\mu_{1}}D^{\mu_{2}}\dots
D^{\mu_{N}}\frac{\lambda_{r}}{2}\psi(z)\right]-\text{trace terms}\,,$ (8)

where $D^{\mu}$ is the covariant derivative, $\psi$ is the quark field
operator, $\lambda_{r}$ denotes the Gell-Mann matrices of
$\mathrm{SU}(3)_{\text{flavour}}$,
$\sigma^{\mu\nu}=\tfrac{\mathrm{i}}{2}(\gamma^{\mu}\gamma^{\nu}-\gamma^{\nu}\gamma^{\mu})$
and $S[\dots]$ denotes the symmetrisation of the Lorentz indices.

The operator in Eq. (8) is related to the transversity distribution
$\Delta_{T}f_{k}$ which appears in the transversity structure function
$h_{1}(x,Q^{2})$ [143]. It can be measured in semi-inclusive deep-inelastic
scattering [144, *Jaffe:1991kp, *Jaffe:1991ra, *Cortes:1991ja] and polarised
Drell-Yan processes [148, *Collins:1992kk, *Jaffe:1993xb, *Tangerman:1994bb,
*Boer:1997nt]. The massive OME of this operator is denoted by
$A_{qq,Q}^{\mathrm{NS,TR},(3)}$ and involves the same diagrams as the vector
non-singlet operator $A_{qq,Q}^{\mathrm{NS},(3)}$. Of course, different
Feynman rules and a different projector have to be used, but the required
master integrals turn out to be the same. Our calculation allows to extract
the $N_{F}$-dependent parts of the anomalous dimensions of the transversity
operator up to 3-loop order. The result for the anomalous dimensions and the
complete OME $A_{qq,Q}^{\mathrm{NS,TR},(3)}$ are given in [124]. In principle
the heavy flavour Wilson coefficients could be constructed in a similar
fashion as for the other structure functions, but at this point the light
flavour Wilson coefficients have not been calculated to a sufficient order
yet.

##  5 Conclusions

Up to now, seven out of eight 3-loop massive operator matrix elements have
been calculated analytically for general values of the Mellin variable $N$.
Here, we discussed both the framework of calculation as well as a selected
number of applications. The analytic calculation of the required 3-loop master
integrals with additional local operator insertions led to the development and
improvement of computer-algebraic and mathematical methods and tools. The
packages Sigma, HarmonicSums, EvaluateMultiSums, SumProduction and
SolveCoupledSystem both enabled the calculation and have benefited greatly
from the challenges posed by the tasks arising from the project.

The largest contribution to the heavy flavour corrections to the structure
function $F_{2}(x,Q^{2})$ is expected to come from the Wilson coefficient
$H_{g,2}^{\text{S}}$, which is not yet known due to the missing OME
$A_{Qg}^{(3)}$. Nonetheless, a first impression of the impact of the heavy
flavour corrections can be seen from our illustrations in Section 3. We also
discussed the influence of the massive OME $A_{qq,Q}^{\mathrm{NS},(3)}$ on the
non-singlet polarised and charged-current structure functions $g_{1}(x,Q^{2})$
and $xF_{3}(x,Q^{2})$, which is of the order of a few percent in both cases.
Due to the vanishing first moment of the non-singlet OME, the polarised
Bjorken and Gross-Llewellyn-Smith sum rules are only modified by shifting
$N_{F}$ to $N_{F}+1$ if compared to the purely massless case with $N_{F}$
quark flavours in the region $Q^{2}\gg m^{2}$. Finally, we also calculated the
massive OME for the non-singlet transversity operator, which enabled us to
obtain the $N_{F}$-dependent parts of the 3-loop anomalous dimension of this
operator.

## References

  * [1] S. Bethke et al., _Workshop on Precision Measurements of $\alpha_{s}$_, 1110.0016
  * [2] S. Moch et al., _High precision fundamental constants at the TeV scale_ , 1405.4781
  * [3] S. Alekhin, J. Bümlein and S. O. Moch, _$\alpha_{s}$ from global fits of parton distribution functions_, _Mod. Phys. Lett._ A31 (2016) 1630023
  * [4] S. Alekhin, J. Blümlein and S. Moch, _The ABM parton distributions tuned to LHC data_ , _Phys. Rev._ D89 (2014) 054028, [1310.3059]
  * [5] P. Jimenez-Delgado and E. Reya, _Delineating parton distributions and the strong coupling_ , _Phys. Rev._ D89 (2014) 074049, [1403.1852]
  * [6] S. Dulat, T.-J. Hou, J. Gao, M. Guzzi, J. Huston, P. Nadolsky et al., _New parton distribution functions from a global analysis of quantum chromodynamics_ , _Phys. Rev._ D93 (2016) 033006, [1506.07443]
  * [7] L. A. Harland-Lang, A. D. Martin, P. Motylinski and R. S. Thorne, _Parton distributions in the LHC era: MMHT 2014 PDFs_ , _Eur. Phys. J._ C75 (2015) 204, [1412.3989]
  * [8] NNPDF collaboration, R. D. Ball et al., _Parton distributions for the LHC Run II_ , _JHEP_ 04 (2015) 040, [1410.8849]
  * [9] ZEUS, H1 collaboration, H. Abramowicz et al., _Combination of measurements of inclusive deep inelastic ${e^{\pm}p}$ scattering cross sections and QCD analysis of HERA data_, _Eur. Phys. J._ C75 (2015) 580, [1506.06042]
  * [10] A. Accardi et al., _A Critical Appraisal and Evaluation of Modern PDFs_ , _Eur. Phys. J._ C76 (2016) 471, [1603.08906]
  * [11] S. Alekhin, J. Blümlein, K. Daum, K. Lipka and S. Moch, _Precise charm-quark mass from deep-inelastic scattering_ , _Phys. Lett._ B720 (2013) 172–176, [1212.2355]
  * [12] S. Alekhin, J. Blümlein, S.-O. Moch and R. Plačakytė, _The new ABMP16 PDF_ , _PoS_ DIS2016 (2016) 016, [1609.03327]
  * [13] M. Buza, Y. Matiounine, J. Smith, R. Migneron and W. L. van Neerven, _Heavy quark coefficient functions at asymptotic values $Q^{2}\gg m^{2}$_, _Nucl. Phys._ B472 (1996) 611–658, [hep-ph/9601302]
  * [14] A. J. Buras, _Asymptotic Freedom in Deep Inelastic Processes in the Leading Order and Beyond_ , _Rev. Mod. Phys._ 52 (1980) 199
  * [15] E. Reya, _Perturbative Quantum Chromodynamics_ , _Phys. Rept._ 69 (1981) 195
  * [16] J. Blümlein, _The Theory of Deeply Inelastic Scattering_ , _Prog. Part. Nucl. Phys._ 69 (2013) 28–84, [1208.6087]
  * [17] J. Ablinger, J. Blümlein, S. Klein, C. Schneider and F. Wißbrock, _3-Loop Heavy Flavor Corrections to DIS with two Massive Fermion Lines_ , in _19th International Workshop on Deep-Inelastic Scattering and Related Subjects (DIS 2011) Newport News, Virginia, April 11-15, 2011_ , 2011\.  1106.5937
  * [18] J. Ablinger, J. Blümlein, A. Hasselhuhn, S. Klein, C. Schneider and F. Wißbrock, _New Heavy Flavor Contributions to the DIS Structure Function $F_{2}(x,Q^{2})$ at $\mathcal{O}(\alpha_{s}^{3})$_, _PoS_ RADCOR2011 (2011) 031, [1202.2700]
  * [19] J. Ablinger, J. Blümlein, A. De Freitas, A. Hasselhuhn, C. Schneider and F. Wißbrock. DESY-14-019
  * [20] I. Bierenbaum, J. Blümlein and S. Klein, _Mellin Moments of the $O(\alpha_{s}^{3})$ Heavy Flavor Contributions to unpolarized Deep-Inelastic Scattering at $Q^{2}\gg m^{2}$ and Anomalous Dimensions_, _Nucl. Phys._ B820 (2009) 417–482, [0904.3563]
  * [21] A. Behring, I. Bierenbaum, J. Blümlein, A. De Freitas, S. Klein and F. Wißbrock, _The logarithmic contributions to the $O(\alpha^{3}_{s})$ asymptotic massive Wilson coefficients and operator matrix elements in deeply inelastic scattering_, _Eur. Phys. J._ C74 (2014) 3033, [1403.6356]
  * [22] 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–182, [hep-ph/0504242]
  * [23] E. Laenen, S. Riemersma, J. Smith and W. L. van Neerven, _Complete $O(\alpha_{s})$ corrections to heavy flavor structure functions in electroproduction_, _Nucl. Phys._ B392 (1993) 162–228
  * [24] E. Laenen, S. Riemersma, J. Smith and W. L. van Neerven, _$O(\alpha_{s})$ corrections to heavy flavor inclusive distributions in electroproduction_, _Nucl. Phys._ B392 (1993) 229–250
  * [25] S. Riemersma, J. Smith and W. L. van Neerven, _Rates for inclusive deep inelastic electroproduction of charm quarks at HERA_ , _Phys. Lett._ B347 (1995) 143–151, [hep-ph/9411431]
  * [26] S. I. Alekhin and J. Blümlein, _Mellin representation for the heavy flavor contributions to deep inelastic structure functions_ , _Phys. Lett._ B594 (2004) 299–307, [hep-ph/0404034]
  * [27] M. Buza, Y. Matiounine, J. Smith and W. L. van Neerven, _Charm electroproduction viewed in the variable flavor number scheme versus fixed order perturbation theory_ , _Eur. Phys. J._ C1 (1998) 301–320, [hep-ph/9612398]
  * [28] M. Buza, Y. Matiounine, J. Smith and W. L. van Neerven, _$O(\alpha_{s}^{2})$ corrections to polarized heavy flavor production at $Q^{2}\gg m^{2}$_, _Nucl. Phys._ B485 (1997) 420–456, [hep-ph/9608342]
  * [29] I. Bierenbaum, J. Blümlein and S. Klein, _Two-Loop Massive Operator Matrix Elements and Unpolarized Heavy Flavor Production at Asymptotic Values $Q^{2}\gg m^{2}$_, _Nucl. Phys._ B780 (2007) 40–75, [hep-ph/0703285]
  * [30] I. Bierenbaum, J. Blümlein, S. Klein and C. Schneider, _Two-Loop Massive Operator Matrix Elements for Unpolarized Heavy Flavor Production to $O(\varepsilon)$_, _Nucl. Phys._ B803 (2008) 1–41, [0803.0273]
  * [31] I. Bierenbaum, J. Blümlein and S. Klein, _The Gluonic Operator Matrix Elements at $O(\alpha_{s}^{2})$ for DIS Heavy Flavor Production_, _Phys. Lett._ B672 (2009) 401–406, [0901.0669]
  * [32] G. ’t Hooft and M. J. G. Veltman, _Regularization and Renormalization of Gauge Fields_ , _Nucl. Phys._ B44 (1972) 189–213
  * [33] J. F. Ashmore, _A Method of Gauge Invariant Regularization_ , _Lett. Nuovo Cim._ 4 (1972) 289–290
  * [34] G. M. Cicuta and E. Montaldi, _Analytic renormalization via continuous space dimension_ , _Lett. Nuovo Cim._ 4 (1972) 329–332
  * [35] C. G. Bollini and J. J. Giambiagi, _Dimensional Renormalization: The Number of Dimensions as a Regularizing Parameter_ , _Nuovo Cim._ B12 (1972) 20–26
  * [36] P. Nogueira, _Automatic Feynman graph generation_ , _J. Comput. Phys._ 105 (1993) 279–289
  * [37] J. A. M. Vermaseren, _New features of FORM_ , math-ph/0010025
  * [38] M. Tentyukov and J. A. M. Vermaseren, _The Multithreaded version of FORM_ , _Comput. Phys. Commun._ 181 (2010) 1419–1427, [hep-ph/0702279]
  * [39] J. Kuipers, T. Ueda, J. A. M. Vermaseren and J. Vollinga, _FORM version 4.0_ , _Comput. Phys. Commun._ 184 (2013) 1453–1467, [1203.6543]
  * [40] T. van Ritbergen, A. N. Schellekens and J. A. M. Vermaseren, _Group theory factors for Feynman diagrams_ , _Int. J. Mod. Phys._ A14 (1999) 41–96, [hep-ph/9802376]
  * [41] J. L. Lagrange, _Nouvelle Recherches sur la nature et la Propagation du son_ , _Miscellanea Taurinensia_ II (1760-61) 11–172, reprinted in Serret, J. A. (ed.), _Oeuvres de Lagrange_ , vol. 1, pp. 151–316. Gauthier-Villars, Paris, 1867
  * [42] C. F. Gauß, _Theoria attactionis corporum sphaeroidicorum ellipticorum homogeneorum methodo novo tractate_ , _Commentationes scietas scientiarum Gottingensis recentiores_ II (1813) , reprinted in _Werke_ , vol. 5, pp. 3–22. Göttingen, 1867
  * [43] G. Green, _An Essay on the Application of Mathematical Analysis to the Theories of Electricity and Magnetism_ , reprinted in Ferrers, N. M. (ed.), _Mathematical papers of the late George Green_ , pp. 1–115. Macmillan, London, 1871
  * [44] M. Ostrogradski, _Note sur une intégrale qui se rencontre dans le calcul de l’attraction des sphéroïdes_ , _Mem. Ac. Sci St. Peters._ 1 (1831) 39–53
  * [45] K. G. Chetyrkin, A. L. Kataev and F. V. Tkachov, _New Approach to Evaluation of Multiloop Feynman Integrals: The Gegenbauer Polynomial x Space Technique_ , _Nucl. Phys._ B174 (1980) 345–377
  * [46] K. G. Chetyrkin and F. V. Tkachov, _Integration by Parts: The Algorithm to Calculate beta Functions in 4 Loops_ , _Nucl. Phys._ B192 (1981) 159–204
  * [47] F. V. Tkachov, _A Theorem on Analytical Calculability of Four Loop Renormalization Group Functions_ , _Phys. Lett._ B100 (1981) 65–68
  * [48] C. Studerus, _Reduze – Feynman Integral Reduction in C++_ , _Comput. Phys. Commun._ 181 (2010) 1293–1300, [0912.2546]
  * [49] A. von Manteuffel and C. Studerus, _Reduze 2 – Distributed Feynman Integral Reduction_ , 1201.4330
  * [50] C. W. Bauer, A. Frink and R. Kreckel, _Introduction to the GiNaC framework for symbolic computation within the C++ programming language_ , _J. Symb. Comput._ 33 (2000) 1, [cs/0004015]
  * [51] R. H. Lewis, “Computer Algebra System Fermat.”  http://home.bway.net/lewis
  * [52] C. F. Gauß, _Disquisitiones generales circa seriam infinitam $1+\tfrac{\alpha\beta}{1\cdot\gamma}x+\text{}\tfrac{\alpha(\alpha+1)\beta(\beta+1)}{1\cdot 2\cdot\gamma(\gamma+1)}x~{}x+\text{etc.}$_, _Commentationes societatis regiae scientarum Gottingensis recentiores_ (1812) 3–48
  * [53] E. E. Kummer, _Über die hypergeometrische Reihe $1+\tfrac{\alpha\cdot\beta}{1\cdot\gamma}x+\tfrac{\alpha(\alpha+1)\beta(\beta+1)}{1\cdot 2\cdot\gamma(\gamma+1)}x^{2}$ $+\text{}\tfrac{\alpha(\alpha+1)(\alpha+2)\beta(\beta+1)(\beta+2)}{1\cdot 2\cdot 3\cdot\gamma(\gamma+1)(\gamma+2)}x^{3}+\dots$_, _J. reine angew. Math._ 15 (1836) 39–83, 127–172
  * [54] W. N. Bailey, _Generalized Hypergeometric Series_.  No. 32 in Cambridge Tracts in Mathematics and Mathematical Physics. Cambridge University Press, London, 1935
  * [55] L. J. Slater, _Generalized hypergeometric functions_.  Cambridge University Press, Cambridge, 1966
  * [56] I. S. Gradstein and I. M. Ryshik, _Tables of series, products and integrals_.  Verlag Harri Deutsch, Thun, 1981
  * [57] P. Appell and J. Kampé de Fériet, _Fonctions hypergéométriques et hypersphériques: polynômes d’Hermite_.  Gauthier-Villars, Paris, 1926
  * [58] P. Appell, _Sur les fonctions hypergéométriques de plusieurs variables, les polynômes d’Hermite et autres fonctions sphériques dans l’hyperespace_.  Gauthier-Villars, Paris, 1925
  * [59] J. Kampé De Fériet, _La fonction hypergéométrique_.  Gauthier-Villars, Paris, 1937
  * [60] H. Exton, _Multiple Hypergeometric Functions and Applications_.  Ellis Horwood, Chichester, 1976
  * [61] H. Exton, _Handbook of hypergeometric integrals_.  Ellis Horwood, Chichester, 1978
  * [62] H. M. Srivastava and P. W. Karlsson, _Multiple Gaussian hypergeometric series_.  Ellis Horwood, Chichester, 1985
  * [63] C. Schneider, _Symbolic Summation in Difference Fields_.  PhD thesis, RISC, J. Kepler University Linz, 2001
  * [64] M. Karr, _Summation in Finite Terms_ , _J. ACM_ 28 (1981) 305–350
  * [65] C. Schneider, _A Collection of Denominator Bounds to Solve Parameterized Linear Difference Equations in $\Pi\Sigma$-Extensions_, _An. Univ. Timisoara Ser. Mat.-Inform._ 42 (2004) 163–179
  * [66] C. Schneider, _Solving Parameterized Linear Difference Equations in Terms of Indefinite Nested Sums and Products_ , _J. Differ. Equations Appl._ 11 (2005) 799–821
  * [67] C. Schneider, _Degree Bounds to Find Polynomial Solutions of Parameterized Linear Difference Equations in $\Pi\Sigma$-Fields_, _Appl. Algebra Engrg. Comm. Comput._ 16 (2005) 1–32
  * [68] C. Schneider, _Simplifying Sums in $\Pi\Sigma^{*}$-Extensions_, _J. Algebra Appl._ 6 (2007) 415–441
  * [69] C. Schneider, _A Symbolic Summation Approach to Find Optimal Nested Sum Representations_ , in _Motives, Quantum Field Theory, and Pseudodifferential Operators_ (A. Carey, D. Ellwood, S. Paycha and S. Rosenberg, eds.), vol. 12 of _Clay Mathematics Proceedings_ , pp. 285–308, Amer. Math. Soc, 2010.  0904.2323
  * [70] C. Schneider, _Parameterized Telescoping Proves Algebraic Independence of Sums_ , _Ann. Comb._ 14 (2010) 533–552, [0808.2596]
  * [71] C. Schneider, _Fast Algorithms for Refined Parameterized Telescoping in Difference Fields_ , in _Computer Algebra and Polynomials, Applications of Algebra and Number Theory_ (J. Gutierrez, J. Schicho and M. Weimann, eds.), no. 8942 in Lecture Notes in Computer Science, pp. 157–191.  Springer, 2015.  1307.7887
  * [72] C. Schneider, _A refined difference field theory for symbolic summation_ , _J. Symbolic Comput._ 43 (2008) 611–644, [0808.2543]
  * [73] C. Schneider, _A difference ring theory for symbolic summation_ , _J. Symbolic Comput._ 72 (2016) 82–127, [1408.2776]
  * [74] C. Schneider, _Summation Theory II: Characterizations of $R\Pi\Sigma$-extensions and algorithmic aspects_, _to appear in J. Symbolic Comput._ (2016) , [1603.04285]
  * [75] C. Schneider, _Symbolic Summation Assists Combinatorics_ , _Sem. Lothar. Combin._ 56 (2007) article 56b
  * [76] C. Schneider, _Simplifying Multiple Sums in Difference Fields_ , in _Computer Algebra in Quantum Field Theory: Integration, Summation and Special Functions_ (C. Schneider and J. Blümlein, eds.), pp. 325–360.  Springer, Vienna, 2013.  1304.4134.  DOI
  * [77] J. Ablinger, J. Blümlein, S. Klein and C. Schneider, _Modern Summation Methods and the Computation of 2- and 3-loop Feynman Diagrams_ , _Nucl. Phys. Proc. Suppl._ 205-206 (2010) 110–115, [1006.4797]
  * [78] J. Blümlein, A. Hasselhuhn and C. Schneider, _Evaluation of Multi-Sums for Large Scale Problems_ , _PoS_ RADCOR2011 (2011) 032, [1202.4303]
  * [79] C. Schneider, _Symbolic Summation in Difference Fields and Its Application in Particle Physics_ , _Computer Algebra Rundbrief_ 53 (2013) 8–12
  * [80] C. Schneider, _Modern Summation Methods for Loop Integrals in Quantum Field Theory: The Packages Sigma, EvaluateMultiSums and SumProduction_ , _J. Phys. Conf. Ser._ 523 (2014) 012037, [1310.0160]
  * [81] J. Ablinger, _A Computer Algebra Toolbox for Harmonic Sums Related to Particle Physics_ , Diploma thesis, J. Kepler University Linz, 2009
  * [82] J. Ablinger, _Computer Algebra Algorithms for Special Functions in Particle Physics_.  PhD thesis, J. Kepler University Linz, 2012.  1305.0687
  * [83] J. Ablinger, J. Blümlein and C. Schneider, _Harmonic Sums and Polylogarithms Generated by Cyclotomic Polynomials_ , _J. Math. Phys._ 52 (2011) 102301, [1105.6063]
  * [84] J. Ablinger, J. Blümlein and C. Schneider, _Analytic and Algorithmic Aspects of Generalized Harmonic Sums and Polylogarithms_ , _J. Math. Phys._ 54 (2013) 082301, [1302.0378]
  * [85] J. Ablinger, J. Blümlein, C. G. Raab and C. Schneider, _Iterated Binomial Sums and their Associated Iterated Integrals_ , _J. Math. Phys._ 55 (2014) 112301, [1407.1822]
  * [86] J. Ablinger, _The package HarmonicSums: Computer Algebra and Analytic aspects of Nested Sums_ , _PoS_ LL2014 (2014) 019, [1407.6180]
  * [87] H. Mellin, _Om definita integraler, hvilka för obegränsadt växande värden af vissa heltaliga parametrar hafva till gränser hypergeometriska funktioner af särskilda ordningar_ , _Acta Societatis Scientiarum Fennicae_ XX (1885) 1–39
  * [88] E. W. Barnes, _A New Development of the Theory of the Hypergeometric Functions_ , _Proc. London Math. Soc._ s2-6 (1908) 141–177
  * [89] E. W. Barnes, _A Transformation of Generalised Hypergeometric Series_ , _Quart. J. Math._ 41 (1910) 136–140
  * [90] V. A. Smirnov, _Analytical result for dimensionally regularized massless on shell double box_ , _Phys. Lett._ B460 (1999) 397–404, [hep-ph/9905323]
  * [91] J. B. Tausk, _Nonplanar massless two loop Feynman diagrams with four on-shell legs_ , _Phys. Lett._ B469 (1999) 225–234, [hep-ph/9909506]
  * [92] V. A. Smirnov, _Feynman integral calculus_.  Springer, Berlin, 2006
  * [93] M. Czakon, _Automatized analytic continuation of Mellin-Barnes integrals_ , _Comput. Phys. Commun._ 175 (2006) 559–571, [hep-ph/0511200]
  * [94] A. V. Smirnov and V. A. Smirnov, _On the Resolution of Singularities of Multiple Mellin-Barnes Integrals_ , _Eur. Phys. J._ C62 (2009) 445–449, [0901.0386]
  * [95] G. Almkvist and D. Zeilberger, _The method of differentiating under the integral sign_ , _J. Symb. Comp._ 10 (1990) 571 – 591
  * [96] M. Apagodu and D. Zeilberger, _Multi-variable Zeilberger and Almkvist–Zeilberger algorithms and the sharpening of Wilf–Zeilberger theory_ , _Adv. Appl. Math._ 37 (2006) 139 – 152
  * [97] A. V. Kotikov, _Differential equations method: New technique for massive Feynman diagrams calculation_ , _Phys. Lett._ B254 (1991) 158–164
  * [98] M. Caffo, H. Czyż, S. Laporta and E. Remiddi, _Master equations for master amplitudes_ , _Acta Phys. Polon._ B29 (1998) 2627–2635, [hep-th/9807119]
  * [99] M. Caffo, H. Czyż, S. Laporta and E. Remiddi, _The Master differential equations for the two loop sunrise selfmass amplitudes_ , _Nuovo Cim._ A111 (1998) 365–389, [hep-th/9805118]
  * [100] T. Gehrmann and E. Remiddi, _Differential equations for two loop four point functions_ , _Nucl. Phys._ B580 (2000) 485–518, [hep-ph/9912329]
  * [101] J. Ablinger, A. Behring, J. Blümlein, A. De Freitas, A. von Manteuffel and C. Schneider, _Calculating Three Loop Ladder and V-Topologies for Massive Operator Matrix Elements by Computer Algebra_ , _Comput. Phys. Commun._ 202 (2016) 33–112, [1509.08324]
  * [102] B. Zürcher, _Rationale Normalformen von pseudo-linearen Abbildungen_ , Diploma thesis, ETH Zürich, 1994
  * [103] S. Gerhold, _Uncoupling Systems of Linear Ore Operator Equations_ , Diploma thesis, RISC, J. Kepler University Linz, 2002
  * [104] J. Ablinger, J. Bluemlein, A. de Freitas and C. Schneider, _A toolbox to solve coupled systems of differential and difference equations_ , _PoS_ RADCOR2015 (2015) 060, [1601.01856]
  * [105] J. Ablinger, A. Behring, J. Bluemlein, A. de Freitas and C. Schneider, _Algorithms to solve coupled systems of differential equations in terms of power series_ , _PoS_ LL2016 (2016) 005, [1608.05376]
  * [106] J. Ablinger, J. Blümlein, S. Klein, C. Schneider and F. Wißbrock, _The $O(\alpha_{s}^{3})$ Massive Operator Matrix Elements of $O(n_{f})$ for the Structure Function $F_{2}(x,Q^{2})$ and Transversity_, _Nucl. Phys._ B844 (2011) 26–54, [1008.3347]
  * [107] J. Blümlein, A. Hasselhuhn, S. Klein and C. Schneider, _The $O(\alpha_{s}^{3}n_{f}T_{F}^{2}C_{A,F})$ Contributions to the Gluonic Massive Operator Matrix Elements_, _Nucl. Phys._ B866 (2013) 196–211, [1205.4184]
  * [108] J. Ablinger, J. Blümlein, A. Hasselhuhn, S. Klein, C. Schneider and F. Wißbrock, _Massive 3-loop Ladder Diagrams for Quarkonic Local Operator Matrix Elements_ , _Nucl. Phys._ B864 (2012) 52–84, [1206.2252]
  * [109] J. Ablinger, J. Blümlein, A. De Freitas, A. Hasselhuhn, A. von Manteuffel, M. Round et al., _The Transition Matrix Element $A_{gq}(N)$ of the Variable Flavor Number Scheme at $O(\alpha_{s}^{3})$_, _Nucl. Phys._ B882 (2014) 263–288, [1402.0359]
  * [110] J. Ablinger, J. Blümlein, A. De Freitas, A. Hasselhuhn, A. von Manteuffel, M. Round et al., _The $O(\alpha_{s}^{3}T_{F}^{2})$ Contributions to the Gluonic Operator Matrix Element_, _Nucl. Phys._ B885 (2014) 280–317, [1405.4259]
  * [111] J. Ablinger, J. Blümlein, C. Raab, C. Schneider and F. Wißbrock, _Calculating Massive 3-loop Graphs for Operator Matrix Elements by the Method of Hyperlogarithms_ , _Nucl. Phys._ B885 (2014) 409–447, [1403.1137]
  * [112] J. A. M. Vermaseren, _Harmonic sums, Mellin transforms and integrals_ , _Int. J. Mod. Phys._ A14 (1999) 2037–2076, [hep-ph/9806280]
  * [113] J. Blümlein and S. Kurth, _Harmonic sums and Mellin transforms up to two loop order_ , _Phys. Rev._ D60 (1999) 014018, [hep-ph/9810241]
  * [114] J. Blümlein, _Analytic continuation of Mellin transforms up to two loop order_ , _Comput. Phys. Commun._ 133 (2000) 76–104, [hep-ph/0003100]
  * [115] J. Blümlein, _Algebraic relations between harmonic sums and associated quantities_ , _Comput. Phys. Commun._ 159 (2004) 19–54, [hep-ph/0311046]
  * [116] J. Blümlein, _Structural Relations of Harmonic Sums and Mellin Transforms up to Weight $w=5$_, _Comput. Phys. Commun._ 180 (2009) 2218–2249, [0901.3106]
  * [117] J. Blümlein, _Structural Relations of Harmonic Sums and Mellin Transforms at Weight $w=6$_, in _Motives, Quantum Field Theory, and Pseudodifferential Operators_ (A. Carey, D. Ellwood, S. Paycha and S. Rosenberg, eds.), vol. 12 of _Clay Mathematics Proceedings_ , p. 167, Amer. Math. Soc, 2010.  0901.0837
  * [118] S. Moch, P. Uwer and S. Weinzierl, _Nested sums, expansion of transcendental functions and multiscale multiloop integrals_ , _J. Math. Phys._ 43 (2002) 3363–3386, [hep-ph/0110083]
  * [119] J. Fleischer, A. V. Kotikov and O. L. Veretin, _Analytic two loop results for selfenergy type and vertex type diagrams with one nonzero mass_ , _Nucl. Phys._ B547 (1999) 343–374, [hep-ph/9808242]
  * [120] A. I. Davydychev and M. Y. Kalmykov, _Massive Feynman diagrams and inverse binomial sums_ , _Nucl. Phys._ B699 (2004) 3–64, [hep-th/0303162]
  * [121] S. Weinzierl, _Expansion around half integer values, binomial sums and inverse binomial sums_ , _J. Math. Phys._ 45 (2004) 2656–2673, [hep-ph/0402131]
  * [122] E. Remiddi and J. A. M. Vermaseren, _Harmonic polylogarithms_ , _Int. J. Mod. Phys._ A15 (2000) 725–754, [hep-ph/9905237]
  * [123] J. Ablinger, A. Behring, J. Blümlein, A. De Freitas, A. von Manteuffel and C. Schneider, _The 3-loop pure singlet heavy flavor contributions to the structure function $F_{2}(x,Q^{2})$ and the anomalous dimension_, _Nucl. Phys._ B890 (2014) 48–151, [1409.1135]
  * [124] J. Ablinger, A. Behring, J. Blümlein, A. De Freitas, A. Hasselhuhn, A. von Manteuffel et al., _The 3-Loop Non-Singlet Heavy Flavor Contributions and Anomalous Dimensions for the Structure Function $F_{2}(x,Q^{2})$ and Transversity_, _Nucl. Phys._ B886 (2014) 733–823, [1406.4654]
  * [125] J. Ablinger et al. DESY-15-112
  * [126] J. Blümlein and H. Böttcher, _QCD Analysis of Polarized Deep Inelastic Scattering Data_ , _Nucl. Phys._ B841 (2010) 205–230, [1005.3113]
  * [127] A. Behring, J. Blümlein, A. De Freitas, A. von Manteuffel and C. Schneider, _The 3-Loop Non-Singlet Heavy Flavor Contributions to the Structure Function $g_{1}(x,Q^{2})$ at Large Momentum Transfer_, _Nucl. Phys._ B897 (2015) 612–644, [1504.08217]
  * [128] D. Boer et al., _Gluons and the quark sea at high energies: Distributions, polarization, tomography_ , 1108.1713
  * [129] A. Accardi et al., _Electron Ion Collider: The Next QCD Frontier - Understanding the glue that binds us all_ , 1212.1701
  * [130] LHeC Study Group collaboration, J. L. Abelleira Fernandez et al., _A Large Hadron Electron Collider at CERN: Report on the Physics and Design Concepts for Machine and Detector_ , _J. Phys._ G39 (2012) 075001, [1206.2913]
  * [131] S. Wandzura and F. Wilczek, _Sum Rules for Spin Dependent Electroproduction: Test of Relativistic Constituent Quarks_ , _Phys. Lett._ B72 (1977) 195
  * [132] J. D. Bjorken, _Inelastic Scattering of Polarized Leptons from Polarized Nucleons_ , _Phys. Rev._ D1 (1970) 1376–1379
  * [133] P. A. Baikov, K. G. Chetyrkin and J. H. Kühn, _Adler Function, Bjorken Sum Rule, and the Crewther Relation to Order $\alpha_{s}^{4}$ in a General Gauge Theory_, _Phys. Rev. Lett._ 104 (2010) 132004, [1001.3606]
  * [134] S. A. Larin, _The singlet contribution to the Bjorken sum rule for polarized deep inelastic scattering_ , _Phys. Lett._ B723 (2013) 348–350, [1303.4021]
  * [135] P. A. Baikov, K. G. Chetyrkin and J. H. Kühn, _Massless Propagators, $R(s)$ and Multiloop QCD_, _Nucl. Part. Phys. Proc._ 261-262 (2015) 3–18, [1501.06739]
  * [136] J. Blümlein, G. Falcioni and A. De Freitas, _The Complete $O(\alpha_{s}^{2})$ Non-Singlet Heavy Flavor Corrections to the Structure Functions $g_{1,2}^{ep}(x,Q^{2})$, $F_{1,2,L}^{ep}(x,Q^{2})$, $F_{1,2,3}^{\nu(\bar{\nu})}(x,Q^{2})$ and the Associated Sum Rules_, _Nucl. Phys._ B910 (2016) 568–617, [1605.05541]
  * [137] A. Behring, J. Blümlein, A. De Freitas, A. Hasselhuhn, A. von Manteuffel and C. Schneider, _$O(\alpha_{s}^{3})$ heavy flavor contributions to the charged current structure function $xF_{3}(x,Q^{2})$ at large momentum transfer_, _Phys. Rev._ D92 (2015) 114005, [1508.01449]
  * [138] J. Blümlein, A. Hasselhuhn, P. Kovačíková and S. Moch, _$O(\alpha_{s})$ Heavy Flavor Corrections to Charged Current Deep-Inelastic Scattering in Mellin Space_, _Phys. Lett._ B700 (2011) 294–304, [1104.3449]
  * [139] J. Blümlein, A. Hasselhuhn and T. Pfoh, _The $O(\alpha_{s}^{2})$ heavy quark corrections to charged current deep-inelastic scattering at large virtualities_, _Nucl. Phys._ B881 (2014) 1–41, [1401.4352]
  * [140] D. J. Gross and C. H. Llewellyn Smith, _High-energy neutrino–nucleon scattering, current algebra and partons_ , _Nucl. Phys._ B14 (1969) 337–347
  * [141] P. A. Baikov, K. G. Chetyrkin and J. H. Kühn, _Adler Function, DIS sum rules and Crewther Relations_ , _Nucl. Phys. Proc. Suppl._ 205-206 (2010) 237–241, [1007.0478]
  * [142] P. A. Baikov, K. G. Chetyrkin, J. H. Kühn and J. Rittinger, _Adler Function, Sum Rules and Crewther Relation of Order $O(\alpha_{s}^{4})$: the Singlet Case_, _Phys. Lett._ B714 (2012) 62–65, [1206.1288]
  * [143] V. Barone, A. Drago and P. G. Ratcliffe, _Transverse polarisation of quarks in hadrons_ , _Phys. Rept._ 359 (2002) 1–168, [hep-ph/0104283]
  * [144] J. P. Ralston and D. E. Soper, _Production of Dimuons from High-Energy Polarized Proton Proton Collisions_ , _Nucl. Phys._ B152 (1979) 109
  * [145] R. L. Jaffe and X.-D. Ji, _Chiral odd parton distributions and polarized Drell-Yan_ , _Phys. Rev. Lett._ 67 (1991) 552–555
  * [146] R. L. Jaffe and X.-D. Ji, _Chiral odd parton distributions and Drell-Yan processes_ , _Nucl. Phys._ B375 (1992) 527–560
  * [147] J. L. Cortes, B. Pire and J. P. Ralston, _Measuring the transverse polarization of quarks in the proton_ , _Z. Phys._ C55 (1992) 409–416
  * [148] X. Artru and M. Mekhfi, _Transversely Polarized Parton Densities, their Evolution and their Measurement_ , _Z. Phys._ C45 (1990) 669
  * [149] J. C. Collins, _Fragmentation of transversely polarized quarks probed in transverse momentum distributions_ , _Nucl. Phys._ B396 (1993) 161–182, [hep-ph/9208213]
  * [150] R. L. Jaffe and X.-D. Ji, _Novel quark fragmentation functions and the nucleon’s transversity distribution_ , _Phys. Rev. Lett._ 71 (1993) 2547–2550, [hep-ph/9307329]
  * [151] R. D. Tangerman and P. J. Mulders, _Polarized twist-three distributions $g_{T}$ and $h_{L}$ and the role of intrinsic transverse momentum_, hep-ph/9408305
  * [152] D. Boer and P. J. Mulders, _Time reversal odd distribution functions in leptoproduction_ , _Phys. Rev._ D57 (1998) 5780–5786, [hep-ph/9711485]


# Universality class of epidemic percolation transitions driven by random
walks

Mohadeseh Feshanjerdi<EMAIL_ADDRESS>Department of Physics,
University of Tehran, P. O. Box 14395-547, Tehran, Iran  Abbas Ali Saberi
(corresponding author<EMAIL_ADDRESS>Department of Physics, University of
Tehran, P. O. Box 14395-547, Tehran, Iran  Institut für Theoretische Physik,
Universität zu Köln, Zülpicher Str. 77, 50937 Köln, Germany

###### Abstract

Inspired by the recent viral epidemic outbreak and its consequent worldwide
pandemic, we devise a model to capture the dynamics and the universality of
the spread of such infectious diseases. The transition from a pre-critical to
the post-critical phase is modeled by a percolation problem driven by random
walks on a two-dimensional lattice with an extra average number $\rho$ of
nonlocal links per site. Using the finite-size scaling analysis, we find that
the effective exponents of the percolation transitions as well as the
corresponding time thresholds, extrapolated to the infinite system size, are
$\rho$-dependent. We argue that the $\rho$-dependence of our estimated
exponents represents a crossover-type behavior caused by the finite-size
effects between the two limiting regimes of the system. We also find that the
universal scaling functions governing the critical behavior in every single
realization of the model, can be well described by the theory of extreme
values for the maximum jumps in the order parameter and by the central limit
theorem for the transition threshold.

##  I Introduction

Percolation is one of the simplest models in probability theory which provides
a suitable platform to formulate and model various natural phenomena that
exhibit geometric phase transitions with universal characteristics r1 ; r2 ;
r2-1 ; r2-2 . Percolation has been applied to describe a wide range of
critical behavior such as, among many others, flow through porous media for
connectivity r3 ; r4 , networks r5 ; r6 ; r7 , magnetic models r8 ; r9 ; r10 ;
r11 ; r11-1 ; r11-2 ; r11-3 ; r11-4 , colloids r12 ; r13 , growth models r14 ,
topography of planets r14-1 ; r14-2 ; r14-3  and epidemic models r14-4 . The
concept of universality lets many microscopically different physical systems
exhibit the same critical behavior with quantitatively identical features,
assigned by a set of critical exponents. Thanks to the emergence of the
powerful conformal field theories (CFTs) in two dimensions (2D), many exact
results have been obtained for the percolation model in 2D r14-5 .

On the other hand, the theory of random walks, also known as Brownian motion
in the scaling limit, plays an essential role in many stochastic processes and
challenging problems in probability and statistical physics r14-6 ; r14-7 ;
r14-8 . In two dimensions, random walks and percolation at its critical point
share various interesting properties other than their conformal invariant
property. For example, the external perimeters of a critical percolation
cluster and a planar random walk in the continuum limit are believed to be
described by the scaling limit of planar self-avoiding random walks (SAWs)
r14-9 . Moreover, the trace of a smart kinetic self-avoiding random walk r22
(i.e., the trace of a random walker that never crosses itself and never gets
trapped) produces statistically the same fractal path in accord with the
perimeter (or the hull) of a critical percolation cluster in 2D.

Here we combine these two fundamental theories with some modifications to
develop a model to study the universal behavior and dynamics of the spread of
infectious diseases over a 2D square lattice possessing a variable percentage
$\rho$ of nonlocal extra links. In the absence of the nonlocal links i.e., the
trivial case with $\rho=0$, there will be a random walker (infected person)
which randomly visits the lattice sites (the people prone to disease) and the
average linear size of the visited (infected) sites grows with the number $N$
of walks as $\sim\sqrt{N}$. For a square lattice of linear size $L$ and
periodic boundary conditions, this means that the percolation of the infected
sites is expected to happen after $t_{c}=N\sim L^{2}$ steps which diverges
when the system size $L$ goes to infinity. In the following, we will indeed
study the dynamics of the model on square lattices with nonzero $\rho$.
We use finite-size scaling (FSS) hypothesis r22-1  to precisely locate the
critical thresholds in the limit of the infinite system size
$L\rightarrow\infty$ and extract various critical exponents of the model for
$\rho>0$. Instead of performing a traditional FSS analysis at the critical
phase transition point, we rather follow the universal framework of a FSS
analysis recently developed in r26  based on the statistics and critical
scaling of the size of the largest gap in the order parameter. The key-point
of r26  is the following: Since the size of the largest gap (within a time-
series of the evolving largest cluster of the infected sites) can be regarded
as an extreme value of the percolation process, the corresponding scaling
function is believed to be governed by extreme-value statistics. Such
intuition has led to unify continuous and discontinuous percolation
transitions by identifying the universal critical scaling functions as the
extreme-value Gumbel distribution r26 . This framework also provides an
alternative approach to estimate various critical exponents and fractal
dimensions which determine the universality class and the type (continuous or
discontinuous) of the transition.

This paper is structured as follows: The following Section II gives the
details of our devised model and the dynamics. Section III goes over the
theoretical backgrounds and definition of various percolation observables that
we are going to measure. Section IV presents our numerical results and their
comparison with our analytical arguments. A full list of the effective
critical exponents estimated for finite-size systems for various choices of
$\rho$ is also summarized. The final Section V concludes and provides
suggestions for future research.

##  II The Model

Consider a square lattice of size $L\times L$ with periodic boundary
conditions in both directions. Initially, each of the $N=L^{2}$ sites (or
susceptible individuals) has four nearest neighbors which are connected by
local links of a lattice constant $a=1$ length. For a given $\rho=$(number of
extra nonlocal links)/$N$, we select a pair of randomly chosen sites and
connect them together by an extra nonlocal link which effectively makes them
nearest neighbors. This procedure continues until the desired value $\rho$ is
reached. The lattice with the only parameter $\rho$ is now ready to study our
stochastic epidemic modeling.

To initiate the dynamics of our model, a random walker (or the infected
individual) starts its journey from a randomly chosen site on the lattice at
$t=0$. For the ease of notation in the remainder of this paper, we refer to
time $t$ as the number of random walk steps over $N$, i.e., $t:=t/N$. At each
time step, the random walker moves to one of its randomly chosen nearest
neighboring sites (either the local or nonlocal ones). Every visited site by
the random walker is marked as an infected site upon visiting (no matter how
many times it is visited) and remains infected forever. A cluster of infected
sites is identified as a set of nearest neighboring infected sites which are
connected only by local links of length $a=1$. At early times $t<t_{c}$, there
will only appear small-size disjoint clusters all over the lattice (epidemic
outbreak phase) and they gradually merge into a single super cluster at
$t>t_{c}$ which spans the entire lattice (pandemic phase). Figure 1
illustrates these two phases (left panels) and steps taken in each phase
(right panels) on a lattice of size $L=40$ and $\rho=0.8$. Characterizing the
nature of this phase transition at $t=t_{c}$ and its universality class by
determining the corresponding critical exponents as well as the critical time
$t_{c}$ where the epidemic prevalence peaks, are the main subjects of our
present study.

Figure 1: Illustration of the infected clusters in the epidemic outbreak
(left-up panel) and pandemic (left-down panel) phases at time $t<t_{c}$ and
$t>t_{c}$, respectively, on a square lattice of linear size $L=40$ and
$\rho=0.8$ with periodic boundary conditions in both directions. Every
disjoint infected cluster is shown by a different color. Uninfected sites are
shown in black. The corresponding steps taken by the random walker in the time
interval [0, $t$] for each phase are shown in the right panels.

##  III theoretical background and Definition of the observables

The relative size $s(t)$ of the largest cluster at time $t$ is known as the
standard order parameter in percolation transitions. Once the random walker
starts moving, the infected clusters form and evolve in time. We monitor the
evolution of $s(t)$ at every time step $\Delta t$, and record its maximum
jump,

$\Delta\equiv\max_{t}[s(t+\Delta t)-s(t)],$ (1)

and the time of incidence $t_{c}$, for every single realization. We also
record the size of the largest cluster just at and after the critical
threshold $t_{c}$ and denote them by $S_{c}^{-}=Ns(t_{c})$ and
$S_{c}^{+}=Ns(t_{c}+\Delta t)$, respectively. All these defined quantities
i.e., $\Delta$, $t_{c}$, $S_{c}^{-}$ and $S_{c}^{+}$, are in fact random
variables whose values are sample dependent. Therefore, our observables in
this paper are the probability distribution function, the average and the root
mean square of these variables. We consider $M=10^{4}$ number of independent
realizations for every given system size $L$. We will then start our analysis
with a set of measured variables as

$Y=\\{Y_{1},Y_{2},...,Y_{M}\\}_{Y=\,\Delta,\,t_{c},\,S_{c}^{-}\mathrm{and}\,S_{c}^{+}},$
(2)

with the average

$\bar{Y}(L)=\langle Y\rangle=\dfrac{1}{M}\sum_{i=1}^{M}\,\,Y_{i},$ (3)

and the root mean square

$\chi_{Y}(L)=\sqrt{\langle Y_{i}^{2}\rangle}.$ (4)

Figure 2: Plots for the computation of critical gap exponents and the fractal
dimensions for a wide range of nonlocal link densities
$\rho=10^{-3},10^{-2},5\times 10^{-2},10^{-1},8\times 10^{-1},2,10$, and
$10^{2}$. The average size of the largest gap $\bar{\Delta}$ (a) and its
fluctuation $\chi_{\Delta}$ (b), $\bar{t}_{c}(L)-t_{c}(\infty)$ (c) and
$P_{t}(t_{c}^{*},L)$ (d) as functions of $L$ give the exponents $\beta_{1}$,
$\beta_{2}$, $\nu_{1}$ and $\nu_{2}$, respectively. The average size of the
largest clusters $\bar{S}^{-}$ (e) and $\bar{S}^{+}$ (g) and their
corresponding fluctuations $\chi_{S^{-}}$ (f) and $\chi_{S^{+}}$ (h) as
functions of $L$ to extract the fractal dimensions. One of the solid lines
show the best fit to our data for $\rho=10^{-3}$ (cross symbols) and the other
shows comparison with the exact predictions available for the standard 2D
percolation which is believed to explain the limit $\rho\rightarrow\infty$ of
our model.

###  III.1 Critical gap exponents

Finite-size scaling hypothesis provides a scaling framework for our
observables as a function of the system size $L$, giving rise to a number of
scaling and critical exponents which are able to classify the model r25 ; r26
. Various gap exponents and fractal dimensions can then be defined by the
following proposed scaling behaviors,

$\bar{\Delta}\,(L)\sim L^{-\beta_{1}},\,\,\,\,\,\chi_{\Delta}\sim
L^{-\beta_{2}},$ (5) $\bar{S}{{}_{c}^{-}}(L)\sim
L^{d_{f_{1}}^{-}},\,\,\,\,\,\chi_{{S}{{}_{c}^{-}}}\sim L^{d_{f_{2}}^{-}},$ (6)
$\bar{S}{{}_{c}^{+}}(L)\sim
L^{d_{f_{1}}^{+}},\,\,\,\,\,\chi_{{S}{{}_{c}^{+}}}\sim L^{d_{f_{2}}^{+}},$ (7)
$\bar{t}{{}_{c}}(L)-\bar{t}{{}_{c}}(\infty)\sim
L^{-1/\nu_{1}},\,\,\,\,\,\chi_{t_{c}}\sim L^{-1/\nu_{2}}.$ (8)

The vanishing exponent $\beta_{1}=0$ can be used to determine if a percolation
transition is discontinuous. However, models with $0<\beta_{1}<1$ can be
either continuous or discontinuous r26-1 . It can be also shown that
$\beta_{1}=\beta_{2}$, and the critical largest clusters share the same
fractal dimensions: $d_{f_{1}}^{-}=d_{f_{1}}^{+}=d_{f_{2}}^{-}=d_{f_{2}}^{+}$.
Moreover, the scaling relation $\beta_{1}=d-d_{f_{1}}$ holds for the gap
exponents, which is similar to the known relation $\beta/\nu=d-d_{f}$ in terms
of the standard order parameter and correlation exponents i.e., $\beta$ and
$\nu$, respectively r26 . $\bar{t}_{c}(\infty)$ in Eq. (8) denotes for the
average critical threshold of the model in the thermodynamic limit where the
system size goes to infinity i.e., in the limit $1/L\rightarrow 0$.

Table 1: The summarized asymptotic critical thresholds $\bar{t}_{c}(\infty)$, the effective critical gap exponents $\nu_{1}$, $\nu_{2}$, $\beta_{1}$, $\beta_{2}$, the fractal dimensions $d_{f_{1}}^{\pm}$, $d_{f_{2}}^{\pm}$, and the effective Fisher exponent $\tau$, estimated for various percentage $\rho$ of nonlocal links. The last row presents the same quantities known for the standard 2D site-percolation model. $\,\,\,\,\,\,\,\,\,\,\rho\,\,\,\,\,\,\,\,\,\,$ | $\,\,\,\,\,\,\,\,\,\,\bar{t}_{c}(\infty)$ |  $1/\nu_{1}$ |  $1/\nu_{2}$ | $\beta_{1}$ | $\beta_{2}$ | $d_{f_{1}}^{-}$ | $d_{f_{1}}^{+}$ | $d_{f_{2}}^{-}$ | $d_{f_{2}}^{+}$ | $\tau$ |
---|---|---|---|---|---|---|---|---|---|---|---
0.001 | 1.50(5) | 0.42(3) | 0.16(3) | 0.07(2) | 0.09(2) | 1.83(1) | 1.86(1) | 1.83(1) | 1.86(1) | 1.03(5) |
0.01 | 1.35(5) | 0.58(3) | 0.43(3) | 0.07(2) | 0.09(2) | 1.83(1) | 1.86(1) | 1.83(1) | 1.86(1) | 1.83(5) |
0.05 | 1.230(5) | 0.66(3) | 0.61(3) | 0.11(2) | 0.12(2) | 1.84(1) | 1.86(1) | 1.84(1) | 1.86(1) | 1.95(3) |
0.1 | 1.170(5) | 0.66(3) | 0.65(3) | 0.12(2) | 0.13(2) | 1.85(1) | 1.86(1) | 1.85(1) | 1.86(1) | 1.99(3) |
0.8 | 1.025(5) | 0.74(2) | 0.74(2) | 0.13(2) | 0.13(2) | 1.86(1) | 1.86(1) | 1.86(1) | 1.87(1) | 2.03(2) |
2 | 0.977(5) | 0.75(1) | 0.75(1) | 0.13(2) | 0.13(2) | 1.87(1) | 1.87(1) | 1.87(1) | 1.87(1) | 2.04(2) |
10 | 0.928(5) | 0.75(1) | 0.75(1) | 0.12(2) | 0.12(2) | 1.88(1) | 1.88(1) | 1.88(1) | 1.88(1) | 2.052(5) |
100 | 0.901(5) | 0.75(1) | 0.75(1) | 0.12(2) | 0.12(2) | 1.89(1) | 1.89(1) | 1.89(1) | 1.89(1) | 2.054(5) |
$2$D site-percolation | 0.592746 | 3/4 | 3/4 | 5/48 | 5/48 | 91/48 | 91/48 | 91/48 | 91/48 | 187/91 |
Figure 3: The average critical threshold $\bar{t}_{c}$ as a function of
inverse size $1/L$ for $\rho=10^{-3},10^{-2},5\times 10^{-2},10^{-1},8\times
10^{-1},2,10$, and $10^{2}$. The dotted-dashed lines show our extrapolations
to the infinite system size at $1/L\rightarrow 0$ to estimate the asymptotic
$\bar{t}_{c}(\infty)$, listed in the first column of Table 1.

###  III.2 Universal scaling functions

One may hypothesize that the distributions $P_{Y}(Y,L)$ of the introduced
random observables $Y=\Delta,S_{c}^{\pm}$ and $t_{c}$, for a system of size
$L$, have to take the following finite-size scaling forms r26-2

$P_{\Delta}(\Delta,L)=L^{\beta_{2}}\,f_{\Delta}(\Delta\,L^{\beta_{2}}),$ (9)

$P_{S^{\pm}_{c}}(S_{c}^{\pm},L)=L^{-d_{f^{\pm}_{2}}}\,f_{S_{c}^{\pm}}(S_{c}^{\pm}\,L^{-d^{\pm}_{f_{2}}}),$
(10)

$P_{t_{c}}(t_{c},L)=L^{1/\nu_{2}}\,f_{t_{c}}(\delta t_{c}\,L^{1/\nu_{2}}),$
(11)

where $f_{Y}(\cdot)$ are the universal scaling functions, and $\delta
t_{c}=t_{c}-\bar{t}_{c}(\infty)$. The size of the largest gap $\Delta$ in the
order parameter, as the maximum of random variables drawn from independent
realizations of the entire percolation process, and the corresponding size of
the largest cluster $S^{\pm}_{c}$ which dominates the system at $t_{c}$, can
be viewed as the extreme events whose distributions are predicted by the
extreme-value theory r26-3 . The critical threshold $t_{c}$, on the other
hand, is a non-extremal variable drawn and averaged from independent
realizations. Therefore, according to the central limit theorem, we expect a
Gaussian distribution for $t_{c}$.

The scaling forms (9)-(11) also suggest an alternative approach to estimate
the involved critical exponents by data collapse onto the universal scaling
curves. In particular, in order to estimate the exponent $\nu_{2}$, we will
use the following finite-size scaling law

$P_{t_{c}}(t_{c}^{*},L)\sim L^{1/\nu_{2}},$ (12)

where $t_{c}^{*}$ is the mode of the distribution $P_{t_{c}}(t_{c},L)$ which
is also very close to $\bar{t}_{c}(\infty)$ due to the Gaussian form of the
distribution.

Figure 4: Estimation of the Fisher exponent $\tau$ at the asymptotic critical
threshold for $\rho=10^{-3},10^{-2},5\times 10^{-2},10^{-1},8\times
10^{-1},2,10$, and $10^{2}$, extracted from the scaling region of the average
number $\bar{n}(s,t_{c}(\infty))$ of clusters of size $s$. The averages are
taken over $10^{6}$ independent realizations for a system of linear size
$L=200$. The lower solid line shows the best fit to our data for
$\rho=10^{-3}$, and the upper solid line compares our data for the largest
$\rho$, with the predicted behavior for the standard 2D percolation model.
Figure 5: The difference time threshold
$\delta\bar{t}_{c}(\infty)=\bar{t}_{c}(\infty)-\bar{t}^{\text{site}}_{c}(\infty)$
as a function of $\rho$. The solid line shows the best power law fit
$\propto\rho^{-\zeta}$ to our data with the exponent $\zeta=0.10(2)$ estimated
from the log-log plot of the data shown in the Inset.  Figure 6: Universal
scaling functions for $\rho=2$. (a–d) Distributions of the observables
$\Delta$, $t_{c}$, $S_{c}^{-}$ and $S_{c}^{+}$ for different system sizes
$L=40,80,100,120,150$ and $200$. (e–h) Corresponding universal scaling
functions of rescaled variables $\Delta L^{\beta_{2}}$, $\delta
t_{c}L^{1/\nu_{2}}$, $S_{c}^{-}L^{-d^{-}_{f_{2}}}$ and
$S_{c}^{+}L^{-d^{+}_{f_{2}}}$. The green solid lines are the best fit to our
data compatible with Frechet extreme value in Eq. 15 with $\xi=+0.064(5)$,
$\mu=0.181(5)$ and $\sigma=0.069(5)$ (e), Gaussian distribution (f), and
Weibull extreme value with $\xi=-166(5)$, $\mu=0.0.390(5)$ and
$\sigma=0.163(5)$ in (g), and $\xi=-0.360(5)$, $\mu=0.616(5)$ and
$\sigma=0.212(5)$ in (h).

##  IV Numerical results

In this Section we are going to present the results of our numerical
simulations of our model previously described in Section II. We have carried
out extensive simulations for different linear sizes $L$ and various choices
of $\rho$ over a range of $5$ orders of magnitude $\rho\in[10^{-3},10^{2}]$ to
best estimate the critical exponents and threshold for every $\rho$ in the
asymptotic limit $L\rightarrow\infty$. Figures 2a,b demonstrate the plots of
$\bar{\Delta}(L)$ and $\chi_{\Delta}$ as a function of $L$ for various choices
of $\rho$. We use logarithmic scales on both the horizontal and vertical axes
to extract the critical exponents in all graphs. The exponents $\beta_{1}$ and
$\beta_{2}$ are measured by examining the scaling relations defined in Eq.
(5). As shown in the figures (and summarized in Table 1) we find a spectrum of
exponents which vary with $\rho$. The lower solid lines show the best fit to
our data for $\rho=10^{-3}$ giving $\beta_{1}\sim 0.07(1)$ and $\beta_{2}\sim
0.09(1)$. The upper solid lines in the graphs show the corresponding exact
exponents predicted for 2D standard percolation model.
In the limit $\rho\rightarrow\infty$, there are infinitely many nonlocal links
which makes the RW’s movement totally nonlocal and the model effectively
becomes equivalent to a 2D site-percolation model in which every site is
randomly chosen and gets infected. Therefore, we expect that in the large
$\rho$ limit, our model falls into the 2D percolation universality class with
exactly known critical exponents. As can be seen in the figures (and Table 1),
the agreement with our expectation is evident for $\rho\geq 1$.
In order to measure the exponents $\nu_{1}$ and $\nu_{2}$ by using Eqs. (8)
and (12), we first estimate the asymptotic average critical threshold
$\bar{t}_{c}(\infty)$ for every $\rho$ by plotting our data for the average
value of the threshold $\bar{t}_{c}(L)$ obtained for a given system size $L$
as a function of $1/L$ and extrapolating to the infinite system size limit
$1/L\rightarrow 0$ (see Fig. 3 and second column in Table 1). Figures 2c,d
show the plots $\bar{t}_{c}(L)-\bar{t}_{c}(\infty)$ and $P_{t}(t^{*}_{c},L)$
as a function of $L$, yielding the exponents $\nu_{1}$ and $\nu_{2}$,
respectively, also summarized in Table 1. We find $1/\nu_{1}=1/\nu_{2}\sim
0.75(1)$ for $\rho>1$, which agrees very well with the exact exponent
$1/\nu=3/4$ for 2D percolation. However, for smaller values of $\rho$ we find
that $\nu_{1}\neq\nu_{2}$ whose values increase as $\rho$ decreases (see the
third and forth columns in Table 1). This finding is intuitively
understandable, since it is consistent with the known observation that the
correlation exponent increases with increase in the underlying dimensionality,
where, in our case, the decrease in $\rho$ decreases the effective
dimensionality from infinity (for $\rho\rightarrow\infty$) to a pure 2D (for
$\rho=0$).
We have also measured the fractal dimensions defined in Eqs. (6) and (7) for
the giant cluster of infected individuals at the critical threshold which are
shown in Figs. 2e-h. For higher values of $\rho$, the sites all over the
lattice become more accessible for the random walker to be visited and, thus,
leads to formation of more compact clusters with higher fractal dimension
(Note that despite the local dependencies introduced by consecutive steps of
an ordinary RW, the overall distribution of the sites visited by the RW in the
presence of the nonlocal extra links is well spread across the system.)

The other geometric exponent that characterizes the universality class of
critical systems is provided by the distribution of cluster sizes which
appears to follow a power law

$\bar{n}(s,\bar{t}_{c}(\infty))\sim s^{-\tau},$ (13)

for large $s$, where $\bar{n}$ gives the average number of clusters of $s$
connected sites, per lattice site, measured at the asymptotic critical point
$\bar{t}_{c}(\infty)$ for each $\rho$. The Fisher exponent $\tau$ is a
universal quantity whose value is the same for all systems of a given class.
Figure 4 shows our computations for $\bar{n}(s,\bar{t}_{c}(\infty))$ for
various values of $\rho$ on a system of linear size $L=200$. The estimated
values of Fisher exponents in the scaling region are presented in the Inset of
Fig. 4 as a function of $\rho$. The distribution is proportional to the
inverse cluster size $\bar{n}(s,\bar{t}_{c}(\infty))\propto 1/s$ with the
exponent $\tau\sim 1$ for the smallest $\rho(=10^{-3})$, while the exponent
converges to the exact value $\tau=187/91$, known for the standard 2D
percolation, in the limit $\rho\rightarrow\infty$. Our estimated Fisher
exponents are summarized in the last column of Table 1.

A closer look at Figure 4 shows that for each intermediate $\rho$, the
corresponding graph shows a crossover behavior between two regimes with small
($s<s_{c}$) and large ($s>s_{c}$) cluster sizes, where the crossover point
$s_{c}$ vanishes as $\rho$ increases. The behavior for $s<s_{c}$ is governed
by the RWs on a pure square lattice with $\tau\sim 1$ and, for $s>s_{c}$, it
is given by $\tau=187/91$ for the ordinary percolation in two dimensions. This
observation strongly suggests that the dependence on $\rho$ in our estimated
exponents listed in Table 1, is just a reminiscence of a crossover behavior
that system undergoes between the 2D percolation and 2D RWs universality
classes. This is however, due to the finite-size of the system, which due to
its high computational cost, we are not able to perform simulations for much
larger systems. This is why we call our estimated exponents ”effective
exponents”.

A remarkable aspect of our model is that it allows us to define and predict
the emergence of a pandemic in a more precise form as the state of an
epidemics where the percolation threshold $\bar{t}_{c}(\infty)$ is reached.
This gives the time at which an epidemic outbreak becomes a pandemic whose
prediction is of great importance in epidemiology. In our model, the temporal
and spatial progress of a global epidemic is totally governed by the
percentage of available nonlocal links in the sample space controlled by the
parameter $\rho$. In the limit $\rho\rightarrow\infty$, the threshold is given
by the critical point of a 2D site-percolation at
$\bar{t}^{\text{site}}_{c}(\infty)\sim 0.5927$. We find that the difference
time threshold
$\delta\bar{t}_{c}(\infty)=\bar{t}_{c}(\infty)-\bar{t}^{\text{site}}_{c}(\infty)$
exhibits a scaling relation with $\rho$ as follows

$\delta\bar{t}_{c}(\infty)\sim\rho^{-\zeta},$ (14)

with the exponent $\zeta\sim 0.10(2)$ (see Fig. 5). This relationship suggests
that the onset of a global pandemic can be algebraically delayed by reducing
the percentage of nonlocal links in our setting.

Figure 7: Shape parameters $\xi$ defined in the Generalized Extreme Value
(GEV) distribution (Eq. 15), as functions of $\rho$ obtained from the best
fits to our data for the universal scaling functions of the observables
$\Delta$ (blue squares), $S_{c}^{-}$ (red circles) and $S_{c}^{+}$ (green
triangles). $f_{\Delta}$ is compatible with Frechet distribution for all
$\rho$, $f_{S_{c}^{-}}$ with Weibull distribution for all $\rho$ expect for
$\rho=10^{-3}$ which is more close to the Gumbel distribution with $\xi=0$.
$f_{S_{c}^{+}}$ is in agreement with the Weibull distribution for all $\rho$
with negative shape parameters $\xi<0$.

In the rest of the present Section, we are going to establish a connection
between the universal scaling functions $f_{Y}(\cdot)$ mentioned in Eqs.
(9)-(11) with the extreme value theory. This connection was first introduced
in r26  for a wide class of percolation problems in different dimensions where
$f_{Y}(\cdot)$ for $Y=\Delta-\bar{\Delta}(L)$ and
$S_{c}^{-}-\bar{S}_{c}^{-}(L)$ were shown to be consistent with the Gumbel
distribution r29 . Here by $Y$ we do not mean the fluctuation of the
corresponding random variables about their average, but we simply mean the
variables themselves i.e., $Y$ denotes for either $\Delta$, $S_{c}^{-}$ or
$S^{+}_{c}$ directly. This introduces computational simplicity, in particular
for problems in which the average of the variables is not in a priori known.

In order to obtain the universal scaling functions $f_{Y}(\cdot)$ defined in
Eqs. (9)-(11), we first compute the distributions of corresponding variables
i.e., $Y=\Delta,S_{c}^{\pm}$ and $t_{c}$, for various system sizes. Figures
6a-d. show an example of the distributions for $\rho=2$ and system sizes
$L=40,80,100,120,150$ and $200$. Strong size-dependence is evident in these
figures. As shown in subsequent figures 6e-h, when the horizontal and vertical
axes in Figs. 6a-d are suitably rescaled, all data collapse onto a single
universal curve. We find excellent data collapse over the entire range of
scaling variables. The best fits to our data are shown with the solid (green)
lines in Figs. 6e-h. We find that the universal function for $Y=t_{c}$, is in
perfect agreement with a Gaussian function with R-square
$R^{2}_{f_{t_{c}}}=0.9977(10)$ (see Fig. 6f). For other variables
$Y=\Delta,S_{c}^{\pm}$, we have examined the following Generalized Extreme
Value (GEV) distribution which unifies the types I, II, and III extreme value
distributions into a single family, by allowing a continuous range of possible
shapes

$\displaystyle f(z)=\frac{1}{\sigma}\,\,t(z)^{\xi+1}\,\,e^{-t(z)}$ (15)

with

$\displaystyle
t(z)=\left\\{\begin{array}[]{rcl}\Big{(}1+\xi\,\,(\frac{z-\mu}{\sigma})\Big{)}^{-1/\xi}&\text{if}&\xi\neq
0,\\\ \exp\big{(}-\frac{(z-\mu)}{\sigma}\big{)}&\text{if}&\xi=0.\\\
\end{array}\right.$ (18)

It is parameterized with location $\mu$, scale parameter $\sigma$, and a shape
parameter $\xi$. When $\xi<0$, the GEV is equivalent to the Weibull (type III)
extreme value. When $\xi>0$, the GEV is equivalent to the Frechet (type II),
and in the limit $\xi\rightarrow 0$, it becomes the Gumbel (type I)
distribution. The best fits in MATLAB to the universal scaling functions
$f_{\Delta}(\Delta\,L^{\beta_{2}})$ (Fig. 6e),
$f_{S^{-}_{c}}(S_{c}^{-}\,L^{-d_{f2}^{-}})$ (Fig. 6g) and
$f_{S^{+}_{c}}(S_{c}^{+}\,L^{-d_{f2}^{+}})$ (Fig. 6h), are attained with the
shape parameters $\xi=+0.064(5)$ ($R^{2}_{f_{\Delta}}=0.9958(10)$) ,
$-0.166(5)$ ($R^{2}_{f_{S_{c}^{-}}}=0.9927(10)$) and $-0.360(5)$
($R^{2}_{f_{S^{+}}}=0.9827(10)$), respectively. Our results indicate that the
universal function for the maximum jumps $\Delta$ is compatible with the
Frechet distribution and the ones for $S_{c}^{\pm}$ can be well described by
Weibull distributions. The shape parameters as a function of $\rho$ for
$f_{\Delta}$, $f_{S_{c}^{-}}$ and $f_{S_{c}^{+}}$ are presented in figure 7.
Despite a slight change, almost the same behavior can be observed for the
whole range of studied $\rho$.

##  V Concluding remarks

We have introduced a variant of percolation model to study the effect of the
underlying network topology, controlled by the density $\rho$ of nonlocal
links, on the spread of an infection/information when an infected/informed
person performs independent random walks on a graph embedded in two
dimensions. Our model reduces to the ordinary nonequilibrium diffusion problem
in 2D in the limit $\rho\rightarrow 0$ and, to the standard site-percolation
model in pure 2D in the limit $\rho\rightarrow\infty$ with a known equilibrium
geometric phase transition. Due to finite-size effects, we find a spectrum of
effective universality classes in terms of the single control parameter
$\rho$, crossing over between these two important prototype models in the
statistical physics. We have carried out extensive numerical simulations to
estimate the effective exponents from a finite-size scaling analysis of the
extreme event records in the temporal evolution of the largest cluster in an
ensemble of independent realizations. This approach facilitates the practical
applications of our model into the realistic situations. We find that the time
threshold at which the onset of a pandemic occurs, exhibits a scaling relation
with the control parameter $\rho$ that is a quantity of great interest in the
epidemiology. We have also established a connection between the extreme value
theory and the universal scaling functions describing the probability
distribution of the extreme observables. Our study calls for further future
research to mathematical understanding of our model to realize the nature of
the transitions as well as exact description of the crossover behavior in
terms of the density of the nonlocal links.

## Acknowledgments

A.A.S. acknowledges the supports from the Alexander von Humboldt Foundation
(DE) and the research council of the University of Tehran. We also thank the
High Performance Computing (HPC) center in the University of Cologne, Germany,
where a part of computations have been carried out. M.F. would like to thank
M. Ebrahimi for her helps on programming.

## References

  * (1) D. Stauffer, and A. Aharony, Introduction To Percolation Theory (Taylor and Francis, 1991).
  * (2) B. Bollobas, and O. Riordan, Percolation (Cambridge University Press, 2006).
  * (3) A.A Saberi, Physics Reports 578 1-32 (2015).
  * (4) M. Sahimi, Applications of percolation theory (CRC Press, 2003).
  * (5) B. Vigolo, C. Coulon, M. Maugey, C. Zakri, and P. Poulin, Science. 309, 920 (2005).
  * (6) C. Grimaldi, and I. Balberg, Phys. Rev. Lett. 96, 066602 (2006).
  * (7) I. Derenyi, G. Palla, and T. Vicsek, Phys. Rev. Lett. 94, 160202 (2005).
  * (8) D. S. Callaway, M. E. J. Newman, S. H. Strogatz, and D.J. Watts, Phys. Rev. Lett. 85 5468 (2010).
  * (9) T. Kalisky, and R. Cohen, Phys. Rev. E. 73 035101 (2006).
  * (10) P. W. Kasteleyn, and C. M. Fortuin, J. Phys. Soc. Jpn. 26 11 (1969).
  * (11) C. M. Fortuin, and P. W. Kasteleyn, Physica. 57 536 (1972).
  * (12) V. S. Dotsenko, M. Piccoa, P. Windeya, G. Harris, E. Martinec, and E. Marinari, Nucl. Phys. B. 448 577 (1995).
  * (13) V. S. Dotsenko, P. Windey, G. Harris, E. Marinari, E. Martinec, and Marco Picco, Phys. Rev. Lett. 71 811 (1993).
  * (14) A.A. Saberi, H. Dashti-Naserabadi, EPL (Europhysics letters), 92(6), 67005 (2011).
  * (15) A.A. Saberi, Journal of Statistical Mechanics: Theory and Experiment, 2009(07), P07030 (2009).
  * (16) H. Dashti-Naserabadi, A.A. Saberi, S.H.E. Rahbari, and H. Park, Physical Review E(R) 100, no. 6 060101 (2019).
  * (17) S. Davatolhagh, M. Moshfeghian, A.A. Saberi, Journal of Statistical Mechanics: Theory and Experiment no. 02 P02015 (2012).
  * (18) S. G. Anekal, P. Bahukudumbi, and M. A. Bevan, Phys. Rev. E. 73 020403 (2006).
  * (19) N. Gnan, E. Zaccarelli, and F. Sciortino, Nature Communications. 5 3267 (2014).
  * (20) A.A. Saberi, Appl. Phys. Lett. 97 154102 (2010).
  * (21) J. Fan, J. Meng, and A.A. Saberi, Physical Review E, 99(2), 022304 (2019).
  * (22) A.A. Saberi, Physical review letters, 110(17), 178501 (2013).
  * (23) A.A. Saberi, The Astrophysical Journal Letters (ApJL) vol. 896, No. 2, L25 (2020).
  * (24) J.L. Cardy, and P. Grassberger, Journal of Physics A: Mathematical and General, 18(6), L267 (1985).
  * (25) J.L. Cardy, Scaling and renormalization in statistical physics (Vol. 5. Cambridge university press, 1996.)
  * (26) W. Feller, An Introduction to Probability Theory and its Applications, 3rd ed. (John Wiley & Sons, Hoboken, NJ, 1968), Vol. I.
  * (27) W. Feller, An Introduction to Probability Theory and its Applications, 2nd ed. (John Wiley & Sons, Hoboken, NJ, 1971), Vol. II.
  * (28) B. D. Hughes, Random Walks and Randon Environments, Volume 1: Random Walks (Oxford University Press, Oxford, UK, 1995).
  * (29) G. Lawler, O. Schramm, W. Werner, On the scaling limit of planar self-avoiding walk, in: Fractal Geometry and Applications, a Jubilee of B. Mandelbrot, in: Proc. Symp. Pure Math., vol. 72, (2004).
  * (30) R. F. Voss, J. Phys. A. 17 L373 (1984).
  * (31) J.L. Cardy, ed. Finite-size scaling (Elsevier, 2012).
  * (32) J. Fan, J. Meng, Y. Liu, A. A. Saberi, J. Kurths, and J. Nagler, Nature. Phys. 16 495 (2020).
  * (33) R. M. D’Souza, and J. Nagler, Nat. Phys. 11, 531–538 (2015).
  * (34) V. Privman, and M.E. Fisher, Phys. Rev. B 30, 322–327 (1984).
  * (35) E. J. Gumbel, Statistics of Extremes (Courier Corporation, 2012).
  * (36) J. Nagler, A. Levina, and M. Timme, Nature. Phys. 7 265 (2011).
  * (37) V. Privman, Finite-size Scaling and Numerical Simulation of Statistical Systems, (World Scientific, Singapore,1900).
  * (38) J. Fan, and X. Chen, Europhysics Letters. 107 28005 (2014).
  * (39) E.J. Gumbel, Statistics of Extremes (Courier Corporation, 2012).