efficiency). The PDEs (8.40) are, however, written in a redundant way due to
the residue matrix constraint. We thus seek to "mod out" by this constraint,
in order to write the system in a minimal basis. To set this up, begin by
writing $\vec{I}\,^{(0)}$ in terms of some dummy symbols, say

$\displaystyle\vec{I}\,^{(0)}=\left[\begin{array}[]{c}I_{1}\\\ \vdots\\\
I_{R}\end{array}\right]\,.$ (8.45)

Recall the definition of the $\mathsf{R}$-matrix from (8.28). Because of the
second equation in (8.40), we have that

$\displaystyle\mathsf{R}\cdot\vec{I}\,^{(0)}=0.$ (8.46)

This equation can be solved to determine which of the $I_{i}$ symbols are
linearly dependent on each other. Removing those linearly dependent entries
from $\vec{I}\,^{(0)}$, the remaining ones are collected into a new
$R^{\prime}$-dimensional basis

$\displaystyle\vec{J}\quad=\quad\left[\begin{array}[]{c}I_{i_{1}}\\\ \vdots\\\
I_{i_{R^{\prime}}}\end{array}\right]\quad=\quad\mathsf{B}\cdot\vec{I}\,^{(0)}\,.$
(8.50)

The newly introduced $\mathsf{B}$-matrix has dimension $R^{\prime}\times R$
and entries $\mathsf{B}_{ij}\in\\{0,1\\}$. In fact, this is precisely the
$\mathsf{B}$-matrix from (8.27). Therefore, both blocks of the
$\mathsf{M}$-matrix (8.32) are now known. By construction, it holds that

$\displaystyle\mathsf{M}\cdot\vec{I}\,^{(0)}=\left[\begin{array}[]{c}\mathsf{B}\\\
\hline\cr\mathsf{R}\end{array}\right]\cdot\vec{I}\,^{(0)}=\left[\begin{array}[]{c}\vec{J}\\\
\hline\cr\mathbf{0}\end{array}\right]\,.$ (8.55)

The Pfaffian matrices $Q_{i}(z^{\prime})$ for the basis $\vec{J}$ now
immediately follow from (8.35), giving the desired, restricted Pfaffian system

$\displaystyle\partial_{i}\bullet\vec{J}(z^{\prime})=Q_{i}(z^{\prime})\cdot\vec{J}(z^{\prime})\quad\text{for}\quad
i=2,\ldots,N\,.$ (8.56)

This system is _simpler_ to solve than the original Pfaffian system (8.36) for
$\vec{I}(z)$ because 1) it has a smaller rank $R^{\prime}<R$, and 2) it
contains one less variable. Once $\vec{J}$ is solved for, then
$\vec{I}\,^{(0)}$ is determined by multiplying equation (8.55) with
$\mathsf{M}^{-1}$.

Algorithm (3) summarizes the preceding discussion. Note that an initial Moser
reduction step is required if the input is not in normal form.

Input:

  * •

Pfaffian system in normal form of rank $R$:
$\partial_{i}\bullet\vec{I}(z)=P_{i}(z)\cdot\vec{I}(z)$ for $i=1,\ldots,N.$

  * •

Integer $k>0$.

Output: The first $k$ coefficient functions of
$\vec{I}(z)=\sum_{n=0}^{\infty}\vec{I}\,^{(n)}(z^{\prime})\,z_{1}^{n}.$

1:Set $\mathsf{R}=\texttt{RowReduce}[P_{1,-1}]$ and
$R^{\prime}=R-\mathrm{rank}\\!\left(\mathsf{R}\right)$.

2:Write $\vec{I}\,^{(0)}=\big{[}I_{1},\ldots,I_{R}\big{]}^{T}$. Eliminate
linearly dependent dummy symbols $I_{i}$ by solving
$\mathsf{R}\cdot\vec{I}\,^{(0)}=0$. Collect linearly independent ones into
$\vec{J}$. Fix $\mathsf{B}$ via $\vec{J}=\mathsf{B}\cdot\vec{I}\,^{(0)}.$

3:Set $\mathsf{M}=\left[\begin{array}[]{c}\mathsf{B}\\\
\hline\cr\mathsf{R}\end{array}\right].$

4:For each $i=2,\ldots,N$, compute
$\big{(}\partial_{i}\bullet\mathsf{M}+\mathsf{M}\cdot
P_{i,0}\big{)}\cdot\mathsf{M}^{-1}$ and save the upper left $R^{\prime}\times
R^{\prime}$ block matrices dubbed $Q_{i}$.

5:Solve $\partial_{i}\bullet\vec{J}=Q_{i}\cdot\vec{J}$ for $\vec{J}$ and set
$\vec{I}\,^{(0)}=\mathsf{M}^{-1}\cdot\left[\begin{array}[]{c}\vec{J}\\\
\hline\cr\mathbf{0}\end{array}\right].$

6:For each $n=1,\ldots,k$, calculate $\vec{I}\,^{(n)}$ using the recursion
(8.41).

7: return $\\{\vec{I}\,^{(0)},\ldots,\vec{I}\,^{(k)}\\}$.

Algorithm 3 : Holomorphic series expansion

A geometric interpretation of the holomorphic restriction is illustrated in
Figure (8.1). There are three axes, one for the variable $z_{2}$ and two for
the vector components of
$\vec{I}\,^{(0)}(z^{\prime})=[I_{1}(0,z_{2}),\,I_{2}(0,z_{2})]^{T},$ where
$I_{1}(z_{1},z_{2})$ and $I_{2}(z_{1},z_{2})$ represent a basis before
restriction.

The  red surface depicts all possible solutions to the rank jump constraint,
i.e. the nullspace of the residue matrix $P_{1,-1}(z_{2})$. For fixed $z_{2}$,
this nullspace spans a 1-dimensional vector space shown as the  red lines. As
$z_{2}$ varies, the lines trace out the whole  red surface555Such a family of
vector spaces is naturally formalized in terms of vector bundles.. The
solutions living on the surface constitute a basis, determined by the row
space of the $\mathsf{B}$-matrix (8.27), for the restricted Pfaffian system.
Since the nullspace is 1-dimensional, the size of this basis is
$R^{\prime}=R-1=1$.

For fixed $z_{2}$, each  purple line represents the 1-dimensional span of the
row vectors in $P_{1,-1}(z_{2})$. Varying $z_{2}$, then these lines trace out
the  purple surface. This surface is thus a depiction of the
$\mathsf{R}$-matrix (8.28).

Since the  red and  purple surfaces are orthogonal to each other, together
they span a 2-dimensional vector space. This implies that the
$\mathsf{M}$-matrix (8.32) is full rank and therefore invertible.

The  yellow curve illustrates a _particular_ solution vector
$\vec{I}\,^{(0)}(z_{2})$, given some initial condition $\vec{I}\,^{(0)}(0)$.
Observe that this solution is constrained to lie on the  red surface - this is
analogous to the situation in classical mechanics where one finds the solution
to some equations of motion given a holonomic constraint. According to (8.40),
the "flow" of the solution vector is governed by the Pfaffian matrix
$P_{2,0}(z_{2})$: given some point $\vec{I}\,^{(0)}(\overset{\circ}{z_{2}})$
on the curve, then the solution is evolved a small step in $z_{2}$ by

$\displaystyle\vec{I}\,^{(0)}(\overset{\circ}{z_{2}}+\Delta z_{2})$
$\displaystyle\ \simeq\ \vec{I}\,^{(0)}(\overset{\circ}{z_{2}})\ +\ \Delta
z_{2}\Big{[}\partial_{2}\bullet\vec{I}\,^{(0)}(z_{2})\Big{]}\Big{|}_{z_{2}=\overset{\circ}{z_{2}}}$
$\displaystyle\ =\ \vec{I}\,^{(0)}(\overset{\circ}{z_{2}})\ +\ \Delta
z_{2}\Big{[}P_{2,0}(\overset{\circ}{z_{2}})\cdot\vec{I}\,^{(0)}(\overset{\circ}{z_{2}})\Big{]}\,.$
(8.57)

Figure 8.1:  The holomorphic restriction of a rank $R=2$ solution vector
$\vec{I}(z_{1},z_{2})=\left[\begin{array}[]{c}I_{1}(z_{1},z_{2})\\\
I_{2}(z_{1},z_{2})\end{array}\right]$ as $z_{1}\to 0$. The  purple and  red
surfaces depict $\mathsf{R}$ and $\mathsf{B}$ respectively.

##### 8.2.4 Logarithmic restriction of solutions

Now suppose that the solution vector $\vec{I}(z)$ to the Pfaffian system
(8.36) is logarithmically singular rather than holomorphic on
$Y^{\prime}=\\{z_{1}=0\\}$. Our goal is to state an algorithm for calculating
the logarithmic series expansion of $\vec{I}(z)$.

The first step is to write down the most general ansatz to a Pfaffian system
with regular singularities. When the Pfaffian system is in normal form, it is
known that the ansatz decomposes into a sum over unique eigenvalues of the
residue matrix [157, Chapter 4]. Write the spectrum of the residue matrix as

$\displaystyle\mathrm{spec}\left(P_{1,-1}(z^{\prime})\right)=\Big{\\{}\underbrace{\lambda_{1},\ldots,\lambda_{1}}_{\Lambda_{1}}\
,\ \underbrace{\lambda_{2},\ldots,\lambda_{2}}_{\Lambda_{2}}\
,\ldots\Big{\\}}\,,$ (8.58)

where $\Lambda_{i}$ is the multiplicity of eigenvalue $\lambda_{i}$. Even
though $P_{1,-1}(z^{\prime})$ depends on $z^{\prime}$, the integrability
relations (4.34) actually dictate that the eigenvalues $\lambda_{i}$ are
independent of $z^{\prime}$ [157, Theorem 12.1]. For FIs in DR, it means that
$\lambda_{i}$ can only depend on $\epsilon$. In all examples known to us, the
dependence is even linear in $\epsilon$:

$\displaystyle\lambda_{i}=a_{i}+b_{i}\epsilon\quad\text{for}\quad
a_{i},b_{i}\in\mathbb{Q}\,.$ (8.59)

We are not aware of a proof of the above statement.

The ansatz for a logarithmically divergent solution vector can now be
represented as

$\displaystyle\vec{I}(z)\quad=\quad\sum_{\lambda\ \in!\
\mathrm{spec}\left(P_{1,-1}\right)}z_{1}^{\lambda}\ \times\
\sum_{n=0}^{\infty}\ z_{1}^{n}\ \times\ \sum_{m=0}^{M_{\lambda}}\
\vec{I}\,^{(\lambda,n,m)}(z^{\prime})\ \times\ \log^{m}(z_{1})\,,$ (8.60)

where "$\in!$" means "unique element of". In DR, $\vec{I}\,^{(\lambda,n,m)}$
can be further expanded in $\epsilon$. Note that there are two sources of
logarithms in that case: from the $\epsilon$-expansion of
$z_{1}^{\lambda_{i}}=z_{1}^{a_{i}+b_{i}\epsilon}$ and from the finite sum over
$m$.

The integer $M_{\lambda}\in\mathbb{Z}_{\geq 0}$ denotes the maximal power of
$\log^{m}(z_{1})$. Its value depends on the Jordan form of the residue matrix.
Decomposing this matrix into blocks labeled by unique eigenvalues

$\displaystyle\texttt{JordanDecomposition}\\!\left[P_{1,-1}\right]=\bigoplus_{\lambda\
\in!\ \mathrm{spec}\left(P_{1,-1}\right)}\,\mathsf{J}^{(\lambda)}\,,$ (8.61)

then666 Recall that a Jordan block has the eigenvalue $\lambda$ on the
diagonal, and $1$’s on the superdiagonal. $0$’s on the superdiagonal
correspond to Jordan blocks of size $1$.

$\displaystyle M_{\lambda}=\big{[}\text{size of the largest Jordan block of}\
\mathsf{J}^{(\lambda)}\big{]}-1\,.$ (8.62)

The ansatz (8.60) can now be inserted into the Pfaffian system (8.36) together
with the $z_{1}$-expansions of $P_{i}(z)$ from (8.14)-(8.15). One thereby
obtains relations among the coefficient functions $\vec{I}\,^{(\lambda,n,m)}$
by comparing the coefficients of $z_{1}^{n}\,\log^{m}(z_{1})$ on both sides of
the equality [157, Chapter 4]. For fixed $\lambda$, among these relations one
finds

$\displaystyle\begin{array}[]{ll}&\partial_{i}\bullet\vec{I}\,^{(\lambda,0,0)}(z^{\prime})\
=\ P_{i,0}(z^{\prime})\ \cdot\
\vec{I}\,^{(\lambda,0,0)}(z^{\prime})\qquad\text{for}\qquad
i=2,\ldots,N\\\\[5.0pt] &\text{Rank jump
constraint}:\quad\big{[}P_{1,-1}(z^{\prime})-\lambda\mathbf{1}\big{]}^{M_{\lambda}+1}\cdot\
\vec{I}\,^{(\lambda,0,0)}(z^{\prime})\ =\ 0\,.\end{array}$ (8.65)

We call this the _logarithmic restriction_ of the solution vector. If
$\lambda=0$, which is the unique integral eigenvalue for a system in normal
form, then (8.65) reduces to the holomorphic restriction (8.40).

The second equation of (8.65) represents a constraint on the vector
$\vec{I}\,^{(\lambda,0,0)}$; it must be a generalized eigenvector for
$P_{1,-1}(z^{\prime})$. Recall that we modded out by the rows of
$P_{1,-1}(z^{\prime})$ in the holomorphic case because of the constraint
$P_{1,-1}(z^{\prime})\cdot\vec{I}\,^{(0)}=0$ (cf. (8.21)). For fixed
$\lambda$, in the logarithmic case we instead mod out the rows of
$\big{[}P_{1,-1}(z^{\prime})-\lambda\mathbf{1}\big{]}^{M_{\lambda}+1}$ (this
would therefore supplant the $\mathcal{N}^{\prime}$ in (8.21)).

The vector $\vec{I}\,^{(\lambda,0,0)}$ hence plays the role of the basis $e$
in (8.26), and we can repeat the construction of Section (8.2.1) by replacing
$\mathsf{R}$ with

$\displaystyle\mathsf{R}^{(\lambda)}=\texttt{RowReduce}\\!\left[\big{(}P_{1,-1}(z^{\prime})-\lambda\mathbf{1}\big{)}^{M_{\lambda}+1}\right]\,.$
(8.66)

This matrix can be inserted into the formula (8.35) together with some basis
matrix $\mathsf{B}^{(\lambda)}$. The result is a restricted Pfaffian system of
size ${R\,^{\prime}}^{(\lambda)}$ for a basis
$\vec{J}\,^{(\lambda)}(z^{\prime})$:

$\displaystyle\partial_{i}\bullet\vec{J}\,^{(\lambda)}(z^{\prime})\ =\
Q_{i}^{(\lambda)}(z^{\prime})\ \cdot\
\vec{J}\,^{(\lambda)}(z^{\prime})\quad\text{for}\quad i=2,\ldots,N\,.$ (8.67)

The rank of this system is precisely the eigenvalue multiplicity of $\lambda$,
which also equals the size of the block matrix $\mathsf{J}^{(\lambda)}$ in
(8.61). That is,

$\displaystyle{R\,^{\prime}}^{(\lambda)}=\Lambda\,.$ (8.68)

The higher-order terms $\vec{I}\,^{(\lambda,n,m)}$ follow from recursion
relations given $\vec{I}\,^{(\lambda,0,0)}$ as input.

  * •

Shifting powers of $\log^{m}(z_{1})$. Set

$\displaystyle\Phi^{(\lambda,m)}=\frac{P_{1,-1}-\lambda\mathbf{1}}{m}\,.$
(8.69)

Then

$\displaystyle\vec{I}\,^{(\lambda,0,m)}\quad=\quad\Phi^{(\lambda,m)}\ \cdot\
\vec{I}\,^{(\lambda,0,m-1)}\quad\text{for}\quad 0<m\leq M_{\lambda}\,.$ (A)

The terminating condition is $\vec{I}\,^{(\lambda,0,M_{\lambda}+1)}=0.$

  * •

Shifting powers of $z_{1}^{n}$. Define

$\displaystyle\Psi^{(\lambda,n)}=\big{[}(\lambda+n)\mathbf{1}-P_{1,-1}\big{]}^{-1}\,.$
(8.70)

It holds that

$\displaystyle\vec{I}\,^{(\lambda,n,m)}\quad=\quad\Psi^{(\lambda,n)}\ \cdot\
\left(\sum_{\nu=0}^{n-1}P_{1,n-\nu-1}\ \cdot\
\vec{I}\,^{(\lambda,\nu,m)}\hskip 5.69046pt-\hskip
5.69046pt(m+1)\vec{I}\,^{(\lambda,n,m+1)}\right)\,.$ (B)

The last term drops out when $m=M_{\lambda}$, leading to777 (B) in fact also
reduces to (A) when $n=0$.

$\displaystyle\vec{I}\,^{(\lambda,n,M_{\lambda})}\quad=\quad\Psi^{(\lambda,n)}\
\cdot\ \sum_{\nu=0}^{n-1}P_{1,n-\nu-1}\ \cdot\
\vec{I}\,^{(\lambda,\nu,M_{\lambda})}\,.$ (C)

The index $n$ runs up to $\infty$, so this recursion ought to be truncated
based on a numerical accuracy condition for the full solution vector
$\vec{I}(z)$.

Starting from $\vec{I}\,^{(\lambda,0,0)}$, we present a flowchart showing how
to use the recursion relations (A), (B) and (C) to compute
$\vec{I}\,^{(\lambda,n,m)}$ for all $0\leq n\leq N_{\lambda}$ and $0\leq m\leq
M_{\lambda}.$ The flowchart is adapted from [157, Page 42], and the tuple
$(n,m)$ represents $\vec{I}\,^{(\lambda,n,m)}$.

$\displaystyle\begin{array}[]{ccccccccc}(0,0)&\overset{\eqref{m_shift}}{\longrightarrow}&(0,1)&\overset{\eqref{m_shift}}{\longrightarrow}&\cdots&\overset{\eqref{m_shift}}{\longrightarrow}&(0,M_{\lambda}-1)&\overset{\eqref{m_shift}}{\longrightarrow}&(0,M_{\lambda})\\\\[10.0pt]
&&&&&&&&\hskip
8.5359pt\rotatebox[origin={c}]{270.0}{$\overset{\eqref{n_shift_simple}}{\longrightarrow}$}\\\\[10.0pt]
(1,0)&\overset{\eqref{n_shift}}{\longleftarrow}&(1,1)&\overset{\eqref{n_shift}}{\longleftarrow}&\cdots&\overset{\eqref{n_shift}}{\longleftarrow}&(1,M_{\lambda}-1)&\overset{\eqref{n_shift}}{\longleftarrow}&(1,M_{\lambda})\\\\[10.0pt]
&&&&&&&&\hskip
8.5359pt\rotatebox[origin={c}]{270.0}{$\overset{\eqref{n_shift_simple}}{\longrightarrow}$}\\\
&&&&&&&&\vdots\\\\[10.0pt] &&&&&&&&\hskip
8.5359pt\rotatebox[origin={c}]{270.0}{$\overset{\eqref{n_shift_simple}}{\longrightarrow}$}\\\\[10.0pt]
(N_{\lambda},0)&\overset{\eqref{n_shift}}{\longleftarrow}&(N_{\lambda},1)&\overset{\eqref{n_shift}}{\longleftarrow}&\cdots&\overset{\eqref{n_shift}}{\longleftarrow}&(N_{\lambda},M_{\lambda}-1)&\overset{\eqref{n_shift}}{\longleftarrow}&(N_{\lambda},M_{\lambda})\end{array}$
(8.78)

Example 8.1. Fix $\lambda$ and drop the $\lambda$-index from
$\vec{I}\,^{(\lambda,n,m)}$, $\Phi^{(\lambda,m)}$ and $\Psi^{(\lambda,n)}$. We
showcase the recursion for $0\leq n\leq 2$ and $0\leq m\leq 2$.

  1. 1.

$n=0$: Using (A) twice,

