%% LyX 2.3.4.2 created this file.  For more info, see http://www.lyx.org/.
%% Do not edit unless you really know what you are doing.
\documentclass[english]{beamer}
\usepackage{amstext}
\usepackage{amssymb}
\usepackage{fontspec}
\setcounter{secnumdepth}{3}
\setcounter{tocdepth}{3}
\usepackage{url}
\usepackage{cancel}
\usepackage{microtype}
\ifx\hypersetup\undefined
  \AtBeginDocument{%
    \hypersetup{unicode=true}
  }
\else
  \hypersetup{unicode=true}
\fi

\makeatletter

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% LyX specific LaTeX commands.
\newcommand{\noun}[1]{\textsc{#1}}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Textclass specific LaTeX commands.
% this default might be overridden by plain title style
\newcommand\makebeamertitle{\frame{\maketitle}}%
% (ERT) argument for the TOC
\AtBeginDocument{%
  \let\origtableofcontents=\tableofcontents
  \def\tableofcontents{\@ifnextchar[{\origtableofcontents}{\gobbletableofcontents}}
  \def\gobbletableofcontents#1{\origtableofcontents}
}
\providecommand*{\code}[1]{\texttt{#1}}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% User specified LaTeX commands.
\input{preamble}

\makeatother

\usepackage{polyglossia}
\setdefaultlanguage[variant=american]{english}
\begin{document}
\title{Loop amplitude calculations\\
on a computer}
\subtitle{Do-It-Yourself guide}
\author{\href{mailto:vitaly.magerya@tx97.net}{Vitaly Magerya}}
\institute{Institut für Theoretische Physik,\\
Karlsruher Institut für Technologie\\
(ITP KIT)}
\date{Hamburg, CAPP 2021}

\makebeamertitle
\global\long\def\V#1{\mathbf{#1}}%

\global\long\def\M#1{\mathbb{#1}}%

\global\long\def\d{\mathrm{d}}%

\begin{frame}{Goals}

Goal of the lecture: take a \emph{loop amplitude}, calculate it \emph{on a computer}:
analytically (and then numerically), from start to finish, with our
own code.
\[
\fig{eetohadrons}
\]

Target quantity: the \emph{total cross-section of $e^{+}e^{-}$ annihilation to hadrons},
\[
\sigma(e^{+}e^{-}\to\text{hadrons})=\sigma(e^{+}e^{-}\to\text{partons});
\]
we shall calculate $\mathcal{O}(\alpha_{s}^{2})$ corrections to it
(and then $\mathcal{O}(\alpha_{s}^{3})$ too).

Model: QCD with \emph{$N_{f}$ massless quarks}, and \emph{$N_{t}$ massive quarks}
of mass $m$.
\end{frame}
%
\begin{frame}{Motivation: the $R$ ratio}

\begin{figure}
\includegraphics[trim=5mm 0mm 0mm 67mm,clip,height=60mm]{fig/rpp2020-sigma_R_ee.pdf}

\caption{$R\equiv\sigma_{\mathrm{tot}}(e^{+}e^{-}\to\text{hadrons})/\sigma_{\mathrm{tot,leading}}(e^{+}e^{-}\to\gamma^{*}\to\text{muons})$,
from \protect\href{https://arxiv.org/abs/hep-ph/0312114}{Ezhela,~Lugovsky,~Zenin~'03},
available at \protect\href{https://pdg.lbl.gov/2020/hadronic-xsections/}{pdg.lbl.gov/2020/hadronic-xsections}.}

\end{figure}
\end{frame}
%
\begin{frame}{Calculating \textit{e}\textsuperscript{+}\textit{e}\textsuperscript{–}
annihilation to hadrons}

The square matrix element of the total cross-section is
\[
\left|M(e^{+}e^{-}\to\mathrm{partons})\right|^{2}=\sum_{n}\int\mathrm{d}\mathrm{PS}_{n}\left|\sum\!\sfig{eetopartons}\right|^{2}.
\]
Often calculating it is easier via the optical theorem:
\[
\left|M\right|^{2}=-2\,\mathrm{Re}\left(\sfig{eeoptical}\right).
\]
Further, we can split this into the leptonic ($L$) and the hadronic
($H$) tensors:
\[
\left|M\right|^{2}=\frac{2}{q^{4}}\,\mathrm{Re}\left(L_{\mu\nu}H^{\mu\nu}\right),\ L_{\mu\nu}\equiv\sfig{eeleptonic},\ H^{\mu\nu}\equiv\sfig{eehadronic}.
\]
All the loop integration and $\alpha_{s}$ corrections are in $H^{\mu\nu}$,
and $L_{\mu\nu}$ is simple.
\end{frame}
%
\begin{frame}{Tensor structures and projectors}

Because $H^{\mu\nu}(q)$ is a Lorentz-covariant tensor, it can only
have this general structure:
\[
H^{\mu\nu}(q)=g^{\mu\nu}F_{1}(q)+q^{\mu}q^{\nu}F_{2}(q).
\]
Because of Ward identities, $q_{\mu}H^{\mu\nu}=H^{\mu\nu}q_{\nu}=0$,
the structure is further restricted to
\[
H^{\mu\nu}(q)=F_{1}(q)\,\left(g^{\mu\nu}-\frac{q^{\mu}q^{\nu}}{q^{2}}\right).
\]
To invert this and get $F_{1}$ from $H^{\mu\nu}$ we must construct
a projector:
\[
F_{1}=H^{\mu\nu}\frac{g_{\mu\nu}}{d-1}.
\]
Calculting scalar values like $F_{1}$ is simpler than tensor values
like $H^{\mu\nu}$, and a tensor decomposition like this along with
projector constuction is almost always needed.
\end{frame}
%
\begin{frame}{Calculation plan}

\begin{enumerate}
\item Generate Feynman diagrams for the process.
\begin{itemize}
\item $F_{1}=\sfig{dia0}+\sfig{dia1}+\sfig{dia8}+\dots$
\item \noun{Qgraf} with \noun{Mathematica} output.
\end{itemize}
\item Apply the Feynman rules.
\begin{itemize}
\item $F_{1}=\int\mathrm{d}^{d}l_{1}\,\mathrm{Tr}\!\left(\gamma^{\mu}\slashed{k}_{1}\gamma^{\nu}\slashed{k}_{2}\right)\cdots+\dots$
\item Custom \noun{Mathematica} code.
\end{itemize}
\item Resolve Dirac and color tensor summation, convert to the scalar integral
families.
\begin{itemize}
\item $F_{1}=N_{f}C_{a}\cdots I_{123}+\dots$
\item \noun{Mathematica} $\to$ \noun{Form}$\to$ \noun{Mathematica}.
\end{itemize}
\item Use IBP relations to reduce to smaller set of “master integrals”.
\begin{itemize}
\item $F_{1}=\left(N_{f}C_{a}\cdots+\dots\right)I_{111}+\dots$
\item \noun{Mathematica} $\to$ \noun{Kira} $\to$ \noun{Mathematica}.
\end{itemize}
\item Evaluate the master integrals.
\begin{itemize}
\item Numerically: sector decomposition with \noun{pySecDec}.
\item Analytically: differential equations with \noun{Fuchsia}.
\end{itemize}
\end{enumerate}
General idea: use \noun{Mathematica} to glue everything together.
\end{frame}
%

\section{Diagram generation}
\begin{frame}{Feynman diagrams with \noun{Qgraf}}

\noun{Qgraf} is a widely used program for Feynman diagram generation
available at \url{http://cfif.ist.utl.pt/~paulo/qgraf.html}.

To generate diagrams with \noun{Qgraf}:
\begin{enumerate}
\item Create \code{qgraf.dat}, defining the process (incoming particles,
outgoing particles, loop count, names of the momenta, vetoed diagrams),
and listing the \emph{model} and \emph{style files}.

See \code{qgraf.dat.example}.
\item Create a \emph{model file} with a list of fields and vertices.

See \code{qgraf-modfile} for the QCD model we'll use.
\item Create a \emph{style file} (output template) defining the output
format.

See \code{qgraf-stylefile} for the Mathematica output style definition
we'll use.

\end{enumerate}

Also introducing: \code{qgraf.sh} and \code{mkdia.py} programs that
will make our \noun{Qgraf} usage simpler.

Use the \noun{Mathematica} notebook \code{show-diagrams.nb} to view
the diagrams.
\end{frame}
%
\begin{frame}{\noun{Qgraf} diagram structure}

Short summary:
\begin{itemize}
\item In \noun{Qgraf} diagrams consist of vertices and propagators (edges).
\item Each vertex has a list of “fields” (legs).
\item Vertex fields (legs) and their momenta are presented as incoming into
the vertex.
\item Both vertices and their legs are identified by their indices.
\item Vertices and internal legs have positive indices.
\item External legs have negative indices: \code{-1}, \code{-3}, \code{-5},
... for incoming particles, \code{-2}, \code{-4}, \code{-6}, ...
for outgoing.
\end{itemize}

\[
\fig{qgraf}
\]
\end{frame}
%

\section{Feynman rules}
\begin{frame}{Notation on a computer}

Feynman rules \emph{in a book}:
\[
\fig{AA1l}=
\]
\[
\int\frac{\mathrm{d}^{d}l}{\left(2\pi\right)^{d}}\,ig_{e}Q_{f_{1}}\delta_{i_{1}i_{2}}\delta_{f_{1}f_{2}}\,\mathrm{Tr}\!\left(\gamma^{\mu}\frac{\slashed{l}-\slashed{q}}{\left(q-l\right)^{2}+i0}\gamma^{\nu}\frac{\slashed{l}}{l^{2}+i0}\right)ig_{e}Q_{f_{2}}\delta_{i_{2}i_{1}}\delta_{f_{2}f_{1}}.
\]

\medskip{}
Feynman rules \emph{on a computer}: 
\begin{align*}
 & \text{\code{In[1]\ :=\ One\ Feynman\ diagram,\ please???}}\\
\\
 & \text{\code{Syntax::sntxf:\ "One\ Feynman\ diagram"}}\\
 & \text{\code{cannot\ be\ followed\ by\ ",\ please???".}}
\end{align*}
\end{frame}
%
\begin{frame}{\noun{Mathematica} notation, I}

To operate on Feynman rules in \noun{Mathematica}, we need to choose
a notation. Any will work; we'll use the following.

Index names:
\begin{itemize}
\item $\mu_{i}$, Lorentz indices, $1\dots d$: \code{lor{[}}$i$\code{{]}};
\item $i_{i}$, fundamental color indices, $1\dots N_{c}$: \code{fun{[}}$i$\code{{]}};
\item $a_{i}$, adjoint color indices, $1\dots N_{a}=\frac{1}{2}\left(N_{c}^{2}-1\right)$:
\code{adj{[}}$i$\code{{]}};
\item $f_{i}$, light quark flavors (up, down, etc), $1\dots N_{f}$: \code{flv{[}}$i$\code{{]}};
\item $t_{i}$, heavy quark flavors (top, maybe bottom), $1\dots N_{t}$:
\code{flvt{[}}$i$\code{{]}};
\item $s_{i}$, Dirac (spinor) indices, $1\dots4$: \code{spn{[}}$i$\code{{]}}.
\end{itemize}

Tensors:
\begin{itemize}
\item $f^{a_{1}a_{2}a_{3}}$ , $SU(N_{c})$ structure constants: \code{colorf{[}}$a_{1}$\code{,}$a_{2}$\code{,}$a_{3}$\code{{]}};
\item $T_{i_{1}i_{2}}^{a}$, $SU(N_{c})$ generators, \code{colorT{[}}$a$\code{,}$i_{1}$\code{,}$i_{2}$\code{{]}};
\item $\left(\gamma^{\mu}\right)_{s_{1}s_{2}}$, Dirac matrices: \code{gammachain{[}gamma{[}}$\mu$\code{{]},}$s_{1}$\code{,}$s_{2}$\code{{]}};
\item $\left(\slashed{p}\right)_{s_{1}s_{2}}$, Dirac slash notation: \code{gammachain{[}slash{[}p{]},}$s_{1}$\code{,}$s_{2}$\code{{]}};
\end{itemize}
\end{frame}
%
\begin{frame}{\noun{Mathematica} notation, II}
\begin{itemize}
\item Propagator denominators:
\begin{itemize}
\item Massless: $\frac{1}{p^{2}+i0}\to$ \code{den{[}p{]}}.
\item Massive: $\frac{1}{p^{2}-m^{2}+i0}\to$~\code{den{[}p,m2{]}}.
\end{itemize}
\item Momenta components: $p^{\mu}\to$ \code{momentum{[}p,}$\mu$\code{{]}}.
\item Scalar products: $p\cdot q\to$ \code{sp{[}p,q{]}}.
\item Delta functions:
\begin{itemize}
\item Quark flavor: $\delta_{f_{1}f_{2}}\to$ \code{deltaf{[}}$f_{1}$\code{,}$f_{2}$\code{{]}}.
\item Heavy quark flavor: $\delta_{t_{1}t_{2}}\to$~\code{deltaft{[}}$t_{1}$\code{,}$t_{2}$\code{{]}}.
\item Generic: $\delta_{xy}\to$ \code{delta{[}}$x$\code{,}$y$\code{{]}}.
\end{itemize}
\item Quark electric charges:
\begin{itemize}
\item Light: $Q_{f_{i}}\to$ \code{chargeQ{[}}$f_{i}$\code{{]}}.
\item Heavy: $Q_{t_{i}}\to$ \code{chargeQt{[}}$t_{i}$\code{{]}}.
\end{itemize}
\end{itemize}
Everything else (products, sums, powers, etc): the normal \noun{Mathematica}
expressions.
\end{frame}
%
\begin{frame}{Feynman rules: propagators}
\begin{itemize}
\item Quark propagator: $i\delta_{f_{1}f_{2}}\delta_{i_{1}i_{2}}\frac{\left(\slashed{p}\right)_{s_{2}s_{1}}}{p^{2}}$.
\item Heavy quark propagator: $i\delta_{t_{1}t_{2}}\delta_{i_{1}i_{2}}\frac{\left(\slashed{p}+\mathds{1}m_{t}\right)_{s_{2}s_{1}}}{p^{2}-m^{2}}$.
\item Gluon propagator: $-i\delta_{a_{1}a_{2}}\left(\frac{g_{\mu_{1}\mu_{2}}}{p^{2}}-\left(\xi-1\right)\frac{p^{\mu_{1}}p^{\mu_{2}}}{\left(p^{2}\right)^{2}}\right)$.
\item Ghost propagator: $i\delta_{a_{1}a_{2}}\frac{1}{p^{2}}$.
\item Photon propagator: not needed, photons are external for $H^{\mu\nu}$.
\end{itemize}
See \code{Amplitude{[}{]}} in \code{allthecode.m} for the implementation.
\end{frame}
%
\begin{frame}{Feynman rules: vertices}
\begin{itemize}
\item Anti-quark/quark/gluon vertex: $ig_{s}\delta_{f_{1}f_{2}}T_{i_{1}i_{2}}^{a_{3}}\left(\gamma^{\mu_{3}}\right)_{s_{1}s_{2}}$.
\item Anti-quark/quark/photon vertex: $ig_{e}Q_{f_{1}}\delta_{f_{1}f_{2}}\delta_{i_{1}i_{2}}\left(\gamma^{\mu_{1}}\right)_{s_{1}s_{2}}$.
\item Anti-ghost/ghost/gluon vertex: $g_{s}f^{a_{3}a_{2}a_{1}}p_{1}^{\mu_{3}}$.
\item Three-gluon vertex: 
\[
g_{s}f^{a_{1}a_{2}a_{3}}\left(g^{\mu_{1}\mu_{2}}\left(p_{1}-p_{2}\right)^{\mu_{3}}+\left(123\to231\right)+\left(123\to312\right)\right).
\]
\item Four-gluon vertex:
\[
-ig_{s}\left(f^{\nu\mu_{1}\mu_{2}}f^{\nu\mu_{3}\mu_{4}}\left(g^{\mu_{1}\mu_{3}}g^{\mu_{2}\mu_{4}}-g^{\mu_{1}\mu_{4}}g^{\mu_{2}\mu_{3}}\right)+\left(1342\right)+\left(1423\right)\right).
\]
Note the fresh summation index $\nu$.
\end{itemize}
\end{frame}
%

\section{Tensor summation}
\begin{frame}{Dirac tensor summation with \noun{Form}}

\noun{Form} can expand traces of gamma matrices:
\begin{align*}
\mathrm{Tr}\left(\gamma^{\mu}\slashed{p}\right)=\, & \left(\gamma^{\mu}\right)_{s_{1}s_{2}}\left(\slashed{p}\right)_{s_{2}s_{1}}\\
=\, & \text{\code{gammachain[gamma[lor[mu]],spn[1],spn[2]]*}}\\
 & \text{\code{*gammachain[slash[p],spn[2],spn[1]]}}\\
\to\, & \text{\code{gammachain[gamma[mu],slash[p],spn[1],spn[1]]}}\\
\to\, & \text{\code{g\_(}}i\text{\code{,mu,p)}}\\
 & \text{\code{traceN} }i;\\
\to\, & \text{\code{p(mu)}}\\
\to\, & \text{\code{momentum[p,lor[mu]]}}
\end{align*}
On the \noun{Form} side:
\begin{itemize}
\item $i$ can be any arbitrary (but unique) integer;
\item \code{mu} must be declared as an index: \code{index mu = d};
\item \code{p} must be declared as a vector: \code{vector p};
\item \code{p} must be a single variable, not an expression.
\end{itemize}
See: \code{example.dirac-trace.frm}.

\end{frame}
%
\begin{frame}{Color tensor summation with \noun{Color.h}}

\noun{Color.h} is a \noun{Form} package for $SU(N)$ tensor summation,
available at \url{https://www.nikhef.nl/~form/maindir/packages/color/}

Usage in summary:
\begin{align*}
\mathrm{Tr}\left(T^{a_{1}}T^{a_{1}}\right)=\, & T_{i_{1}i_{2}}^{a_{1}}T_{i_{2}i_{1}}^{a_{1}}\\
=\, & \text{\code{colorT[ajd[1],fun[1],fun[1]]*}}\\
 & \text{\code{*colorT[ajd[2],fun[2],fun[1]]}}\\
\to\, & \text{\code{T(fun1,fun2,adj1)*T(fun2,fun1,adj1)}}\\
 & \text{\code{\#call\ docolor}}\\
\to\, & \text{\code{NA*I2R}}\\
\to\, & N_{a}\,T_{f}
\end{align*}
See: \code{example.color-trace.frm}.
\end{frame}
%
\begin{frame}{Quark flavor summation}

Flavor-related factors are complicated by the dependence of charge
on flavor ($Q_{f}$). Three cases are relevant:
\begin{align*}
\sfig{flavor.nf} & \sim N_{f},\\
\sfig{flavor.sumqf} & \sim{\textstyle \sum}Q_{f},\\
\sfig{flavor.sumqf2} & \sim{\textstyle \sum}Q_{f}^{2},
\end{align*}
and the same three cases for the heavy quarks ($N_{t}$, $\sum Q_{t}$,
$\sum Q_{t}^{2}$).

No library for this; just some \noun{Form} code that:
\begin{enumerate}
\item Apply each $\delta_{f_{1}f_{2}}$ factor by renaming $f_{1}$ into
$f_{2}$ (if $f_{1}\ne f_{2}$).
\item Recognize the possible remaining cases ($Q_{f}^{2}\delta_{ff}$, $Q_{f}\delta_{ff}$,
$\delta_{ff}$).
\end{enumerate}
See: \code{example.flavor-trace.frm}.
\end{frame}
%

\section{IBP relations}
\begin{frame}{Integral family construction (for IBP)}

An IBP integral family with $L$ loop momenta $l_{i}$, and $E$ external
momenta $p_{i}$, is the set of integrals
\[
I_{\nu_{1},\nu_{2},\dots,\nu_{N}}\equiv\int\frac{\mathrm{d}^{d}l_{1}}{\left(2\pi\right)^{d}}\dots\frac{\mathrm{d}^{d}l_{L}}{\left(2\pi\right)^{d}}\frac{1}{D_{1}^{\nu_{1}}\dots D_{N}^{\nu_{N}}},
\]
where $N=L\left(L+1\right)/2+LE$, and the denominators,
\[
D_{i}=\left(c_{ij}\,l_{j}+c'_{ik}\,p_{k}\right)^{2}-m_{i}^{2}+i0,
\]
are linearly independent when viewed as polynomials in the scalar
products
\[
s\equiv\left\{ l_{i}\!\cdot\!l_{j},\ l_{i}\!\cdot\!p_{j}\right\} .
\]
This ensures that each $s_{i}$ can be uniquely expressed in terms
of $D_{j}$:
\[
D_{i}=M_{ij}\,s_{j}+C_{i}\quad\Rightarrow\quad s_{i}=\left(M^{-1}\right)_{ij}\left(D_{j}-C_{j}\right).
\]
To express an integral in $I_{\nu\cdots}$ notation: replace all $s$
by $D$, count $D_{i}$ powers.
\end{frame}
%
\begin{frame}{Integral symmetries}

Compare these two integrals:
\[
I_{1}=\int\frac{\mathrm{d}^{d}l}{\left(p_{1}-l\right)^{2}\left(p_{1}+p_{2}-l\right)^{2}l^{2}}\text{ and }I_{2}=\int\frac{\mathrm{d}^{d}l'}{\left(p_{1}+l'\right)^{2}\left(p_{2}-l'\right)^{2}l'^{2}}.
\]
To see that they are equal, use Feynman parameters:
\[
I_{i}=\Gamma\!\left(3-\frac{d}{2}\right)\int\mathrm{d}x_{1}\mathrm{d}x_{1}\mathrm{d}x_{1}\,x_{1}x_{2}x_{3}\,\delta(1-x_{1}-x_{2}-x_{3})\,\mathcal{U}_{i}^{3-d}\mathcal{F}_{i}^{d/2-3},
\]
\[
\mathcal{U}_{1}=x_{1}+x_{2}+x_{3},\ \mathcal{F}_{1}=\left(x_{1}+x_{2}\right)x_{3}p_{1}^{2}+\left(x_{1}+x_{3}\right)x_{2}p_{2}^{2}+2x_{2}x_{3}p_{1}\!\cdot\!p_{2},
\]
\[
\mathcal{U}_{2}=x_{1}+x_{2}+x_{3},\ \mathcal{F}_{2}=\left(x_{2}+x_{3}\right)x_{1}p_{1}^{2}+\left(x_{1}+x_{3}\right)x_{2}p_{2}^{2}+2x_{1}x_{2}p_{1}\!\cdot\!p_{2}.
\]
Both expressions become identical under $x_{1,2,3}\leftrightarrow x_{3,2,1}$.
In momenta space this corresponds to $l'=l-p_{1}$.\sidenote{\href{https://arxiv.org/abs/1111.0868}{Pak~'11};~\href{https://github.com/magv/feynson}{\noun{Feynson}}}

Automated implementation: \noun{Feynson} (\href{https://github.com/magv/feynson}{github.com/magv/feynson}).

Example: \code{example.feynson.symmetrize.in}.
\end{frame}
%
\begin{frame}{Scaleless (zero) integral detection}

For efficiency, scaleless integrals should be put to zero early.

Consider the triangle family with one off-shell leg:
\[
I_{a,b,c}\equiv\sfig{triangle}
\]
Two subsectors of this family are zero:
\[
I_{0,b,c}\equiv\sfig{triangle-noa}=0,\qquad\text{and}\qquad I_{a,b,0}\equiv\sfig{triangle-noc}=0.
\]

\emph{Sufficient criteria} (\href{https://arxiv.org/abs/1310.1145}{Lee~'13}):
a family (or a subsector) is zero if there are such $x$-independent
$k_{i}$, that
\[
\sum_{i}k_{i}\,x_{i}\frac{\partial}{\partial x_{i}}\left(\mathcal{F}(x)+\mathcal{U}(x)\right)=\mathcal{F}(x)+\mathcal{U}(x),
\]
where $\mathcal{F}$, $\mathcal{U}$, and $x$ give the corresponding
Feynman parameterization.

Implementation: any IBP solver, also \noun{Feynson}.

Example: \code{example.feynson.zero-sectors.in}.
\end{frame}
%
\begin{frame}{IBP software}

\noun{Fire6} (\href{https://bitbucket.org/feynmanIntegrals/fire}{bitbucket.org/feynmanIntegrals/fire})\sidenote{\href{https://arxiv.org/abs/1901.07808}{Smirnov,~Chukharev~'19}}
\begin{itemize}
\item \emph{Fast and parallel} Laporta-style IBP reduction implementation.
\item Has a (terrible) \noun{Mathematica} interface, but C++ core.
\item Can use modular arithmetic to control intermediate expression swell,
and for much greater parallelizibility (thousands of cores).
\item Requires \noun{LiteRed} to discover symmetries within integral families.
\end{itemize}
\noun{Kira} (\href{https://kira.hepforge.org/}{kira.hepforge.org}):\sidenote{\href{http://arxiv.org/abs/1705.05610}{Maierhöfer,~Usovitsch,~Uwer~'18}}
\begin{itemize}
\item \emph{Fast and parallel} Laporta-style IBP reduction implementation.
\item Automatically finds symmetries within and between families.
\item Optionally uses modular arithmetic methods.\sidenote{\href{https://arxiv.org/abs/2008.06494}{Klappert,~Lange,~et~al~'20}}
\item Main drawback: excessive memory use (e.g. 200GB for a 4-loop problem).
\end{itemize}
\end{frame}
%
\begin{frame}{IBP software, contd.}

\noun{LiteRed} (\href{https://inp.nsk.su/~lee/programs/LiteRed/}{inp.nsk.su/$\sim$lee/programs/LiteRed})\sidenote{\href{https://arxiv.org/abs/1310.1145}{Lee~'13}}
\begin{itemize}
\item Heuristic-driven IBP relation solution for general indices.
\item Written in \noun{Mathematica}, \emph{easy to use}.
\item Not parallelizable, and \emph{too slow} for larger problems.
\item Guarantees the minimal number of masters.
\item Contains auxiliary functions for integral differentiation, Feynman
parameterization, and dimentional recurrence construction.
\end{itemize}
\noun{Forcer} (\href{https://github.com/benruijl/forcer}{github.com/benruijl/forcer})\sidenote{\href{https://arxiv.org/abs/1704.06650}{Ruijl,~Ueda,~Vermaseren~'17}}
\begin{itemize}
\item Hand-crafted reduction for massless 2-point functions up to 4~loops.
\item Written in \noun{Form}, parallelizable.
\item The \emph{fastest thing for 2-point functions}.
\end{itemize}
Many others: \href{https://arxiv.org/abs/1201.4330}{\noun{Reduze}},
\href{https://arxiv.org/abs/hep-ph/0404258}{\noun{Air}}, \href{https://arxiv.org/abs/1707.01710}{\noun{Fmft}}/\href{https://arxiv.org/abs/hep-ph/0009029}{\noun{Matad}},
private implementations, etc.
\end{frame}
%
\begin{frame}{IBP reduction with \noun{Kira}}

Usage in short:\sidenote{\href{https://kira.hepforge.org/}{kira.hepforge.org}}
\begin{itemize}
\item Define kinematics (\code{config/kinematics.yaml}).
\item List integral families (\code{config/integralfamilies.yaml}).
\item Create a jobs file (e.g. \code{jobs.yaml}), defining
\begin{itemize}
\item for which integrals to \emph{write down} IBP relations ($r$ and
$s$ bounds);
\item for which integrals to \emph{solve} IBP relations ($r$, $s$, and
$d$ bounds, or a list).
\end{itemize}
\item Get the results as \noun{Mathematica} (or \noun{Form}) substitution
tables.
\end{itemize}
Guide to the notation:
\begin{itemize}
\item \emph{$r$} is the sum of denominator powers (positive indices);
\item \emph{$s$} is the sum of numerator powers (negative indices);
\item \emph{$d$} is the sum of dots (indices $\ge2$);
\item \emph{$t$} is the number of denominators.
\end{itemize}
\[
\begin{array}{l}
\left.\begin{array}{lcccc}
\hphantom{I_{3,3,3,-3,0}}\, & r & s & d & t\end{array}\right.\\
\left.\begin{array}{lcccc}
I_{1,2,1,0,0} & 4 & 0 & 1 & 3\\
I_{1,1,3,-1,0} & 5 & 1 & 2 & 3
\end{array}\right\} 11100_{2l}=\text{Sector\ id }7\\
\left.\begin{array}{lcccc}
I_{0,2,2,0,0} & 4 & 0 & 2 & 2\\
I_{0,1,2,-2,0} & 3 & 2 & 1 & 2
\end{array}\right\} 01100_{2l}=\text{Sector\ id }6
\end{array}
\]
\end{frame}
%

\section{Sector decomposition}
\begin{frame}{Sector decomposition software}

\noun{pySecDec} (\href{https://github.com/gudrunhe/secdec}{github.com/gudrunhe/secdec}):
\begin{itemize}
\item \noun{Python} core that generates C++ files.
\item Can use different Monte-Carlo integration algorithms: Randomized Quasi-Monte
Carlo~(\href{https://github.com/mppmu/qmc}{\noun{Qmc}}), \noun{Vegas}/\noun{Suave}/\noun{Divonne}/\noun{Cuhre}~(\href{http://www.feynarts.de/cuba/}{\noun{Cuba}}),
\noun{Cquad}~(\href{http://www.gnu.org/software/gsl/}{\noun{Gsl}}).
\item Optional GPU support (with \noun{Qmc}).
\end{itemize}
\noun{Fiesta} (\href{https://bitbucket.org/feynmanIntegrals/fiesta/}{bitbucket.org/feynmanIntegrals/fiesta}):
\begin{itemize}
\item \noun{Mathematica} core, generates/compiles to C++ behind the scenes.
\item Custom Monte-Carlo integration on both CPU and GPU.
\end{itemize}
\end{frame}
%
\begin{frame}{Sector decomposition with \noun{pySecDec}}

Usage overview (details: \href{https://secdec.readthedocs.io/}{secdec.readthedocs.io}):
\begin{enumerate}
\item Setup \noun{pySecDec} by installing dependencies and running \code{make}
in its release archive.
\item Use the \noun{Python} module \code{pySecDec} to define the integral
and generate the code for the integration library.
\item Compile the integration library.
\item Import the integration library from \noun{Python} (or \noun{C++}),
call it to perform integration.
\end{enumerate}
\end{frame}
%

\section{Differential equations}
\begin{frame}{Method of differential equations}

Consider a family of integrals depending on external momenta $p_{i}$
and masses $m_{i}$:
\[
I_{\nu_{1},\nu_{2},\dots,\nu_{N}}\equiv\int\frac{\mathrm{d}^{d}l_{1}}{\left(2\pi\right)^{d}}\dots\frac{\mathrm{d}^{d}l_{L}}{\left(2\pi\right)^{d}}\frac{1}{D_{1}^{\nu_{1}}\dots D_{N}^{\nu_{N}}}=I_{\nu_{1},\nu_{2},\dots,\nu_{N}}\!\left(\left\{ p_{i}\!\cdot\!p_{j},m_{i}^{2}\right\} \right).
\]
With dimensional analysis one of the parameters can be scaled out,
making remaining arguments dimensionless:
\[
I\left(\left\{ p_{i}\!\cdot\!p_{j},m_{i}\right\} \right)=\left(p_{1}^{2}\right)^{\frac{d}{2}L-\sum_{i}\nu_{i}}\,I\left(\left\{ x_{i}\right\} \right),\qquad\left\{ x_{i}\right\} =\left\{ \frac{p_{i}\!\cdot\!p_{j}}{p_{1}^{2}},\frac{m_{i}^{2}}{p_{1}^{2}}\right\} .
\]
Idea: if $\left\{ x_{i}\right\} \ne\emptyset$, we can construct differential
equations in $x_{i}$, and solve them instead of performing the loop
integration directly.

(In practice: just set $p_{1}^{2}=1$, restore it in the end by dimensionality).
\end{frame}
%
\begin{frame}{Constructing differential equations}

Suppose IBP relations were solved, and we have a set of master integrals
$I_{i}$.
\begin{enumerate}
\item Differentiate each $I_{i}$ by one of the parameters, $x=m_{a}^{2}\text{ or }p_{a}\cdot p_{b}$:
\[
\partial_{m_{a}^{2}}I=\int\frac{\mathrm{d}^{d}l_{1}}{\left(2\pi\right)^{d}}\dots\frac{\mathrm{d}^{d}l_{L}}{\left(2\pi\right)^{d}}\,\partial_{m_{a}}\frac{1}{D_{1}^{\nu_{1}}\dots D_{N}^{\nu_{N}}},
\]
\[
\partial_{p_{a}\cdot p_{b}}I=\left(G^{-1}\right)_{ia}p_{i}\!\cdot\!\left(\partial_{p_{b}}I\right),\quad\partial_{p_{a}\cdot p_{a}}I=\frac{1}{2}\left(G^{-1}\right)_{ia}p_{i}\!\cdot\!\left(\partial_{p_{a}}I\right),
\]
where $G$ is the Gram matrix: $G_{ij}\equiv p_{i}\cdot p_{j}$.
\item Express the derivatives as integrals in the same family:
\[
\partial_{x}I_{i}=\sum_{k}C_{k}\,I_{\nu_{1}^{(k)},\nu_{2}^{(k)},\dots,\nu_{N}^{(k)}}.
\]
\item Use the IBP tables to reduce those integrals back to the master integrals,
thus obtaining a linear differential equation system for $I_{i}$:
\begin{align*}
\partial_{x}I_{i} & \overset{\text{IBP}}{=}\sum_{j}M_{ij}\,I_{j}.
\end{align*}
\end{enumerate}
\end{frame}
%
\begin{frame}{Example: self-energy with one mass}

Consider a family of self-energy with one mass integrals
\[
I_{a,b}\equiv\sfig{selfenergy}=\int\frac{\mathrm{d}^{d}l}{\left(2\pi\right)^{d}}\frac{1}{\left(q-l\right)^{2a}\left(l^{2}-m^{2}\right)^{b}}.
\]
This family has two master integrals:
\begin{align*}
I_{1} & \equiv I_{0,1}=\sfig{selfenergy-01}=\int\frac{\mathrm{d}^{d}l}{\left(2\pi\right)^{d}}\frac{1}{l^{2}-m^{2}},\\
I_{2} & \equiv I_{1,1}=\sfig{selfenergy-11}=\int\frac{\mathrm{d}^{d}l}{\left(2\pi\right)^{d}}\frac{1}{\left(q-l\right)^{2}\left(l^{2}-m^{2}\right)}.
\end{align*}
Constructing differential equation system:
\[
\partial_{m^{2}}\underbrace{\begin{pmatrix}\vphantom{\frac{2}{\left(m^{2}\right)}}I_{1}\\
\vphantom{\frac{2}{\left(m^{2}\right)}}I_{2}
\end{pmatrix}}_{\V I}=\begin{pmatrix}\vphantom{\frac{2}{\left(m^{2}\right)}}I_{0,2}\\
\vphantom{\vphantom{\frac{2}{\left(m^{2}\right)}}}I_{1,2}
\end{pmatrix}\overset{\text{IBP}}{=}\underbrace{\begin{pmatrix}\vphantom{\frac{2}{\left(m^{2}\right)}}\frac{-2+d}{2m^{2}} & 0\\
\vphantom{\frac{2}{\left(m^{2}\right)}}\frac{2-d}{2m^{2}\left(m^{2}-q^{2}\right)} & \frac{-3+d}{m^{2}-q^{2}}
\end{pmatrix}}_{\M M}\underbrace{\begin{pmatrix}\vphantom{\frac{2}{\left(m^{2}\right)}}I_{1}\\
\vphantom{\frac{2}{\left(m^{2}\right)}}I_{2}
\end{pmatrix}}_{\V I}.
\]
\end{frame}
%
\begin{frame}{Solution via the epsilon-form}

If we have a differential equation system for $\V I$,
\[
\partial_{x}\V I(x,d)=\M M(x,d)\,\V I(x,d),\qquad\text{where }d=4-2\varepsilon,
\]
we can transform to a new basis $\V J$ via the transformation matrix
$\M T$,
\[
\V I(x,\varepsilon)=\M T(x,\varepsilon)\,\V J(x,\varepsilon),
\]
and for $\V J$ the same differential equation system will look like
\[
\partial_{x}\V J=\M T^{-1}\left(\M M\M T-\partial_{x}\M T\right)\V J\equiv\M M'\,\V J.
\]
Idea: if this system is in an \emph{$\varepsilon$-form} (has the
dependence on $\varepsilon$ factorized),
\[
\partial_{x}\V J(x,\varepsilon)=\varepsilon\,\M S(x)\,\V J(x,\varepsilon),
\]
then the solution for $\V J$ as a series in $\varepsilon$ becomes
trivial:\sidenote{\href{https://arxiv.org/abs/1304.1806}{Henn~'13}}
\[
\V J(x,\varepsilon)\equiv\varepsilon^{k_{0}}\sum_{k=0}^{\infty}\varepsilon^{k}\V J^{(k)}(x),
\]
\[
\sum\varepsilon^{k}\V J^{(k)}=\sum\varepsilon^{k+1}\M S\,\V J^{(k)}\quad\Rightarrow\quad\V J^{(k)}(x)=\int^{x}\M S(x')\,\V J^{(k-1)}(x')\mathrm{d}x'+\V C^{(k)}.
\]
\end{frame}
%
\begin{frame}{Example, cont.: the epsilon-form}

The $\varepsilon$-form can be achieved with a transformation to the
following basis $\V J$:
\[
\V I=\begin{pmatrix}\left(-1+2\varepsilon\right)m^{2} & 0\\
-1+\varepsilon & \left(-1+\varepsilon\right)\left(m^{2}-q^{2}\right)
\end{pmatrix}\V J.
\]
Differential equation system in $\V J$ then has the form
\[
\partial_{m^{2}}\V J=\varepsilon\begin{pmatrix}-\frac{1}{m^{2}} & 0\\
-\frac{1}{q^{2}}\frac{1}{m^{2}}-\frac{1}{q^{2}}\frac{1}{m^{2}-q^{2}} & -\frac{2}{m^{2}-q^{2}}
\end{pmatrix}\V J.
\]
Accordingly, the solution is
\begin{align*}
\V J^{\left(0\right)} & =\V C^{(0)},\\
\V J^{(1)} & =\V C^{(1)}+\begin{pmatrix}-C_{1}^{(0)}\int^{m^{2}}\frac{\mathrm{d}m'^{2}}{m'^{2}}\\
\dots
\end{pmatrix},\\
\V J^{(2)} & =\V C^{(2)}+\begin{pmatrix}-C_{1}^{(1)}\int^{m^{2}}\frac{\mathrm{d}m'^{2}}{m'^{2}}+C_{1}^{(0)}\int^{m^{2}}\frac{\mathrm{d}m'^{2}}{m'^{2}}\int^{m'^{2}}\frac{\mathrm{d}m''^{2}}{m''^{2}}\\
\dots
\end{pmatrix}.
\end{align*}
\end{frame}
%
\begin{frame}{Iterated integrals}

The solution to differential equation in an $\varepsilon$-form always
come as iterated integrals of the form of \emph{multiple polylogarithms}
(a.k.a Goncharov polylogarithms):\sidenote{\href{https://arxiv.org/abs/1105.2076}{Goncharov~'98}}
\[
G(w_{1},w_{2},\dots,w_{n};x)\equiv\int_{0}^{x}\frac{\mathrm{d}t_{1}}{t_{1}-w_{1}}\int_{0}^{t_{1}}\frac{\mathrm{d}t_{2}}{t_{2}-w_{2}}\cdots\int_{0}^{t_{n}}\frac{\mathrm{d}t_{n}}{t_{n}-w_{n}}.
\]
Special case for the trailing zeros:
\[
G(\underbrace{0,\dots,0}_{n};x)\equiv\frac{1}{n!}\log^{n}(x).
\]
Integration and differention, the simple case:
\[
\int\mathrm{d}x\frac{1}{x-a}G(\vec{w};x)=G(a,\vec{w};x)+C,\quad\frac{\mathrm{d}}{\mathrm{d}x}G(w_{1},\vec{w};x)=\frac{1}{x-w_{1}}G(\vec{w};x).
\]
Software: \href{https://www.ginac.de/}{\noun{GiNaC}}, \href{https://www.physik.uzh.ch/data/HPL/}{\noun{Hpl}},
\href{https://www.nikhef.nl/~form/maindir/packages/harmpol/harmpol.html}{\noun{Harmpol}},
\href{https://bitbucket.org/PanzerErik/hyperint/}{\noun{HyperInt}},
\href{https://gitlab.com/pltteam/plt}{\noun{PolyLogTools}}, etc.
\end{frame}
%
\begin{frame}{Example, cont.: the solution in GPLs}

Rewriting the iterated integrals in terms of $G$:
\begin{align*}
\V J^{\left(0\right)} & =\V C^{(0)},\\
\V J^{(1)} & =\V C^{(1)}+\begin{pmatrix}-C_{1}^{(0)}G(0;m^{2})\\
\dots
\end{pmatrix},\\
\V J^{(2)} & =\V C^{(2)}+\begin{pmatrix}-C_{1}^{(1)}G(0;m^{2})+C_{1}^{(0)}G(0,0;m^{2})\\
\dots
\end{pmatrix}.
\end{align*}
It is trivial to generate this answer to any required order. Note
that the answer will only have $w_{i}\in\{0,q^{2}\}$. By rescaling
the weights (or setting $q^{2}$ to~1), this can be turned into $w_{i}\in\left\{ 0,1\right\} $,
the \emph{harmonic polylogarithms}.
\end{frame}
%
\begin{frame}{Multiple polylogarithms, more properties}

\emph{Differentiation in the general case} when $w_{i}$ may depend
on $x$:
\[
\frac{\d}{\d x}G(w_{1},\dots,w_{n};y)=\sum_{i}G(w_{1},\dots,\xcancel{w_{i}},\dots,w_{n};y)\frac{\d}{\d x}\log\frac{w_{i}-w_{i-1}}{w_{i}-w_{i+1}},
\]
where $w_{0}\equiv y$, and $w_{i+1}\equiv0$.

\emph{Integration in the general case} when $w_{i}$ may depend on
$x$ is nontrivial, one must first represent $G$ in a way that the
integration variable only appears in the last argument.\sidenote{\href{https://arxiv.org/abs/0804.1660}{Brown~'08}}

Automated via \noun{HyperInt} (\href{https://bitbucket.org/PanzerErik/hyperint/}{bitbucket.org/PanzerErik/hyperint}):\sidenote{\href{https://arxiv.org/abs/1403.3385}{Panzer~'14}}
\begin{align*}
 & \text{\code{\$\ maple}}\\
 & \text{\code{>\ read\ "HyperInt.mpl";}}\\
 & \text{\code{>\ fibrationBasis(Hlog(x,[y,x,y]),[x,y]);}}\\
 & \text{\code{Hlog(x,[y,y,y])-Hlog(x,[y,0,y])}}
\end{align*}
\end{frame}
%
\begin{frame}{Multiple polylogarithms and their relatives}

The integral representation of multiple polylogarithms ($G$) is equivalent
to the infinite sum representation ($\mathrm{Li}$):
\begin{align*}
\mathrm{Li}_{m_{1},\dots,m_{n}}\!\left(x_{1},\dots,x_{n}\right) & =\sum_{i_{1}>\dots>i_{n}>0}\frac{x_{1}^{i_{1}}}{i_{1}^{m_{1}}}\cdots\frac{x_{n}^{i_{n}}}{i_{n}^{m_{n}}}=\\
 & =\left(-1\right)^{n}G\bigl(\underbrace{0,\dots,0}_{m_{1}-1},\frac{1}{x_{i}},\dots,\underbrace{0,\dots,0}_{m_{2}-1},\frac{1}{x_{1}x_{2}\cdots x_{n}};1\bigr).
\end{align*}
\begin{itemize}
\item Note: there are conflicting conventions for the order of the indices
in the $\mathrm{Li}$ summation. Above is the “physicist” notation
(used in e.g. \noun{Hpl}, \noun{GiNaC}, and the MZV datamine).
\item The “mathematician” notation is reverse: $0<i_{1}<\dots<i_{n}$;
it was used by Goncharov, and is also used in \noun{HyperInt}. The
order of indices in the Multiple Zeta Values is also reversed there.
\end{itemize}
\end{frame}
%
\begin{frame}{Multiple polylogarithms and their relatives, II}
\begin{itemize}
\item \emph{Logarithms} are GPLs with a single weight: 
\[
\log x=G(0;x),\qquad\log\left(\frac{a-x}{a}\right)=G(a;x).
\]
\item \emph{Nielsen’s generalized polylogarithms} are GPLs with $w_{i}\in\{0,1\}$:
\[
S_{n,p}(x)=\left(-1\right)^{p}G(\underbrace{0,...,0}_{n},\underbrace{1,...,1}_{p};x).
\]
\item \emph{Harmonic polylogarithms} (HPLs) are GPLs with $w_{i}\in\{0,\pm1\}$:
\[
H_{\dots,+m,-n,0}(x)=(-1)^{m}G(\dots,\underbrace{0,\dots,0,1}_{m},\underbrace{0,\dots,0,-1}_{n},0;x).
\]
\item \emph{Two-dimensional HPLs} are GPLs with $w_{i}\in\left\{ 0,1,1-z,-z\right\} $.
\item \emph{Multiple Zeta Values} $\zeta_{\vec{w}}$ are just $H_{\vec{w}}(1)$
(in the “physicist” notation).
\end{itemize}
\end{frame}
%
\begin{frame}{Fixing the integration constants}

Differential equations only give the solution up to the integration
constants. Finding these constants is the essential difficulty of
the method. There are many ways.
\begin{itemize}
\item By evaluating the integrals in a limit where they simplify.
\begin{itemize}
\item Large mass limit. Small mass limit.\sidenote{\href{https://www.springer.com/us/book/9783540423348}{Smirnov~'02}}
\item Even massless integrals can be evaluated by adding masses to them,
and connecting the large mass limit to the massless limit via the
differential equations.
\end{itemize}
\item Using the knowledge of the analytic properties of the integrals.
\begin{itemize}
\item E.g.: enforcing regularity in the kinematic limits.\sidenote{\href{https://arxiv.org/abs/hep-ph/0008287}{Gehrmann,~Remiddi~'00}}
\end{itemize}
\item From partial knowledge of the integrals values, such as a Mellin moment.
\begin{itemize}
\item One can integrate over the semi-inclusive integrals to obtain the
fully inclusive ones.\sidenote{\href{https://arxiv.org/abs/1512.02045}{Gituliar~'15}}
\end{itemize}
\end{itemize}
\end{frame}
%
\begin{frame}{Example, cont.: integration constants}

First, we consider $I_{1}$ (the vacuum bubble) “simple”, and
look it up in a book:
\[
\sfig{selfenergy-01}=-\frac{i\pi^{\frac{d}{2}}}{\left(2\pi\right)^{d}}\Gamma\left(1-\frac{d}{2}\right)\left(m^{2}-i0\right)^{\frac{d}{2}-1}.
\]
For $I_{2}$: in the limit of zero mass the massive diagram becomes
massless,
\[
\lim_{m^{2}\to0}\sfig{selfenergy-11}=\sfig{selfenergy-11-nomass}.
\]
This massless integral we also know from a book:
\[
\sfig{selfenergy-11-nomass}=\frac{i\pi^{\frac{d}{2}}}{\left(2\pi\right)^{d}}\frac{\Gamma^{2}\!\left(\frac{d}{2}-1\right)\Gamma\!\left(2-\frac{d}{2}\right)}{\Gamma\!\left(d-2\right)}\left(-q^{2}-i0\right)^{\frac{d}{2}-2}.
\]
Then, we can expand these in series' in $\varepsilon$, compare with
the series from $\varepsilon$-form solution, and fix all $\V C^{(k)}$.

See: \code{example.massive-self-energy.nb}.
\end{frame}
%
\begin{frame}{The fundamental solution}

A \emph{fundamental solution} to $\partial_{x}\V I=\M M\,\V I$ is
an $n$ by $n$ matrix of independent solutions, such that any solution
can be expressed as a linear combination of its columns.

A fundametal solution for a system in an $\varepsilon$-form, $\partial_{x}\V J=\varepsilon\,\M S\,\V J$,
can be constructed as a series in $\varepsilon$:
\[
\M W=\mathds{1}+\varepsilon\int_{x_{0}}^{x}\d x'\,\M S(x')+\varepsilon^{2}\int_{x_{0}}^{x}\d x'\,\M S(x')\int_{x_{0}}^{x'}\d x''\,\M S(x'')+\dots.
\]
The general solution is then just $\M W$ multiplied by a vector of
integration constants $\V C$:
\[
\V J(x,\varepsilon)=\M W(x,\varepsilon)\,\V C(\varepsilon),
\]
where the constants themselves are a series in $\varepsilon$,
\[
\V C(\varepsilon)\equiv\varepsilon^{k_{0}}\sum_{k=0}^{\infty}\varepsilon^{k}\V C^{(k)}.
\]
This is an alternative way to write down a solution for $\V J$, with
the benefit of being immediately extendable to the multivariate case.
\end{frame}
%
\begin{frame}{The multivariate case}

If differential equations in multiple variables are considered, and
a combined $\varepsilon$-form is achieved,
\[
\partial_{x_{i}}\V J(\vec{x},\varepsilon)=\varepsilon\,\M S_{i}(\vec{x})\,\V J(\vec{x},\varepsilon),
\]
then writing down the solution in this case can be made easy by:
\begin{enumerate}
\item Choosing an integration contour along the axes $x_{i}$ in an some
order, for example from $(0,0,\dots)$ to $(x_{1},0,\dots)$, then
to $(x_{1},x_{2},0,\dots)$, etc.
\item Writing down the fundamental solutions along each segment, $\M W_{i}$.

Because segments are chosen such that only $x_{i}$ changes along
each, $\M W_{i}$ can be calculated the same as in single-variate
case.

\end{enumerate}
The general solution for $\V J$ is then
\[
\V J(\vec{x},\varepsilon)=\M W_{n}(x_{1},\dots,x_{n},\varepsilon)\cdots\M W_{n}(x_{1},\varepsilon)\,\V C(\varepsilon).
\]
Overall this result will be a sum of terms of this form:
\[
G({\text{arguments depending on }x_{1},\dots,x_{n-1}};x_{n})\cdots G({\text{constants}};x_{1}).
\]
\end{frame}
%

\section{Lee algorithm}
\begin{frame}{Lee algorithm for finding the epsilon-form}

Overall idea: to go from a differential equation system
\[
\partial_{x}\V I=\M M\,\V I,\quad\text{where }\M M(x,\varepsilon)=\sum_{i}\frac{\M A_{i}(\varepsilon)}{\left(x-x_{i}\right)^{k_{i}}},
\]
to an $\varepsilon$-form
\[
\partial_{x}\V J=\varepsilon\,\M S\,\V J,\quad\text{where }\M S(x)=\sum_{i}\frac{\M S_{i}}{x-x_{i}},
\]
apply a series of simple basis transformation $\V I=\M T\,\V J$,
such that each brings the system a bit closer to an $\varepsilon$-form.\sidenote{\href{https://arxiv.org/abs/1411.0911}{Lee~'14}}
\begin{enumerate}
\item If a higher pole is present (i.e. $k_{i}\ne1$) then use a transformation
that reduces the rank of $\M A_{i}$, eventually eliminating it. (“Fuchsification”).
\item Else, for the eigenvalues of $\M A_{i}$ of the form $n+k\varepsilon$,
use a transformation that shifts $n$ by $\pm1$, eventually setting
it to zero. (“Normalization”).
\item If all $\M A_{i}$ eigenvalues are proportional to $\varepsilon$,
use a transformation that makes the whole $\M A_{i}$ proportional
to $\varepsilon$. (“Factorization”).
\end{enumerate}

\end{frame}
%
\begin{frame}{Lee algorithm, fuchsification}

Consider \emph{a “balance” transformation between $x_{1}$ and $x_{2}$}:
\[
\M T(x,\varepsilon)=\overline{\M P}(\varepsilon)+\frac{x-x_{2}}{x-x_{1}}\M P(\varepsilon),\quad\text{with }\M P^{2}=\M P,\text{ and }\M P+\overline{\M P}=\mathds{1}.
\]
If the matrix $\M M$ has a pole at $x=x_{1}$ of power $n>1$,
\[
\M M=\M A_{-n}\left(x-x_{1}\right)^{-n}+\M A_{-n+1}\left(x-x_{1}\right)^{-n+1}+\dots,
\]

then either 
\begin{enumerate}
\item there is such $\M P$ and $x_{2}$ , that the transformed $\M A_{-n}$
is of lower rank than $\M A_{-n}$ , while the leading expansion order
around $x=x_{2}$ does increase beyound $n=1$; or
\item the $\varepsilon$-form cannot be reached by any rational transformation.
\end{enumerate}
Then, a series of balance transformations can decrease the rank of
all $\M A_{-n}$ with $n\ne1$ to zero, transforming $\M M$ into
the \emph{Fuchsian form}:
\[
\M M(x,\varepsilon)=\sum\frac{\M A_{i}(\varepsilon)}{x-x_{i}}.
\]
\end{frame}
%
\begin{frame}{Lee algorithm, normalization}

Differential equations for master integrals are observed to have a
special feature: when transformed into Fuchsian form,
\[
\M M=\sum\frac{\M A_{i}}{x-x_{i}},
\]
the eigenvalues of $\M A_{i}$ often have the form of $n+k\varepsilon$,
where $n\in\mathbb{N}$.

Now, a balance can be found between $x_{1}$ and $x_{2}$ that does
not spoil the Fuchsian form. Such a balance will shift one of the
eigenvalues of $\M A_{1}$ by $+1$, one of $\M A_{2}$ by $-1$.

Then, a series of such balances can transform $\M M$ into a \emph{normalized Fuchsian form}:
where all $\M A_{i}$ have eigenvalues proportional to $\varepsilon$.
\begin{itemize}
\item Sometimes eigenvalues of the form $\frac{1}{2}+n+k\varepsilon$ are
encountered. If they are present at up to three different $\M A_{i}$
, then it is possible to change the integration variable from $x$
to such $y$, that the differential equation in $y$ has all the eigenvalues
in the form $n+k\varepsilon$.
\item Sometimes eigenvalues are not linear in $\varepsilon$. This often
means the master integral basis is linearly dependent.
\end{itemize}
\end{frame}
%
\begin{frame}{Lee algorithm, factorization}

Finally, once all the eigenvalues of $\M M$ residues are proportional
to $\varepsilon$, then either
\begin{enumerate}
\item the whole matrix can be made proportional to $\varepsilon$ by a transformation
that does not depend on $x$; or
\item the $\varepsilon$-form can not be reached.
\end{enumerate}
Such transformation is searched for via an ansatz.

Note that on this step, and on all previous ones too, there exists
multiple transformations that can achieve the desired form. As a result,
the $\varepsilon$-form is not unique.
\end{frame}
%
\begin{frame}{Lee algorithm, the multivariate case}

If a system of differential equations in multiple variables is considered:
\[
\partial_{x_{i}}\V I(\vec{x},\varepsilon)=\M M_{i}(\vec{x},\varepsilon)\,\V I(\vec{x},\varepsilon),
\]
then it is useful to have a single transformation $\M T(x,\varepsilon)$
that transforms all of the differential equation systems into $\varepsilon$-form
simultaneously,
\[
\partial_{x_{i}}\V J(\vec{x},\varepsilon)=\varepsilon\,\M S_{i}(\vec{x})\,\V J(\vec{x},\varepsilon),\qquad\text{with }\V I=\M T\,\V J.
\]
Single-variable Lee algorithm can be reused for this by:
\begin{enumerate}
\item Reducing $\M M_{1}$ to an $\varepsilon$-form with $\M T_{1}(x_{1},x_{2},\dots,\varepsilon)$.
\item Transforming all the equations with $\M T_{1}$.
\item Reducing $\M M_{2}$ to an $\varepsilon$-form with a such a $\M T_{2}(x_{2},\dots,\varepsilon)$
that is independent of $x_{1}$ and $\varepsilon$.
\item Transforming all the equations with $\M T_{2}$. This will not spoil
the $\varepsilon$-form in $x_{1}$ because $\M T_{2}$ does not depend
on $x_{1}$ or $\varepsilon$.
\item Repeating similarly for the rest of $\M M_{i}.$
\end{enumerate}
\end{frame}
%
\begin{frame}{Epsilon-form software}

\noun{Fuchsia} (\href{https://github.com/magv/fuchsia.cpp}{github.com/magv/fuchsia.cpp})\sidenote{\href{https://arxiv.org/abs/1701.04269}{Gituliar,~Magerya~'17}}
\begin{itemize}
\item Implements the Lee algorithm.
\item Initial version in \noun{Python}/\noun{SageMath}, newer version
in C++/\noun{GiNaC}.
\item Differential equation systems in multiple variables are supported.
\item Can certify if a system is irreducible (with exceptions).
\item Can suggest variable changes that will make it reducible.
\end{itemize}
\noun{Canonica} (\href{https://github.com/christophmeyer/CANONICA}{github.com/christophmeyer/CANONICA})\sidenote{\href{https://arxiv.org/abs/1701.04269}{Meyer~'17}}
\begin{itemize}
\item Uses the rational ansatz method.

Constructs an ansatz for the transformation matrix up to some powers
of the variables, then solves for the coefficients.
\item Written in \noun{Mathematica}.
\item Supports systems with multiple variables.
\end{itemize}
\noun{Epsilon} (\href{\%20https://github.com/mprausa/epsilon}{github.com/mprausa/epsilon})\sidenote{\href{https://arxiv.org/abs/1701.00725}{Prausa~'17}}
\begin{itemize}
\item Implements the Lee algorithm. 
\item Doesn't support systems with multiple variables.
\end{itemize}
\end{frame}
%

\section{Summary}
\begin{frame}{Summary}

We have talked about:
\begin{itemize}
\item Computing loop amplitudes with \noun{Qgraf}, \noun{Form}, \noun{Color.h},
\noun{Feynson}.
\item \noun{Mathematica} as the hot glue.
\item IBP reduction with \noun{Kira}.
\item Numerical evaluation with \noun{pySecDec}.
\item Differential equations, the $\varepsilon$-form, \noun{Fuchsia}.
\end{itemize}
\vspace{5mm}

\begin{center}
\textit{\only<2>{\emph{Thank you for your attention.}}}
\par\end{center}

\end{frame}

\end{document}
