# Evaluation of Feynman integrals with arbitrary complex masses via series
expansions

Tommaso Armadillo<EMAIL_ADDRESS>Roberto Bonciani
<EMAIL_ADDRESS>Simone Devoto<EMAIL_ADDRESS>Narayan
Rana<EMAIL_ADDRESS>Alessandro Vicini<EMAIL_ADDRESS>Centre for Cosmology, Particle Physics and Phenomenology (CP3), Université
Catholique de Louvain, Chemin du Cyclotron, B-1348 Louvain la Neuve, Belgium
Dipartimento di Fisica “Aldo Pontremoli”, University of Milano, Via Celoria
16, I-20133 Milano, Italy  Dipartimento di Fisica, Università di Roma “La
Sapienza”, Piazzale Aldo Moro 5, I-00185 Roma, Italy  INFN Sezione di Roma,
Piazzale Aldo Moro 5, I-00185 Roma, Italy  INFN Sezione di Milano, Via Celoria
16, I-20133 Milano, Italy  Department of Physics, Indian Institute of
Technology Kanpur, 208016 Kanpur, India

###### Abstract

We present an algorithm to evaluate multiloop Feynman integrals with an
arbitrary number of internal massive lines, with the masses being in general
complex-valued, and its implementation in the Mathematica package SeaSyde. The
implementation solves by series expansions the system of differential
equations satisfied by the Master Integrals. At variance with respect to other
existing codes, the analytical continuation of the solution is performed in
the complex plane associated to each kinematical invariant. We present the
results of the evaluation of the Master Integrals relevant for the NNLO QCD-EW
corrections to the neutral-current Drell-Yan processes.

††journal: CPC

Program Title: SeaSyde
CPC Library link to program files: (to be added by Technical Editor)
Developer’s repository link: https://github.com/TommasoArmadillo/SeaSyde
Code Ocean capsule: (to be added by Technical Editor)
Licensing provisions(please choose one): GPLv3
Programming language: Mathematica
Supplementary material:
Nature of problem: Solution of multi-loop Feynman integrals with an arbitrary
number of internal complex-valued massive lines.
Solution method: Solution via series expansions of the systems of differential
equations satisfied by the loop integrals. Analytic continuation of the
solution to any point of the kinematic variables. The solution is transported
from the initial to the final point, by moving in the complex plane associated
to each kinematic variable. In multi-variable problems, the transport is
computed one variable at the time.

##  1 Introduction

The high quality and precision of the large amount of data collected by the
experiments at the CERN LHC calls for a corresponding theoretical effort, to
interpret the results in a meaningful way. The evaluation of the radiative
corrections at second and higher orders in the perturbative coupling constant
expansion has become a need. These challenging calculations require the
solution of multiloop Feynman integrals, which are in several cases still
unknown. The complexity of these problems grows with the number of energy
scales: independent kinematical invariants and masses. The latter might be
complex valued, as it is the case for unstable particles like the gauge and
the Higgs bosons or the top quark. The gauge symmetry of the scattering
amplitude leads to non-trivial cancellations among the different terms
associated to the Feynman diagrams and possible problems of precision loss
must be kept under control. The analytical knowledge of the solution of the
Feynman integrals is thus desirable, to have an arbitrary number of digits
available. When this is not the case, then different approaches might be
useful. Comprehensive reviews on the Feynman integrals, covering the various
strategies can be found in Refs.[1, 2, 3, 4, 5].

It is difficult to design a general algorithm which allows the evaluation of
an arbitrary Feynman integral, because of the presence of singular points and
branch cuts, connected to the physical thresholds and pseudo-thresholds of the
diagram. For this reason a purely numerical approach [6, 7] can encounter
problems whenever an integrable singularity is present. Isolating the singular
points and working out a dedicated analytical rearrangement of the integrand
function may yield the successful evaluation of the integral, decomposed in
the sum of several regular terms. Whenever the identification of the singular
structure or the following rearrangement can not be performed, then numerical
instabilities are expected. The analytical solution is clearly preferable in
all such cases, because the singular behaviour can be exposed and arbitrary
precision can be achieved everywhere else.

We have analytical control over the solution of an integral, if it can be
expanded as a power series at every point of its domain. We dub analytical
solutions those which can be expressed in closed form as a combination of
elementary and special functions, whose power expansion is known. We call
semi-analytical solutions those which can only be represented as power
expansions, without additional functional relations (the latter are typically
available in the closed form case).

Beside the different strategies to solve the problem of integration of the
Feynman integrals, the differential equations approach [8, 9, 10, 11, 12, 13,
14, 15] has become very popular, covering a large number of cases relevant in
multi-loop calculations. The Feynman integrals can be written as a combination
of basic integrals called Master Integrals (MIs) and the latter satisfy, in
fact, systems of first order linear differential equations. These systems can
be solved, in several relevant cases, in closed form: the system is expanded
in the dimensional regularisation parameter $\varepsilon=(4-d)/2$, where $d$
is the dimension of the space-time and, order by order in $\varepsilon$, the
solutions are expressed in terms of a suitable functional basis. A first class
of problems is represented by the systems that admit solutions expressed in
terms of generalised polylogarithmic functions [16, 17, 18, 19]. The algebraic
properties of these functions allow to simplify the final expression of the
scattering amplitude. Moreover, their numerical evaluation is under control
and provided by general public routines [20, 21, 22, 23, 24, 25]. Recently,
another class of problems became relevant for phenomenological applications,
namely cases in which the systems admit solutions expressed in terms of more
general special functions, e.g. of elliptic kind [26, 27, 28, 29, 30, 31, 3].
Also in this case, the properties of the functions under consideration allow
for a simplification in the final expressions and a control on the complexity
of the analytic formulae. The numerical evaluation of such functions is less
general, but nevertheless possible using power expansion representations [32].
Cases in which a closed form instead is not known or not suitable for power
expansions need to be treated using a different strategy.

We can thus summarise that the study of a multi-loop Feynman integral leads to
two distinct cases: $i)$ the solution is achievable in terms of special
functions; $ii)$ the closed form solution is not available. We focus in this
paper on the second case and consider the recent developments of semi-
analytical algorithms to solve via series expansion the systems of
differential equations satisfied by a Feynman integral 111 Different numerical
approaches for the solution of differential equations have been proposed in
Refs. [33, 34]. . The solution is first written as a power series with unknown
coefficients and is inserted in the differential equations; an infinite number
of algebraic equations for the unknown coefficients can be written and
recursively solved. The solution is thus known inside the convergence radius
of the series expansion. The extension to different regions can be achieved
via analytic continuation. This method was applied to the evaluation of one-
dimensional problems (sunrise and vertex corrections, in which the MIs depend
upon a single dimensionless variable) [35, 36, 37, 38, 39, 40, 41]. Then it
was generalised to multi-dimensional problems in [42]. Recently a Mathematica
package, called DiffExp, has been presented [43], allowing the solution of an
arbitrary system of differential equations, provided that the boundary
conditions and the necessary prescriptions for the analytic continuation are
known. This package has been successfully applied to several different
problems [44, 45, 46, 47, 48, 49, 50], always considering the case of real-
valued masses for the particles running inside the loops.

The electroweak (EW) precision physics program at the LHC and future colliders
requires the evaluation of EW radiative corrections which involve the exchange
of $W$, $Z$, and the Higgs bosons. The latter are unstable particles and the
gauge invariant definition of their masses can be achieved in the Complex-
Mass-Scheme (CMS)[51, 52], identifying the complex-valued pole of their
propagator. Precise predictions for the production of the top quark require
the inclusion of off-shellness effects including the top decay width. In this
paper we present a package which is able to deal with complex-valued masses in
the internal lines of the Feynman integral. We develop an original independent
algorithm to implement the analytic continuation.

The paper is organised as follows. In Section 2 we provide a pedagogical
introduction to the solution of a system of first order linear differential
equations by series expansions, with a particular focus on how the analytic
continuation of the result can be safely obtained to an arbitrary point in the
complex plane. We outline the implementation in the Mathematica package
SeaSyde (Series Expansion Approach for SYstems of Differen- tial Equations) of
the algorithm that solves a system of differential equations, for arbitrary
complex-valued kinematical variables and we illustrate the different
computational strategies available. In Section 3, as a first practical
application of the algorithm, we discuss the solution of the MIs needed to
evaluate the two-loop QCD-EW virtual corrections to the neutral-current Drell-
Yan processes [48, 53]. The detailed expressions and evaluation properties of
the integrals used in Ref.[53] constitute an original result of this paper.
Finally, in Section 4 we draw our conclusions.

The latest version of the Mathematica package SeaSyde can be downloaded from
https://github.com/TommasoArmadillo/SeaSyde, while its full documentation for
Version 1.0 is provided in A.

##  2 The solution algorithm

The solution of a system of linear differential equations by series expansions
is well known in the mathematical literature. In this Section we present some
basic definitions and a pedagogical introduction to the procedure, and we show
how it can be turned into an algorithm and applied to the specific case of
Feynman integrals.

###  2.1 Solving differential equations by series expansion

The general method to obtain the solution as a series expansion can be
illustrated with a simple example. Let us consider the differential equation