$\displaystyle\vec{I}\,^{(0,1)}$
$\displaystyle=\Phi^{(1)}\cdot\vec{I}\,^{(0,0)}$ (8.79)
$\displaystyle\vec{I}\,^{(0,2)}$
$\displaystyle=\Phi^{(2)}\cdot\vec{I}\,^{(0,1)}\,.$ (8.80)

  2. 2.

$n=1$: Using (C) once,

$\displaystyle\vec{I}\,^{(1,2)}=\Psi^{(1)}\cdot
P_{1,0}\cdot\vec{I}\,^{(0,2)}\,.$ (8.81)

Using (B) twice,

$\displaystyle\vec{I}\,^{(1,1)}$
$\displaystyle=\Psi^{(1)}\cdot\Big{(}P_{1,0}\cdot\vec{I}\,^{(0,1)}-2\vec{I}\,^{(1,2)}\Big{)}$
(8.82) $\displaystyle\vec{I}\,^{(1,0)}$
$\displaystyle=\Psi^{(1)}\cdot\Big{(}P_{1,0}\cdot\vec{I}\,^{(0,0)}-1\vec{I}\,^{(1,1)}\Big{)}\,.$
(8.83)

  3. 3.

$n=2$: Using (C) once,

$\displaystyle\vec{I}\,^{(2,2)}=\Psi^{(2)}\cdot\Big{(}P_{1,1}\cdot\vec{I}\,^{(0,2)}+P_{1,0}\cdot\vec{I}\,^{(1,2)}\Big{)}\,.$
(8.84)

Using (B) twice,

$\displaystyle\vec{I}\,^{(2,1)}$
$\displaystyle=\Psi^{(2)}\cdot\Big{(}P_{1,1}\cdot\vec{I}\,^{(0,1)}+P_{1,0}\cdot\vec{I}\,^{(1,1)}-2\vec{I}\,^{(2,2)}\Big{)}$
(8.85) $\displaystyle\vec{I}\,^{(2,0)}$
$\displaystyle=\Psi^{(2)}\cdot\Big{(}P_{1,1}\cdot\vec{I}\,^{(0,0)}+P_{1,0}\cdot\vec{I}\,^{(1,0)}-1\vec{I}\,^{(2,1)}\Big{)}\,.$
(8.86)

Note that this recursion can be performed over a finite field for the sake of
efficiency. $\blacksquare$

The discussion above is summarized in the following algorithm.

Input:

  * •

Pfaffian system in normal form:
$\partial_{i}\bullet\vec{I}(z)=P_{i}(z)\cdot\vec{I}(z)$ for $i=1,\ldots,N.$

  * •

Integer $k>0$.

Output: The series expansion coefficients $\vec{I}\,^{(\lambda,n,m)}$ from
(8.60) up to order $\mathcal{O}(z_{1}^{k})$.

1:Compute $\texttt{JordanDecomposition}\\!\left[P_{1,-1}(z^{\prime})\right]$
with $z^{\prime}$ fixed to generic numbers. Save

  * •

unique eigenvalues $\\{\lambda_{1},\lambda_{2},\ldots\\}$

  * •

eigenvalue multiplicities $\\{\Lambda_{1},\Lambda_{2},\ldots\\}$

  * •

integers $\\{M_{\lambda_{1}},M_{\lambda_{2}},\ldots\\}$ from (8.62).

2:for each unique $\lambda$ do

3: Repeat steps 1. to 5. in Algorithm (3) with the replacements

  * •

$\mathsf{R}\leftarrow\mathsf{R}^{(\lambda)}=\texttt{RowReduce}\\!\left[\big{(}P_{1,-1}(z^{\prime})-\lambda\mathbf{1}\big{)}^{M_{\lambda}+1}\right]$

  * •

$R^{\prime}\leftarrow\Lambda$

  * •

$\vec{I}\,^{(0)}\leftarrow\vec{I}\,^{(\lambda,0,0)}$

and save the result for $\vec{I}\,^{(\lambda,0,0)}$.

4: for $0\leq n\leq k$ and $0\leq m\leq M_{\lambda}$ do

5: Input $\vec{I}\,^{(\lambda,0,0)}$ into the recursion relations (A), (B),
(C) to compute $\vec{I}\,^{(\lambda,n,m)}$.

6: end for

7:end for

8: return The coefficient functions $\bigcup_{\lambda}\
\bigcup_{0\,\leq\,n\,\leq k}\ \bigcup_{0\,\leq\,m\,\leq M_{\lambda}}\
\vec{I}\,^{(\lambda,n,m)}$

Algorithm 4 : Logarithmic series expansion

We expect this to be an efficient method for computing series expansions of
MIs. Its efficiency stems from trading a "big" Pfaffian system (satisfied by
the full solution vector $\vec{I}(z)$) for several "small" Pfaffian systems
(satisfied by the restricted bases $\vec{J}\,^{(\lambda)}(z^{\prime})$), each
one containing one less variable. For instance, if the "big" system includes
unrationalizable square roots and elliptic functions, then it is likely that
the "small" systems only contain simple GPLs.

The catch is that the algorithm does not produce results valid in the entire
phase space; there must exist a variable $z_{1}\ll 1$ that can be used for the
expansion. Luckily, there are many such situations in high-energy physics:
$z_{1}$ could be a small mass, a threshold variable (e.g. $s-4m^{2}$), an
angle between two collinear momenta (e.g. in Regge theory), a small or large
transverse momentum ($z_{1}=p_{T}$ or $z_{1}=1/p_{T}$), a soft momentum
variable in Post-Newtonian/Minkowskian gravity, and so forth.

The _method of regions_ is the classic approach to asymptotic expansions [159,
160, 161]. More recent approaches include [162, 163, 164, 165, 166]. It would
be interesting to import some of the wisdom from these studies in order to
enhance the logarithmic restriction approach presented here.

####  Section 8.3Restriction at the level of Macaulay matrices

While the protocol of the previous section is applicable to any regularly
singular Pfaffian system, it does require an unrestricted system as input,
which might be hard to obtain. It should also be mentioned that this method is
only conjectured to yield the rational restriction of the twisted cohomology
group (6.59) when it is viewed as a $\mathcal{D}$-module [83, Remark 3.5].

Building on the Macaulay matrix method introduced in Chapter (5), this section
presents a holomorphic restriction protocol which works entirely at the level
of $\mathcal{D}$-modules. Theorems from [83, Section 4] ensure that the method
yields correct results. Some steps of the construction are technical, so we
refer to the previous reference for details.

##### 8.3.1 Restriction to a hyperplane

Let $\mathcal{I}$ be a holonomic ideal in $\mathcal{D}_{Y}$. This could e.g.
be the annihilating ideal for the GKZ system. The goal is to compute a
Pfaffian system for the rational restriction module onto the hyperplane given
by $Y^{\prime}=\\{z_{1}=\cdots=z_{N^{\prime}}=0\\}.$

The first question to answer is that of a basis for the restriction module. In
[83, Algorithm 4.7], a systematic search algorithm is formulated which
succeeds when the dimension of the basis matches the expected holonomic rank
$R^{\prime}$ of the restriction module. The algorithm requires $R^{\prime}$ as
input, but this is easily obtained for the case of Euler integrals by
computing the Euler characteristic $\chi$ via the likelihood equations (6.60).
Moreover, this algorithm relies entirely on numerical linear algebra, with all
parameters fixed to generic numbers, meaning that it will terminate with
modest computing power even for relatively large Pfaffian systems.

The resulting standard monomials for the restriction module, call them
$\mathrm{RStd}$, now serve as input for the Macaulay matrix $M$ in Algorithm
(1). Recall that $M$ is built by acting with "seed" monomials on the
generators of $\mathcal{I}$. Since $\mathcal{I}$ contains unrestricted
variables defined on $Y$, so will $M$. The key observation is this: all the
entries of $M$ are polynomial in $z\in Y$, so we can simply fix
$z_{1}=\ldots=z_{N^{\prime}}=0$ in $M$, and then proceed with the usual
Macaulay matrix algorithm.

Input:

  * •

The restricted basis $\mathrm{RStd}$.

  * •

Generators $\\{d_{1},\ldots,d_{D}\\}$ of a holonomic ideal $\mathcal{I}$ in
$\mathcal{D}_{Y}$.

  * •

Direction $i$.

Output: Pfaffian matrix $Q_{i}(z^{\prime})$ for the restriction

$\displaystyle\mathcal{R}_{Y^{\prime}}\otimes_{\mathcal{D}_{Y^{\prime}}}\mathcal{D}_{Y}/\big{(}\mathcal{I}+z_{1}\mathcal{D}_{Y}+\ldots+z_{N^{\prime}}\mathcal{D}_{Y}\big{)}\,.$
(8.87)

1:Call Algorithm (1) with the replacements

  * •

$\mathrm{Std}\leftarrow\mathrm{RStd}$

  * •

$M\leftarrow\left[\begin{array}[]{c|c}M_{\mathrm{Ext}}&M_{\mathrm{RStd}}\end{array}\right]\Big{|}_{z_{1}=\ldots=z_{N^{\prime}}=0}$

2:return $Q_{i}(z^{\prime})=C_{\mathrm{RStd}}^{(i)}-C^{(i)}\cdot
M_{\mathrm{RStd}}$

Algorithm 5 : Rationally restricted Pfaffian system by the Macaulay matrix

Let us remark that the basis $\mathrm{RStd}$ can also be guessed, assuming
$R^{\prime}$ is known, by executing Algorithm (5) several times with all
parameters set to numbers until it succeeds.

Algorithm (5) has been implemented in the asir package mt_mm by Nobuki
Takayama, as showcased in the following example.

Example 8.1. Let $\mathcal{I}=\langle d_{1},d_{2}\rangle$ be the ideal
generated by

$\displaystyle d_{1}$
$\displaystyle=z_{1}(1-z_{1})\partial_{1}^{2}-z_{1}z_{2}\partial_{1}\partial_{2}+(c_{1}-z_{1}(a+b_{1}+1))\partial_{1}-b_{1}z_{2}\partial_{2}-ab_{1}$
(8.88) $\displaystyle d_{2}$
$\displaystyle=z_{2}(1-z_{2})\partial_{2}^{2}-z_{1}z_{2}\partial_{1}\partial_{2}+(c_{2}-z_{2}(a+b_{2}+1))\partial_{2}-b_{2}z_{1}\partial_{1}-ab_{2}\,,$
(8.89)

which annihilates Appell’s $F_{2}$ function

$\displaystyle F_{2}(a,b_{1},$ $\displaystyle
b_{2},c_{1},c_{2}|z_{1},z_{2})=\frac{\Gamma(c_{1})\Gamma(c_{2})}{\Gamma(b_{1})\Gamma(b_{2})\Gamma(c_{1}-b_{1})\Gamma(c_{2}-b_{2})}\
\times$
$\displaystyle\int_{(0,1)^{2}}(1-x_{1})^{c_{1}-b_{1}-1}(1-x_{2})^{c_{2}-b_{2}-1}(1-z_{1}x_{1}-z_{2}x_{2})^{-a}x_{1}^{b_{1}}x_{2}^{b_{2}}\frac{\mathrm{d}x}{x}\,.$
(8.90)

This function appears often in the FI calculus [167, 168].

The following asir script computes the holomorphic restriction to
$Y^{\prime}=\\{z_{1}=0\\}$.

⬇

import("mt_mm.rr")$

Vars = [z1, z2]$

DVars = [dz1, dz2]$

D1 = z1*(1-z1)*dz1^2 - z1*z2*dz1*dz2 + (c1-z1*(a+b1+1))*dz1 - b1*z2*dz2 -
a*b1;

D2 = z2*(1-z2)*dz2^2 - z1*z2*dz1*dz2 + (c2-z2*(a+b2+1))*dz2 - b2*z1*dz1 -
a*b2;

Ideal = [D1, D2];

NumericRule = [[z2,z2+1/3], [a,1/2], [b1,1/3], [b2,1/5], [c1,1/7], [c2,1/11]];

IdealNumeric = base_replace(Ideal, NumericRule);

RStd = mt_mm.restriction_to_pt_(IdealNumeric, Gamma = 2, KK = 4, Vars | p = 10^8);

RStd = reverse(map(dp_ptod, RStd[0], DVars));

Ideal = map(dp_ptod, Ideal, DVars);

MMData = mt_mm.find_macaulay(Ideal, RStd, Vars | restriction_var = [z1]);

Q2 = mt_mm.find_pfaffian(MMData, Vars, 2 | use_orig = 1);

The inputs of mt_mm.restriction_to_pt_() merit explanation. The list
IdealNumeric contains generators with $z_{2}$ and exponent parameters fixed to
generic numbers. The integers Gamma and KK refer to the $\gamma$ and $k$ of
the basis search algorithm [83, Algorithm 4.7]. The value p = 10^8 means that
the function is executed over a finite field of order $p=100000007$, which is
more efficient than a symbolic calculation.

The result of this script is the Pfaffian system

$\displaystyle\partial_{2}\,\mathrm{RStd}=Q_{2}\cdot\mathrm{RStd}$ (8.91)

given by

$\displaystyle\partial_{2}\left[\begin{array}[]{c}1\\\
\partial_{2}\end{array}\right]=\left[\begin{array}[]{cc}0&1\\\
\frac{-ab_{2}}{z_{2}(z_{2}-1)}&\frac{c_{2}-(a+b_{2}+1)}{z_{2}(z_{2}-1)}\end{array}\right]\cdot\left[\begin{array}[]{c}1\\\
\partial_{2}\end{array}\right]\,.$ (8.98)

$\blacksquare$

##### 8.3.2 Restriction to a hypersurface

This chapter has, up until now, only dealt with restrictions onto hyperplanes.
Namely, we have assumed that $Y^{\prime}=\\{z_{1}=\ldots=z_{N^{\prime}}=0\\}$
is defined by an affine _linear_ coordinate change. To date, there is no known
algorithm for restrictions of $\mathcal{D}$-modules onto general
hypersurfaces. Surprisingly, the Macaulay matrix restriction can be
generalized to this scenario.

Suppose that $Y^{\prime}=\\{f=0\\}$ is the vanishing locus of an irreducible
polynomial $f=f(z)\in\mathbb{C}[z_{1},\ldots,z_{N}]$. The idea is,
essentially, just to solve the Macaulay matrix equation (5.13) for $C^{(i)}$
modulo $f$:

$\displaystyle C_{\mathrm{Ext}}^{(i)}\ -\ C^{(i)}\cdot M_{\mathrm{Ext}}\equiv
0\quad\text{mod}\ f\,.$ (8.99)

We give one example below to illustrate this procedure. See [169] for further
results on restrictions onto hypersurface singularities of Horn’s
hypergeometric functions $\\{H_{1},\ldots,H_{7}\\}$.

Example 8.2. We compute the restriction of Appell’s $F_{4}$ function

$\displaystyle
F_{4}(a,b,c_{1},c_{2}|z_{1},z_{2})=\frac{\Gamma(c_{1})\Gamma(c_{2})}{\Gamma(a)\Gamma(b)\Gamma(c_{1}-a)\Gamma(c_{2}-b)}\
\times$
$\displaystyle\int_{(0,1)^{2}}(1-x_{1})^{c_{1}-a-1}(1-x_{2})^{c_{2}-b-1}(1-z_{1}x_{1})^{a-c_{1}-c_{2}+1}\
\times$ (8.100) $\displaystyle\hskip
31.2982pt(1-z_{2}x_{2})^{b-c_{1}-c_{2}+1}(1-z_{1}x_{1}-z_{2}x_{1})^{c_{1}+c_{2}-a-b-1}x_{1}^{a}x_{2}^{b}\frac{\mathrm{d}x}{x}$

onto its singular locus defined by

$\displaystyle L=z_{1}z_{2}f=0\quad\text{with}\quad
f=(z_{1}-z_{2})^{2}-2(z_{1}+z_{2})+1\,.$ (8.101)

In terms of Euler operators $\theta_{i}=z_{i}\,\partial_{i}$, the annihilating
ideal $\mathcal{I}=\langle d_{1},d_{2}\rangle$ is generated by

$\displaystyle d_{1}$
$\displaystyle=\theta_{1}(\theta_{1}+c_{1}-1)-z_{1}(\theta_{1}+\theta_{2}+a)(\theta_{1}+\theta_{2}+b)$
(8.102) $\displaystyle d_{2}$
$\displaystyle=\theta_{2}(\theta_{2}+c_{2}-1)-z_{2}(\theta_{1}+\theta_{2}+a)(\theta_{1}+\theta_{2}+b)\,.$
(8.103)

The holonomic rank outside the singular locus is $R=4$. On
$Y^{\prime}=\\{f=0\\}$, the rank drops to $R^{\prime}=3$. Both of these ranks
assume that the exponent parameters $(a,b,c_{1},c_{2})$ are generic.

We search for a restricted basis by executing Algorithm (5) with different
guesses for $\mathrm{RStd}$ until it succeeds. The parameters are set to
generic numbers in this search. The result is

$\displaystyle\mathrm{RStd}=\big{[}1,\partial_{1},\partial_{2}\big{]}^{T}\,.$
(8.104)

Next we build a $6\times 10$ Macaulay matrix
$M_{1}=\big{[}M_{\mathrm{Ext}}\big{|}M_{\mathrm{RStd}}\big{]}$ of degree $1$.
The left block contains $7$ columns labeled by
$\mathrm{Ext}=\big{[}\partial_{1}^{3},\,\partial_{1}^{2}\partial_{2},\,\partial_{1}\partial_{2}^{2},\,\partial_{2}^{3},\,\partial_{1}^{2},\,\partial_{1}\,\partial_{2},\,\partial_{2}^{2}\big{]}^{T}\,.$

The goal is now to solve (8.99) modulo $f$. Since $f$ is quadratic in $z_{1}$,
we can eliminate all instances of $z_{1}^{p}$ for $p\geq 2$ via polynomial
division by $f$. It is therefore enough to consider the ansatz

$\displaystyle C^{(i)}=C_{0}^{(i)}+C_{1}^{(i)}z_{1}\,,$ (8.105)

where $C_{1,2}^{(i)}$ are independent of $z_{1}$. The reduction modulo $f$ can
be performed via a sequence of replacements $z_{1}^{p}\to
g_{p}(z_{2})+h_{p}(z_{2})z_{1}.$ For instance, $f=0$ immediately gives

$\displaystyle z_{1}^{2}\to(-1+2z_{2}-z_{2}^{2})+2(1+z_{2})z_{1}\,.$ (8.106)

Inserting this into $z_{1}f=0$, we further have

$\displaystyle z_{1}^{3}\to
2(-1+z_{2}+z_{2}^{2}-z_{2}^{3})+(3+10z_{2}+3z_{2}^{2})z_{1}\,,$ (8.107)

and so forth.

After solving for $C_{1,2}^{(i)}$, the Pfaffian matrices are computed as usual
by

$\displaystyle Q_{i}=C_{\mathrm{RStd}}^{(i)}-C^{(i)}\cdot
M_{\mathrm{RStd}}\,.$ (8.108)

Setting $(a,b,c_{1},c_{2})=(-2/3,\,1/3,\,1/3,\,1/3)$ to simplify the output,
we finally get

$\displaystyle Q_{1}$ $\displaystyle=\left[\begin{array}[]{ccc}0&1&0\\\
-\frac{5(3z_{1}z_{2}+z_{1}-5z_{2}(z_{2}+2)-1)}{162(z_{2}-1)^{3}z_{2}}&\frac{z_{1}(41z_{2}-5)-z_{2}(87z_{2}+62)+5}{108(z_{2}-1)^{2}z_{2}}&\frac{5(3z_{1}z_{2}+z_{1}-5z_{2}(z_{2}+2)-1)}{108(z_{2}-1)^{3}}\\\
\frac{4(-z_{1}+2z_{2}+2)}{81(z_{2}-1)^{2}z_{2}}&-\frac{2}{27z_{2}}&\frac{2(z_{1}-2(z_{2}+1))}{27(z_{2}-1)^{2}}\end{array}\right]$
(8.112) $\displaystyle Q_{2}$
$\displaystyle=\left[\begin{array}[]{ccc}0&0&1\\\
\frac{4(-z_{1}+2z_{2}+2)}{81(z_{2}-1)^{2}z_{2}}&-\frac{2}{27z_{2}}&\frac{2(z_{1}-2(z_{2}+1))}{27(z_{2}-1)^{2}}\\\
\frac{5(-z_{1}+3z_{2}+1)}{162(z_{2}-1)z_{2}^{2}}&-\frac{5(z_{1}+z_{2}-1)}{108z_{2}^{2}}&\frac{5z_{1}-51z_{2}+31}{108(z_{2}-1)z_{2}}\end{array}\right]$
(8.116)

See the script [169, 2023-01-26-F4-on-sing.py] for comparisons of known
analytic results with numerical solutions to this Pfaffian system.
$\blacksquare$

### Chapter 9Restriction Examples

This chapter gives examples of how to compute Pfaffian systems for restriction
modules related to FIs. Sections (9.2) and (9.3) also include partial results
for logarithmic series expansions related to 1- and 2-loop Bhabha scattering.

####  Section 9.1Holomorphic restriction: 2-loop N-box

This is a continuation of Example (8). The goal is to apply the holomorphic
restriction protocol of Section (8.2.3) to that GKZ system. Schematically, we
restrict

$\displaystyle A=\left[\begin{array}[]{ccccc ccccc}1&1&1&1&1&1&1&1&1&1\\\
1&1&1&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&1&\scalebox{0.7}{$\cdot$}\\\
1&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&1&1&1&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&1&1\\\
\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&1&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&1&1&\scalebox{0.7}{$\cdot$}&1\\\
\scalebox{0.7}{$\cdot$}&1&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&1&\scalebox{0.7}{$\cdot$}&1&\scalebox{0.7}{$\cdot$}&1&\scalebox{0.7}{$\cdot$}\\\
\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&1&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&1&\scalebox{0.7}{$\cdot$}&1&\scalebox{0.7}{$\cdot$}&1\end{array}\right]\quad\xrightarrow[i\,\neq\,10]{z_{i}=1}\quad\parbox[c]{0.0pt}{\includegraphics{Figs/2L_nbox.pdf}}\,.$
(9.7)

We begin with a rank $R=9$ Pfaffian system, obtained via the Macaulay matrix
method, in terms of variables $z\in Y=\mathbb{A}^{10}.$ After a coordinate
shift
$(z_{1},\ldots,z_{9},z_{10})\to(z_{1},\ldots,z_{9},z_{10})+(1,\ldots,1,0),$
the restriction surface is $Y^{\prime}=\\{z_{1}=\ldots=z_{9}=0\\}.$ Six
variables can immediately be restricted using homogeneity (cf. Section (6.5)),
say $\\{z_{1},z_{2},z_{3},z_{4},z_{5},z_{9}\\}.$ So the restriction protocol
need only be applied to $\\{z_{6},z_{7},z_{8}\\}$. By calculating the Euler
characteristic on $Y^{\prime}$, we expect a rank $R^{\prime}=3$ Pfaffian
system in the end.

##### 9.1.1 Basis

The method of [118] (implemented in the asir command mt_gkz.cbase_by_euler(A))
gives a basis of 9 standard monomials

$\displaystyle\mathrm{Std}=\left[\begin{array}[]{c}\partial_{8}^{2}\\\
\partial_{6}\partial_{10}\\\ \partial_{7}\partial_{10}\\\ \partial_{10}^{2}\\\
\partial_{6}\\\ \partial_{7}\\\ \partial_{8}\\\ \partial_{10}\\\
1\end{array}\right]\,.$ (9.17)

The Macaulay matrix method then outputs a Pfaffian system

$\displaystyle\partial_{i}\,\mathrm{Std}=P_{i}\cdot\mathrm{Std}\quad\text{for}\quad
i=6,7,8,10\,.$ (9.18)

For the sake of verifying our results, an independent IBP calculation is
performed for the following MIs in $\mathrm{D}=4-2\epsilon$ dimensions:

$\displaystyle\left[\begin{array}[]{c}I_{02202}\\\ I_{22020}\\\
I_{11111}\end{array}\right]=G_{0}\cdot e^{(I)}\,.$ (9.22)

The subscripts on the integrals label powers of the inverse propagators

$\displaystyle D_{1}=\ell_{1}^{2}\,,\ D_{2}=\ell_{2}^{2}\,,\
D_{3}=(p_{1}-\ell_{1})^{2}\,,\ D_{4}=(p_{12}-\ell_{12})^{2}\,,\
D_{5}=(-p_{123}+\ell_{12})^{2}\,,$ (9.23)

and $G_{0}$ is a matrix containing powers of $s=2p_{1}\cdot p_{2}$ that render
$e^{(I)}$ dimensionless.

We would like to land on the basis (9.22) after computing the restriction of
(9.18), in order to compare our Pfaffian matrices with those coming from IBPs.
By the formulas of Section (6.3.3), we therefore translate the integrals
(9.18) into $\mathcal{D}$-module elements:

$\displaystyle e^{(I)}\to
e^{(\mathcal{D})}=\mathrm{c}\cdot\left[\begin{array}[]{c}\partial_{10}^{2}\\\
\epsilon(2+5\delta)(1+2\epsilon+5\epsilon\delta)+2z_{10}(1+2\epsilon+5\epsilon\delta)\partial_{10}+z_{10}^{2}\partial_{10}\\\
(1+\epsilon+2\epsilon\delta)\partial_{10}+z_{10}\partial_{10}^{2}-z_{7}\partial_{7}\partial_{10}+z_{6}\partial_{6}\partial_{10}\end{array}\right]\,.$
(9.27)

The $3\times 3$ diagonal matrix $\mathrm{c}$ contains prefactors depending on
$\\{\epsilon,\delta,s\\}$. Thus, there is a $3\times 9$ $\mathsf{B}$-matrix
(8.27) with rows and columns labeled by $e^{(\mathcal{D})}$ and $\mathrm{Std}$
respectively:

$\displaystyle\mathsf{B}=\left[\begin{array}[]{ccccc
ccccc}\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&1&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}\\\\[4.0pt]
\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&z_{10}^{2}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&2z_{10}(1+2\epsilon+5\epsilon\delta)&\epsilon(2+5\delta)(1+2\epsilon+5\epsilon\delta)\\\\[4.0pt]
\scalebox{0.7}{$\cdot$}&z_{6}&-z_{7}&z_{10}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&1+\epsilon+2\epsilon\delta&\scalebox{0.7}{$\cdot$}\end{array}\right]\,.$
(9.31)

This matrix satisfies111 One may employ (A.10) in cases where
$e^{(\mathcal{D})}$ contains some monomials that are not contained in
$\mathrm{Std}$.
$e^{(\mathcal{D})}=\mathrm{c}\cdot\mathsf{B}\cdot\mathrm{Std}.$

##### 9.1.2 Normal form

The Pfaffian system (9.18) must be brought to normal form before the
restriction procedure can be applied (cf. Section (8.1.3)). The system turns
out to be in normal form w.r.t. the variables $z_{6}$ and $z_{7}$, but $P_{8}$
has a second-order pole in $z_{8}$:

$\displaystyle P_{8}\big{|}_{z_{6},\,z_{7}\to
1}=(z_{8}-1)^{2}\,P_{8,-2}\big{|}_{z_{6},\,z_{7}\to 1}+\ldots\,.$ (9.32)

The second-order pole can be cured by Moser reduction, namely a gauge
transformation by

$\displaystyle
G_{1}=\texttt{Diagonal}\\!\left[(z_{8}-1)^{-1},\,1,\ldots,1\right]\,.$ (9.33)

The notation

$\displaystyle G[P_{i}]=G^{-1}\cdot(P_{i}\cdot G-\partial_{i}\bullet G)$
(9.34)

is used throughout this chapter to denote gauge transformations. In addition
to curing the second-order pole, the new residue matrices
$G[P_{i}]_{-1},i\in\\{6,7,8\\},$ also obtain non-resonant spectra.
Consequently, the system is in normal form. To simplify notation, let us

$\displaystyle\text{replace}\ G_{1}[P_{i}]\ \to\ P_{i}\ \text{in the rest of
this example}.$ (9.35)

##### 9.1.3 Restriction

The Pfaffian-level restriction protocol implores us to build the matrix

$\displaystyle\mathsf{M}=\left[\begin{array}[]{c}\mathsf{B}\\\
\hline\cr\mathsf{R}\end{array}\right]$ (9.38)

coming from (8.32). The $\mathsf{B}$-matrix was already found in (9.31). The
$\mathsf{R}$-matrix is constructed from the independent rows of all the
(gauged transformed) residue matrices:

$\displaystyle\mathsf{R}=\texttt{RowReduce}\\!\left[\begin{array}[]{c}P_{6,-1}\\\
\hline\cr P_{7,-1}\\\ \hline\cr P_{8,-1}\end{array}\right]\,.$ (9.42)

We remind that the $\texttt{RowReduce}\\!\left[...\right]$ operation includes
the deletion of zero-rows.

The $\mathsf{M}$-matrix can now be inserted into (8.35), giving

$\displaystyle\big{(}\partial_{10}\bullet\mathsf{M}+\mathsf{M}\cdot
P_{10,0}\big{)}\cdot\mathsf{M}^{-1}=\left[\begin{array}[]{c|c}Q_{10}&\star\\\
\hline\cr\mathbf{0}_{6\times 3}&\star\end{array}\right]\,,$ (9.45)

where

$\displaystyle P_{10,0}=P_{10}\big{|}_{z_{6},\,z_{7},\,z_{8}\to 0}\,,$ (9.46)

which is well-defined because the system is in normal form. Sending the
regulator $\delta\to 0$ in $Q_{10}$, the final result is

$\displaystyle
Q_{10}=\left[\begin{array}[]{ccc}-\frac{2(1+\epsilon)}{z_{10}}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}\\\
\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}\\\
-\frac{z_{10}\epsilon}{6(1+z_{10})}&\frac{\epsilon}{6z_{10}(1+z_{10})}&-\frac{z_{10}+2\epsilon}{z_{10}(1+z_{10})}\end{array}\right]\,.$
(9.50)

This Pfaffian matrix agrees with an independent IBP calculation for the basis
$e^{(I)}$ in (9.22), where $z_{10}=t/s$ is a ratio of Mandelstam variables.

In conclusion, the $\mathcal{D}$-module associated to the massless 2-loop
N-box diagram is isomorphic to the rational restriction of a GKZ system.

####  Section 9.2Logarithmic restriction: 1-loop Bhabha scattering

This example uses Algorithm (4) to compute a logarithmic restriction. Our goal
is the small-mass expansion of the 1-loop box topology for Bhabha scattering:

$\displaystyle\parbox[c]{0.0pt}{\includegraphics{Figs/1L_bhabha_box.pdf}}\quad\xrightarrow{m\to
0}\quad\parbox[c]{0.0pt}{\includegraphics{Figs/1L_massless_box.pdf}}\,.$

The momentum space integral family in $\mathrm{D}=4-2\epsilon$ dimensions is

$\displaystyle
I_{\nu_{1}\ldots\nu_{5}}(m)=\int\frac{\mathrm{d}^{\mathrm{D}}\ell}{\big{[}\ell^{2}-m^{2}]^{\nu_{1}}\big{[}(\ell+p_{1})^{2}]^{\nu_{2}}\big{[}(\ell+p_{12})^{2}-m^{2}]^{\nu_{3}}\big{[}(\ell+p_{123})^{2}]^{\nu_{4}}}\,,$
(9.51)

and the kinematic variables are

$\displaystyle p_{1}^{2}=p_{2}^{2}=p_{3}^{2}=p_{4}^{2}=m^{2},\quad
s=p_{12}^{2}\,,\quad t=p_{23}^{2}\,.$ (9.52)

IBP reduction reveals a basis of 5 MIs:

$\displaystyle\left[\begin{array}[]{c}I_{0010}(m)\\\ I_{0101}(m)\\\
I_{1010}(m)\\\ I_{0111}(m)\\\
I_{1111}(m)\end{array}\right]=G_{0}\cdot\vec{I}\,,\quad
G_{0}=(-s)^{-\epsilon}\,\texttt{Diagonal}\\!\left[(-s),\,1,\,1,\,(-s)^{-1},\,(-s)^{-2}\right]\,.$
(9.58)

The purpose of $G_{0}$ is to factor out the mass dimension of the integrals.
Using IBPs, the MIs $\vec{I}$ are found to obey a Pfaffian system

$\displaystyle\partial_{i}\bullet\vec{I}=P_{i}\cdot\vec{I}\quad\text{with}\quad
z_{1}=\frac{m^{2}}{-s}\quad\text{and}\quad z_{2}=\frac{t}{s}\,.$ (9.59)

Using the method of Section (8.2.4), we now compute the logarithmic
restriction of solutions onto $Y^{\prime}=\\{z_{1}=0\\},$ i.e. for $m\sim 0$.

##### 9.2.1 Normal form

The matrices $P_{1,2}$ in (9.59) are not in normal form w.r.t. $Y^{\prime}$.
In particular, $P_{1}$ has a second-order pole there, and $P_{2}$ a first-
order pole. This is cured via Moser reduction:

$\displaystyle\vec{I}\to G_{1}\cdot\vec{I}\quad\text{with}\quad
G_{1}=\texttt{Diagonal}\\!\left[z_{1},\,1,\,1,\,1,\,1\right]\,.$ (9.60)

To simplify notation, we

$\displaystyle\text{replace}\ G_{1}[P_{i}]\ \to\ P_{i}\ \text{in this
example}.$ (9.61)

The gauge transformed Pfaffian matrices now have the desired pole structure:

$\displaystyle\begin{array}[]{lllllll}P_{1}&=&z_{1}^{-1}&\times&P_{1,0}(z_{2})&+&\mathcal{O}(z_{1}^{0})\\\\[2.0pt]
P_{2}&=&z_{1}^{0}&\times&P_{2,0}(z_{2})&+&\mathcal{O}(z_{1})\,.\end{array}$
(9.64)

In addition, the spectrum of the residue matrix $P_{1,-1}$ is non-resonant,
meaning that the system is now in normal form.

From

$\displaystyle\texttt{JordanDecomposition}\\!\left[P_{1,-1}\right]=\left[\begin{array}[]{ccccc}0&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}\\\
\scalebox{0.7}{$\cdot$}&0&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}\\\
\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&0&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}\\\
\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&-\epsilon&1\\\
\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&-\epsilon\end{array}\right]$
(9.70)

we read off the data

$\displaystyle\begin{array}[]{ccl}\\{\lambda_{1},\,\lambda_{2}\\}&=&\\{0,\,-\epsilon\\}\\\
\\{\Lambda_{1},\,\Lambda_{2}\\}&=&\\{3,\,2\\}\\\
\\{M_{\lambda_{1}},\,M_{\lambda_{2}}\\}&=&\\{0,\,1\\}\,,\end{array}$ (9.74)

in the notation of step 1. of Algorithm (4). The ansatz for the logarithmic
expansion (8.60) thus takes the form

$\displaystyle\vec{I}\quad$ $\displaystyle=\quad
z_{1}^{\lambda_{1}}\sum_{n=0}^{\infty}\vec{I}\,^{(\lambda_{1},n)}(z_{2})z_{1}^{n}$
(9.75) $\displaystyle+\quad
z_{1}^{\lambda_{2}}\sum_{n=0}^{\infty}\big{[}\vec{I}\,^{(\lambda_{2},n,0)}(z_{2})\hskip
5.69046pt+\hskip
5.69046pt\vec{I}\,^{(\lambda_{2},n,1)}(z_{2})\times\log(z_{1})\big{]}z_{1}^{n}\,.$

The logarithm appears due to the "$1$" in the $\lambda_{2}$-eigenvalue block
of (9.70). Recall that the vector $\vec{I}\,^{(\lambda_{2},0,1)}$ is
immediately obtained from $\vec{I}\,^{(\lambda_{2},0,0)}$ via the recursion
relation (A). Higher-order terms in $z_{1}^{n}$ are obtained from the
recursions (B) and (C). If what follows, we determine
$\vec{I}\,^{(\lambda_{1},0)}$ from the solution to a Pfaffian system of rank
$\Lambda_{1}=3$, and $\vec{I}\,^{(\lambda_{2},0,0)}$ from the solution to a
system of rank $\Lambda_{2}=2$. Both systems are univariate.

##### 9.2.2 Eigenvalue $\lambda_{1}=0$

There is no logarithm for the eigenvalue $\lambda_{1}$, so this case can in
fact be treated with the holomorphic restriction of Section (8.2.3).

###### Rank jump

Recall the constraint from the second equation of (8.40):

$\displaystyle\mathsf{R}^{(\lambda_{1})}\cdot\vec{I}\,^{(\lambda_{1},0)}$
$\displaystyle=\texttt{RowReduce}\\!\left[P_{1,-1}\right]\cdot\vec{I}\,^{(\lambda_{1},0)}$
(9.76)
$\displaystyle=\left[\begin{array}[]{ccccc}1&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}\\\
\scalebox{0.7}{$\cdot$}&1&\scalebox{0.7}{$\cdot$}&\frac{z_{2}\epsilon}{1-2\epsilon}&\scalebox{0.7}{$\cdot$}\end{array}\right]\cdot\vec{I}\,^{(\lambda_{1},0)}$
(9.79) $\displaystyle=0\,.$ (9.80)

Two constraints are here imposed on the entries of the $5$-dimensional vector
$\vec{I}\,^{(\lambda_{1},0)}$. This is the reason why the rank drops to
$5-2=3$.

The entries of $\vec{I}\,^{(\lambda_{1},0)}$ apparently correspond to the MIs
(9.58) with the limit $z_{1}\to 0$ taken at the _integrand_ level222 This is
called the _hard region_ in the method of regions [159, 160]. , i.e. setting
$m^{2}=0$ inside the propagators of (9.51). The interpretation of the two
constraints from (9.79) is now clear:

  1. 1.

The first constraint $I^{(\lambda_{1},0)}_{1}=0$ means that the massless
tadpole $I_{0010}(m=0)$ should vanish, which is indeed true in DR.

  2. 2.

The second constraint
$I^{(\lambda_{1},0)}_{2}+\frac{z_{2}\epsilon}{1-2\epsilon}I^{(\lambda_{1},0)}_{4}=0$
is an IBP relation between the massless $t$-channel bubble $I_{0101}(m=0)$ and
the massless triangle $I_{0111}(m=0)$. We have verified this relation using
LiteRed.

###### Basis

Following step 2. of Algorithm (3), a basis for the restricted Pfaffian system
can be found by simply eliminating linearly dependent entries of
$\vec{I}\,^{(\lambda_{1},0)}$. We should eliminate two integrals according to
the above discussion, so we have333 We choose the bubble instead of the
triangle for the first basis element since it is preferred to have basis
elements with the fewest number of distinct propagators.

$\displaystyle\vec{J}\,^{(\lambda_{1})}=\left[\begin{array}[]{c}I^{(\lambda_{1},0)}_{2}\\\
I^{(\lambda_{1},0)}_{3}\\\
I^{(\lambda_{1},0)}_{5}\end{array}\right]=(-s)^{\epsilon}\left[\begin{array}[]{c}I_{0101}(0)\\\
I_{1010}(0)\\\ (-s)^{2}I_{1111}(0)\end{array}\right]\,.$ (9.87)

The factors of $(-s)$ stem from the matrix $G_{0}$ in (9.58). The
$\mathsf{B}$-matrix (8.27) is therefore

$\displaystyle\mathsf{B}^{(\lambda_{1})}=\left[\begin{array}[]{ccccc}\scalebox{0.7}{$\cdot$}&1&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}\\\
\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&1&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}\\\
\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&1\end{array}\right]\quad\text{such
that}\quad\vec{J}\,^{(\lambda_{1})}=\mathsf{B}^{(\lambda_{1})}\cdot\vec{I}\,^{(\lambda_{1},0)}\,.$
(9.91)

###### Restriction

The $\mathsf{R}$\- and $\mathsf{B}$-matrices from (9.79) and (9.91) combine to

$\displaystyle\mathsf{M}^{(\lambda_{1})}=\left[\begin{array}[]{c}\mathsf{B}^{(\lambda_{1})}\\\
\hline\cr\mathsf{R}^{(\lambda_{1})}\end{array}\right]\quad\text{such
that}\quad\mathsf{M}^{(\lambda_{1})}\cdot\vec{I}\,^{(\lambda_{1},0)}=\left[\begin{array}[]{c}\vec{J}\,^{(\lambda_{1})}\\\
\hline\cr 0\\\ 0\end{array}\right]\,.$ (9.97)