$\left\\{\begin{array}[]{l}f^{\prime}(x)+\frac{1}{x^{2}-4x+5}f(x)=\frac{1}{x+2}\\\
f(0)=1\end{array}\;;\right.$ (1)

with $x$ being a real variable.

In order to solve the homogeneous equation we use the well-known Frobenius
method. To this aim we introduce
$f_{hom}(x)=x^{r}\sum_{k=0}^{\infty}c_{k}x^{k}$ as the expansion of the
homogeneous solution around $x=0$, with $c_{k}$ being arbitrary coefficients
that we need to determine. By replacing $f_{hom}(x)$ with its power series in
the associated homogeneous differential equation and by collecting all the
terms with the same power of $x$, we obtain an infinite set of algebraic
equations in the unknowns $c_{k}$. In our example, Eq. (1) leads to:

$\left\\{\begin{array}[]{l}{r}{c_{0}}=0\\\ \frac{1}{5}c_{0}+c_{1}(r+1)=0\\\
\frac{4}{25}c_{0}+\frac{1}{5}c_{1}+c_{2}(2+r)=0\\\
\frac{11}{125}c_{0}+\frac{4}{25}c_{1}+\frac{1}{5}c_{2}+c_{3}(3+r)=0\\\
\dots\end{array}\right.$ (2)

The equation associated to the lowest power of $x$ is called indicial
equation, and it assigns the value of the exponent $r$. In our case, we get
${r}{c_{0}}=0$, which corresponds to $r=0$. The other equations determine all
the $c_{k}$ but one; the latter can be chosen arbitrarily, and we set e.g.
$c_{0}=5$. We thus obtain the solution of the homogeneous differential
equation associated to the one in Eq. (1), expressed as the following series
expansion:

$f_{hom}(x)=5-x-\frac{3}{10}x^{2}-\frac{11}{150}x^{3}+...$ (3)

A particular solution for the original problem can now be obtained by applying
the variation of the constant method, where the inverse of the homogeneous
solution, multiplied by the inhomogeneous term, is expanded about $x=0$ and
easily integrated:

$f_{part}(x)=f_{hom}(x)\int_{0}^{x}dx^{\prime}\,\frac{1}{(x^{\prime}+2)}\,f^{-1}_{hom}(x^{\prime})=\frac{1}{2}x-\frac{7}{40}x^{2}+\frac{2}{75}x^{3}+...$
(4)

The general solution is finally given by

$f(x)=f_{part}(x)+Cf_{hom}(x)\;.$ (5)

In order to satisfy the boundary condition $f(0)=1$ we set $C=1/5$.

###  2.2 Singularities and branch cuts

The validity of a solution obtained as a series expansion, with the procedure
outlined in Section 2.1, is limited by its convergence radius. We now discuss
how the latter is defined and how the solution can be extended, via analytic
continuation, from its initial region of convergence to an external arbitrary
point.

A solution written in series representation converges in the complex-$z$ plane
inside a disc centered around the expansion point with a convergence radius
limited by the closest singularity. As an example, we consider again the
first-order linear differential equation presented in Eq. (1), now as a
function of a complex variable $z$. By analysing the differential equation, we
can identify three singularities: the factor $(z+2)$ in the denominator of the
inhomogeneous part generates a pole in $z=w_{0}=-2$, while the factor
$(z^{2}-4z+5)$ in the denominator of the homogeneous part implies the presence
of two poles in $z=w_{\pm}=2\pm i$. The solution presented in Eq. (4) thus
converges in the complex plane within a disc $\Gamma_{0}$ centered in
$z_{0}=0$ with radius $2$, as illustrated by the blue circle in Figure 1.

We can analytically continue the solution into a new disc, centered at an
arbitrary point $z_{1}$ internal to $\Gamma_{0}$. Also in this case, the
convergence radius is defined by the closest singular point, as illustrated in
Figure 1 by the orange circle $\Gamma_{1}$. It is possible to demonstrate that
this procedure is unique.

Figure 1:  Example of analytic continuation.

In the case of simple poles in the homogeneous coefficient, we obtain a
logarithmic behaviour of the solution, that requires some additional care in
the analytic continuation. The presence of logarithmic functions makes the
solution multi-valued: it can be made single-valued by adding cuts in the
complex plane, thus specifying the Riemann sheet where the function is
evaluated. This feature is clearly not present in the power series
representation, but must be introduced to allow for a physical interpretation
of the results. To take it into account, we need to define in an unambiguous
way the cut associated to each singularity that shows a logarithmic behaviour:
their presence will then affect the process of analytic expansion.

In our example, we now consider a logarithmic behaviour in the points $w_{0}$,
$w_{\pm}$, and we choose as the branch-cuts the horizontal lines parallel to
the real axis that go from the singular point to $-\infty$. In the following
discussion we will always assume this convention, which represents the
standard choice for logarithmic branch cuts222This conventions is used, for
instance, in Mathematica by the function $\tt{Log[z]}$, for complex values of
the variable $\tt{z}$..

Figure 2:  Example of the effect of branch-cuts on the convergence of the
expanded solution: reduced convergence area (left) and different path for the
analytic continuation (right).

The first consequence of the presence of branch-cuts in the procedure of
analytic continuation is on the area of convergence of the series expansion.
While the proximity with a branch-cut does not modify the convergence radius
of the series itself, it limits the area of the disc where the expansion
converges to the desired value: once the cut is crossed, the series is
evaluated on a Riemann sheet different than the one that we have chosen with
the cuts, and thus the result can not be considered as the solution of the
problem in the assigned domain. This effect is illustrated on the left panel
of Figure 2, that shows how the convergence discs $\Gamma_{0}$, $\Gamma_{1}$
introduced in Figure 1 are modified once the branch cuts are inserted: while
the discs themselves do not shrink, only the reduced highlighted area provides
the correct evaluation of the solution.

The second consequence is on the path that is necessary to follow in order to
extend the solution to an arbitrary point in the complex plane. The path has
to be chosen in such a way that it does not cross any branch cut. This is
illustrated on the right panel of Figure 2, where the dotted path from $z=0$
to $z_{1}$ is now forbidden by the presence of the branch cut corresponding to
the pole $w_{+}$. The cut needs to be avoided with a path that goes around to
the right of the singularity, as shown by the solid line. For the sake of
readability, we have not drawn all the discs along the path, representing the
intermediate steps of the analytic continuation.

###  2.3 Choice of the path

In the Mathematica package SeaSyde, we use the same convention for the branch-
cuts as the one we introduced in the previous section, with branch-cuts
parallel to the real axis and going from the singularities to $-\infty$. As it
has been shown in Section 2.2, with this choice the value of the solution in
each point of the complex plane is unambiguously defined, as long as the path
we use for the analytic continuation does not cross any branch-cut. In the
following, we briefly describe how such path can be defined when connecting
two arbitrary points on the complex plane and the algorithm for its
determination implemented within the package. The analytic continuation
described in the previous Section, moving in the complex plane of each
kinematical invariant, allows to treat exactly in the same way the cases of
real- and complex-valued masses. This represents the main remark of this
paper.

Figure 3:  The two possible approaches in the definition of the path: the one
requiring a logarithmic expansion (left) and the one only relying on the
Taylor expansion (right).

There are two main approaches that can be used while defining the path: we can
avoid the singularities or we can use them as an expansion point. This two
possibilities are illustrated in Figure 3, where we consider once again the
pole structure given by Eq. (1), and we aim to move the boundary conditions
from the point $z_{0}$ to the point $z_{1}$. By choosing the singularity in
$w_{+}$ as one of the possible expansion points, as shown in the left panel of
Figure 3, we obtain a logarithmic expansion around $w_{+}$, while if we avoid
the singularity, as shown in the right panel of Figure 3, we only rely on the
ordinary Taylor expansion. While leading exactly to the same result, the two
approaches have different advantages and disadvantages. By using the
logarithmic expansion, we can expect a larger radius of convergence of the
series, since the latter is not affected anymore by the presence of the pole
upon which we are expanding: this can effectively reduce the number of steps
required to reach the point of interest. Furthermore, the presence of the
branch-cut, starting from the singularity we are expanding upon, is
automatically taken into account by the explicit logarithms appearing in the
series. On the other hand, in the current implementation of SeaSyde, the
expansion on a singular point requires a longer evaluation time with respect
to the ordinary Taylor expansion on a regular point of the complex plane. An
optimised choice between the two different approaches thus requires to
carefully balance these two effects. In the physical applications we have
dealt with while testing the code, we have observed how the presence of
multiple singularities close to each other does not allow for a significant
improvement of the convergence radius while using the logarithmic expansion.
This fact has led us to prefer, for our final implementation, the more regular
behaviour of the ordinary Taylor expansion.

With this choice, the algorithm to determine the path is straightforwardly
implemented. Since all the branch-cuts are parallel to each other and in the
same direction, it is always possible to connect two arbitrary points on the
complex plane with a path that goes around to the right of the rightmost
singularity that lays between them, as it is shown in the right panel of
Figure 3.

There is only one edge case that we need to discuss in more detail, that is
the one in which the starting and ending point lay down on the same branch-
cut, which is usually the real axis. In this case, indeed, there is an
ambiguity, since it is not clear where the points are located with respect to
the cut. This ambiguity can be discussed and solved in the light of the
Feynman prescription for the particle propagators, which enforces the
causality requirement of the theory. If we consider, for the sake of
simplicity, the propagator $\Pi(s)=i/(s-m^{2}+i\delta)$ of a massive scalar
field, despite the kinematic invariant $s$ and the mass $m$ being real, the
Feynman prescription shifts the first one to $s+i\delta$, i.e. in the upper
half of the complex-$s$ plane. Having a prescription which pushes the solution
of the differential equation in one specific half of the complex-$z$ plane,
removes any ambiguity associated with the evaluation of the function on the
branching cut and allows us to uniquely determine the solution for any given
Boundary Condition (BC). From a practical point of view, we might encounter
two different scenarios: without or with poles between the start and the end
points. In the first one the choice of the path is trivial and we can move
directly from one point to the other, because, by construction, we are never
crossing the cut. We can see this in the left panel of Figure 4.

Figure 4:  Examples of two possible paths for linking two points laying on the
same branch-cut: no singularities between them are present (left) or at least
one is (right).

In the second one, instead, we could again either expand on the singularity or
avoid it, as previously discussed. The choice we have made to avoid the
singularities, however, forces us to adopt an horseshoe path, as depicted in
the right panel of Figure 4. A prescription of $-i\delta$, on the contrary,
would have implied an horseshoe in the lower half of the complex plane. The
prescription is not hard-coded within SeaSyde and must be passed as an input
parameter. As a consequence, the package can automatically perform the
analytic continuation of the solution even if the user adopts different
conventions.

###  2.4 Solving a system of differential equations

The evaluation of a scattering amplitude requires the numerical calculation of
a large set of Feynman integrals. Their number can range from a few, for one-
loop corrections, to hundreds or thousands, for multi-loop ones. Fortunately
not all of them are independent from each other. Integration-by-Parts
Identities [54, 55, 56] allow us to express these integrals as a linear
combination of MIs333Several computer implementations of the procedure are
available [57, 58, 59, 60, 61, 62, 63, 64, 65]..

The computation of the MIs can be faced in different ways. One of these
possibilities is by writing the system of first order linear differential
equations with respect to the kinematical invariants $\mathbf{s}$, satisfied
by the MIs. If $\vec{I}(\varepsilon,\mathbf{s})$ is a vector with the MIs of
the problem under consideration as components, the systems can be expressed as
follows:

$\frac{\partial}{\partial
s_{\alpha}}\vec{I}(\varepsilon,\mathbf{s})=\mathbf{A}_{\alpha}(\varepsilon,\mathbf{s})\,\vec{I}(\varepsilon,\mathbf{s})\,,$
(6)

where $s_{\alpha}$ is one of the invariants, and
$\mathbf{A}_{\alpha}(\varepsilon,\mathbf{s})$ is the matrix of the
coefficients of the system, containing rational functions of all the
invariants and depending upon the parameter $\varepsilon$. We look for a
solution in form of a Laurent series in $\varepsilon$,
$\vec{I}=\sum_{k=-N}^{\infty}\varepsilon^{k}\,\vec{I}^{(k)}$, where $N$ is the
order of the maximal pole of the MIs. We substitute this form of the solution
in Eq. (6) and, order-by-order in $\varepsilon$, we study the resulting system
of equations, starting from the one associated to the lowest power of the
expansion in $\varepsilon$.

A particularly suitable situation is the one in which the choice of a basis of
MIs yields systems of equations that, order-by-order in $\varepsilon$, are in
triangular form. This allows for a straightforward solution of the problem.
The solution of the system can proceed in cascade: we solve a first-order
linear differential equation at a time, starting from the one that involves a
single MI. The solution of this equation is plugged in the following one, that
becomes a first-order non-homogeneous differential equation. Also the solution
of the second equation is plugged, together with the solution of the first
one, into the third equation, that in turn becomes a first-order non-
homogeneous differential equation as well. We proceed in this way up to the
complete solution of the system444This is the class of problems that can admit
solutions in closed form, in terms of polylogarithmic functions. The system
simplifies further, if it is possible, by an appropriate change of basis, to
cast it in “canonical form” [15].. Each first-order differential equation can
be solved using the procedure outlined in section 2.1, i.e. in a semi-analytic
way.

However, the general situation is the one in which we face problems that are
described by systems that include subsets of equations that are coupled555This
means that we cannot solve just first-order linear differential equations, but
we must solve, in principle, also higher-order differential equations,
equivalent to the coupled subset under consideration. at each order in
$\varepsilon$. In these cases, the strategy that SeaSyde uses to solve the
systems is the following. The equations that are in triangular form are solved
in cascade (as described above), while the subsets of coupled equations are
solved expanding the corresponding masters in power series, with unknown
coefficients. These coefficients are then fixed by the solution of linear
systems of equations, constrained by the solution of the system of
differential equations.

The Mathematica package SeaSyde, in its current implementation, can solve
systems both in triangular and in coupled forms, covering the entire (current)
set of possibilities encountered so far in the evaluation of Feynman diagrams
for collider processes.

The systems of differential equations for a fixed order in $\varepsilon$
depend on multiple kinematic variables. We solve the problem by analysing only
one of them at a time. In doing so, the dependence on all the other variables
is parametric, and we reduce the problem to a system of differential equations
in one variable. The strategy adopted in the package SeaSyde is the following.
Suppose that we are dealing with a problem depending upon two dimensionless
variables, $x$ and $y$ . We ask SeaSyde to evolve the solution from the point
in which we fixed the boundary conditions, say $(x_{0},y_{0})$, to the point
$(x_{1},y_{1})$. The solution of the problem will proceed solving the system
in $x$, with a fixed value of $y=y_{0}$, from $x_{0}$ to $x_{1}$ and then
solving the system in $y$ with $x=x_{1}$ from $y=y_{0}$ to $y=y_{1}$.
Integrating in $x$, at fixed $y=y_{0}$, the eventual cuts will depend on
$y_{0}$. However, this value is fixed and then the cuts in the complex $x$
plane are given. To overcome the cuts, SeaSyde uses the procedure outlined in
section 2.1. The same situation holds when we integrate in $y$ at fixed
$x=x_{1}$. Now the eventual cuts will be in the complex $y$ plane and they
will depend upon $x_{1}$, which is a fixed value. This means that also the
cuts in $y$ are now given and SeaSyde is able to overcome them with a path
that moves at the right of the right-most branching point in $y$ (section
2.1).

We can start to solve the system from the lowest order in the $\varepsilon$
expansion and work our way up to the desired order.

##  3 Two-loop MIs for the neutral-current Drell-Yan

As a first practical application of the algorithm discussed in section 2, we
provide the results for the evaluation of the MIs relevant for the mixed QCD-
EW corrections to the neutral-current Drell-Yan process, recently presented in
Ref [53]. We focus on the integrals numbered 32-36 in Ref. [66], because their
complexity offers an interesting case to study the performances of our new
package. The solution proposed in Ref. [66] in terms of Chen-Goncharov
repeated integrals has severe limitations in the numerical evaluation in the
physical region, and a semi-analytical approach represents an interesting
alternative.

We first introduce the problem, then we show a comparison between the solution
obtained with real and complex internal masses, and, finally, we discuss the
numerical and time performances of our code.

###  3.1 Master integrals for the mixed EW-QCD corrections to the neutral-
current Drell-Yan

The NNLO mixed QCD-EW corrections to the neutral current Drell-Yan have been
calculated recently in Ref. [53]. In the complete amplitude 204 MIs appear and
they can be evaluated using the method of differential equations.

These integrals can be classified into three different groups [66], according
to the number of their internal massive lines: zero, one or two. The four-
point functions corresponding to the first two sets are known in closed
analytic form in terms of generalized polylogarithms (GPLs). The last group is
the most challenging due to the presence of different kinematic scales. In
particular 36 MIs belong to this set. For 31 of them, a solution in terms of
GPLs is known. The remaining 5 have been formally solved in [66] using a
formal expression in terms of Chen-Goncharov integrals [67]. These MIs belong
to the following integral family (see the classification in [53]):

$\\{{\cal D}_{1},{\cal D}_{2}-m_{V}^{2},{\cal D}_{12},{\cal D}_{1;1},{\cal
D}_{2;1},{\cal D}_{1;12},{\cal D}_{2;12}-m_{V}^{2},{\cal D}_{1;3},{\cal
D}_{2;3}\\}$ (7)

where $V$ can be either $W$ or $Z$ and $m^{2}_{V}$ is the mass-squared of the
$V$ boson. We defined ${\cal D}$s to be:

${\cal D}_{i}=k_{i}^{2},\;{\cal D}_{ij}=(k_{i}-k_{j})^{2},\;{\cal
D}_{i;j}=(k_{i}-p_{j})^{2},\;{\cal D}_{i;jl}=(k_{i}-p_{j}-p_{l})^{2},$ (8)

where $p_{i}$ are the external momenta and $k_{j}$ the loop ones. In
particular the 5 MIs are:

$\displaystyle\text{Master 32}\;:\>\\{0,1,1,1,0,1,1,0,1\\}\quad$
$\displaystyle\text{Master 33}\;:\>\\{1,1,1,1,0,1,1,0,1\\}$
$\displaystyle\text{Master 34}\;:\>\\{1,1,1,1,-1,1,1,0,1\\}\quad$
$\displaystyle\text{Master 35}\;:\>\\{1,1,1,1,0,1,1,-1,1\\}$
$\displaystyle\text{Master 36}\;:\>\\{1,1,1,1,-1,1,1,-1,1\\}\;,$ (9)

where we keep the numbering from Ref. [53]. It is possible to numerically
evaluate the Chen iterated integrals in the Euclidean region, where the
boundary conditions are imposed. Their analytic continuation in the physical
region, however, is non trivial and no standard technique to perform it is
available. For this reason their computation in Ref. [53] has been performed
by using a semi-analytical approach, that is by solving the relevant system of
differential equations using a series expansion method as described in Section
2. By approaching the numerical evaluation of these MIs in such a semi-
analytical way, we have to solve the complete system (36x36) of differential
equations which includes the integrals from Eq. (9).

These MIs depend on two dimensionless kinematic invariants $x$ and $y$:

$s=(p_{1}+p_{2})^{2}\;,\qquad t=(p_{1}-p_{3})^{2}\;,\qquad
x=-\frac{s}{M^{2}}\;,\qquad y=-\frac{t}{M^{2}}\;,$ (10)

where $p_{1,2}$ are incoming momenta, $p_{3}$ is an outgoing momentum, $s,t$
are the Mandlestam invariant, and $M$ is the mass of the vector boson. Let us
consider first the case of real-valued masses $M=m_{V}$, with $V=W,Z$.

Once we have chosen the point $(x,y)$ in which we would like to evaluate the
solution, the latter can be obtained in a straightforward way using SeaSyde:

ConfigurationNCDY={EpsilonOrder->4,ExpansionOrder->50};UpdateConfiguration[ConfigurationNCDY];SetSystemOfDifferentialEquation[SystemOfEquations,
BoundaryConditions, MIs, {x-I,y+I}, pointBC];TransportBoundaryConditions[ {
tValue , sValue } ]SolutionValue[]

We use in this case a form of the MIs in which they are re-scaled by a power
of $\varepsilon$, so that they do not contain explicit divergences in
$\varepsilon$. Due to this re-scaling, in order to obtain the finite part for
MI32-36 listed in Eq. (9) we need to keep 5 terms in the $\varepsilon$
expansion. We have chosen to keep 50 terms in the expansion with respect to
the kinematical variable, because it is a good compromise between execution
time and precision; we will comment later how the number of terms affects the
numerical precision. In SetSystemOfDifferentialEquation we are preparing the
package for solving the system and, finally, using TransportBoundaryConditions
we are both solving the system and performing the analytic continuation, which
let us extend the solution from pointBC, where the boundary conditions are
imposed, to x=sValue and y=tValue. At this point, by calling again the routine
TransportBoundaryConditions we could transport the boundary conditions from
{x=sValue, y=tValue} to {x=sValueNEW, y=tValueNEW} and so on. By repeating
this procedure multiple times we can easily obtain a numerical grid for all
the MIs. In A we provide the full documentation of the package.

###  3.2 Feynman integrals with complex-valued masses

When dealing with intermediate unstable particles, such as $W$s and $Z$s, it
is useful to perform the calculations in the CMS, in order to regularise the
behaviour at the resonance while preserving the gauge invariance of the
scattering amplitude. To this end, we introduce the complex mass, defined as:

$\mu^{2}_{V}=m_{V}^{2}-i\Gamma_{V}m_{V},$ (11)

where $\Gamma_{V}$, a real parameter, labels the decay rate of the boson. The
complex mass $\mu_{V}$ then replaces the real mass $m_{V}$ in all the steps of
the computation. From Eq. (10), in particular, we can observe that the
dimensionless kinematic variables $x$ and $y$ become complex-valued, once we
set $M=\mu_{V}$.

The complex mass regularises the integrand function in the threshold region:
indeed, all the propagators

$\frac{1}{s-\mu_{V}^{2}+i\delta}$ (12)

do not diverge for any real value of $s$, and, hence, in the kinematic region
near the resonances, the MIs are smooth functions.

This fact let us really appreciate the discussion about the analytic
continuation in the complex plane presented in Section 2: thanks to the
algorithm for performing the analytic continuation of the solution included in
SeaSyde, we are able to evaluate the 5 MIs of interest at every point of the
physical region and for arbitrary complex variables, thus implementing the CMS
in a straightforward way.

Figure 5: Comparison between real- and complex-valued masses for MI10 and MIs
from 32 to 36. The plots show the value of the $\varepsilon^{0}$ order for
$130\;\text{GeV}\leq\sqrt{s}\leq 190\;\text{GeV}$, $\cos\theta=0$ and
$M^{2}=m_{W}^{2}$ or $M^{2}=\mu_{W}^{2}$. In the MI10 plot, the red crosses
represent the analytical value obtained with complex masses with GiNaC.

In Fig. 5 we report the comparison between real- and complex-valued internal
masses for six selected MIs, evaluated with SeaSyde. The plots present the MIs
as a function of $\sqrt{s}$ at fixed $\cos\theta=0$, where we label with
$\theta$ the scattering angle in the centre of mass reference frame. The
latter can be written in terms of the kinematic invariant $t$ introduced in
Eq. (10) with the following relation:

$t=-\frac{s}{2}\left(1-\cos\theta\right).$ (13)

For each plot in Fig. 5, the blue lines represent the masters for real-valued
masses, while the orange the complex-valued ones. The solid and dashed lines
represent the real and imaginary part of the solution, respectively.

The first plot shows the result for MI 10, defined as follows:

$\text{Master 10}\;:\>\\{1,1,2,0,0,0,1,0,0\\}\;.$ (14)

It provides a validation of our code: for the first 31 MIs, indeed, we have an
analytic expression in terms of GPLs that can be checked against SeaSyde, also
for complex values of the internal mass. In the plot we can observe an
excellent agreement between our result and the red crosses, that represent the
value of the analytic result, evaluated with GiNaC [22, 23] for a few selected
points.

The remaining 5 plots show the comparison between real- and complex-valued
internal masses for MI32-36, and represent an original result of this paper.
We can observe that the masters with real-valued masses exhibit a non
differentiable behaviour around the internal threshold $\sqrt{s}=2m_{W}$,
while they become smooth when moving to the CMS, making the latter
particularly important for phenomenological studies in the threshold region.

Figure 6: Comparison between different values of the decay rate $\Gamma$ for
MI35. The plot shows the value of the $\varepsilon^{0}$ order for
$150\;\text{GeV}\leq\sqrt{s}\leq 190\;\text{GeV}$ and $\cos{\theta}=0$ with
the mass of the $W$ boson. The blue line corresponds to the real-valued mass
case, while the other colours are for the complex-valued case with different
$\Gamma$ values. The solid and dashed lines represent the real and imaginary
part of the solution.

The definition of the complex mass given in Eq. (11) shows that the case of
real masses can be thought also as the limit for $\Gamma\to 0$. Hence, by
considering values of $\Gamma$ gradually smaller, all the MIs should converge
to the one obtained with real-valued internal masses. This is shown explicitly
for MI35 in Fig. 6, although it holds for all the other MIs. In particular we
have compared the results for real- and complex-valued masses, where the
latter have been calculated using different values of the decay-width
$\Gamma$. In the plot we observe that all the lines approach the blue one,
representing the real mass case, as expected. We note that in the case of
complex internal masses the dimensionless kinematic variables $(x,y)$ become
complex. Hence, a Feynman prescription is not necessary, since we can always
link two points in the complex plane in a consistent way, as described in
Section 2.3. On the contrary, when working with real-valued masses, $(x,y)$
may lay on a branch-cut, namely the real axis. In this case providing a
Feynman prescription is mandatory for performing the analytic continuation in
an unambiguous way. In conclusion, this example besides showing that SeaSyde
can robustly handle arbitrary internal complex masses, indicates also that the
complex mass scheme is consistent with the Feynman prescription, and,
moreover, it can be used to automatically implement the causality of the
theory.

###  3.3 Checks and performances

As anticipated in the previous section, we have checked our numerical results
against other publicly available tools, namely GiNaC, DiffExp, AMFlow [68]. In
particular, we have used DiffExp and/or AMFlow to check the value of the 36
MIs considering internal real-valued masses in all the regions of the phase
space. We also have compared our results against pySecDec [7], but only in few
points of the Euclidean region, with a limited number of significant digits.
We have used GiNaC to check the first 31 masters with internal complex-valued
masses and arbitrary values of the Mandelstam invariants. The value of the
last 5 MIs with internal complex masses is, instead, a new prediction.

The strength of the series expansion approach is that we can achieve an
arbitrary level of precision, by simply adding more terms in the expansion.
The increase in precision, however, corresponds to an increase in the amount
of time and computational resources necessary to carry out the calculation.
The two aspects, hence, have to be carefully balanced.

As a check, we evaluated the MIs in different points of the phase space, below
and above the thresholds at $\sqrt{s}=m_{V}$ and $\sqrt{s}=2m_{V}$ and for
different values of $\cos\theta$, using a different number of terms in the
series. In particular, using $150$ terms, we have obtained perfect agreement
between SeaSyde, DiffExp and GiNaC up to the 40th decimal digit, i.e.
quadruple floating-point precision.

Regarding the time performances, they strongly depend not only on the number
of terms in the series, but also on the presence of singularities between the
starting and ending points. In Table 1 we report the time necessary for
transporting the boundary conditions from the Euclidean region to a test
point, namely $\sqrt{s}=155$ GeV and $\cos\theta=0$, for different number of
terms, along with the precision of the result.

Number of terms | Precision | Execution time
---|---|---
50 terms | $10^{-14}$ |  $\sim 14$ min
75 terms | $10^{-19}$ |  $\sim 26$ min
100 terms | $10^{-25}$ |  $\sim 50$ min
125 terms | $10^{-33}$ |  $\sim 75$ min
150 terms | $10^{-40}$ |  $\sim 90$ min
Table 1: Precision achieved and execution time for transporting the boundary
conditions from the Euclidean region to $\sqrt{s}=155$ GeV and $\cos\theta=0$,
for different number of terms in the series.

This estimate is obtained on a MacBook Pro (2015) equipped with a 2.7 GHz
Dual-Core Intel i5. The fact that the execution time increases so much with
the number of terms, is mainly due to the poor-performances of the Mathematica
function Series when dealing with series with a big number of terms. Anyhow,
for most current precision studies, a $10^{-19}$ precision is sufficient. The
execution time depends also on the distance between the starting point, where
the boundary conditions are imposed, and the final one, in which we would like
to evaluate the solution. For reaching the physical region we usually need a
number between $15$ and $20$ steps. Given the results in Table 1 we can
estimate a total of $50$, $160$ and $300$ seconds per step for $50$, $100$ and
$150$ terms, respectively. Once in the physical region, less singularities
need to be avoided and, hence, more direct path can be chosen. In this case a
number between $1$ and $10$ steps is sufficient to reach every point in the
phase-space.

Given the time performances of our package, and in general of packages
implementing the series expansion approach, a direct implementation in Monte
Carlo generators is out of question. However it is possible to reformulate the
problem in such a way that an evaluation on-the-fly is not necessary: we can
use SeaSyde to obtain a numerical grid for the MIs first; a grid for the total
correction can be computed later, covering the whole phase space of the
process; once this latter grid is ready, we can obtain the value of the
correction at any point, simply by a linear interpolation, thanks to the
smoothness of the total correction. In Ref. [53] we presented a numerical grid
for the MI32-36. This grid is constituted by $(130\times 25)$ points which
samples the phase-space for different values of $(\sqrt{s},\cos\theta)$ in the
range $\sqrt{s}\in[40,13000]$ GeV and $\cos\theta\in[-1,1]$. The full
evaluation was performed in parallel on 26 cores on an Intel Xeon Silver 4110
2.1 GHz, considering $75$ terms, and took 11 hours in total. The development
of interpolation techniques is an active field of research (see e.g. [69, 70,
71] ), in part because of the increasing evaluation time per phase-space point
of the two-loop virtual amplitudes. A detailed discussion of these
alternatives is however beyond the scope of this paper.

##  4 Conclusions

The solution via series expansion of multi-loop Feynman integrals has received
a renewed attention after the work of Ref. [42] and the publication of the
package DiffExp [43]. The definition of the masses of the gauge bosons in the
Complex-Mass-Scheme implies the need to evaluate multi-loop Feynman integrals
with complex-valued internal massive lines and the latter could not be studied
with the public version of DiffExp. In this paper we have presented a new
implementation of the series expansion method in the package SeaSyde, with an
original algorithm to perform the analytic continuation of the results, in the
whole complex plane of the kinematical invariants. The usage of this package
allows the consistent evaluation of two-loop and higher-order electroweak
corrections to scattering processes of interest at the LHC and future
colliders.

## Acknowledgments

S.D. and A.V. are supported by the Italian Ministero della Università e della
Ricerca (grant PRIN201719AVICI_01). R. B. is partly supported by the italian
Ministero della Università e della Ricerca (MIUR) under grant PRIN 20172LNEEZ.
T.A. acknowledges the support of the F.R.S.-FNRS under the grant
F.R.S.-FNRS–PDR–T.0142.18.

##  Appendix A Package documentation

The SeaSyde package (Series Expansion Approach for SYstems of Differential
Equations) can be imported in Mathematica using the command Get[…], i.e.
through << SeaSyde.m. Note that it has been developed and tested on
Mathematica 12.0, no previous version of Mathematica has been tested.

Below we present all different functions and their functionalities.

  * 1.

CurrentConfiguration[]

It returns the current configuration of the package. The configuration
parameters can be modified using the function UpdateConfiguration.

  * 2.

UpdateConfiguration[NewConfig_]

NewConfig must be a List whose elements are replacement rules NameParameter ->
NewValue. See Table 2 for a complete overview on all the parameters that the
user can modify.

NameParameter Value type Default Description EpsilonOrder Integer 2 The
maximum order in the dimensional regulator $\varepsilon$ at which the system
is expanded. Note that the minimum order is determined by the boundary
conditions, i.e. if they contain a term $1/\varepsilon^{2}$ the minimum order
will be $-2$. ExpansionOrder Integer 50 The maximum order in the kinematic
variable at which the solution is expanded. InternalWorkingPrecision Integer
250 Specifies the number of digits that are used in internal calculations. If
it is too high, the execution will require more time and space, if it is too
low, we may face rounding errors. ChopPrecision Integer 150 Replaces number
which are smaller than $10^{-\texttt{ChopPrecision}}$ by exactly $0$ in
internal computations. LogarithmicExpansion Bool False It specifies which
method to use for transporting boundary conditions. If it is set to True,
SeaSyde will expand also on top of singularities. RadiusOfConvergence Integer
2 It controls how fast we move at every expansion. If RadiusOfConvergence is
$n$, and the maximum radius of convergence of the solution is $r$, then the
new point will be distant $r/n$ from the center of the series. Note that $r$
is determined internally at every step, based on the position of
singularities. LogarithmicSingularities List {} The user can explicitly state
which singularities are of Logarithmic type, i.e. develop a branch-cut, and
which ones are not. This might speed up the evaluation of numerical grids
since it allows more direct paths in the phase-space. The format in which the
singular points are passed must the one returned by CheckSingularities[].
SafeSingularities List {} Same as LogarithmicSingularities.

Table 2: All the parameters that can be modified by the user.

  * 3.

ReadFrom[FilePath_]

Utility function which reads from the specified path and returns the content
of the file.

  * 4.

SetSystemOfDifferentialEquation[System_, BCs_, MIs_, Variables_, PointBC_,
Param_:{}]

It sets all the internal variables of the package and prepare the system of
differential equations. It receives in input

    * (a)

System_: the system of differential equations. The system equations must be
given in triangular form and it must be ordered so that, order by order in
$\varepsilon$, every equation contains only MIs from previous equations. If
there are multiple kinematic variables, for example x and y, the first $n$
equations in System must be the ones with respect to x, while the last $n$ the
ones with respect to y, where $n$ is the number of MIs. The equations can be
given expanded in $\varepsilon$ or in a closed form in $\varepsilon$.
Example:

{ B1(1,0)[x,y] == 0 , B2(1,0)[x,y] == (- 1/x - /x)B2[x,y] - B1[x,y] / x,
B1(0,1)[x,y] == 0 , B2(0,1)[x,y] == 0 }

    * (b)

BCs_: the boundary conditions for the given equation. They can be given in a
closed form in $\varepsilon$ or as a series. They can also be given as an
asymptotic limit. They can be exact or floating numbers, in this case make
sure that the precision of the boundary condition is sufficient for your final
precision goal. Note that the argument of the function is ignored, hence, the
point in which the boundary conditions are imposed is only fixed by the
PointBC parameter.
Example:

{ B1[1,1] == - 1/32 + /32 + 2( - 1/32 + 2/96 ), B2[1,1] == - Log[2]/16 +
2(-2/192 + Log[2]2/8) }

    * (c)

MIs_: the list of master integrals, as they appear in the equations.
Example:

{ B1[x,y], B2[x,y] }

    * (d)

Variables_: the list of variables appear in the equations, together with their
Feynman prescriptions.
Example:

{ x - I , y + I }

    * (e)

PointBC_: the point in the phase-space in which the boundary conditions are
imposed.
Example:

{ 1, 1 }

    * (f)

Param_: this is an optional parameter. Some equations might contain some
external parameters, for example some masses Mw, Mz. This substitutions are
performed before solving the system.
Example:

{ MW -> 80.38, MZ -> 91.19 }

  * 5.

GetSystemOfDifferentialEquation[] and
GetSystemOfDifferentialEquationExpanded[]

They return the system of differential equations before and after it has been
expanded in $\varepsilon$. These functions can be used to check if everything
has been set correctly.

  * 6.

SolveSystem[Variable_]

It solves the system of differential equations with respect to the kinematic
variable Variable. The series solution in centred in the point where the
boundary conditions are imposed. After solving the system of differential
equations, it is possible to obtain the solution through Solution[],
SolutionValue[] or SolutionTable[].

  * 7.

GetPoint[]

It returns the current point in which the boundary conditions are imposed.

  * 8.

TransportVariable[Variable_, Destination_, Line_:{}]

It transports the boundary conditions for the variable Variable from the
current point to Destination. After transporting the boundary conditions, the
point in which the boundary conditions are imposed is updated to Destination.
The Line parameter is optional. If the user is not satisfied by the path
automatically chosen by the package can use their own. The Line object must be
created with CreateLine

  * 9.

CreateLine[Points_]

It returns a line object that can be used in TransportVariable

  * 10.

TransportBoundaryConditions[PhaseSpacePoint_]

PhaseSpacePoint must be a List whose length is given by the number of
kinematic variables. Its first element must be the final value for the first
variable, its second element the value for the second variable, and so on. The
order of the kinematic variables is the same passed as an input
inSetSystemOfDifferentialEquation. After transporting the boundary conditions,
the point in which they are imposed is updated to PhaseSpacePoint.

  * 11.

Solution[]

It returns the series solution in the current point. The coefficients of the
series are given with InternalWorkingPrecision digits. The result is given as
a List and every MI as a Laurent series in $\varepsilon$.

  * 12.

SolutionValue[]

It returns the value of the MIs, in the centre of the series, as a Laurent
expansion in $\varepsilon$. The coefficients of the $\varepsilon$-series are
given with InternalWorkingPrecision digits.

  * 13.

SolutionTable[]

It returns the value of the MIs, in the centre of the series, as a List of
List going from the first to the last master and from the minimum to the
maximum order in $\varepsilon$. The coefficients of the $\varepsilon$-series
are given with InternalWorkingPrecision digits.

  * 14.

CheckSingularities[]

It checks whether the singularities are logarithmic, i.e. if they develop a
branch-cut. The check is done by performing a path round the singularity and
checking if the solution has developed, or not, an imaginary part for at least
one of the coefficients. If it is not, it means that in doing so we crossed a
branch-cut and, hence, the singularity is logarithmic. If this function is not
called, SeaSyde consider every singularity as logarithmic, and the analytic
continuation is still possible. However, if the user is planning to make an
intensive use of the function TransportBoundaryConditions, e.g. create a
numerical grid, knowing the position of non-logarithmic singularities might
allow for more direct paths in the complex plane. The output of
CheckSingularities can be passed back in the SafeSingularities and
LogarithmicSingularities parameters for future runs.

  * 15.

CreateGraph[MI_, EpsOrder_, Left_, Right_, OtherFunctions_:{}]

Draws a ReImPlot of the solution of order EpsOrder for the master MI. The
graph runs from Left to Right. The argument OtherFunctions may contains other
functions to be plotted in the same graph.

In the Example/ folder of the GitHub repository of SeaSyde, the user can find
working examples to play with.

##  Appendix B Exact boundary conditions

We report the exact boundary conditions we used for solving the 36x36 system
of differential equations for the neutral-current Drell-Yan. The MIs are
labelled $\mathcal{T}_{i}$ and numbered as in Ref. [66]. The boundary
conditions are taken in the asymptotic limit

$x\to 0^{+}\;,\qquad y=1,$ (15)

where $x$ and $y$ are the kinematic variables introduced in Section 3.

$\displaystyle\mathcal{T}_{1}$
$\displaystyle=\frac{1}{16\varepsilon^{2}}\bigg{[}-\frac{1}{2}-\frac{\varepsilon}{2}-\varepsilon^{2}\Big{(}\frac{1}{2}+\zeta_{2}\Big{)}+\varepsilon^{3}\Big{(}-\frac{1}{2}-\zeta_{2}+\zeta_{3}\Big{)}$
$\displaystyle\qquad\quad+\varepsilon^{4}\Big{(}-\frac{1}{2}-\zeta_{2}-\frac{9}{5}\zeta_{2}^{2}+\zeta_{3}\Big{)}\bigg{]},$
$\displaystyle\mathcal{T}_{2}$
$\displaystyle=\frac{1}{16\varepsilon^{2}}\bigg{[}\varepsilon+\varepsilon^{2}+\varepsilon^{3}\Big{(}1+2\zeta_{2}\Big{)}+\varepsilon^{4}\Big{(}1+2\zeta_{2}-2\zeta_{3}\Big{)}\bigg{]},$
$\displaystyle\mathcal{T}_{3}$
$\displaystyle=\frac{1}{16\varepsilon^{2}}\bigg{[}-1-2\varepsilon+\varepsilon^{2}\Big{(}-2-2\zeta_{2}\Big{)}+\varepsilon^{3}\Big{(}-2-4\zeta_{2}+2\zeta_{3}\Big{)}$
$\displaystyle\qquad\quad+\varepsilon^{4}\Big{(}-2-4\zeta_{2}-\frac{18}{5}\zeta_{2}^{2}+4\zeta_{3}\Big{)}\bigg{]},$
$\displaystyle\mathcal{T}_{4}$
$\displaystyle=\frac{1}{16x\varepsilon^{2}}\bigg{[}1-\varepsilon\log(x)+\varepsilon^{2}\Big{(}\frac{\log^{2}(x)}{2}-\zeta_{2}\Big{)}+\varepsilon^{3}\Big{(}-2\zeta_{3}+\zeta_{2}\log(x)-\frac{1}{6}\log^{3}(x)\Big{)}$
$\displaystyle\qquad\quad+\varepsilon^{4}\Big{(}-\frac{9}{10}\zeta_{2}^{2}+2\zeta_{3}\log(x)-\frac{1}{2}\zeta_{2}\log^{2}(x)+\frac{\log^{4}(x)}{24}\Big{)}\bigg{]},$
$\displaystyle\mathcal{T}_{5}$
$\displaystyle=\frac{1}{16\varepsilon^{2}}\bigg{[}-1+2\varepsilon^{2}\zeta_{2}+10\varepsilon^{3}\zeta_{3}+\frac{22}{5}\varepsilon^{4}\zeta_{2}^{2}\bigg{]},$
$\displaystyle\mathcal{T}_{6}$
$\displaystyle=\frac{1}{16x\varepsilon^{3}}\bigg{[}\frac{1}{4}-\frac{1}{2}\varepsilon\log(x)+\frac{1}{2}\varepsilon^{2}\log^{2}(x)+\varepsilon^{3}\Big{(}-\frac{1}{3}\log^{3}(x)-2\zeta_{3}\Big{)}$
$\displaystyle\qquad\quad+\varepsilon^{4}\Big{(}-\frac{6}{5}\zeta_{2}^{2}+4\zeta_{3}\log(x)+\frac{\log^{4}(x)}{6}\Big{)}\bigg{]},$
$\displaystyle\mathcal{T}_{7}$
$\displaystyle=\frac{1}{16\varepsilon^{3}}\bigg{[}\frac{1}{4}-2\varepsilon^{3}\zeta_{3}-\frac{6}{5}\varepsilon^{4}\zeta_{2}^{2}\bigg{]},$
$\displaystyle\mathcal{T}_{8}$
$\displaystyle=\frac{1}{16\varepsilon^{3}}\bigg{[}\varepsilon^{2}+\varepsilon^{3}+\varepsilon^{4}\Big{(}1+2\zeta_{2}\Big{)}\bigg{]},$
$\displaystyle\mathcal{T}_{9}$
$\displaystyle=\frac{1}{16\varepsilon^{2}}\bigg{[}-\frac{\varepsilon}{2}-\frac{3\varepsilon^{2}}{2}+\varepsilon^{3}\Big{(}-\frac{3}{2}-\zeta_{2}\Big{)}+\varepsilon^{4}\Big{(}-\frac{3}{2}-3\zeta_{2}+\zeta_{3}\Big{)}\bigg{]},$
$\displaystyle\mathcal{T}_{10}$
$\displaystyle=\frac{1}{16\varepsilon^{2}}\bigg{[}\varepsilon+\varepsilon^{2}+\varepsilon^{3}\Big{(}1+2\zeta_{2}\Big{)}+\varepsilon^{4}\Big{(}1+2\zeta_{2}-2\zeta_{3}\Big{)}\bigg{]},$
$\displaystyle\mathcal{T}_{11}$
$\displaystyle=\frac{1}{16x\varepsilon^{2}}\bigg{[}-\frac{\varepsilon}{2}+\frac{1}{2}\varepsilon^{2}\log(x)+\varepsilon^{3}\Big{(}-\frac{1}{4}\log^{2}(x)+\frac{\zeta_{2}}{2}\Big{)}$
$\displaystyle\qquad\quad+\varepsilon^{4}\Big{(}-\frac{1}{2}\zeta_{2}\log(x)+\frac{\log^{3}(x)}{12}+\zeta_{3}\Big{)}\bigg{]},$
$\displaystyle\mathcal{T}_{12}$
$\displaystyle=\frac{1}{16\varepsilon^{3}}\bigg{[}\varepsilon^{2}\zeta_{2}+2\varepsilon^{3}\zeta_{3}+2\varepsilon^{4}\zeta_{2}^{2}\bigg{]},$
$\displaystyle\mathcal{T}_{13}$
$\displaystyle=\frac{1}{16\varepsilon^{4}}\bigg{[}\varepsilon^{4}\Big{(}-3-2\zeta_{2}+2\log(x)-\frac{1}{2}\log^{2}(x)\Big{)}\bigg{]},$
$\displaystyle\mathcal{T}_{14}$
$\displaystyle=\frac{1}{16x\varepsilon^{3}}\bigg{[}\frac{\varepsilon}{2}+\varepsilon^{3}\Big{(}-1-\zeta_{2}+\log(x)-\frac{1}{2}\log^{2}(x)\Big{)}$
$\displaystyle\qquad\quad+\varepsilon^{4}\Big{(}-3-\zeta_{2}+2\zeta_{3}+\Big{(}3+\zeta_{2}\Big{)}\log(x)-\frac{3}{2}\log^{2}(x)+\frac{\log^{3}(x)}{2}\Big{)}\bigg{]},$
$\displaystyle\mathcal{T}_{15}$
$\displaystyle=\frac{1}{16\varepsilon^{4}}\bigg{[}\varepsilon^{3}+\varepsilon^{4}\bigg{]},$
$\displaystyle\mathcal{T}_{16}$
$\displaystyle=\frac{1}{16\varepsilon^{3}}\bigg{[}\frac{\varepsilon}{2}+\frac{3\varepsilon^{2}}{2}+\varepsilon^{4}\Big{(}\frac{3}{2}+3\zeta_{2}-\zeta_{3}\Big{)}+\varepsilon^{3}\Big{(}\frac{3}{2}+\zeta_{2}\Big{)}\bigg{]},$
$\displaystyle\mathcal{T}_{17}$
$\displaystyle=\frac{1}{16\varepsilon^{3}}\bigg{[}\varepsilon^{3}(-2+\log(x))+\varepsilon^{4}\Big{(}-6+3\zeta_{2}+3\log(x)-\frac{1}{2}\log^{2}(x)\Big{)}\bigg{]},$
$\displaystyle\mathcal{T}_{18}$
$\displaystyle=\frac{1}{16x\varepsilon^{3}}\bigg{[}-\varepsilon^{2}+\varepsilon^{4}\Big{(}-1+\log(x)-\frac{1}{2}\log^{2}(x)+\zeta_{2}\Big{)}+\varepsilon^{3}(-1+\log(x))\bigg{]},$
$\displaystyle\mathcal{T}_{19}$
$\displaystyle=\frac{1}{16x\varepsilon^{3}}\bigg{[}-\frac{9}{4}+\varepsilon^{3}\Big{(}24\zeta_{3}-12\zeta_{2}\log(x)-\frac{1}{2}\log^{3}(x)\Big{)}$
$\displaystyle\qquad\quad+\varepsilon^{4}\Big{(}-18\zeta_{3}\log(x)+6\zeta_{2}\log^{2}(x)+\frac{3\log^{4}(x)}{8}\Big{)}+12\varepsilon^{2}\zeta_{2}+\frac{3}{2}\varepsilon\log(x)\bigg{]},$
$\displaystyle\mathcal{T}_{20}$
$\displaystyle=\frac{1}{16\varepsilon^{4}}\bigg{[}2\varepsilon^{3}\zeta_{3}+\frac{13}{5}\varepsilon^{4}\zeta_{2}^{2}\bigg{]},$
$\displaystyle\mathcal{T}_{21}$
$\displaystyle=\frac{1}{16\varepsilon^{3}}\bigg{[}-\frac{1}{4}+\varepsilon^{2}\zeta_{2}+5\varepsilon^{3}\zeta_{3}+\frac{27}{5}\varepsilon^{4}\zeta_{2}^{2}\bigg{]},$
$\displaystyle\mathcal{T}_{22}$
$\displaystyle=\frac{1}{16\varepsilon^{4}}\bigg{[}2\varepsilon^{3}\zeta_{3}+\frac{13}{5}\varepsilon^{4}\zeta_{2}^{2}\bigg{]},$
$\displaystyle\mathcal{T}_{23}$
$\displaystyle=\frac{1}{16\varepsilon^{3}}\bigg{[}-\frac{1}{2}+3\varepsilon^{2}\zeta_{2}+13\varepsilon^{3}\zeta_{3}+\frac{57}{5}\varepsilon^{4}\zeta_{2}^{2}\bigg{]},$
$\displaystyle\mathcal{T}_{24}$
$\displaystyle=\frac{1}{16\varepsilon^{3}}\bigg{[}-\varepsilon^{2}\zeta_{2}+\varepsilon^{3}\Big{(}-2\zeta_{2}-3\zeta_{3}\Big{)}+\varepsilon^{4}\Big{(}-4\zeta_{2}-\frac{21}{5}\zeta_{2}^{2}-6\zeta_{3}\Big{)}\bigg{]},$
$\displaystyle\mathcal{T}_{25}$
$\displaystyle=\frac{1}{16\varepsilon^{3}}\bigg{[}-\frac{1}{4}+\varepsilon^{2}\zeta_{2}+5\varepsilon^{3}\zeta_{3}+\frac{27}{5}\varepsilon^{4}\zeta_{2}^{2}\bigg{]},$
$\displaystyle\mathcal{T}_{26}$
$\displaystyle=\frac{1}{16\varepsilon^{3}}\bigg{[}-\varepsilon^{2}+\varepsilon^{3}\Big{(}-4+\zeta_{2}\Big{)}+\varepsilon^{4}\Big{(}-12+2\zeta_{2}+2\zeta_{3}\Big{)}\bigg{]},$
$\displaystyle\mathcal{T}_{27}$
$\displaystyle=\frac{1}{16\varepsilon^{3}}\bigg{[}-\varepsilon^{2}+\varepsilon^{3}\Big{(}2-3\zeta_{2}\Big{)}+\varepsilon^{4}\Big{(}-4+4\zeta_{2}-6\zeta_{3}\Big{)}\bigg{]},$
$\displaystyle\mathcal{T}_{28}$
$\displaystyle=\frac{1}{16\varepsilon^{3}}\bigg{[}\frac{1}{4}+\varepsilon^{2}\Big{(}-1-\zeta_{2}\Big{)}+\varepsilon^{3}\Big{(}1-3\zeta_{2}-5\zeta_{3}\Big{)}+\varepsilon^{4}\Big{(}-1-\frac{27}{5}\zeta_{2}^{2}-9\zeta_{3}+\zeta_{2}\Big{)}\bigg{]},$
$\displaystyle\mathcal{T}_{29}$
$\displaystyle=\frac{1}{16\varepsilon^{4}}\bigg{[}\varepsilon^{4}\Big{(}1+2\zeta_{2}-\log(x)+\frac{\log^{2}(x)}{2}\Big{)}\bigg{]},$
$\displaystyle\mathcal{T}_{30}$
$\displaystyle=\frac{1}{16x\varepsilon^{4}}\bigg{[}\varepsilon^{2}\zeta_{2}-\varepsilon^{3}\zeta_{2}\log(x)+\varepsilon^{4}\Big{(}-\frac{43}{10}\zeta_{2}^{2}+\frac{1}{2}\zeta_{2}\log^{2}(x)\Big{)}\bigg{]},$
$\displaystyle\mathcal{T}_{31}$
$\displaystyle=\frac{1}{16x\varepsilon^{3}}\bigg{[}-\varepsilon+\varepsilon^{2}\Big{(}3+\log(x)-2\zeta_{2}\Big{)}+\varepsilon^{3}\Big{(}-9+11\zeta_{2}-\Big{(}3-2\zeta_{2}\Big{)}\log(x)-\frac{1}{2}\log^{2}(x)\Big{)}$
$\displaystyle+\varepsilon^{4}\Big{(}27-33\zeta_{2}+\frac{43}{5}\zeta_{2}^{2}+26\zeta_{3}+\Big{(}9-11\zeta_{2}\Big{)}\log(x)+\Big{(}\frac{3}{2}-\zeta_{2}\Big{)}\log^{2}(x)+\frac{\log^{3}(x)}{6}\Big{)}\bigg{]},$
$\displaystyle\mathcal{T}_{32}$
$\displaystyle=\frac{1}{16\varepsilon^{4}}\bigg{[}-\varepsilon^{3}\zeta_{2}-6\varepsilon^{4}\zeta_{3}\bigg{]},$
$\displaystyle\mathcal{T}_{33}$
$\displaystyle=\frac{1}{16x\varepsilon^{4}}\bigg{[}-\varepsilon^{2}+\varepsilon^{3}\Big{(}3+\log(x)-2\zeta_{2}\Big{)}+\varepsilon^{4}\Big{(}-9+11\zeta_{2}-\Big{(}3-2\zeta_{2}\Big{)}\log(x)$
$\displaystyle\qquad\quad-\frac{1}{2}\log^{2}(x)\Big{)}\bigg{]},$
$\displaystyle\mathcal{T}_{34}$
$\displaystyle=\frac{1}{16x\varepsilon^{4}}\bigg{[}\varepsilon^{2}+\varepsilon^{3}\Big{(}3-\log(x)-2\zeta_{2}\Big{)}+\varepsilon^{4}\Big{(}5-7\zeta_{2}+\Big{(}-3+2\zeta_{2}\Big{)}\log(x)+\frac{\log^{2}(x)}{2}\Big{)}\bigg{]},$
$\displaystyle\mathcal{T}_{35}$
$\displaystyle=\frac{1}{16x\varepsilon^{4}}\bigg{[}\varepsilon^{2}+\varepsilon^{3}\Big{(}-3-\log(x)+4\zeta_{2}\Big{)}+\varepsilon^{4}\Big{(}9-11\zeta_{2}+12\zeta_{3}+\Big{(}3-4\zeta_{2}\Big{)}\log(x)$
$\displaystyle\qquad\quad+\frac{\log^{2}(x)}{2}\Big{)}\bigg{]},$
$\displaystyle\mathcal{T}_{36}$
$\displaystyle=\frac{1}{16x\varepsilon^{4}}\bigg{[}-\varepsilon^{2}+\varepsilon^{3}\Big{(}-1+\log(x)+2\zeta_{2}-4\zeta_{3}\Big{)}+\varepsilon^{4}\Big{(}-1-\frac{33}{5}\zeta_{2}^{2}+12\zeta_{3}$
$\displaystyle\qquad\quad+\Big{(}1-2\zeta_{2}+4\zeta_{3}\Big{)}\log(x)-\frac{1}{2}\log^{2}(x)+\zeta_{2}\Big{)}\bigg{]},$

where $\zeta_{n}=\zeta(n)$ denotes the Riemann zeta function.

## References

  * [1] G. Heinrich, Collider Physics at the Precision Frontier, Phys. Rept. 922 (2021) 1–69.  arXiv:2009.00516, doi:10.1016/j.physrep.2021.03.006.
  * [2] S. Weinzierl, Feynman Integrals (1 2022).  arXiv:2201.03593.
  * [3] J. L. Bourjaily, et al., Functions Beyond Multiple Polylogarithms for Precision Collider Physics, in: 2022 Snowmass Summer Study, 2022.  arXiv:2203.07088.
  * [4] S. Abreu, R. Britto, C. Duhr, The SAGEX Review on Scattering Amplitudes, Chapter 3: Mathematical structures in Feynman integrals (3 2022).  arXiv:2203.13014.
  * [5] J. Blümlein, C. Schneider, The SAGEX Review on Scattering Amplitudes, Chapter 4: Multi-loop Feynman Integrals (3 2022).  arXiv:2203.13015.
  * [6] A. V. Smirnov, FIESTA4: Optimized Feynman integral calculations with GPU support, Comput. Phys. Commun. 204 (2016) 189–199.  arXiv:1511.03614, doi:10.1016/j.cpc.2016.03.013.
  * [7] S. Borowka, G. Heinrich, S. Jahn, S. P. Jones, M. Kerner, J. Schlenk, T. Zirke, pySecDec: a toolbox for the numerical evaluation of multi-scale integrals, Comput. Phys. Commun. 222 (2018) 313–326.  arXiv:1703.09692, doi:10.1016/j.cpc.2017.09.015.
  * [8] A. Kotikov, Differential equations method: New technique for massive Feynman diagrams calculation, Phys.Lett. B254 (1991) 158–164.  doi:10.1016/0370-2693(91)90413-K.
  * [9] A. V. Kotikov, Differential equation method: The Calculation of N point Feynman diagrams, Phys. Lett. B 267 (1991) 123–127, [Erratum: Phys.Lett.B 295, 409–409 (1992)].  doi:10.1016/0370-2693(91)90536-Y.
  * [10] Z. Bern, L. J. Dixon, D. A. Kosower, Dimensionally regulated pentagon integrals, Nucl. Phys. B 412 (1994) 751–816.  arXiv:hep-ph/9306240, doi:10.1016/0550-3213(94)90398-0.
  * [11] E. Remiddi, Differential equations for Feynman graph amplitudes, Nuovo Cim. A110 (1997) 1435–1452.  arXiv:hep-th/9711188.
  * [12] T. Gehrmann, E. Remiddi, Differential equations for two loop four point functions, Nucl.Phys. B580 (2000) 485–518.  arXiv:hep-ph/9912329, doi:10.1016/S0550-3213(00)00223-6.
  * [13] M. Argeri, P. Mastrolia, Feynman Diagrams and Differential Equations, Int.J.Mod.Phys. A22 (2007) 4375–4436.  arXiv:0707.4037, doi:10.1142/S0217751X07037147.
  * [14] J. M. Henn, Multiloop integrals in dimensional regularization made simple, Phys.Rev.Lett. 110 (2013) 251601.  arXiv:1304.1806, doi:10.1103/PhysRevLett.110.251601.
  * [15] J. M. Henn, Lectures on differential equations for Feynman integrals, J. Phys. A48 (2015) 153001.  arXiv:1412.2296, doi:10.1088/1751-8113/48/15/153001.
  * [16] A. Goncharov, Polylogarithms in arithmetic and geometry, Proceedings of the International Congress of Mathematicians 1,2 (1995) 374–387.
  * [17] A. B. Goncharov, Multiple polylogarithms, cyclotomy and modular complexes, Math. Res. Lett. 5 (1998) 497–516.  arXiv:1105.2076, doi:10.4310/MRL.1998.v5.n4.a7.
  * [18] A. B. Goncharov, Multiple polylogarithms and mixed Tate motives (3 2001).  arXiv:math/0103059.
  * [19] E. Remiddi, J. Vermaseren, Harmonic polylogarithms, Int.J.Mod.Phys. A15 (2000) 725–754.  arXiv:hep-ph/9905237, doi:10.1142/S0217751X00000367.
  * [20] T. Gehrmann, E. Remiddi, Numerical evaluation of harmonic polylogarithms, Comput.Phys.Commun. 141 (2001) 296–312.  arXiv:hep-ph/0107173, doi:10.1016/S0010-4655(01)00411-8.
  * [21] T. Gehrmann, E. Remiddi, Numerical evaluation of two-dimensional harmonic polylogarithms, Comput.Phys.Commun. 144 (2002) 200–223.  arXiv:hep-ph/0111255, doi:10.1016/S0010-4655(02)00139-X.
  * [22] C. W. Bauer, A. Frink, R. Kreckel, Introduction to the GiNaC framework for symbolic computation within the C++ programming language (2000).  arXiv:cs/0004015.
  * [23] J. Vollinga, S. Weinzierl, Numerical evaluation of multiple polylogarithms, Comput.Phys.Commun. 167 (2005) 177.  arXiv:hep-ph/0410259, doi:10.1016/j.cpc.2004.12.009.
  * [24] R. Bonciani, G. Degrassi, A. Vicini, On the Generalized Harmonic Polylogarithms of One Complex Variable, Comput. Phys. Commun. 182 (2011) 1253–1264.  arXiv:1007.1891, doi:10.1016/j.cpc.2011.02.011.
  * [25] L. Naterop, A. Signer, Y. Ulrich, handyG —Rapid numerical evaluation of generalised polylogarithms in Fortran, Comput. Phys. Commun. 253 (2020) 107165.  arXiv:1909.01656, doi:10.1016/j.cpc.2020.107165.
  * [26] L. Adams, C. Bogner, A. Schweitzer, S. Weinzierl, The kite integral to all orders in terms of elliptic polylogarithms, J. Math. Phys. 57 (12) (2016) 122302\.  arXiv:1607.01571, doi:10.1063/1.4969060.
  * [27] E. Remiddi, L. Tancredi, An Elliptic Generalization of Multiple Polylogarithms, Nucl. Phys. B 925 (2017) 212–251.  arXiv:1709.03622, doi:10.1016/j.nuclphysb.2017.10.007.
  * [28] J. Broedel, C. Duhr, F. Dulat, L. Tancredi, Elliptic polylogarithms and iterated integrals on elliptic curves. Part I: general formalism, JHEP 05 (2018) 093.  arXiv:1712.07089, doi:10.1007/JHEP05(2018)093.
  * [29] J. Broedel, C. Duhr, F. Dulat, L. Tancredi, Elliptic polylogarithms and iterated integrals on elliptic curves II: an application to the sunrise integral, Phys. Rev. D 97 (11) (2018) 116009.  arXiv:1712.07095, doi:10.1103/PhysRevD.97.116009.
  * [30] J. Broedel, C. Duhr, F. Dulat, B. Penante, L. Tancredi, Elliptic symbol calculus: from elliptic polylogarithms to iterated integrals of Eisenstein series, JHEP 08 (2018) 014.  arXiv:1803.10256, doi:10.1007/JHEP08(2018)014.
  * [31] J. Ablinger, J. Blümlein, A. De Freitas, M. van Hoeij, E. Imamoglu, C. G. Raab, C. S. Radu, C. Schneider, Iterated Elliptic and Hypergeometric Integrals for Feynman Diagrams, J. Math. Phys. 59 (6) (2018) 062305.  arXiv:1706.01299, doi:10.1063/1.4986417.
  * [32] M. Walden, S. Weinzierl, Numerical evaluation of iterated integrals related to elliptic Feynman integrals, Comput. Phys. Commun. 265 (2021) 108020.  arXiv:2010.05271, doi:10.1016/j.cpc.2021.108020.
  * [33] M. Czakon, Tops from Light Quarks: Full Mass Dependence at Two-Loops in QCD, Phys. Lett. B 664 (2008) 307–314.  arXiv:0803.1400, doi:10.1016/j.physletb.2008.05.028.
  * [34] M. K. Mandal, X. Zhao, Evaluating multi-loop Feynman integrals numerically through differential equations, JHEP 03 (2019) 190.  arXiv:1812.03060, doi:10.1007/JHEP03(2019)190.
  * [35] S. Pozzorini, E. Remiddi, Precise numerical evaluation of the two loop sunrise graph master integrals in the equal mass case, Comput. Phys. Commun. 175 (2006) 381–387.  arXiv:hep-ph/0505041, doi:10.1016/j.cpc.2006.05.005.
  * [36] U. Aglietti, R. Bonciani, L. Grassi, E. Remiddi, The Two loop crossed ladder vertex diagram with two massive exchanges, Nucl. Phys. B789 (2008) 45–83.  arXiv:0705.2616, doi:10.1016/j.nuclphysb.2007.07.019.
  * [37] R. N. Lee, A. V. Smirnov, V. A. Smirnov, Solving differential equations for Feynman integrals by expansions near singular points, JHEP 03 (2018) 008.  arXiv:1709.07525, doi:10.1007/JHEP03(2018)008.
  * [38] R. N. Lee, A. V. Smirnov, V. A. Smirnov, Evaluating elliptic master integrals at special kinematic values: using differential equations and their solutions via expansions near singular points, JHEP 07 (2018) 102.  arXiv:1805.00227, doi:10.1007/JHEP07(2018)102.
  * [39] R. Bonciani, G. Degrassi, P. P. Giardino, R. Groeber, A Numerical Routine for the Crossed Vertex Diagram with a Massive-Particle Loop, Comput. Phys. Commun. 241 (2019) 122–131.  arXiv:1812.02698, doi:10.1016/j.cpc.2019.03.014.
  * [40] M. Fael, F. Lange, K. Schönwald, M. Steinhauser, A semi-analytic method to compute Feynman integrals applied to four-loop corrections to the $\overline{\mathrm{MS}}$-pole quark mass relation, JHEP 09 (2021) 152.  arXiv:2106.05296, doi:10.1007/JHEP09(2021)152.
  * [41] M. Fael, F. Lange, K. Schönwald, M. Steinhauser, Massive Vector Form Factors to Three Loops, Phys. Rev. Lett. 128 (17) (2022) 172003.  arXiv:2202.05276, doi:10.1103/PhysRevLett.128.172003.
  * [42] F. Moriello, Generalised power series expansions for the elliptic planar families of Higgs + jet production at two loops, JHEP 01 (2020) 150.  arXiv:1907.13234, doi:10.1007/JHEP01(2020)150.
  * [43] M. Hidding, DiffExp, a Mathematica package for computing Feynman integrals in terms of one-dimensional series expansions, Comput. Phys. Commun. 269 (2021) 108125\.  arXiv:2006.05510, doi:10.1016/j.cpc.2021.108125.
  * [44] R. Bonciani, V. Del Duca, H. Frellesvig, J. M. Henn, M. Hidding, L. Maestri, F. Moriello, G. Salvatori, V. A. Smirnov, Evaluating a family of two-loop non-planar master integrals for Higgs + jet production with full heavy-quark mass dependence, JHEP 01 (2020) 132.  arXiv:1907.13156, doi:10.1007/JHEP01(2020)132.
  * [45] H. Frellesvig, M. Hidding, L. Maestri, F. Moriello, G. Salvatori, The complete set of two-loop master integrals for Higgs + jet production in QCD, JHEP 06 (2020) 093.  arXiv:1911.06308, doi:10.1007/JHEP06(2020)093.
  * [46] S. Abreu, H. Ita, F. Moriello, B. Page, W. Tschernow, M. Zeng, Two-Loop Integrals for Planar Five-Point One-Mass Processes, JHEP 11 (2020) 117.  arXiv:2005.04195, doi:10.1007/JHEP11(2020)117.
  * [47] I. Dubovyk, A. Freitas, J. Gluza, K. Grzanka, M. Hidding, J. Usovitsch, Evaluation of multi-loop multi-scale Feynman integrals for precision physics (1 2022).  arXiv:2201.02576.
  * [48] R. Bonciani, L. Buonocore, M. Grazzini, S. Kallweit, N. Rana, F. Tramontano, A. Vicini, Mixed strong$-$electroweak corrections to the Drell$-$Yan process (6 2021).  arXiv:2106.11953.
  * [49] M. Becchetti, R. Bonciani, V. Del Duca, V. Hirschi, F. Moriello, A. Schweitzer, Next-to-leading order corrections to light-quark mixed QCD-EW contributions to Higgs boson production, Phys. Rev. D 103 (5) (2021) 054037.  arXiv:2010.09451, doi:10.1103/PhysRevD.103.054037.
  * [50] M. Becchetti, F. Moriello, A. Schweitzer, Two-loop amplitude for mixed QCD-EW corrections to $gg\to Hg$ (12 2021).  arXiv:2112.07578.
  * [51] A. Denner, S. Dittmaier, M. Roth, D. Wackeroth, Predictions for all processes e+ e- —$>$ 4 fermions + gamma, Nucl. Phys. B 560 (1999) 33–65.  arXiv:hep-ph/9904472, doi:10.1016/S0550-3213(99)00437-X.
  * [52] A. Denner, S. Dittmaier, M. Roth, L. H. Wieders, Electroweak corrections to charged-current e+ e- —$>$ 4 fermion processes: Technical details and further results, Nucl. Phys. B 724 (2005) 247–294, [Erratum: Nucl.Phys.B 854, 504–507 (2012)].  arXiv:hep-ph/0505042, doi:10.1016/j.nuclphysb.2011.09.001.
  * [53] T. Armadillo, R. Bonciani, S. Devoto, N. Rana, A. Vicini, Two-loop mixed QCD-EW corrections to neutral current Drell-Yan, JHEP 05 (2022) 072.  arXiv:2201.01754, doi:10.1007/JHEP05(2022)072.
  * [54] F. Tkachov, A Theorem on Analytical Calculability of Four Loop Renormalization Group Functions, Phys.Lett. B100 (1981) 65–68.  doi:10.1016/0370-2693(81)90288-4.
  * [55] K. Chetyrkin, F. Tkachov, Integration by Parts: The Algorithm to Calculate beta Functions in 4 Loops, Nucl.Phys. B192 (1981) 159–204.  doi:10.1016/0550-3213(81)90199-1.
  * [56] S. Laporta, High precision calculation of multiloop Feynman integrals by difference equations, Int.J.Mod.Phys. A15 (2000) 5087–5159.  arXiv:hep-ph/0102033, doi:10.1016/S0217-751X(00)00215-7.
  * [57] C. Anastasiou, A. Lazopoulos, Automatic integral reduction for higher order perturbative calculations, JHEP 07 (2004) 046.  arXiv:hep-ph/0404258, doi:10.1088/1126-6708/2004/07/046.
  * [58] C. Studerus, Reduze-Feynman Integral Reduction in C++, Comput.Phys.Commun. 181 (2010) 1293–1300.  arXiv:0912.2546, doi:10.1016/j.cpc.2010.03.012.
  * [59] A. von Manteuffel, C. Studerus, Reduze 2 - Distributed Feynman Integral Reduction (2012).  arXiv:1201.4330.
  * [60] R. N. Lee, Presenting LiteRed: a tool for the Loop InTEgrals REDuction (2012).  arXiv:1212.2685.
  * [61] R. N. Lee, LiteRed 1.4: a powerful tool for reduction of multiloop integrals, J. Phys. Conf. Ser. 523 (2014) 012059.  arXiv:1310.1145, doi:10.1088/1742-6596/523/1/012059.
  * [62] A. V. Smirnov, Algorithm FIRE – Feynman Integral REduction, JHEP 10 (2008) 107\.  arXiv:0807.3243, doi:10.1088/1126-6708/2008/10/107.
  * [63] A. V. Smirnov, FIRE5: a C++ implementation of Feynman Integral REduction, Comput. Phys. Commun. 189 (2014) 182–191.  arXiv:1408.2372, doi:10.1016/j.cpc.2014.11.024.
  * [64] P. Maierhoefer, J. Usovitsch, P. Uwer, Kira - A Feynman Integral Reduction Program (2017).  arXiv:1705.05610.
  * [65] J. Klappert, F. Lange, P. Maierhöfer, J. Usovitsch, Integral reduction with Kira 2.0 and finite field methods, Comput. Phys. Commun. 266 (2021) 108024.  arXiv:2008.06494, doi:10.1016/j.cpc.2021.108024.
  * [66] R. Bonciani, S. Di Vita, P. Mastrolia, U. Schubert, Two-Loop Master Integrals for the mixed EW-QCD virtual corrections to Drell-Yan scattering, JHEP 09 (2016) 091.  arXiv:1604.08581, doi:10.1007/JHEP09(2016)091.
  * [67] K.-T. Chen, Iterated path integrals, Bull. Am. Math. Soc. 83 (1977) 831–879.  doi:10.1090/S0002-9904-1977-14320-6.
  * [68] X. Liu, Y.-Q. Ma, AMFlow: a Mathematica package for Feynman integrals computation via Auxiliary Mass Flow (1 2022).  arXiv:2201.11669.
  * [69] H. A. Chawdhry, M. L. Czakon, A. Mitov, R. Poncelet, NNLO QCD corrections to three-photon production at the LHC, JHEP 02 (2020) 057.  arXiv:1911.00479, doi:10.1007/JHEP02(2020)057.
  * [70] F. Bishara, M. Montull, (Machine) Learning amplitudes for faster event generation (12 2019).  arXiv:1912.11055.
  * [71] R. Winterhalder, V. Magerya, E. Villa, S. P. Jones, M. Kerner, A. Butter, G. Heinrich, T. Plehn, Targeting multi-loop integrals with neural networks, SciPost Phys. 12 (4) (2022) 129.  arXiv:2112.09145, doi:10.21468/SciPostPhys.12.4.129.