The $3\times 3$ Pfaffian matrix $Q_{2}^{(\lambda_{1})}$ associated to
$\vec{J}\,^{(\lambda_{1})}$ is now computed with (8.35). Namely, by gauge
transforming $P_{2,0}$ from (9.64) with
$\left(\mathsf{M}^{(\lambda_{1})}\right)^{-1},$ we get

$\displaystyle\left(\partial_{2}\bullet\mathsf{M}^{(\lambda_{1})}+\mathsf{M}^{(\lambda_{1})}\cdot
P_{2,0}\right)\cdot\left(\mathsf{M}^{(\lambda_{1})}\right)^{-1}=\left[\begin{array}[]{c|c}Q_{2}^{(\lambda_{1})}&\star\\\
\hline\cr\mathbf{0}_{2\times 3}&\star\end{array}\right]\,.$ (9.100)

Explicitly, we have derived the Pfaffian system

$\displaystyle\partial_{2}\bullet\vec{J}\,^{(\lambda_{1})}=Q_{2}^{(\lambda_{1})}\cdot\vec{J}\,^{(\lambda_{1})}\quad\text{with}\quad
Q_{2}^{(\lambda_{1})}=\left[\begin{array}[]{ccc}-\frac{\epsilon}{z_{2}}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}\\\
\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}\\\
\frac{2(1-2\epsilon)}{z_{2}^{2}(z_{2}+1)}&\frac{2(2\epsilon-1)}{z_{2}(z_{2}+1)}&\frac{-z_{2}-\epsilon-1}{z_{2}(z_{2}+1)}\end{array}\right]\,.$
(9.104)

This system agrees with an independent IBP calculation performed directly on
the basis of massless integrals (9.87). See [146] for a pedagogical exposition
on how to solve these equations by passing to a canonical basis.

##### 9.2.3 Eigenvalue $\lambda_{2}=-\epsilon$

The Pfaffian system associated to the eigenvalue $\lambda_{2}=-\epsilon$ turns
out to be even simpler than the one for $\lambda_{1}=0$. Let’s build it.

###### Rank jump

The rank jump constraint from (8.65) given $M_{\lambda_{2}}=1$ reads

$\displaystyle\mathsf{R}^{(\lambda_{2})}\cdot\vec{I}\,^{(\lambda_{2},0,0)}$
$\displaystyle=\texttt{RowReduce}\\!\left[(P_{1,-1}-\lambda_{2}\mathbf{1})^{2}\right]\cdot\vec{I}\,^{(\lambda_{2},0,0)}$
(9.105)
$\displaystyle=\left[\begin{array}[]{ccccc}1&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\frac{z_{2}\epsilon}{1-\epsilon}&\frac{z_{2}\epsilon}{2(1-\epsilon)}\\\
\scalebox{0.7}{$\cdot$}&1&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}\\\
\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&1&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}\end{array}\right]\cdot\vec{I}\,^{(\lambda_{2},0,0)}$
(9.109) $\displaystyle=0\,.$ (9.110)

There are 3 equations imposed on the entries of the $5$-dimensional vector
$\vec{I}\,^{(\lambda_{2},0,0)}$, wherefore the rank drops to $2$.

In this case, it is less clear what the relations (9.109) mean! The present
author does not know what the momentum space representations looks like for
other eigenvalues besides $\lambda=0$. This is worth clarifying, since it
would help with the determination of boundary constants via direct
integration.

###### Basis

The constraints (9.109) fix
$I_{2}^{(\lambda_{2},0,0)}=I_{3}^{(\lambda_{2},0,0)}=0,$ and
$I_{5}^{(\lambda_{2},0,0)}$ becomes a linear combination of
$I_{1}^{(\lambda_{2},0,0)}$ and $I_{4}^{(\lambda_{2},0,0)}$. So

$\displaystyle\vec{J}\,^{(\lambda_{2})}=\left[\begin{array}[]{c}I^{(\lambda_{2},0,0)}_{1}\\\
I^{(\lambda_{2},0,0)}_{4}\end{array}\right]$ (9.113)

is a linearly independent basis. The $\mathsf{B}$-matrix is therefore given by

$\displaystyle\mathsf{B}^{(\lambda_{2})}=\left[\begin{array}[]{ccccc}1&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}\\\
\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&1&\scalebox{0.7}{$\cdot$}\end{array}\right]\quad\text{such
that}\quad\vec{J}\,^{(\lambda_{2})}=\mathsf{B}^{(\lambda_{2})}\cdot\vec{I}\,^{(\lambda_{2},0,0)}\,.$
(9.116)

###### Restriction

Stacking the $\mathsf{R}$\- and $\mathsf{B}$-matrices from (9.109) and (9.116)
on top of each other, we get

$\displaystyle\mathsf{M}^{(\lambda_{2})}=\left[\begin{array}[]{c}\mathsf{B}^{(\lambda_{2})}\\\
\hline\cr\mathsf{R}^{(\lambda_{2})}\end{array}\right]\quad\text{such
that}\quad\mathsf{M}^{(\lambda_{2})}\cdot\vec{I}\,^{(\lambda_{2},0,0)}=\left[\begin{array}[]{c}\vec{J}\,^{(\lambda_{2})}\\\
\hline\cr 0\\\ 0\\\ 0\end{array}\right]\,.$ (9.123)

The Pfaffian system for $\vec{J}\,^{(\lambda_{2})}$ now follows from (8.35):

$\displaystyle\left(\partial_{2}\bullet\mathsf{M}^{(\lambda_{2})}+\mathsf{M}^{(\lambda_{2})}\cdot
P_{2,0}\right)\cdot\left(\mathsf{M}^{(\lambda_{2})}\right)^{-1}=\left[\begin{array}[]{c|c}Q_{2}^{(\lambda_{2})}&\star\\\
\hline\cr\mathbf{0}_{3\times 2}&\star\end{array}\right]\,.$ (9.126)

Explicitly, we have that

$\displaystyle\partial_{2}\bullet\vec{J}\,^{(\lambda_{2})}=Q_{2}^{(\lambda_{2})}\cdot\vec{J}\,^{(\lambda_{2})}\quad\text{with}\quad
Q_{2}^{(\lambda_{2})}=\left[\begin{array}[]{cc}0&0\\\
\frac{\epsilon-1}{z_{2}^{2}}&\frac{-1}{z_{2}}\end{array}\right]\,.$ (9.129)

The general solution is swiftly found to be

$\displaystyle\vec{J}\,^{(\lambda_{2})}=\left[\begin{array}[]{c}c_{1}(\epsilon)\\\
\frac{c_{1}(\epsilon)(1-\epsilon)\log(z_{2})}{z_{2}}+\frac{c_{2}(\epsilon)}{z_{2}}\end{array}\right]\,,$
(9.132)

where $c_{1,2}(\epsilon)$ are two boundary constants.

##### 9.2.4 Result for the logarithmic series

The solution vectors in (9.75) at $0$th order in $z_{1}$ take the form

$\displaystyle\vec{I}\,^{(\lambda_{1},0)}$
$\displaystyle=\left(\mathsf{M}^{(\lambda_{1})}\right)^{-1}\cdot\left[\begin{array}[]{c}\vec{J}\,^{(\lambda_{1})}\\\
\hline\cr 0\\\ 0\end{array}\right]$ (9.136)
$\displaystyle\vec{I}\,^{(\lambda_{2},0,0)}$
$\displaystyle=\left(\mathsf{M}^{(\lambda_{2})}\right)^{-1}\cdot\left[\begin{array}[]{c}\vec{J}\,^{(\lambda_{2})}\\\
\hline\cr 0\\\ 0\\\ 0\end{array}\right]$ (9.141)
$\displaystyle\vec{I}\,^{(\lambda_{2},0,1)}$
$\displaystyle=(P_{1,-1}-\lambda_{2}\mathbf{1})\cdot\vec{I}\,^{(\lambda_{2},0,0)}\,,$
(9.142)

where the last equation comes from the recursion (A). The vector
$\vec{J}\,^{(\lambda_{1})}$ comes from solving the massless 1-loop box system
(9.104), and $\vec{J}\,^{(\lambda_{2})}$ was given in (9.132). We do not
include higher-order terms in $z_{1}^{n}$ (coming from the recursions (B) and
(C)) in order to show a compact, final result.

Inserting the solution vectors above into (9.75) we get444We are neglecting
the DR scale $\mu$ in factors such as $(-s/\mu^{2})^{-\epsilon}$.

$\displaystyle{\scriptstyle\left[\begin{array}[]{c}I_{0010}(m)\\\\[3.0pt]
I_{0101}(m)\\\\[3.0pt] I_{1010}(m)\\\\[3.0pt] I_{0111}(m)\\\\[3.0pt]
I_{1111}(m)\end{array}\right]\overset{m\to
0}{=}\left[\begin{array}[]{c}0\\\\[3.0pt]
(-s)^{-\epsilon}I_{0101}(0)\\\\[3.0pt] (-s)^{-\epsilon}I_{1010}(0)\\\\[3.0pt]
\frac{(-s)^{1-\epsilon}(1-2\epsilon)}{t\epsilon}I_{0101}(0)\\\\[3.0pt]
(-s)^{-2-\epsilon}I_{1111}(0)\end{array}\right]+\left[\begin{array}[]{c}\frac{c_{1}}{\epsilon-1}(m^{2})^{1-\epsilon}\\\\[3.0pt]
0\\\\[3.0pt] 0\\\\[3.0pt]
\frac{-(m^{2})^{-\epsilon}}{t}\big{[}c_{1}\log(-t/m^{2})+c_{2}\big{]}\\\\[3.0pt]
\frac{-2(m^{2})^{-\epsilon}}{st\epsilon}\big{[}c_{1}(\epsilon\log(-t/m^{2})-1)+\epsilon
c_{2}\big{]}\end{array}\right]}$ (9.158)

The first and second terms on the LHS come from $\lambda_{1}$ and
$\lambda_{2}$ respectively. Here we reinstated all prefactors, wrote
$z_{1}=m^{2}/(-s)$ and $z_{2}=t/s$, and lastly rescaled
$c_{1}\to\frac{c_{1}}{\epsilon-1}$ to simplify the expression. Note that only
the triangle and the box integrals receive logarithmic corrections. The
bubbles would naturally receive power corrections in $z_{1}$, which we are
neglecting.

The two boundary constants are determined as follows.

  * •

$c_{1}$ is fixed by comparing with the exact expression
$I_{0010}(m)=-\Gamma(1-\epsilon)(m^{2})^{1-\epsilon}$ for the tadpole on the
LHS. This gives $c_{1}=(1-\epsilon)\Gamma(\epsilon-1).$

  * •

$c_{2}=c_{2}(\epsilon)$ can be fixed order-by-order in $\epsilon$ by
numerically matching with an integral on the LHS, say the 5th element
$I_{1111}(m)$. On the LHS we insert the ansatz
$c_{2}=\sum_{i=-2}^{2}c^{(i)}_{2}\epsilon^{i},$ which starts at $i=-2$ because
the massless box $I_{1111}(0)$ has a double pole in $\epsilon$. On the RHS we
numerically evaluate $I_{1111}(m)$ at the point $(s,t,m^{2})=(-1,-9,10^{-3})$
with AMFlow [34]. By comparing the LHS and RHS we get

$\displaystyle\Big{[}c^{(-2)}_{2},c^{(-1)}_{2},c^{(0)}_{2},c^{(1)}_{2},c^{(2)}_{2}\Big{]}=\big{[}1,-0.588997,-7.12155,-4.69031,-20.4068\big{]}\,.$

AMFlow is easily able to compute $\mathcal{O}(10^{2})$ digits for this
topology. Analytical expressions for these constants could therefore, if
needed, be found with PSLQ [77].

Evaluating at other phase space points with small $m^{2}$, we find that the
LHS of (9.158) is a fair approximation to the $\epsilon$-expansions of
$I_{0111}(m)$ and $I_{1111}(m)$ up to $\mathcal{O}(\epsilon^{3})$. For
instance, at $(s,t,m^{2})=(-1,-8,10^{-3})$ we find agreement within $0.01\%$
accuracy for the triangle, and within $0.001\%$ accuracy for the box. The non-
leading terms in $\epsilon$ for the bubbles $I_{0101}(m)$ and $I_{1010}(m)$
give errors between $5-20\%$ because we have neglected power corrections.

####  Section 9.3Logarithmic restriction: 2-loop Bhabha scattering

This is another example of logarithmic restriction using Algorithm (4). We
provide partial results on the small-mass logarithmic expansion of the 2-loop
Bhabha double-box topology:

$\displaystyle\parbox[c]{0.0pt}{\includegraphics{Figs/2L_bhabha_box.pdf}}\quad\xrightarrow{m\to
0}\quad\parbox[c]{0.0pt}{\includegraphics{Figs/2L_massless_box.pdf}}\,.$
(9.159)

The momentum space integral family in $\mathrm{D}=4-2\epsilon$ dimensions
reads

$\displaystyle
I_{\nu_{1}\ldots\nu_{9}}(m)=\int\frac{\mathrm{d}^{\mathrm{D}}\ell_{1}\mathrm{d}^{\mathrm{D}}\ell_{2}}{D_{1}^{\nu_{1}}\cdots
D_{9}^{\nu_{9}}}\,,$ (9.160)

with

$\displaystyle\begin{array}[]{llllllllll}&D_{1}&=&\ell_{1}^{2}-m^{2}&D_{2}&=&(\ell_{1}+p_{1}+p_{2})^{2}-m^{2}&D_{3}&=&\ell_{2}^{2}-m^{2}\\\
&D_{4}&=&(\ell_{2}+p_{1}+p_{2})^{2}-m^{2}&D_{5}&=&(\ell_{1}+p_{1})^{2}&D_{6}&=&(\ell_{1}-\ell_{2})^{2}\\\
&D_{7}&=&(\ell_{2}-p_{3})^{2}&D_{8}&=&(\ell_{2}+p_{1})^{2}&D_{9}&=&(\ell_{1}-p_{3})^{2}\,,\end{array}$
(9.164)

where $D_{8}$ and $D_{9}$ denote ISPs. Solutions for this integral family were
found in [170] (see also [171]).

Inspired by [170], we choose the following basis for the unrestricted Pfaffian
system:

$\displaystyle\begin{array}[]{lllllll}G_{0}\cdot\vec{I}=\big{[}&I_{202000000},&I_{000022100},&I_{001022000},&I_{122000000},&I_{012002000},&I_{022001000},\\\
&I_{001112000},&I_{012011000},&I_{100011200},&I_{121200000},&I_{011011100},&I_{011012100},\\\
&I_{021011100},&I_{022011100},&I_{101111000},&I_{101211000},&I_{102011100},&I_{101011200},\\\
&I_{110011200},&I_{111011100},&I_{1110111-10},&I_{111111100},&I_{1111111-10}&\big{]}^{T}\,.\end{array}$
(9.169)

The $23\times 23$ diagonal matrix $G_{0}$ factors out the mass dimension of
each integral via powers of $(-s)$. IBP reduction with kira [10] then provides
a Pfaffian system $\partial_{i}\bullet\vec{I}=P_{i}\cdot\vec{I}$ in terms of
the kinematic variables

$\displaystyle z_{1}=\frac{m^{2}}{-s}\quad\text{and}\quad
z_{2}=\frac{t}{s}\,.$ (9.170)

Our goal is to restrict onto the surface $Y^{\prime}=\\{z_{1}=0\\}$.

##### 9.3.1 Normal form

The Pfaffian system associated to the basis (9.169) is not in normal form. The
first issue to resolve is the higher-order singularities on $Y^{\prime}$:

$\displaystyle\begin{array}[]{lllllll}P_{1}&=&z_{1}^{-3}&\times&P_{1,-3}(z_{2})&+&\mathcal{O}(z_{1}^{-2})\\\\[2.0pt]
P_{2}&=&z_{1}^{-2}&\times&P_{2,-2}(z_{2})&+&\mathcal{O}(z_{1}^{-1})\,.\end{array}$
(9.173)

Moser reduction suggests to gauge transform with

$\displaystyle
G_{1}=\texttt{Diagonal}\big{[}\underbrace{z_{1},\ldots,z_{1}}_{8\
\text{times}},\,1,\ldots,\,1\big{]}\quad\text{and}\quad
G_{2}=\texttt{Diagonal}\big{[}\underbrace{z_{1},\ldots,z_{1}}_{12\
\text{times}},\,1,\ldots,\,1\big{]}\,.$

Introducing the shorthand notation $G_{i\ldots j}=G_{i}\cdot\ldots\cdot
G_{j},$ then the new pole structure becomes

$\displaystyle\begin{array}[]{lllllll}G_{21}[P_{1}]&=&z_{1}^{-1}&\times&G_{21}[P_{1}]_{-1}&+&\mathcal{O}(z_{1}^{0})\\\\[2.0pt]
G_{21}[P_{2}]&=&z_{1}^{-1}&\times&G_{21}[P_{2}]_{-1}&+&\mathcal{O}(z_{1}^{0})\,.\end{array}$
(9.176)

So we have cured the higher-order poles in $P_{1}$. But $P_{2}$ is still not
finite on $Y^{\prime}$ because the residue matrix $G_{21}[P_{1}]_{-1}$ has a
resonant spectrum555 See [83, Appendix B] for a proof of why non-resonance of
the residue matrix $P_{1,-1}$ implies finiteness of the remaining Pfaffian
matrices $\\{P_{2},\ldots P_{N}\\}$ on $\\{z_{1}=0\\}$. . In [83, Appendix B],
it is shown how to construct a sequence of gauge transformations (involving
the Jordan decomposition of the residue matrix) to cure the spectrum666 The
Mathematica scripts in [169] contain the gauge transformations for this
particular example. . It takes a matter of seconds on a laptop to compute
these gauge transformations, even without fixing $z$-variables to generic
numbers (recall that the Jordan decomposition is independent of $z$). Writing
$G_{3}$ for the product of said gauge matrices, we get

$\displaystyle G_{321}[P_{2}]=G_{321}[P_{2}]_{0}+\mathcal{O}(z_{1})\,,$
(9.177)

wherefore the system is now in normal form. For brevity, we

$\displaystyle\text{replace}\ G_{321}[P_{i}]\ \to\ P_{i}\ \text{in the rest of
this example}.$ (9.178)

$P_{1,-1}$ denotes the residue matrix in this newfound basis. Its Jordan form
enjoys an eigenvalue decomposition

$\displaystyle\texttt{JordanDecomposition}\\!\left[P_{1,-1}\right]=\mathsf{J}^{(\lambda_{1})}\oplus\mathsf{J}^{(\lambda_{2})}\oplus\mathsf{J}^{(\lambda_{3})}\oplus\mathsf{J}^{(\lambda_{4})}$
(9.179)

in terms of block matrices

$\displaystyle\mathsf{J}^{(\lambda_{1})}$ $\displaystyle=\mathbf{0}_{8\times
8}\,,\quad\mathsf{J}^{(\lambda_{2})}=\left[\begin{array}[]{c}\lambda_{2}\end{array}\right]$
(9.181) $\displaystyle\mathsf{J}^{(\lambda_{3})}$
$\displaystyle=\left[\begin{array}[]{ccccccc}\lambda_{3}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}\\\
\scalebox{0.7}{$\cdot$}&\lambda_{3}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}\\\
\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\lambda_{3}&1&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}\\\
\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\lambda_{3}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}\\\
\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\lambda_{3}&1&\scalebox{0.7}{$\cdot$}\\\
\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\lambda_{3}&\scalebox{0.7}{$\cdot$}\\\
\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\lambda_{3}\end{array}\right]$
(9.189) $\displaystyle\mathsf{J}^{(\lambda_{4})}$
$\displaystyle=\left[\begin{array}[]{ccccccc}\lambda_{4}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}\\\
\scalebox{0.7}{$\cdot$}&\lambda_{4}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}\\\
\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\lambda_{4}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}\\\
\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\lambda_{4}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}\\\
\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\lambda_{4}&1&\scalebox{0.7}{$\cdot$}\\\
\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\lambda_{4}&1\\\
\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\lambda_{4}\end{array}\right]\,,$
(9.197)

where the eigenvalues are

$\displaystyle\\{\lambda_{1},\,\lambda_{2},\,\lambda_{3},\,\lambda_{4}\\}=\left\\{0,\,2(\mathrm{D}-5)\,,\frac{\mathrm{D}-4}{2},\,\mathrm{D}-4\right\\}\,.$
(9.198)

Their associated multiplicities

$\displaystyle\\{\Lambda_{1},\,\Lambda_{2},\,\Lambda_{3},\,\Lambda_{4}\\}=\\{8,\,1,\,7,\,7\\}$
(9.199)

count the sizes of the four different Pfaffian systems that will be obtained
via logarithmic restriction. The Jordan blocks above specify the maximum
logarithm powers

$\displaystyle\\{M_{\lambda_{1}},\,M_{\lambda_{2}},\,M_{\lambda_{3}},\,M_{\lambda_{3}}\\}=\\{0,\,0,\,1,\,2\\}$
(9.200)

for the series ansatz (8.60). The vector of MIs $\vec{I}$ can therefore be
represented as

$\displaystyle\begin{array}[]{llllllll}G_{321}\cdot\vec{I}(z)&=&&z_{1}^{\lambda_{1}}\vec{I}\,^{(\lambda_{1},0)}(z)&+&z_{1}^{\lambda_{2}}\vec{I}\,^{(\lambda_{2},0)}(z)&&\\\\[6.0pt]
&+&&z_{1}^{\lambda_{3}}\big{[}\vec{I}\,^{(\lambda_{3},0)}(z)&+&\vec{I}\,^{(\lambda_{3},1)}(z)\log(z_{1})\big{]}&&\\\\[6.0pt]
&+&&z_{1}^{\lambda_{4}}\big{[}\vec{I}\,^{(\lambda_{4},0)}(z)&+&\vec{I}\,^{(\lambda_{4},1)}(z)\log(z_{1})&+&\vec{I}\,^{(\lambda_{4},2)}(z)\log^{2}(z_{1})\big{]}\,,\end{array}$
(9.204)

where $z=(z_{1},z_{2})$. The logarithm-free terms above all have holomorphic
expansions in $z_{1}$ of the form

$\displaystyle\vec{I}\,^{(\lambda,0)}(z)=\sum_{n=0}^{\infty}\vec{I}\,^{(\lambda,n,0)}(z_{2})\,z_{1}^{n}\,.$
(9.205)

Recall that all higher-order terms $\vec{I}\,^{(\lambda,n,m)}$ in
$z_{1}^{n}\,\log^{m}(z_{1})$ are uniquely specified given
$\vec{I}\,^{(\lambda,0,0)}$ as input via the recursion relations (A), (B) and
(C).

##### 9.3.2 Restriction

We proceed to compute the $\Lambda\times\Lambda$ Pfaffian matrices
$Q_{2}^{(\lambda)}(z_{2})$ associated to each $\vec{I}\,^{(\lambda,0,0)}$. The
full expression (9.204) has not yet been mounted because the boundary
constants are missing, so this step is left for future study.

###### Eigenvalue $\lambda_{1}$

The $\mathsf{R}$-matrix for $\lambda_{1}=0$ splits into a $0$-block and an
identity matrix block:

$\displaystyle\mathsf{R}^{(\lambda_{1})}=\texttt{RowReduce}\\!\left[P_{1,-1}\right]=\left[\begin{array}[]{c|c}\mathbf{0}_{15\times
8}&\mathbf{1}_{15}\end{array}\right]\,.$ (9.207)

Since $\mathsf{R}^{(\lambda_{1})}\cdot\vec{I}\,^{(\lambda_{1},0,0)}=0$, this
means that $I^{(\lambda_{1},0,0)}_{9}=\ldots=I^{(\lambda_{1},0,0)}_{23}=0.$
The $8$ remaining entries of $\vec{I}\,^{(\lambda_{1},0,0)}$ thus make up a
basis for the restricted Pfaffian system:

$\displaystyle\vec{J}\,^{(\lambda_{1})}=\left[\begin{array}[]{c}I^{(\lambda_{1},0,0)}_{1}\\\
\vdots\\\ I^{(\lambda_{1},0,0)}_{8}\end{array}\right]\,.$ (9.211)

This corresponds to the $\mathsf{B}$-matrix

$\displaystyle\mathsf{B}^{(\lambda_{1})}=\left[\begin{array}[]{c|c}\reflectbox{$\mathbf{1}$}_{8}&\mathbf{0}_{8\times
15}\end{array}\right]\quad\text{such
that}\quad\vec{J}\,^{(\lambda_{1})}=\mathsf{B}^{(\lambda_{1})}\cdot\vec{I}\,^{(\lambda_{1},0,0)}\,,$
(9.213)

where $\reflectbox{$\mathbf{1}$}_{8}$ represents the "reversed identity
matrix", i.e. the $8\times 8$ matrix with $1$’s on the anti-diagonal. Setting

$\displaystyle\mathsf{M}^{(\lambda_{1})}=\left[\begin{array}[]{c}\mathsf{B}^{(\lambda_{1})}\\\
\hline\cr\mathsf{R}^{(\lambda_{1})}\end{array}\right]\,,$ (9.216)

then the Pfaffian matrix associated to $\vec{J}\,^{(\lambda_{1})}$ is

$\displaystyle\left(\partial_{2}\bullet\mathsf{M}^{(\lambda_{1})}+\mathsf{M}^{(\lambda_{1})}\cdot
P_{2,0}\right)\cdot\left(\mathsf{M}^{(\lambda_{1})}\right)^{-1}=\left[\begin{array}[]{c|c}Q_{2}^{(\lambda_{1})}&\star\\\
\hline\cr\mathbf{0}_{15\times 8}&\star\end{array}\right]\,.$ (9.219)

We find that

$\displaystyle\partial_{2}\bullet\vec{J}\,^{(\lambda_{1})}=Q_{2}^{(\lambda_{1})}\cdot\vec{J}\,^{(\lambda_{1})}$
(9.220)

is indeed a Pfaffian system for the massless double-box topology, as verified
by an independent IBP computation. Each entry of $\vec{J}\,^{(\lambda_{1})}$
is a linear combination777 It is a linear combination because the gauge
transformation matrix $G_{3}$ from (9.177) is non-diagonal.  of integrals in
the family (9.160) with $m^{2}=0$ at the integrand level. The matrix
$Q_{2}^{(\lambda_{1})}$ has the form
$\frac{Q_{a}(\epsilon)}{z_{1}^{2}}+\frac{Q_{b}(\epsilon)}{z_{1}}+\frac{Q_{c}(\epsilon)}{z_{1}+1}$
where $Q_{a,b,c}(\epsilon)$ contain rational functions in $\epsilon$ with
maximal degree $\epsilon^{2}$. Using the software canonica [172], a canonical
form for $Q_{2}^{(\lambda_{1})}$ is found within seconds.

###### Eigenvalue $\lambda_{2}$

The system for the eigenvalue $\lambda_{2}=2(\mathrm{D}-5)=-2-4\epsilon$ is
peculiar. The rank jump constraint

$\displaystyle\mathsf{R}^{(\lambda_{2})}\cdot\vec{I}\,^{(\lambda_{2},0,0)}=\texttt{RowReduce}\\!\left[P_{1,-1}-\lambda_{2}\mathbf{1}\right]\cdot\vec{I}\,^{(\lambda_{2},0,0)}=0$
(9.221)

sets every entry of $\vec{I}\,^{(\lambda_{2},0,0)}$ to zero except for the
$15$th entry $J^{(\lambda_{2})}=I\,^{(\lambda_{2},0,0)}_{15}$. The Pfaffian
system for this single, non-trivial function works out to

$\displaystyle\partial_{2}\bullet
J^{(\lambda_{2})}=-\frac{J^{(\lambda_{2})}}{z_{2}}\quad\implies\quad
J^{(\lambda_{2})}=\frac{c^{(\lambda_{2})}(\epsilon)}{z_{2}}\,.$ (9.222)

###### Eigenvalue $\lambda_{3}$

The system associated to eigenvalue
$\lambda_{3}=\frac{\mathrm{D}-4}{2}=-\epsilon$ has rank $\Lambda_{3}=7$.
Working out the $\mathsf{R}$\- and $\mathsf{B}$-matrices as usual, one finds a
basis

$\displaystyle\vec{J}\,^{(l_{3})}=\left[\begin{array}[]{c}I^{(\lambda_{3},0,0)}_{16}\\\\[2.0pt]
\vdots\\\\[2.0pt] I^{(\lambda_{3},0,0)}_{22}\end{array}\right]\,.$ (9.226)

The associated Pfaffian matrix $Q_{2}^{(\lambda_{3})}$ has the form
$\frac{Q_{a}(\epsilon)}{z_{2}^{2}}+\frac{Q_{b}(\epsilon)}{z_{2}}+\frac{Q_{c}(\epsilon)}{(z_{2}+1)^{2}}+\frac{Q_{d}(\epsilon)}{z_{2}-1},$
where $Q_{a,b,c,d}(\epsilon)$ contain rational functions in $\epsilon$ of
maximal degree $\epsilon^{3}$. This system is swiftly $\epsilon$-factorized
with canonica.

###### Eigenvalue $\lambda_{4}$

The last eigenvalue $\lambda_{4}=\mathrm{D}-4=-2\epsilon$ also has rank
$\Lambda_{4}=7$. Its associated rank jump constraint
$\mathsf{R}^{(\lambda_{4})}\cdot\vec{I}\,^{(\lambda_{4},0,0)}=0$ contains
linear relations with rather complicated coefficients. This is in contrast to
the rank jump constraints for the other eigenvalues, which were all of the
form $I^{(\lambda,0,0)}_{i}=0.$ We do not know why this is the case.

Regardless, we are easily able to calculate the Pfaffian system associated to
the basis

$\displaystyle\vec{J}\,^{(\lambda_{4})}=\left[\begin{array}[]{c}I^{(\lambda_{4},0,0)}_{5}\\\\[5.0pt]
I^{(\lambda_{4},0,0)}_{17}\\\\[5.0pt] I^{(\lambda_{4},0,0)}_{19}\\\\[5.0pt]
I^{(\lambda_{4},0,0)}_{20}\\\\[5.0pt] I^{(\lambda_{4},0,0)}_{21}\\\\[5.0pt]
I^{(\lambda_{4},0,0)}_{22}\\\\[5.0pt]
I^{(\lambda_{4},0,0)}_{23}\end{array}\right]\,.$ (9.234)

The Pfaffian matrix has the singularity structure
$Q_{2}^{(\lambda_{3})}=\frac{Q_{a}(\epsilon)}{z_{2}^{2}}+\frac{Q_{b}(\epsilon)}{z_{2}}+\frac{Q_{c}(\epsilon)}{(z_{2}+1)^{2}}+\frac{Q_{d}(\epsilon)}{z_{2}-1},$
where $Q_{a,b,c,d}(\epsilon)$ have maximal degree $\epsilon^{3}$. Canonica
brings this system into canonical form with ease.

##### 9.3.3 Discussion

General solutions to the Pfaffian systems above are easily obtained in terms
of HPLs in the variable $z_{2}=t/s$. What is more, every step in the
construction of those systems can be automatized.

It is worth comparing the present method with previous studies. The seminal
calculation carried out in [170] was complicated by an unrationalizable square
root. For another 2-loop Bhabha scattering topology [171], the authors had to
deal with elliptic functions. Because our restriction protocol involves
smaller, univariate Pfaffian systems, such complications are more easily
avoided. The cost is that a smaller part of phase space is covered, since
$m^{2}$ is assumed to be small (it could e.g. be the electron mass $m_{e}\sim
0.511\ \text{MeV}$).

As mentioned, the final expression for the logarithmic series (9.204) still
lacks the determination of boundary constants. These could be fixed by
matching with the full set of two-variable MIs $\vec{I}$ from (9.169) after
imposing physically motivated constraints. For instance, one could impose
real-valuedness in the Euclidean region, or regularity at certain points in
phase space. Another way is just to numerically evaluate the two-variable
integrals $\vec{I}$ and match this with the logarithmic series. Both of these
approaches are, however, unsatisfactory from the point of view of the
restriction protocol. It would be preferable to somehow fix the boundary
constants solely within each univariate system
$\partial_{2}\bullet\vec{J}\,^{(\lambda)}=Q_{2}^{(\lambda)}\cdot\vec{J}\,^{(\lambda)},$
rather than comparing with the two-variate integrals. This is indeed possible
for the eigenvalue $\lambda=0$, since in this case we know the momentum space
representation of $\vec{J}\,^{(\lambda)}$. A speculative suggestion is to
perform an _additional_ restriction on these univariate systems, in which case
the rank jump constraints might provide relations among boundary constants.

####  Section 9.4Macaulay-level restriction: 2-loop diagonal box

This example concerns the rational restriction of the GKZ system associated to
a 2-loop massless diagonal-box diagram. We employ Algorithm (5), which is
based on the Macaulay matrix.

##### 9.4.1 Setup

The LP monomials are given by

$\displaystyle\mathcal{G}(z|x)=\ $ $\displaystyle
z_{1}x_{1}x_{4}+z_{2}x_{1}x_{5}+z_{3}x_{1}x_{6}+z_{4}x_{2}x_{4}+z_{5}x_{2}x_{5}+z_{6}x_{2}x_{6}+$
$\displaystyle
z_{7}x_{3}x_{4}+z_{8}x_{3}x_{5}+z_{9}x_{3}x_{6}+z_{10}x_{4}x_{5}+z_{11}x_{4}x_{6}+z_{12}x_{4}x_{7}+$
$\displaystyle
z_{13}x_{5}x_{7}+z_{14}x_{6}x_{7}+z_{15}x_{1}x_{3}x_{4}+z_{16}x_{1}x_{3}x_{5}+z_{17}x_{1}x_{3}x_{6}+$
(9.235) $\displaystyle
z_{18}x_{1}x_{4}x_{5}+z_{19}x_{2}x_{4}x_{6}+z_{20}x_{2}x_{4}x_{7}+z_{21}x_{2}x_{5}x_{7}+z_{22}x_{2}x_{6}x_{7}\,.$

Recall that this enters the integrand of the GFI (7.7). $\mathcal{G}(z|x)$
gives rise to to a GKZ system with $A$-matrix

$\displaystyle
A=\left[\begin{array}[]{cccccccccccccccccccccc}1&1&1&1&1&1&1&1&1&1&1&1&1&1&1&1&1&1&1&1&1&1\\\
1&1&1&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&1&1&1&1&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}\\\
\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&1&1&1&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&1&1&1&1\\\
\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&1&1&1&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&1&1&1&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}\\\
1&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&1&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&1&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&1&1&1&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&1&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&1&1&1&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}\\\
\scalebox{0.7}{$\cdot$}&1&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&1&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&1&\scalebox{0.7}{$\cdot$}&1&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&1&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&1&\scalebox{0.7}{$\cdot$}&1&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&1&\scalebox{0.7}{$\cdot$}\\\
\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&1&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&1&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&1&\scalebox{0.7}{$\cdot$}&1&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&1&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&1&\scalebox{0.7}{$\cdot$}&1&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&1\\\
\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&1&1&1&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&\scalebox{0.7}{$\cdot$}&1&1&1\\\
\end{array}\right]\,.$ (9.244)

The restriction we compute is

$\displaystyle A\quad\xrightarrow[1\leq i\,\leq\,18\ \text{and}\
19\,\leq\,j\,\leq\,22]{z_{i}=1\ \text{and}\
z_{j}=t/s}\quad\parbox[c]{0.0pt}{\includegraphics{Figs/2L_diagonal_box.pdf}}\,,$
(9.245)

where $s$ and $t$ stand for Mandelstam variables. To be clear, after a linear
shift of coordinates, this is a restriction onto the affine hyperplane

$\displaystyle
Y^{\prime}=\big{\\{}(z_{1},\ldots,z_{22})\in\mathbb{A}^{22}\,\big{|}\,z_{1}=\ldots=z_{18}=0\
\text{and}\ z_{19}=z_{20}=z_{21}=z_{22}=z\ \text{is generic}\big{\\}}\,.$

Using homogeneity (cf. Section (6.5)) we can immediately restrict $z_{i}=0$
for $i\in\\{1,2,5,7,9,11,13,18\\}$. This means that Algorithm (5) should be
applied to restrict $z_{i}=0$ for $i\in\\{3,4,6,8,10,12,14,15,16,17\\}$ and
$z_{i}=z=t/s$ for $i\in\\{19,20,21,22\\}$.

##### 9.4.2 Restriction

For generic values of the analytic regulators $\epsilon$ and $\delta$, the GKZ
system is of rank $R=115$. Computing the Euler characteristic by counting the
number of critical points of (6.60) on $Y^{\prime}$, we find a rank drop to
$R^{\prime}=7$. This number is verified by counting the number of MIs using
IBPs.

We must find a basis of the restriction module before calling Algorithm (5).
Executing [83, Algorithm 4.7] as in Example (5) with parameters
$\gamma=3,\,k=3,\,p=100000007$ and $z=1/7$, we obtain the $7$-dimensional
basis888 This took 40.56s on a machine with the following specs: Intel(R)
Core(TM) i7-10700K CPU @ 3.80GHz, 8G memory.

$\displaystyle\mathrm{RStd}=\left[\begin{array}[]{c}\partial_{21}\partial_{22}\\\
\partial_{22}^{2}\\\ \partial_{6}\\\ \partial_{17}\\\ \partial_{21}\\\
\partial_{22}\\\ 1\end{array}\right]\,.$ (9.253)

A valid Macaulay matrix
$M_{1}=\left[\begin{array}[]{c|c}M_{\mathrm{Ext}}&M_{\mathrm{RStd}}\end{array}\right]$
is then found999 It took 65.76s on the same machine to compute this Macaulay
matrix.  for degree $Q=1$, meaning that the equation (5.13) given by
$C_{\mathrm{Ext}}-C\cdot M_{\mathrm{Ext}}=0$ has a solution for the matrix
$C$. The matrix $M_{\mathrm{Ext}}$ has dimension $1305\times 2763$, so there
are $2763$ exterior monomials $\mathrm{Ext}$. The $1305$ rows are not all
independent; by row reducing $M_{\mathrm{Ext}}$ with $(\epsilon,\delta,z)$
fixed to generic numbers, we identify $912$ linearly independent equations. It
then takes less than $3$ minutes on a laptop to determine $C$ via FiniteFlow.
The $7\times 7$ Pfaffian matrix $P_{z}$ finally results from
$P_{z}=C_{\mathrm{RStd}}-C\cdot M_{\mathrm{RStd}}.$ This becomes a Pfaffian
matrix for the proper 2-loop diagonal-box topology after multiplying by the LP
prefactors from (7.7) and sending $\delta\to 0$.

As an aside, note how large the discrepancy is between the ranks $R=115$ and
$R^{\prime}=7$ for the GKZ system and the proper FI. So even though we benefit
from being able to immediately write down an annihilating ideal and a basis in
the GKZ setting, the additional complexity brought about by unspecified
$z$-variables may quickly become too much to handle.

##  Part III Numerical Integration

### Chapter 10Tropical Integration

This thesis has, so far, predominantly been an _analytic_ study of Euler and
Feynman integrals. In particular, we have focused on the DEQs obeyed by said
integrals. The ideal situation is that 1) these DEQs can be built with
reasonable computer resources, and 2) they can be solved symbolically in terms
of special functions. For sufficiently difficult problems - e.g. the
evaluation of 2- or 3-loop diagrams with several kinematic scales - either or
both of these steps become insurmountable with the presently available
analytical methods. These difficulties are in many cases unavoidable if the
accuracy of theoretical predictions is made to compete with the experimental
precision at the LHC. For instance, the slow convergence rate of the QCD
perturbative series (stemming from the coupling constant being
$\alpha_{s}(M_{Z})\sim 0.1$ near the weak scale) often requires 3-loop results
to get cross sections at $\%$-level accuracy. Furthermore, the study of Higgs
phenomenology typically requires the calculation of diagrams with many
internal masses that do not evaluate to GPLs, which complicates the
determination of canonical bases in the method of DEQs for FIs.

In cases where analytical methods fall short, _numerical_ methods may instead
offer a way forward. The past few years have indeed seen a surge in the
development of public codes tailored to the numerical evaluation of FIs. Using
series of expansions of Pfaffian systems, the codes AMFlow [34], DiffExp [173,
174] and SeaSyde [175] are able to to achieve phenomenal accuracy. They derive
Pfaffian systems by fixing parameters to numbers _before_ solving the IBPs,
which yields a significant speed boost. The pySecDec codebase [176], employing
_Monte Carlo_ integration in combination with _sector decomposition_ , has
been a staple of high-energy physics calculations for many years, and was
recently upgraded in [177]. Sophisticated manipulations of momentum space
integrals via the _Loop-tree duality_ led to efficient and stable numerical
integration schemes in [178, 179]. The authors of [180] presented a high-
precision numerical algorithm by interpreting FIs as _period integrals_.
Further promising studies were carried out in [181, 182, 183], though so far
without publicly available codes.

Part III of this thesis is a contribution to this flourishing research
program. The present chapter describes a Monte Carlo integration scheme called
_tropical integration_. It was established by Borinsky in [40] (see also
[184]), and subsequently generalized by Borinsky, Tellander and the present
author in [185]. A public code named Feyntrop was developed alongside the
aforementioned article. Feyntrop is able to swiftly numerically evaluate FIs
at high loop-order with many mass scales at ‰-level accuracy given modest
computer resources. The next chapter gives a tutorial on how to use the code,
and showcases several state-of-the-art examples. The present chapter
explicates the theoretical background.

####  Section 10.1Tropical Monte Carlo sampling, in brief

Tropical integration is a special instance of Monte Carlo (MC) integration. MC
methods deal with integrals of the form

$\displaystyle I_{f}=\int_{\Gamma}f(x)\mu\,,$ (10.1)

where $f(x)$ is a square-integrable function integrated over an
$n$-dimensional domain $\Gamma$, and $\mu$ is a probability measure that is
positive on $\Gamma$ and normalized to unity, $\int_{\Gamma}\mu=1.$ By
sampling many independent random points $\\{x^{(1)},\ldots,x^{(N)}\\}$ from
$\mu$, it follows from the central limit theorem that

$\displaystyle I_{f}\simeq I_{f}^{(N)}\quad\text{where}\quad
I_{f}^{(N)}=\frac{1}{N}\sum_{i=1}^{N}f\big{(}x^{(i)}\big{)}\,.$ (10.2)

Here $I_{f}^{(N)}$ is simply the average of the integrand $f$ evaluated at the
random points $x^{(i)}$. For sufficiently large $N$, the expected error of
this approximation is

$\displaystyle\sigma_{f}=\sqrt{\frac{I_{f^{2}}-I_{f}^{2}}{N}}\quad\text{where}\quad
I_{f^{2}}=\int_{\Gamma}f(x)^{2}\mu\,.$ (10.3)

This error can itself be estimated by

$\displaystyle\sigma_{f}\simeq\sigma_{f}^{(N)}\quad\text{where}\quad\sigma_{f}^{(N)}=\sqrt{\frac{I_{f^{2}}^{(N)}-\left[I_{f}^{(N)}\right]^{2}}{N-1}}\quad\text{and}\quad
I_{f^{2}}^{(N)}=\frac{1}{N}\sum_{i=1}^{N}f\big{(}x^{(i)}\big{)}^{2}\,.$ (10.4)

A common choice for $\mu$ is just the uniform measure $\mu=\mathrm{d}^{n}x$ on
the $n$-dimensional unit cube $\Gamma=[0,1]^{n}$. Tropical integration,
instead, makes use of polytopal and tropical geometry (see [35, 186] for
seminal textbooks on these topics) to build a more sophisticated _tropical
measure_ $\mu^{\mathrm{tr}}$. This measure is tailored to estimating integrals
of functions $f(x)=\frac{g(x)^{a}}{h(x)^{b}}$ given by homogeneous polynomials
$g,h\in\mathbb{C}[x]$ and rational numbers $a,b\in\mathbb{Q}$. A good example
hereof is the evaluation of the integral
$\int_{0}^{1}\int_{0}^{1}\frac{\mathrm{d}x\wedge\mathrm{d}y}{x+y}=2\log 2$.
"Ordinary" MC fails because the variance of $\frac{1}{x+y}$ is infinite w.r.t.
the uniform measure. It is, however, easily evaluated with the tropical
measure. See [187, Section 3.1.1] for a detailed discussion.

##### 10.1.1 Tropical approximation

As we are ultimately interested in integrating homogeneous polynomials, let us
now define the notion of _real projective space_ $\mathbb{R}\mathbb{P}^{n-1}$;
this which will serve as an integration domain111In practice, this domain can
be parametrized as a simplex polytope. later on. $\mathbb{R}\mathbb{P}^{n-1}$
is defined by the set of points in $\mathbb{R}^{n}\setminus\\{0\\}$ modulo the
equivalence relation

$\displaystyle\begin{array}[]{rcl}(x_{0},\ldots,x_{n-1})&\sim&(y_{0},\ldots,y_{n-1})\\\
&\Updownarrow&\\\
\exists\lambda\in\mathbb{R}\setminus\\{0\\}\,:\,(x_{0},\ldots,x_{n-1})&=&(\lambda
y_{0},\ldots,\lambda y_{n-1})\,.\end{array}$ (10.8)

We denote the coordinates on $\mathbb{R}\mathbb{P}^{n-1}$ by
$[x_{0}:\ldots:x_{n-1}].$ Two points such as $[1:0:2]$ and
$[\tfrac{1}{2}:0:1]$ would be equivalent in $\mathbb{R}\mathbb{P}^{2}$ because
they are related by $(1,\,0,\,2)=(2\cdot\tfrac{1}{2},\,2\cdot 0,\,2\cdot 1).$
_Positive_ projective space is further given by

$\displaystyle\mathbb{P}^{n-1}_{+}=\\{x=[x_{0}:\cdots:x_{n-1}]\in\mathbb{R}\mathbb{P}^{n-1}\,|\,x_{i}>0\\}\,.$
(10.9)

Above we made the unusual choice of starting the index at $0$ rather than $1$.
This convention allows for seamless interoperability with the python
programming language, to be used in the next chapter.

Let $p$ be a homogeneous polynomial. In multi-index notation,

$\displaystyle p(x)=\sum_{k\in\mathrm{supp}\left(p\right)}c_{k}\,x^{k}\ ,$
(10.10)

where $x^{k}=x_{0}^{k_{0}}\cdots x_{n-1}^{k_{n-1}}$ and
$|k|=k_{0}+\ldots+k_{n-1}$ is equal for all multi-indices $k$. The _tropical
approximation_ [39] to $p$ is the central concept in tropical integration:

$\displaystyle
p^{\mathrm{tr}}(x)=\max_{k\in\mathrm{supp}\left(p\right)}x^{k}\,.$ (10.11)

If, for instance, $p(x)=x_{0}x_{1}x_{2}^{2}-2x_{1}^{4}+42ix_{0}^{2}x_{1}^{2},$
then
$p^{\mathrm{tr}}(x)=\max\big{\\{}x_{0}x_{1}x_{2}^{2},\,x_{1}^{4},\,x_{0}^{2}x_{1}^{2}\big{\\}}.$
In other words, the tropical approximation neglects the monomial coefficients,
and only cares about the largest monomial. The set of monomial exponents
$\mathrm{supp}\left(p\right)$ is hence the essential data for
$p^{\mathrm{tr}}$. Interpreting each multi-index in
$\mathrm{supp}\left(p\right)$ as a vector in $\mathbb{Z}^{n}_{\geq 0}$ and
taking the convex hull, we see that $p^{\mathrm{tr}}$ is simply a function
avatar of the _Newton polytope_ $\mathbf{N}[p]$ from (6.75).

Suppose that $\mathbf{N}[p]$ lives in $\mathbb{R}^{n}$ as parametrized by $x$.
By virtue of being a polytope, $\mathbf{N}[p]$ is given by a finite
intersection of _half-spaces_ , each of the form
$\\{v\in\mathbb{R}^{n}\,|\,v\cdot x\leq c\\}$ for some scalar
$c\in\mathbb{R}$. A subset $F\subset\mathbf{N}[p]$ is called a _face_ of
$\mathbf{N}[p]$ iff

  1. 1.

There exist $v\in\mathbb{R}^{n}$ and $c\in\mathbb{R}$ such that
$\mathbf{N}[p]$ is contained in $\\{v\in\mathbb{R}^{n}\,|\,v\cdot x\leq c\\}.$

  2. 2.

$F$ is the intersection of $\mathbf{N}[p]$ with
$\\{v\in\mathbb{R}^{n}\,|\,v\cdot x=c\\}.$

The polynomial $p$ truncated to a face $F$ is given by

$\displaystyle
p_{F}(x)=\sum_{k\,\in\,F\,\cap\,\mathrm{supp}\left(p\right)}c_{k}\,x^{k}\,.$
(10.12)

$p$ is called _completely non-vanishing_ on a domain $X$ if for each face $F$
of the Newton polytope we have $p_{F}(x)\neq 0,\,x\in X.$

Example 10.1. Consider the polynomial

$\displaystyle
p(x)=x_{0}+x_{1}+2x_{0}^{2}+x_{1}^{2}+3x_{0}x_{1}+2x_{0}^{2}x_{1}+x_{0}x_{1}^{2}+x_{0}^{3}\,.$
(10.13)

It is non-homogeneous, but that is inconsequential for the purpose of this
example. Its associated Newton polytope $\mathbf{N}[p]$ is shown in Figure
(10.1). As an intersection of half-spaces, the Newton polytope is represented
by

$\displaystyle\mathbf{N}[p]=\big{\\{}x_{0}\geq
0\big{\\}}\,\cap\,\big{\\{}x_{1}\geq 0\big{\\}}\,\cap\,\big{\\{}x_{1}\leq
2\big{\\}}\,\cap\,\big{\\{}x_{0}+x_{1}\geq
1\big{\\}}\,\cap\,\big{\\{}x_{0}+x_{1}\leq 3\big{\\}}\,.$ (10.14)

Consider a face $F$ of $\mathbf{N}[p]$, for example the bottom line segment
touching the $x_{0}$-axis given by $F=\\{(x_{0},x_{1})\,|\,1\leq x_{0}\leq 3\
\text{and}\ x_{1}=0\\}.$ $F$ is indeed a face since 1) $\mathbf{N}[p]$ is
contained in the half-space $\\{x_{1}\geq 0\\},$ and 2) $F$ is given by an
intersection $F=\mathbf{N}[p]\,\cap\,\\{x_{1}=0\\}.$ The polynomial truncated
to this face is given by

$\displaystyle p_{F}(x)=x_{0}+2x_{0}^{2}+x_{0}^{3}\,.$ (10.15)

Because each monomial coefficient of $p$ is positive, we have that $p$ is
completely non-vanishing on, say, $X=\mathbb{R}^{2}_{>0}$. $\blacksquare$

Figure 10.1:  Newton polytope for the polynomial $p(x)$ in (10.13). The black
dots represent the monomial exponents of $p(x)$, and the pink area is their
convex hull in $\mathbb{R}^{2}$ as parametrized by $(x_{0},x_{1})$.

Borinsky proved that there is really a sense in which $p^{\mathrm{tr}}(x)$
_approximates_ the values of $p(x)$; this is the backbone behind tropical
integration. Take a homogeneous polynomial $p\in\mathbb{C}[x]$ that is
completely non-vanishing on positive projective space $\mathbb{P}^{n-1}_{+}$.
Then there exist two constants $C_{1},C_{2}>0$ such that [40, Theorem 8]

$\displaystyle C_{1}\ \leq\ \frac{|p(x)|}{p^{\mathrm{tr}}(x)}\ \leq\
C_{2}\quad\text{for all}\quad x\in\mathbb{P}^{n-1}_{+}\,.$ (10.16)

In other words, the tropical approximation of $p$ bounds the original
polynomial from above and from below.

We can now describe the intuition behind tropical integration (a precise
algorithm is given in Section (10.3.2)). Consider a convergent integral

$\displaystyle I=\int_{\mathbb{P}^{n-1}_{+}}\frac{g(x)^{a}}{h(x)^{b}}\,\Omega$
(10.17)

given by two homogeneous polynomials $g$ and $h$ that are completely non-
vanishing on the integration domain, and $a,b\in\mathbb{Q}$. We insert a
clever way of writing the number "$1$" into the integrand:

$\displaystyle\begin{array}[]{lllll}I&=&\displaystyle\int_{\mathbb{P}^{n-1}_{+}}\frac{g(x)^{a}}{h(x)^{b}}\,\Omega&\times&\displaystyle\left(\frac{g^{\mathrm{tr}}(x)}{g^{\mathrm{tr}}(x)}\right)^{a}\displaystyle\left(\frac{h^{\mathrm{tr}}(x)}{h^{\mathrm{tr}}(x)}\right)^{b}\vspace{0.3cm}\\\
&=&\displaystyle\underbrace{\int_{\mathbb{P}^{n-1}_{+}}\frac{g^{\mathrm{tr}}(x)^{a}}{h^{\mathrm{tr}}(x)^{b}}\,\Omega}_{\propto\
\text{tropical measure}\
\mu^{\mathrm{tr}}}&\times&\underbrace{\displaystyle\left(\frac{g(x)}{g^{\mathrm{tr}}(x)}\right)^{a}\displaystyle\left(\frac{h^{\mathrm{tr}}(x)}{h(x)}\right)^{b}}_{\text{bounded
due to }\eqref{tropical_approximation_property}}\,.\end{array}$ (10.20)

The integrand has been factored into something resembling a probability
measure (after being suitably normalized) multiplied by a function that is
square-integrable as a consequence of being bounded. This is now a suitable
setup for MC integration: we sample a large number of random points $x^{(i)}$
from $\mu^{\mathrm{tr}}$, and then take the average of the function
$(g/g^{\mathrm{tr}})^{a}\,(h^{\mathrm{tr}}/h)^{b}$ evaluated at these points.

####  Section 10.2Projective Feynman integrals

As we will present a numerical code later on, let us be explicit about our
conventions for FIs. Take $G$ to be a one-particle irreducible Feynman diagram
with edge set $E$ and vertex set $V$. With our index convention, we have
$E=\\{0,1,\ldots,|E|-1\\}$ and $V=\\{0,1,\ldots,|V|-1\\}.$ Every edge $e\in E$
is equipped with a mass $m_{e}$ and an edge weight $\nu_{e}$. Every vertex
$v\in V$ comes with an incoming momentum vector $p_{v}$. Internal vertices are
distinguished by having vanishing incoming momenta; i.e., writing
$V=V_{\text{ext}}\,\sqcup\,V_{\text{int}},$ then $p_{v}=0$ if $v\in
V_{\text{int}}$.

The _incidence matrix_ $\mathcal{E}$ of $G$ is defined by choosing an
arbitrary orientation for the edges, and setting $\mathcal{E}_{v,e}=\pm 1$ if
edge $e$ points to/from vertex $v$ and $\mathcal{E}_{v,e}=0$ if $e$ is not
incident to $v$.

In terms of this data, we write the FI associated to $G$ as

$\displaystyle I_{G}=\int\prod_{e\in
E}\frac{\mathrm{d}^{\mathrm{D}}\ell_{e}}{i\pi^{\mathrm{D}/2}}\left(\frac{-1}{\ell_{e}^{2}-m_{e}^{2}+i\varepsilon}\right)^{\nu_{e}}\prod_{v\in
V\setminus\\{v_{0}\\}}i\pi^{\mathrm{D}/2}\delta^{(\mathrm{D})}\Big{(}p_{v}+\sum_{e\in
E}\mathcal{E}_{v,e}\,\ell_{e}\Big{)}\,.$ (10.21)

The loop momenta $\ell_{e}$ are integrated over $\mathrm{D}$-dimensional
Minkowski space with mostly-minus signature,
$\ell_{e}^{2}=\left(\ell^{0}_{e}\right)^{2}-\left(\ell^{1}_{e}\right)^{2}-\left(\ell^{2}_{e}\right)^{2}-\ldots,$
and an overall momentum-conserving $\delta$-function
$\delta(p_{v_{0}}+\ldots+p_{v_{|V|-1}})$ has been extracted by removing a
vertex $v_{0}\in V$.

To evaluate $I_{G}$ numerically, we use an equivalent parametric
representation (see [44, Section 2.5] for a derivation starting from momentum
space):

$\displaystyle I_{G}$
$\displaystyle=\Gamma(\omega)\int_{\mathbb{P}^{|E|-1}_{+}}\phi\qquad\text{with}$
(10.22) $\displaystyle\phi$ $\displaystyle=\left(\prod_{e\in
E}\frac{x_{e}^{\nu_{e}}}{\Gamma(\nu_{e})}\right)\frac{1}{\mathcal{U}(x)^{\mathrm{D}/2}}\left(\frac{1}{\mathcal{V}(x)-i\varepsilon\sum_{e\in
E}x_{e}}\right)^{\omega}\Omega\,.$

This is an integral over positive projective space $\mathbb{P}^{|E|-1}_{+}$,
defined in (10.9), w.r.t. the volume form

$\displaystyle\Omega=\sum_{e=0}^{|E|-1}(-1)^{|E|-e-1}\,\frac{\mathrm{d}x_{0}\wedge\cdots\wedge\widehat{\mathrm{d}x_{e}}\wedge\cdots\wedge\mathrm{d}x_{|E|-1}}{x_{1}\cdots\widehat{x}_{e}\cdots
x_{|E|-1}}\,,$ (10.23)

where hats denote omission. The superficial degree of divergence of the graph
$G$ is given by $\omega=\sum_{e\in E}\nu_{e}-\mathrm{D}\cdot L/2$ where
$L=|E|-|V|+1$ is the number of loops.

$\mathcal{U}(x)$ is the first Symanzik polynomial, previously defined in
(7.1). Recalling also the second Symanzik polynomial $\mathcal{F}(x)$ from
(7.2), we set

$\displaystyle\mathcal{V}(x)=\frac{\mathcal{F}(x)}{\mathcal{U}(x)}\,.$ (10.24)

Since $\mathcal{U}(\lambda x)=\lambda^{L}\mathcal{U}(x)$ and
$\mathcal{F}(\lambda x)=\lambda^{L+1}\mathcal{F}(x),$ one can quickly verify
that the integrand $\phi$ is invariant under the coordinate rescaling
$x\to\lambda x,$ wherefore it is appropriate to call (10.22) a _projective
integral representation_.

The graph-theoretic formulas (7.1) and (7.2) for $\mathcal{U}$ and
$\mathcal{F}$ are in fact not so efficient for numerically evaluating millions
of points during MC integration. The reason is that the number of monomials
grows quickly as a function of $|E|$ and $L$, making it slow to add up all
terms during each numerical evaluation. Let us therefore present more
efficient formulas in terms of matrices that scale linearly with the
complexity of $G$. To start, define the $(|V|-1)\times(|V|-1)$ _reduced
Laplacian matrix_ component-wise by

$\displaystyle\mathcal{L}_{uv}=\sum_{e\in
E}\frac{\mathcal{E}_{u,e}\,\mathcal{E}_{v,e}}{x_{e}}\quad\text{for}\quad
u,v\,\in\,V\setminus\\{v_{0}\\}\,.$ (10.25)

As this matrix is symmetric, it is only necessary to explicitly calculate the
diagonal and the upper triangular block. The lower triangular block can then
be copied from the upper one. It turns out that $\mathcal{L}$ is positive
definite.

Let $\mathcal{P}_{uv}=p_{u}\cdot p_{v}$ be the $(|V|-1)\times(|V|-1)$ matrix
containing Minkowski scalar products between all independent momenta.

We then have the identities

$\displaystyle\mathcal{U}(x)$
$\displaystyle=\det\big{[}\mathcal{L}(x)\big{]}\prod_{e\in E}x_{e}$ (10.26)
$\displaystyle\mathcal{F}(x)$
$\displaystyle=\mathcal{U}(x)\left[-\sum_{u,v\,\in\,V\setminus\\{v_{0}\\}}\big{(}\mathcal{L}^{-1}(x)\big{)}_{uv}\,\mathcal{P}_{uv}+\sum_{e\in
E}m_{e}^{2}\,x_{e}\right]\,.$ (10.27)

Note that the first sum in $\mathcal{F}$ can be expressed as
$\text{tr}\left[\mathcal{L}^{-1}\cdot\mathcal{P}\right].$ Determinants and
traces are computationally cheap, but the matrix inverse $\mathcal{L}^{-1}$
appears to be expensive at first glance. However, notice that we are only
interested in the combination $\mathsf{X}=\mathcal{L}^{-1}\cdot\mathcal{P},$
not $\mathcal{L}^{-1}$ itself. Owing to the fact that $\mathcal{L}$ is
positive definite, then $\mathsf{X}$ can be swiftly evaluated using _Cholesky
decomposition_. More precisely, in $O(|V|^{3})$ time, it is possible to
decompose $\mathcal{L}=\mathsf{L}\cdot\mathsf{D}\cdot\mathsf{L}^{T}$ in terms
of a lower triangular matrix $\mathsf{L}$ and a diagonal matrix $\mathsf{D}$.
Afterwards, the solution $\mathsf{X}$ to the linear system
$\mathcal{L}\cdot\mathsf{X}=\mathcal{P}$ is immediately obtained by solving

$\displaystyle\mathsf{L}\cdot\mathsf{Y}_{1}=\mathcal{P}\,,\quad\mathsf{D}\cdot\mathsf{Y}_{2}=\mathsf{Y}_{1}\,,\quad\mathsf{L}^{T}\cdot\mathsf{X}=\mathsf{Y}_{2}$
(10.28)

in cascade.

Later on, we will actually require $\mathcal{L}$ to take on complex values
rather than real ones. Then $\mathcal{L}$ is no longer Hermitian, wherefore
Cholesky decomposition fails. In this case, one must use the less efficient LU
(lower-upper) decomposition.

##### 10.2.1 $i\varepsilon$-prescription via contour deformation

For the purposes of this thesis, say that a FI computation is in the
_Euclidean regime_ if $\mathcal{F}(x)\geq 0$ for all
$x\in\mathbb{P}^{|E|-1}_{+}$ (a more refined notion of the Euclidean regime is
given in [185, Section 2.2]). Then it is only the case that $\mathcal{F}(x)=0$
on the boundary $x\to 0$ of the integration domain, which is assumed to be an
integrable singularity. In our convention, being in the Euclidean regime is
analogous to having a FI with scalar products computed in the Euclidean all-
minus metric. Life is simple in the Euclidean regime because one avoids the
singularities in (10.21) coming from setting the time-components equal to
$\ell^{0}=\pm\sqrt{\vec{\ell}\cdot\vec{\ell}+m^{2}}$ in the propagators.
Consequently, Feynman’s $i\varepsilon$-prescription can be dropped, and the
resulting FI will be purely real.

Alas, life is not so simple in Minkowski space. Let us say that we are in the
_Minkowski regime_ when we are not in the Euclidean regime. In this case, the
monomial coefficients of $\mathcal{F}(x)$ can be both positive and negative.
The $i\varepsilon$ regulates singularities coming from having
$\mathcal{F}(x)=0$ inside the integration domain. In momentum space, the
$i\varepsilon$ regulates the $\ell^{0}$-singularities by pushing them slightly
away from the real axis. An alternative description is to deform the
integration contour given by $\ell^{0}\in(-\infty,\infty)$ to one that goes
around the poles222 With two poles there are four possible choices for how to
go around them, but only one is consistent with the notion of causality. This
contour corresponds to a propagator $\ell^{2}-m^{2}+i\varepsilon$. Two non-
casual examples are
$(\ell^{0}+i\varepsilon)^{2}-\vec{\ell}\cdot\vec{\ell}-m^{2}$ and
$\ell^{2}-m^{2}-i\varepsilon$. The correct contour goes below the singularity
at $-\sqrt{\vec{\ell}\cdot\vec{\ell}+m^{2}}$ and above the one at
$+\sqrt{\vec{\ell}\cdot\vec{\ell}+m^{2}}.$ I.e.  in the complex
$\ell^{0}$-plane. .

The $i\varepsilon$-prescription only dictates on which side of the poles we
integrate. But for the purpose of numerical integration we must make an
_explicit_ choice for the contour, preferably with decent numerical stability.
Moreover, the deformation ought to preserve the projective symmetry of the
integral in order to play well with tropical MC sampling. We find that the
following deformation into the complex plane proposed in [188, 189] fulfils
both of these criteria:

$\displaystyle
X_{e}(x)=x_{e}\exp\left[-i\lambda\frac{\partial\mathcal{V}(x)}{\partial
x_{e}}\right]\quad\text{for}\quad\lambda\in\mathbb{R}_{>0}\,.$ (10.29)

Because $\mathcal{U}$ and $\mathcal{F}$ are homogeneous polynomials of degrees
$L$ and $L+1$ respectively, then $\mathcal{V}=\mathcal{F}/\mathcal{U}$ is a
rational function of degree $1$, wherefore
$\partial_{e}\mathcal{V}=\partial\mathcal{V}(x)/\partial x_{e}$ has
homogeneous degree $0$. Thus, (10.29) preserves projective symmetry.

The value of parameter $\lambda$ is finite but small. This is in contrast to
the infinitesimal $\varepsilon$ of the $i\varepsilon$-prescription. As long as
$\lambda$ is is small enough that no poles of the integrand $\phi$ are hit,
then Cauchy’s theorem guarantees that the deformation is valid. Note that
$\lambda$ must have a mass dimension in order to cancel that of
$\partial_{e}\mathcal{V}$ in the exponential. In practice, we shall choose
$\lambda$ to be of order $\mathcal{O}(1/\Lambda^{2})$, where $\Lambda^{2}$ is
the largest kinematic scale in $\mathcal{V}$.

To see why (10.29) is consistent with the $i\varepsilon$-prescription, let us
consider the _Schwinger representation_ of a FI [44, Section 2.5.2]:

$\displaystyle
I=\int_{0}^{\infty}\ldots\exp\left[i\left(-\mathcal{V}(x)+i\varepsilon\sum_{e\in
E}x_{e}\right)\right]\,.$ (10.30)

Here we ignored factors in the integrand which are irrelevant for the present
argument. Observe that the $i\varepsilon$ term cures the divergence at
$x_{e}\to\infty$ by gifting an exponential damping $\sim\exp\left[-\varepsilon
x_{e}\right].$ Notice that a similar damping is achieved by adding a
_positive_ imaginary part to the first term $-\mathcal{V}(x)$. In that case
the $i\varepsilon$ could be dropped. Indeed, changing coordinates according to
(10.29), then a Taylor expansion in $\lambda$ shows that

$\displaystyle-\mathcal{V}(X)=-\mathcal{V}(x)+i\lambda\sum_{e\in
E}x_{e}\left[\frac{\partial\mathcal{V}(x)}{\partial
x_{e}}\right]^{2}+\mathcal{O}(\lambda^{2})\,.$ (10.31)

Assuming that 1) there are no solutions to the _Landau equations_

$\displaystyle x_{e}\frac{\partial\mathcal{V}(x)}{\partial
x_{e}}=0\quad\text{for each}\quad e\in E\quad\text{and for any}\quad
x\in\mathbb{P}^{|E|-1}_{+}\,,$ (10.32)

and 2) $\lambda$ is small enough such that the $\mathcal{O}(\lambda^{2})$
terms to not change the sign of the $\mathcal{O}(\lambda)$ term, then (10.31)
has a positive imaginary part, as desired. For the remaining part of this
thesis, it is assumed that the Landau equations have no solutions.

A short calculation shows that the $|E|\times|E|$ Jacobian matrix associated
to (10.29) is

$\displaystyle(\mathcal{J}_{\lambda}(x))_{eh}=\delta_{eh}-i\lambda
x_{e}\frac{\partial^{2}\mathcal{V}(x)}{\partial x_{e}\partial
x_{h}}\quad\text{for}\quad e,h\in E\,.$ (10.33)

Coming back to the projective FI (10.22), its deformed version now reads

$\displaystyle
I_{G}=\Gamma(\omega)\int_{\mathbb{P}^{|E|-1}_{+}}\left(\prod_{e\in
E}\frac{X_{e}^{\nu_{e}}}{\Gamma(\nu_{e})}\right)\frac{\det\mathcal{J}_{\lambda}(x)}{\mathcal{U}(X)^{\mathrm{D}/2}\
\times\ \mathcal{V}(X)^{\omega}}\,\Omega\,,$ (10.34)

where $X=(X_{0},\ldots,X_{|E|-1})$ with $X_{e}$ given by (10.29).

##### 10.2.2 $\epsilon$-expansion for quasi-finite integrals

We have so far not made any assumptions regarding the finiteness of the FI
(10.34) as the DR parameter $\epsilon$ tends to zero. Recall that $\epsilon$
appears inside $\mathrm{D}=\mathrm{D}_{0}-2\epsilon$ and $\omega=\sum_{e\in
E}\nu_{e}-\mathrm{D}L/2$ in this representation.

We shall restrict our attention to the category of integrals dubbed as _quasi-
finite_ in [190]. This means that $I_{G}$ in (10.34) can have overall
$\epsilon$-poles coming from the prefactor $\Gamma(\omega)$, but the integral
$\int_{\mathbb{P}^{|E|-1}_{+}}(\ldots)$ is itself finite as $\epsilon\to 0$.
Note that this is less restrictive than requiring the whole FI to be finite.

Given this, the $\epsilon$-expansion of $I_{G}$ is derived by merely Taylor
expanding the integrand. Writing $\omega_{0}=\sum_{e\in
E}\nu_{e}-\mathrm{D}_{0}L/2,$ the expansion becomes

$\displaystyle I_{G}$ $\displaystyle=\Gamma(\omega_{0}+\epsilon
L)\sum_{k=0}^{\infty}\frac{\epsilon^{k}}{k!}\ \times$ (10.35)
$\displaystyle\int_{\mathbb{P}^{|E|-1}_{+}}\left(\prod_{e\in
E}\frac{X_{e}^{\nu_{e}}}{\Gamma(\nu_{e})}\right)\frac{\det\mathcal{J}_{\lambda}(x)}{\mathcal{U}(X)^{\mathrm{D}_{0}/2}\
\times\
\mathcal{V}(X)^{\omega_{0}}}\log^{k}\left[\frac{\mathcal{U}(X)}{\mathcal{V}(X)^{L}}\right]\Omega\,.$

If the integral is finite for $k=0$, then the $\log^{k}$ factors cannot spoil
the convergence for higher orders in $k$.

Our ultimate goal is to numerically evaluate the integrals appearing in the
$\epsilon$-expansion above. To this end, it is important that the tropical
approximation property (10.16) holds true, but this is now complicated by the
fact that $X$ depends on $\lambda$. So we make the following _assumption_ :
there exist $\lambda$-dependent constants $C_{1}(\lambda),C_{2}(\lambda)>0$
such that for small enough $\lambda>0$,

$\displaystyle C_{1}(\lambda)\ \leq\
\left|\left(\frac{\mathcal{U}^{\mathrm{tr}}(x)}{\mathcal{U}(X)}\right)^{\mathrm{D}_{0}/2}\left(\frac{\mathcal{V}^{\mathrm{tr}}(x)}{\mathcal{V}(X)}\right)^{\omega_{0}}\right|\
\leq\ C_{2}(\lambda)\quad\text{for all}\quad x\in\mathbb{P}^{|E|-1}_{+}\,.$
(10.36)

We have checked this assumption through extensive numerical testing, but lack
a mathematical argument for its validity.

The statement (10.36) suggest an optimal value of $\lambda$: choose it such
that the bounds are as tight as possible. This is an interesting research
direction, as it would constitute the first _canonical_ contour deformation
without any free parameters.

##### 10.2.3 Convergence of Euler integrals

It was proven by Berkesch, Forsgård and Passare that a quasi-finite
representation can always be found [191, Theorem 2.4]. The proof of this
theorem gives a concrete prescription for factoring out $\epsilon$-poles in
front of finite integrals, but we have not yet implemented this protocol333
Their method requires expanding out derivatives of $\mathcal{U}$ and
$\mathcal{V}$, which can lead to a huge proliferation of terms if done
symbolically. This problem is naturally avoided in the setting of numerical MC
integration. Indeed, we shall provide efficient formulas in Section (10.3.3)
for derivatives of $\mathcal{U}$ and $\mathcal{V}$ that could be used in a
future implementation. .

Assuming a quasi-finite basis of MIs is known, one can alternatively employ a
combination of dimensional recurrence relations and IBPs as in [190] to write
a given FI as linear combination of said basis.

To find a quasi-finite basis, one ought to have test for finiteness of Euler
integrals (in contrast to the GKZ setting studied earlier in this thesis, we
are now allowing the monomial coefficients of the Euler integrands to be fixed
to special values). Convergence theorems for Euler integrals were first
discovered in the mathematics literature by Berkesch, Forsgård, Nilsson and
Passare [192, Theorem 1] [191, Theorem 2.2]. In the physics literature,
convergence criteria were later established by Panzer, Schabinger and von
Manteuffel in [190], followed by works of Arkani-Hamed et al. [193, Claim 1]
and Borinsky [40, Theorem 3].

The most general, proven statement is roughly as follows. Suppose
$\\{f_{1}(x),\ldots,f_{l}(x)\\}$ are completely non-vanishing polynomials on
$\mathbb{R}^{n}_{>0}$ (an equivalent statement can be formulated in projective
space). Let $s=(s_{1},\ldots,s_{l})\in\mathbb{C}^{l}$ and
$\nu=(\nu_{1},\ldots,\nu_{n})\in\mathbb{C}^{n}$ be complex parameters. Assume
$\mathrm{Re}(s_{i})>0$. The integral

$\displaystyle I=\int_{\mathbb{R}^{n}_{>0}}\frac{x_{1}^{\nu_{1}}\cdots
x_{n}^{\nu_{n}}}{f_{1}^{s_{1}}\cdots f_{l}^{s_{l}}}\frac{\mathrm{d}x}{x}$
(10.37)

is convergent iff444 An additional, weak assumption is that the Minkowski sum
of Newton polytopes $\mathbf{N}[f_{1}]+\ldots+\mathbf{N}[f_{l}]$ should have
dimension $n$.

$\displaystyle\left[\begin{array}[]{c}\mathrm{Re}(\nu_{1})\\\ \vdots\\\
\mathrm{Re}(\nu_{n})\end{array}\right]\ \in\
\mathrm{int}(\mathbf{N}(s))\quad\text{where}\quad$ (10.41)
$\displaystyle\mathbf{N}(s)=\mathrm{Re}(s_{1})\cdot\mathbf{N}[f_{1}]+\cdots+\mathrm{Re}(s_{l})\cdot\mathbf{N}[f_{l}]\,.$

  * •

The notation $s\cdot\mathbf{N}[f]$ means that the Newton polytope associated
to $f$ is rescaled by a constant $s\in\mathbb{R}_{\geq 0}$. When
$\mathbf{N}[f]$ is determined by a collection half-spaces $\\{v_{i}\cdot x\leq
c_{i}\\},$ then $s\cdot\mathbf{N}[f]$ has the effect of rescaling the $c$’s,
that is $\\{v_{i}\cdot x\leq s\,c_{i}\\}.$

  * •

The _Minkowski sum_ of two polytopes $\mathbf{N}[f_{A}]$ and
$\mathbf{N}[f_{B}]$ is given by
$\mathbf{N}[f_{A}]+\mathbf{N}[f_{B}]=\big{\\{}a+b\,\big{|}\,a\in\mathbf{N}[f_{A}]\
\text{and}\ b\in\mathbf{N}[f_{B}]\big{\\}}.$

  * •

A point $x$ is in the _interior_ of $\mathbf{N}(s)$,
$x\in\mathrm{int}(\mathbf{N}(s))$, if there exists a small ball (w.r.t. the
Euclidean metric) centered at $x$ that is completely contained in
$\mathbf{N}(s)$.

The condition (10.41) can be checked in practice as follows. First compute the
half-space representation of $\mathbf{N}(s)$ (for example with sage or
polymake). This gives a set of inequalities of the form $\\{v\cdot x\leq c\\}$
for the hyperplanes that cut out $\mathbf{N}(s)$. For FIs we have
$s_{1}=\mathrm{D}/2$ and $s_{2}=\omega$, so these inequalities would depend
parametrically on the DR parameter $\epsilon$. A vector
$v=\left[\mathrm{Re}(\nu_{1}),\ldots,\mathrm{Re}(\nu_{n})\right]^{T}$ is then
in the interior of $\mathbf{N}(s)$ if the inequalities are strict, namely
$\\{v\cdot x<c\\}.$

This suggests the following heuristic algorithm for determining a quasi-finite
basis: check whether the criterion (10.41) is fulfilled during a bottom-up
search in the space $(\nu,\mathrm{D})$ of propagator powers and spacetime
dimensions555 The search space for $\nu$ is cut down drastically in case _any_
other basis is known, so as to pinpoint a list of relevant integral sectors. .

We remark that the alternative quasi-finite basis search algorithm from [190]
has been automated in Reduze2 [12] via the command find_finite_masters.

####  Section 10.3Tropical Monte Carlo sampling, in detail

We have already briefly explained the idea behind tropical MC sampling in
Section (10.1). This algorithm is here described in detail. The goal is to
tropically integrate the terms in the $\epsilon$-expansion (10.35), so let us
factorize the integrals in a way analogous to (10.20):

$\displaystyle I_{G}=I^{\mathrm{tr}}\ \frac{\Gamma(\omega_{0}+\epsilon
L)}{\prod_{e\in E}\Gamma(\nu_{e})}\
\sum_{k=0}^{\infty}\frac{\epsilon^{k}}{k!}I^{(k)}\,,$ (10.42)

where the $k$th order integral reads as

$\displaystyle
I^{(k)}=\int_{\mathbb{P}^{|E|-1}_{+}}\frac{\det\mathcal{J}_{\lambda}(x)\
\times\ \prod_{e\in
E}(X_{e}/x_{e})^{\nu_{e}}}{\big{(}\mathcal{U}(X)/\mathcal{U}^{\mathrm{tr}}(x)\big{)}^{\mathrm{D}_{0}/2}\
\times\
\big{(}\mathcal{V}(X)/\mathcal{V}^{\mathrm{tr}}(x)\big{)}^{\omega_{0}}}\log^{k}\left[\frac{\mathcal{U}(X)}{\mathcal{V}(X)^{L}}\right]\
\times\ \mu^{\mathrm{tr}}\,,$ (10.43)

and the _tropical probability measure_ is

$\displaystyle\mu^{\mathrm{tr}}=\frac{1}{I^{\mathrm{tr}}}\ \frac{\prod_{e\in
E}x_{e}^{\nu_{e}}}{\mathcal{U}^{\mathrm{tr}}(x)^{\mathrm{D}_{0}/2}\ \times\
\mathcal{V}^{\mathrm{tr}}(x)^{\omega_{0}}}\ \Omega\,.$ (10.44)

The normalization factor $I^{\mathrm{tr}}$ ensures that
$\int_{\mathbb{P}^{|E|-1}_{+}}\mu^{\mathrm{tr}}=1$ (a formula for
$I^{\mathrm{tr}}$ is given further down in (10.54)).

The idea of tropical MC is to sample $N$ points $\\{x^{(1)},\ldots,x^{(N)}\\}$
from the measure $\mu^{\mathrm{tr}}$, then evaluate the integrand of (10.43)
at those points, and finally take the average of these evaluations. Borinsky
introduced two different algorithms that sample from $\mu^{\mathrm{tr}}$ in
[40]. The first algorithm [40, Section 5] works for any projective integral
over rational functions. A key step in the algorithm requires the
triangulation of a certain _normal fans_ 666 A normal fan is essentially a
collection of cones that can be associated to the normal directions of each
face of a given polytope [186].  associated to the Newton polytopes of
$\mathcal{U}$ and $\mathcal{F}$. The computation of such triangulations
unfortunately becomes a bottleneck when a Feynman diagram has many edges.

The second algorithm [40, Section 6] takes special features of $\mathcal{U}$
and $\mathcal{F}$ into account, thereby avoiding the triangulation step.
Namely, the Newton polytope of $\mathcal{U}$ is always a so-called
_generalized permutahedron_ , and the same holds true for the Newton polytope
of $\mathcal{F}$ given sufficiently generic kinematics. The triangulation step
of the previous algorithm gets traded for a much simpler preprocessing step,
involving the calculation of a certain recursive function called $J(\gamma)$
below.

##### 10.3.1 Generalized permutahedra

Let $E$ be any finite set. Consider a function777The notation $2^{E}$ stands
for the power set of $E$. $z:2^{E}\to\mathbb{R}$ that assigns a number to each
subset $\gamma$ of $E$. Given such a $z$, the _base polytope_
$\mathbf{P}[z]\subset\mathbb{R}^{|E|}$ is defined by the set of points
$(a_{0},\ldots,a_{|E|-1})\in\mathbb{R}^{|E|}$ satisfying

$\displaystyle\sum_{e\in E}a_{e}$ $\displaystyle=z(E)\quad\text{and}\quad$
(10.45) $\displaystyle\sum_{e\in\gamma}a_{e}$ $\displaystyle\geq
z(\gamma)\quad\text{for all}\quad\gamma\subsetneq E\,.$

The function $z$ is called _supermodular_ if it satisfies the inequalities

$\displaystyle z(\gamma)+z(\delta)\ \leq\
z(\gamma\,\cap\,\delta)+z(\gamma\,\cup\,\delta)\quad\text{for all pairs of
subsets}\quad\gamma,\delta\subset E\,.$ (10.46)

Given a supermodular function $z$, Aguiar and Ardila proved that the base
polytope $\mathbf{P}[z]$ is always a _generalized permutahedron_ (GP) [194,
Theorem 12.3] (as a shorthand, we shall often say that a polytope "has the GP
property", or simply "is GP"). This is a wonderful class of polytopes defined
by Postnikov in [195]. For the purposes of this discussion, we shall simply
_define_ a GP to a supermodular base polytope. To build some intuition for
these polytopes, let us nevertheless give a geometric description of
Postnikov’s original construction. To start, define the _standard
permutahedron_ by

$\displaystyle\Pi_{n}=\texttt{ConvexHull}\\!\left[\text{permutations of}\
(1,\ldots,n)\right]\,\subset\,\mathbb{R}^{n}\,.$ (10.47)

In words, we take the convex hull in $\mathbb{R}^{n}$ of the vector
$\mathrm{v}=(1,\ldots,n)$ and all its $n!$ permutations
$\mathrm{v}_{\sigma}=(\sigma(1),\ldots,\sigma(n)).$ Equivalently, each
permutation of $\mathrm{v}$ corresponds to a vertex. The dimension of
$\Pi_{n}$ is $n-1$ rather than $n$ because its points are all constrained to
lie on the hyperplane $1+2+\ldots+n=n(n+1)/2.$ Postnikov then defines a GP as
a special deformation of a standard permutahedron: one is allowed to move any
facet along its normal ray as long as it does not pass through any vertices.
Equivalently, one can move any vertex as long as the directions are preserved.
Examples of generalized permutahedra are shown in Figure (10.2).

Figure 10.2:  (I) The standard permutahedron $\Pi_{3}$. Each vertex is labeled
by a permutation of $(1,2,3)$. Since the vector components
$(x_{1},x_{2},x_{3})$ of each vertex satisfies $x_{1}+x_{2}+x_{3}=6$, the
dimension of the polytope drops from 3 to 2. (II) Deforming the two red
vertices while preserving edge direction gives a GP. (III) Continuing to
deform the vertices creates another GP. A facet has now degenerated to a
point.

###### GP property of $\mathcal{U}$

Let us now see how this connects to FIs. Take the set $E$ to be the edge set
of a Feynman graph $G$. A subset of $E$ corresponds to a subgraph $\gamma$ of
$G$, so $z$ can be thought of as a collection of numbers for each subgraph of
$G$. The integer $L_{\gamma}$ denotes the number of loops of $\gamma$. It was
shown by Schultka [196] that the function

$\displaystyle z_{\mathcal{U}}(\gamma)=L_{\gamma}$ (10.48)

is supermodular, and Newton polytope $\mathbf{N}[\mathcal{U}]$ of the
$\mathcal{U}$-polynomial equals the base polytope
$\mathbf{P}[z_{\mathcal{U}}]$. According to the theorem by Aguiar and Ardila,
this means that $\mathbf{N}[\mathcal{U}]$ is GP.

###### GP property of $\mathcal{F}$

The situation is much more complicated for the $\mathcal{F}$-polynomial.
Following Brown [197], call a subgraph $\gamma\subset E$ _mass-momentum
spanning_ if the $\mathcal{F}$-polynomial restricted to the cograph
$G\setminus\gamma$ vanishes identically, $\mathcal{F}_{G\setminus\gamma}=0$.
Since the $\mathcal{F}$-polynomial contains all the kinematic scales
($s,t,m^{2}$ etc.) of a given FI, an equivalent definition is to say that
$\mathcal{F}_{G\setminus\gamma}$ is scaleless.

We further define the notion of _generic kinematics_. Consider the
$\mathcal{F}$-polynomial of a 1-loop bubble graph with two different internal
masses:
$\mathcal{F}_{\text{bubble}}(x)=(p^{2}-m_{0}^{2}-m_{1}^{2})x_{0}x_{1}+m_{0}^{2}x_{0}^{2}+m_{1}^{2}x_{1}^{2}.$
Notice that if the kinematics are tuned to enforce
$p^{2}-m_{0}^{2}-m_{1}^{2}=0$, then the first monomial drops out. Special
choices of kinematics can therefore change the Newton polytope
$\mathbf{N}[\mathcal{F}]$, as its vertices stem from the monomial exponents of
$\mathcal{F}$. The question is now how to avoid this degeneration of
$\mathbf{N}[\mathcal{F}]$ by setting a criterion which prevents cancellation
among momentum and mass terms in $\mathcal{F}$. To this end, begin by
recalling the scalar product matrix $\mathcal{P}_{uv}=p_{u}\cdot p_{v}$, where
$u,v\in V$ are vertices of $G$. Call a vertex $v$ _internal_ if
$\mathcal{P}_{uv}=0$ for all $u\in V$ and _external_ otherwise. We then say
that a kinematic configuration is _generic_ if

$\displaystyle\sum_{u,v\,\in\,V^{\prime}}\mathcal{P}_{uv}\ \neq\
\sum_{e\,\in\,E^{\prime}}m_{e}^{2}$ (10.49)

for each proper subset of external vertices $V^{\prime}\subsetneq
V_{\text{ext}}$ and each non-empty subset of edges $E^{\prime}\subset E$. When
the kinematics are not generic, they are called _exceptional_. For example,
kinematics are always generic in the Euclidean regime if 1) $m_{e}>0$ for all
$e\in E$, or 2) $\sum_{u,v\,\in\,V^{\prime}}\mathcal{P}_{uv}<0$ for all
$V^{\prime}\subsetneq V_{\text{ext}}$. See [185, Section 2.2] for additional
details.

Now, assuming generic kinematics, it was shown by Schultka [196]888 Schultka
technically proved this theorem with kinematics in the Euclidean regime. The
same result also holds in the Minkowski regime because the Newton polytopes
coincide in both regimes given the genericity assumption.  that the function

$\displaystyle z_{\mathcal{F}}=\begin{cases}L_{\gamma}+1&\text{if $\gamma$ is
mass-momentum spanning}\\\ L_{\gamma}&\text{else}\,.\end{cases}$ (10.50)

is supermodular, and the Newton polytope $\mathbf{N}[\mathcal{F}]$ equals the
base polytope $\mathbf{P}[z_{\mathcal{F}}]$. Therefore, for generic
kinematics, the Newton polytope $\mathbf{N}[\mathcal{F}]$ is GP. There are
reasons to expect that $\mathbf{N}[\mathcal{F}]$ is also GP even for
exceptional kinematics in the Euclidean regime [185, Conjecture 3.7].

It should be emphasized that $\mathbf{N}[\mathcal{F}]$ is generally _not_ GP
outside of the Euclidean regime. Indeed, the polytopes for three box diagrams
in Figures (10.3(a)), (10.3(b)) and (10.3(c)) are _not_ GP. The GP property
_is_ fulfilled for the polytope associated to the diagram (10.3(d)) though.

The story is further complicated by an interesting observation of Arkani-
Hamed, Hillman and Mizera [36, Section IV]: even for exceptional kinematics,
the Newton polytope of $\mathcal{F}$ is often equal to the base polytope
$\mathbf{P}[z_{\mathcal{F}}]$ with $z_{\mathcal{F}}$ given by (10.50)! In
fact, for all four examples in Figure (10.3) we have that
$\mathbf{N}[\mathcal{F}]=\mathbf{P}[z_{\mathcal{F}}].$ This is expected for
the polytope of the last graph (10.3(d)) since it is GP, but how is this
consistent with the fact that the first three cases (10.3(a)), (10.3(b)) and
(10.3(c)) are not GP? The answer is that $z_{\mathcal{F}}$ is not supermodular
for those graphs. So the function $z_{\mathcal{F}}$ does _not_ fulfil the
inequalities (10.46) for the first three boxes, but it _does_ so for the last
box.

(a) On-shell box.
Not GP.

(b) One off-shell leg.
Not GP.

(c) Two off-shell legs (adjacent). Not GP.

(d) Two off-shell legs (crossed). Is GP.

Figure 10.3:  1-loop box diagrams with massless internal propagators and
different combinations of on- or off-shell external legs. On-shell ($p^{2}=0$)
legs are depicted by dashed lines, and off-shell ($p^{2}\neq 0$) legs are
drawn as solid lines. In all four cases we have
$\mathbf{N}[\mathcal{F}]=\mathbf{P}[z_{\mathcal{F}}]$, but only (10.3(d)) is
GP.

While $\mathbf{N}[\mathcal{F}]=\mathbf{P}[z_{\mathcal{F}}]$ for generic
kinematics, we are not aware of a set of necessary and sufficient conditions
that specify when $\mathbf{N}[\mathcal{F}]$ equals
$\mathbf{P}[z_{\mathcal{F}}]$ for exceptional kinematics. It is at least the
case that $\mathbf{N}[\mathcal{F}]\subset\mathbf{P}[z_{\mathcal{F}}]$, because
we can only lose monomials in the exceptional case. A counter example999We
thank Erik Panzer for sharing this example with us. to equality is given by
the following 1-loop triangle:

$\displaystyle\mathcal{F}_{\text{triangle}}=\mathcal{F}\Big{(}\parbox[c]{0.0pt}{\includegraphics{Figs/1L_exceptional_triangle.pdf}}\Big{)}=m^{2}(x_{0}^{2}+x_{1}^{2})+(2m^{2}-p^{2})x_{0}x_{1}\,.$
(10.51)

The dashed line is massless, the double line has off-shell momentum $p^{2}\neq
0$, and every single solid line has mass $m$. This is an exceptional kinematic
configuration. The Newton polytope of $\mathcal{F}_{\text{triangle}}$ is the
convex hull of the points $(2,0),(0,2)$ and $(1,1)$. Since they all lie on a
line, the dimension of $\mathbf{N}[\mathcal{F}_{\text{triangle}}]$ equals 1.
However, the base polytope $\mathbf{P}[z_{\mathcal{F}_{\text{triangle}}}]$ can
be shown to have dimension 2. Hence
$\mathbf{P}[z_{\mathcal{F}_{\text{triangle}}}]\neq\mathbf{N}[\mathcal{F}_{\text{triangle}}].$

Why is this discussion relevant for tropical MC sampling? The answer is that
the sampling algorithm of Feyntrop assumes that $\mathbf{N}[\mathcal{F}]$ is
GP. Curiously, although we have highlighted the complications with this
assumption above, it appears not to matter that much in practice. Indeed, by
evaluating the examples from this section with Feyntrop (i.e. quasi-finite
versions of the four boxes in Figure (10.3) and the triangle from (10.51)), we
find agreement with other numerical codes in all cases. The catch is that
numerical convergence for Feyntrop becomes slower in the non-GP cases. We
suspect that Feyntrop will fail for more complicated diagrams with exceptional
kinematics, but we did not yet find such an example.

##### 10.3.2 Tropical sampling algorithm

| Euclidean | Minkowski
---|---|---
Generic | no deformation / always GP | deformation / always GP
Exceptional | no deformation / always GP | deformation / not always GP
Table 10.1:  Summary of whether contour deformation is required and whether
the GP property of $\mathbf{N}[\mathcal{F}]$ is fulfilled for each kinematic
regime.

Feyntrop employs a modification of [40, Algorithm 4] to sample from the
tropical measure $\mu^{\mathrm{tr}}$ given in (10.44), under the assumption
that $\mathbf{N}[\mathcal{F}]$ is GP. The algorithm includes a preprocessing
step, to be performed once, followed by a sampling step. The preprocessing
step scales as $\mathcal{O}(n2^{n})$, where $n=|E|$ is the number of edges in
a Feynman diagram. This scaling is determined from the preparation of a
certain table of rational numbers $J(\gamma)$ for each subgraph $\gamma\subset
E$. The sampling step scales as $\mathcal{O}\left(n^{2}\delta^{-2}\right)$,
where $\delta$ is the demanded relative accuracy of the result. This is just
ordinary MC scaling.

Because of the preprocessing step, the time complexity of this algorithm is
not fully polynomial. It is a very interesting question whether a polynomial
time algorithm exists for the numerical evaluation of FIs.

###### Preprocessing

The first task is to determine the kinematic regime, namely Euclidean vs.
Minkowski and generic vs. exceptional.

The regime is Euclidean if the scalar product matrix $P_{uv}$ is negative
semi-definite, and Minkowski otherwise. If the kinematics are in the Minkowski
regime, then contour deformation needs to be switched on.

If the kinematic point is exceptional and Minkowski, then
$\mathbf{N}[\mathcal{F}]$ might not be GP, and it might not even be equal to
$\mathbf{P}[z_{\mathcal{F}}]$. In this case, Feyntrop prints a warning saying
that the integration might give the wrong result (though, as mentioned, we
have not yet witnessed erroneous numerics). The program then bravely continues
under the assumption that
$\mathbf{N}[\mathcal{F}]=\mathbf{P}[z_{\mathcal{F}}]$. For generic kinematics,
regardless of whether we are in the Euclidean or Minkowski regimes, the
sampling algorithm is guaranteed to converge to the correct result [40,
Proposition 31].

A summary of the preceding paragraphs is given in Table (10.1). For brevity,
we have left out a technicality in that table regarding the Euclidean regime:
there is in fact a larger regime, containing the Euclidean one, called the
_pseudo-Euclidean regime_ ; it does not require contour deformation, but might
witness a failure of the GP property for exceptional kinematics. See [185,
Section 2.2] for details.

The next task is to prepare data associated to each of the $2^{|E|}$ subgraphs
$\gamma\subset E$:

  1. 1.

We calculate the loop number $L_{\gamma}=|\gamma|-|V_{\gamma}|+1,$ where
$V_{\gamma}$ is the vertex set associated to $\gamma$.

  2. 2.

We check whether or not $\gamma$ is mass-momentum spanning, i.e. whether
$\mathcal{F}_{G\setminus\gamma}=0$. By definition,
$\delta_{\gamma}^{\mathrm{m.m}}=1$ if $\gamma$ is mass-momentum spanning and
$\delta_{\gamma}^{\mathrm{m.m}}=0$ otherwise.

  3. 3.

We calculate the _generalized degree of divergence_ [40, Section 7.2]

$\displaystyle\omega(\gamma)=\sum_{e\in\gamma}\nu_{e}-\mathrm{D}_{0}L_{\gamma}/2-\omega_{0}\,\delta_{\gamma}^{\mathrm{m.m}}\,,$
(10.52)

where $\mathrm{D}_{0}$ is the integer part of the spacetime dimension and
$\omega_{0}=\sum_{e\in E}\nu_{e}-\mathrm{D}_{0}L/2.$ If $\omega(\gamma)\leq 0$
for any proper subgraph (i.e. not the whole graph $\gamma=E$), then $G$ has a
subdivergence. It implies that every term of the $\epsilon$-expansion (10.35)
is divergent, wherefore Feyntrop prints an error and terminates. In this case,
one has to shift the integrand exponents $(\nu,\mathrm{D})$ to obtain a quasi-
finite integral (cf. Section (10.2.3)).

  4. 4.

If $\omega(\gamma)>0$ for all $\gamma\subset E$, then we _recursively_ compute
a table of numbers $J(\gamma)$ via the formulas

$\displaystyle J(\emptyset)$ $\displaystyle=\omega(\emptyset)=1$
$\displaystyle J(\gamma)$
$\displaystyle=\sum_{e\in\gamma}\frac{J(\gamma\setminus
e)}{\omega(\gamma\setminus e)}\quad\text{for all}\quad\gamma\subset E\,,$
(10.53)

where $\gamma\setminus e$ means that we delete the edge $e$ from $\gamma$. The
last step of the recursion is for $\gamma=E$. In [40, Proposition 29],
Borinsky proved the surprising identity

$\displaystyle J(E)=I^{\mathrm{tr}}\,,$ (10.54)

where $I^{\mathrm{tr}}$ is the tropical normalization factor from (10.44). The
rational numbers $J(\gamma)$ are intimately related to the values
$z_{\mathcal{F}}(\gamma)$ [40, Section 6.1]. Indeed, this is the step during
preprocessing wherein the GP assumption is employed.

The data $(L_{\gamma},\,\delta_{\gamma}^{\mathrm{m.m}})$ feeds into the
computation of $z_{\mathcal{F}}(\gamma)$ via (10.50). For an exceptional
kinematic point in the Minkowski regime, Feyntrop checks whether
$z_{\mathcal{F}}$ is supermodular by testing a (more efficient) variant of the
inequalities (10.46), see [185, Equation (26)]. Feyntrop prints a message if
the inequalities are satisfied to notify that $\mathbf{P}[z_{\mathcal{F}}]$ is
GP, and a warning otherwise.

The output of this preprocessing is the data $\bigcup_{\gamma\subset
E}\big{\\{}L_{\gamma},\,\delta_{\gamma}^{\mathrm{m.m}},\,\omega(\gamma),\,J(\gamma)\big{\\}}\bigcup\big{\\{}I^{\mathrm{tr}}\big{\\}},$
which is stored in computer memory.

###### MC sampling

The validity of this tropical sampling algorithm was proven in [40,
Proposition 31].

Initialize the variables $\gamma=E$ and $\kappa=U=1$.

1:while $\gamma\neq\emptyset$ do

2: Pick a random edge $e\in\gamma$ with probability
$p_{e}^{\gamma}=\frac{1}{J(\gamma)}\frac{J(\gamma\setminus
e)}{\omega(\gamma\setminus e)}.$

3: Set $x_{e}=\kappa$.

4: If $\gamma$ is mass-momentum spanning but $\gamma\setminus e$ is not, set
$V=x_{e}$.

5: If $L_{\gamma\setminus e}<L_{\gamma}$, multiply $U$ with $x_{e}$ and store
the result in $U$, i.e. set $U\leftarrow x_{e}\cdot U$.

6: Remove the edge $e$ from $\gamma$, i.e. set
$\gamma\leftarrow\gamma\setminus e$.

7: Pick a uniformly distributed random number $\xi\in[0,1]$.

8: Multiply $\kappa$ with $\xi^{1/\omega(\gamma)}$ and store the result in
$\kappa$, i.e. set $\kappa\leftarrow\kappa\,\xi^{1/\omega(\gamma)}$.

9:end while

return $x=[x_{0},\ldots,x_{|E|-1}]\in\mathbb{P}^{|E|-1}_{+}$,
$\mathcal{U}^{\mathrm{tr}}(x)=U$ and $\mathcal{V}^{\mathrm{tr}}(x)=V$.

Algorithm 6 Generating a sample $x$ distributed as $\mu^{\mathrm{tr}}$ from
(10.44)

The main idea of this algorithm is to interpret $p_{e}^{\gamma}$ as a
probability distribution over the edges $e$ of a given subgraph $\gamma$.
Indeed, it follows from the definition of $J(\gamma)$ in (10.53) that
$p_{e}^{\gamma}\geq 0$ and $\sum_{e\in\gamma}p_{e}^{\gamma}=1.$ Starting from
$\gamma=E$, the algorithm proceeds to _cut_ an edge $e$ of the graph $G$ with