# Evaluating Feynman Integrals Using $\mathcal{D}$-modules and Tropical
Geometry

Henrik Jessen Munch

###### Abstract

Feynman integrals play a central role in the modern scattering amplitudes
research program. Advancing our methods for evaluating Feynman integrals will,
therefore, strengthen our ability to compare theoretical predictions with data
from particle accelerators such as the Large Hadron Collider. Motivated by
this, the present manuscript purports to study mathematical concepts related
to Feynman integrals. In particular, we present both numerical and analytical
algorithms for the evaluation of Feynman integrals.

The content is divided into three parts.

Part I focuses on the method of DEQs for evaluating Feynman integrals. An
otherwise daunting integral expression is thereby traded for the comparatively
simpler task of solving a system of DEQs. We use this technique to evaluate a
family of two-loop Feynman integrals of relevance for dark matter detection.

Part II situates the study of DEQs for Feynman integrals within the framework
of $\mathcal{D}$-_modules_ , a natural language for studying PDEs
algebraically. Special emphasis is put on a particular $\mathcal{D}$-module
called the _GKZ system_ , a set of higher-order PDEs that annihilate a
generalized version of a Feynman integral. In the course of matching the
generalized integral to a Feynman integral proper, we discover an algorithm
for evaluating the latter in terms of logarithmic series.

Part III develops a numerical integration algorithm. It combines Monte Carlo
sampling with _tropical geometry_ , a particular offspring of algebraic
geometry that studies "piecewise-linear" polynomials. Feynman’s
$i\varepsilon$-prescription is incorporated into the algorithm via contour
deformation. We present a new open-source program named Feyntrop that
implements this algorithm, and use it to numerically evaluate Feynman
integrals between 1-5 loops and 0-5 legs in physical regions of phase space.

monthyeardate[],

See pages - of Frontpage/frontpage.pdf

Author’s Declaration

I hereby declare that this thesis is solely authored by me. The work presented
herein was carried out by the author during his Ph.D. studies at Università
degli Studi di Padova, Dipartimento di Fisica e Astronomia “Galileo Galilei”.
The following manuscripts were co-authored by me during my Ph.D.

###### Articles

  1. 1.

Raghuveer Garani, Federico Gasparotto, Pierpaolo Mastrolia, Henrik J. Munch,
Sergio Palomares-Ruiz and Amedeo Primo, _Two-photon exchange in leptophilic
dark matter scenarios_ , JHEP 12 (2021), DOI: 10.1007/JHEP12(2021)212, arXiv:
2105.12116 [hep-ph].

  2. 2.

Vsevolod Chestnov, Federico Gasparotto, Manoj K. Mandal, Pierpaolo Mastrolia,
Saiei J. Matsubara-Heo, Henrik J. Munch and Nobuki Takayama, _Macaulay Matrix
for Feynman Integrals: Linear Relations and Intersection Numbers_ , JHEP 09
(2022), DOI: 10.1007/JHEP09(2022)187, arXiv: 2204.12983 [hep-th].

  3. 3.

Michael Borinsky, Henrik J. Munch and Felix Tellander, _Tropical Feynman
integration in the Minkowski regime_ , COMPHY 292 (2023), DOI:
https://doi.org/10.1016/
j.cpc.2023.108874, arXiv: 2302.08955 [hep-ph].

  4. 4.

Vsevolod Chestnov, Saiei J. Matsubara-Heo, Henrik J. Munch and Nobuki Takayama
_Restrictions of Pfaffian Systems for Feynman Integrals_ , arXiv: 2305.01585
[hep-th].

###### Conference proceedings

  1. 1.

Henrik J. Munch, _Feynman Integral Relations from GKZ Hypergeometric Systems_
, PoS, LL2022 (2022), DOI: 10.22323/1.416.0042, arXiV: 2207.09780 [hep-th].

  2. 2.

E. Bernardini, M. Carli, M. Y. Elkhashab, A. Ferroglia, M. Fiolhais, L.
Gabelli, H. J. Munch, D. Krym, P. Mastrolia, G. Ossola, O. Pantano, J.
Postiglione, J. S. Poveda Correa, C. Sirignano and F. Soramel, _Using
Smartphones to Innovate Laboratories in Introductory Physics Courses_ , GIREP
2022: Effective Learning in Physics from Contemporary Physics to Remote
Settings, Ljubljana, Slovenia.

Acknowledgements

> It takes a village to raise a child.
>
> African proverb

Let me here express my sincerest gratitude to all those who have so
wonderfully shaped my Ph.D. journey over the past three years.

I first and foremost thank my supervisor Pierpaolo Mastrolia for passionately
introducing me to the subject of scattering amplitudes, and for allowing me to
pursue my research goals.

It is not an exaggeration to say that most of what I learnt throughout my
Ph.D. came from interactions with Federico Gasparotto and Seva Chestnov. I
thoroughly enjoyed collaborating with you, and I cannot thank you enough for
sharing your knowledge with me.

Let me give further thanks to all the people associated to the Padova
Amplitudes group, with whom I had many stimulating discussions on topics in
physics, namely Manoj Kumar Mandal, Luca Mattiazzi, Giulio Crisanti, Giacomo
Brunello, Sid Smith, and Hjalte Frellesvig; I thank Hjalte in particular for
hosting me at the Niels Bohr Institute during the Amplitudes 2021 conference.
I also appreciate all the good times spent with Yousry, Marah, Prajwal and
Marco here in Padova.

I feel extremely fortunate to have worked with Felix Tellander and Michael
Borinsky, and I am so proud of the work we achieved together. I am grateful
for the hospitality shown to us three by Nima Arkani-Hamed, Sebastian Mizera
and Aaron Hillman during our visit to the Institute for Advanced Study in
Princeton.

Though I was officially enrolled in a physics Ph.D. program, my interests lead
me into more mathematical areas of research. In this regard, I had the
pleasure of collaborating with and learning from two brilliant mathematicians,
Saiei-Jaeyeong Matsubara-Heo and Nobuki Takayama. These interactions had a
tremendous impact on my intellectual development, and instilled in me a great
appreciation for the beauty of mathematical structures.

Speaking of mathematics, I want to thank everyone with whom I engaged with at
the Max Planck Institute for Mathematics in the Sciences in Leipzig. I am
especially grateful to Bernd Sturmfels for his hospitality, and for creating
such a unique and engaging environment at the institute. I thank my
collaborators Saiei and Claudia Fevola for all the joyous times spent there
diving into the mathematics of hyperplane arrangements and more.

I thank Jonah Stalknecht, Robert William Moerman and David Damgaard for co-
organizing the Geometry and Scattering Amplitudes Journal Club with me, and I
appreciate all the speakers and audience members for sharing their discoveries
within the field of scattering amplitudes.

Kæmpe tak til mine drenge Skov, Topic og Jeppe for en sindssyg fed tur rundt i
Italien. Det bragte mange gode minder.

Finally, I wish to thank my mother Ljudmila and my partner Laura for their
love and support throughout this whole period.

###### Contents

  1. List of Acronyms
  2. 1 Introduction
  3. I Multi-loop Technology
    1. 2 Feynman Integrals
      1. 2.1 Definition
        1. 2.1.1 Feynman diagram
        2. 2.1.2 The integrand
        3. 2.1.3 General features
      2. 2.2 Integration-by-parts identities
        1. 2.2.1 Irreducible scalar products
      3. 2.3 The method of differential equations
        1. 2.3.1 Canonical bases
        2. 2.3.2 Magnus expansion
      4. 2.4 Generalized polylogarithms
    2. 3 Two-loop Feynman Integrals for Dark Matter Detection
      1. 3.1 The model
      2. 3.2 The scattering amplitude
      3. 3.3 Computing the Feynman integrals
        1. 3.3.1 Two-mass case
        2. 3.3.2 Equal-mass case
        3. 3.3.3 Soft limit
      4. 3.4 Form factor results
        1. 3.4.1 Two-mass case
        2. 3.4.2 Equal-mass case
        3. 3.4.3 Soft limit
  4. II $\mathcal{D}$-modules
    1. 4 Introduction to $\mathcal{D}$-modules
      1. 4.1 What is a $\mathcal{D}$-module?
      2. 4.2 Holonomic $\mathcal{D}$-modules
        1. 4.2.1 Regular holonomic
      3. 4.3 Gröbner bases
      4. 4.4 Pfaffian systems
    2. 5 Macaulay Matrices
      1. 5.1 From Pfaffian to Macaulay matrix
      2. 5.2 From Macaulay to Pfaffian matrix
    3. 6 GKZ Systems
      1. 6.1 Euler integrals
      2. 6.2 GKZ $\mathcal{D}$-module
      3. 6.3 Twisted de Rham cohomology
        1. 6.3.1 "Traditional" de Rham cohomology
        2. 6.3.2 Twisted de Rham cohomology
        3. 6.3.3 Representing Euler integrals inside $\mathcal{D}$-modules
      4. 6.4 Basis of standard monomials
      5. 6.5 Integrand rescaling
      6. 6.6 Recurrence relations
        1. 6.6.1 Shifting $\beta$ by matrix factorials
        2. 6.6.2 General $\beta$ shift algorithm
    4. 7 Generalized Feynman Integrals
      1. 7.1 Lee-Pomeransky representation as a GKZ system
        1. 7.1.1 Generalized Feynman integral
      2. 7.2 Macaulay matrix example: 1-loop massless box
        1. 7.2.1 Basis
        2. 7.2.2 Macaulay matrix
      3. 7.3 Macaulay matrix example: 1-loop pentagon with off-shell leg
        1. 7.3.1 Basis
        2. 7.3.2 Macaulay matrix
      4. 7.4 Recurrence relation example: 1-loop bubble
    5. 8 Restriction of $\mathcal{D}$-modules
      1. 8.1 Preliminaries
        1. 8.1.1 Tensor products
        2. 8.1.2 Definition of restriction
        3. 8.1.3 Normal form
      2. 8.2 Restriction at the level of Pfaffian systems
        1. 8.2.1 Restriction of a logarithmic connection
        2. 8.2.2 Constructing the restricted Pfaffian system
        3. 8.2.3 Holomorphic restriction of solutions
        4. 8.2.4 Logarithmic restriction of solutions
      3. 8.3 Restriction at the level of Macaulay matrices
        1. 8.3.1 Restriction to a hyperplane
        2. 8.3.2 Restriction to a hypersurface
    6. 9 Restriction Examples
      1. 9.1 Holomorphic restriction: 2-loop N-box
        1. 9.1.1 Basis
        2. 9.1.2 Normal form
        3. 9.1.3 Restriction
      2. 9.2 Logarithmic restriction: 1-loop Bhabha scattering
        1. 9.2.1 Normal form
        2. 9.2.2 Eigenvalue $\lambda_{1}=0$
        3. 9.2.3 Eigenvalue $\lambda_{2}=-\epsilon$
        4. 9.2.4 Result for the logarithmic series
      3. 9.3 Logarithmic restriction: 2-loop Bhabha scattering
        1. 9.3.1 Normal form
        2. 9.3.2 Restriction
        3. 9.3.3 Discussion
      4. 9.4 Macaulay-level restriction: 2-loop diagonal box
        1. 9.4.1 Setup
        2. 9.4.2 Restriction
  5. III Numerical Integration
    1. 10 Tropical Integration
      1. 10.1 Tropical Monte Carlo sampling, in brief
        1. 10.1.1 Tropical approximation
      2. 10.2 Projective Feynman integrals
        1. 10.2.1 $i\varepsilon$-prescription via contour deformation
        2. 10.2.2 $\epsilon$-expansion for quasi-finite integrals
        3. 10.2.3 Convergence of Euler integrals
      3. 10.3 Tropical Monte Carlo sampling, in detail
        1. 10.3.1 Generalized permutahedra
        2. 10.3.2 Tropical sampling algorithm
        3. 10.3.3 Evaluating deformed Feynman integrands
    2. 11 The Program feyntrop
      1. 11.1 Basic usage of Feyntrop
        1. 11.1.1 Tutorial
      2. 11.2 Examples of Feynman integral evaluations
        1. 11.2.1 5-loop 2-point zigzag diagram
        2. 11.2.2 3-loop 4-point envelope diagram
        3. 11.2.3 2-loop 4-point $\mu e$-scattering diagram
        4. 11.2.4 2-loop 5-point QCD-like diagram
        5. 11.2.5 2-loop 5-point Higgs production diagram
        6. 11.2.6 4-loop 0-point QED-like diagram
        7. 11.2.7 An elliptic, conformal, 4-point integral
    3. 12 Conclusion and Outlook
    4. A Details on $\mathcal{D}$-modules
      1. A.1 Basis change in a holonomic $\mathcal{D}$-module
      2. A.2 Basis change without derivatives
      3. A.3 General formula for homogeneity rescaling
      4. A.4 GKZ systems are isomorphic to twisted cohomology groups

### List of Acronyms

(CAS): | Computer algebra system
---|---
(DM): | Dark matter
(DR): | Dimensional regularization
(FI): | Feynman integral
(GFI): | Generalized Feynman integral
(GKZ): | Gel’fand Kapranov Zelevinsky
(GP): | Generalized permutahedron
(GPL): | Generalized polylogarithm
(HPL): | Harmonic polylogarithm
(MC): | Monte Carlo
(ISP): | Irreducible scalar product
(LPr): | Lee-Pomeransky representation
(MI): | Master integral
(SM): | Standard model
(QCD): | Quantum chromodynamics
(QED): | Quantum electrodynamics
(QFT): | Quantum field theory

### Chapter 1Introduction

Eugene Wigner famously pondered about the _unreasonable effectiveness of
mathematics in the natural sciences_ [1]. How can it be that esoteric
mathematical structures such as group theory, calculus and geometry - often
developed for entirely independent reasons - are able to describe the "laws of
Nature" so powerfully? A prime example is that of quantum field theory (QFT).
This theory combines two pillars of our modern understanding of Nature:
quantum mechanics and special relativity. Although some intuition can be
gained from describing QFT in plain English, it is only in terms of
mathematics that we can properly (insofar as possible) _define_ the theory and
_do_ something with it!

The standard model (SM) is the crown jewel of the modern QFT program. It
describes three out of four known forces of Nature: electromagnetism, the weak
and the strong nuclear forces. Although it does not incorporate gravity (the
last, known force) at the quantum level, the mathematical techniques that were
developed to probe the SM are now also being used to study gravity in the weak
field limit [2].

Physicists test the SM in the same way Rutherford first investigated the sub-
atomic world a century ago: by scattering particles together and measuring
their cross sections (roughly speaking the effective range of interaction of
the scattered particles). Figure (1.1) compares a host of experimentally
determined total cross sections, from ATLAS data, with corresponding
theoretical predictions made within the SM. The level of agreement is
striking.

A QFT, such as the SM, is typically defined in the abstract language of
Lagrangians, i.e. polynomials in quantum fields and derivatives thereof. How
can this be connected to concrete, observable quantities such as the
aforementioned cross sections? The connection is made through the calculation
of _scattering amplitudes_. These are complex-valued functions (depending on
data such the masses, spins and momenta of the particles involved in a given
scattering experiment) which measure the quantum mechanical probability for
turning an initial state into a final state. Given such an amplitude
$\mathcal{A}$, one obtains a total cross section by integrating
$|\mathcal{A}|^{2}$ over a relevant region of phase space.

Because the SM Lagrangian contains several small parameters in the form of
coupling constants, it is natural to compute scattering amplitudes as
_perturbative_ expansions in said constants. The numerical accuracy of a given
amplitude is then determined by how many orders of this expansion we are able
to compute. Ever since the work of Feynman [3], physicists have organized
these expansions in terms of _Feynman diagrams_.

Figure 1.1:  Summary plot anno 2019 comparing experimental vs. theoretical
values for Standard Model cross sections [4].

Schematically, the perturbative expansion of an amplitude $\mathcal{A}$ is
written as

(1.1)

The _Feynman rules_ , derived from the Lagrangian of a given theory, give
instructions for how to calculate each diagram. The leading-order
contributions, written as sums of _tree_ diagrams, are, arguably, now under
control, and so will not be a focus of this thesis. The subleading-order
contributions are written as sums of _loop_ diagrams. In these cases, the
Feynman rules instruct us to compute integrals over the momenta of virtual
particles running inside the loops. These _Feynman integrals_ (FIs) are the
source of nearly all the computational complexity that goes into the
calculation of perturbative scattering amplitudes beyond leading order. If we
want to probe the SM at increased accuracy in the hopes of detecting signs of
physics _beyond_ the SM (which we know must exist, due to unexplained
phenomena such as neutrino masses, dark matter and the supposedly "unnatural"
value of the Higgs mass), it thus behooves us to strengthen our mathematical
capabilities for computing FIs. Indeed, following Wigner, we should use the
unreasonable effectiveness of mathematics to our benefit!

#### The context of this thesis

There is a very long road from writing down a Lagrangian in, say, QCD to
actually arriving at a cross section at some given order in perturbation
theory (this includes writing down all the relevant Feynman diagrams, taking
care of tensor structure, renormalizing and regularizing UV and IR divergences
respectively, feeding in non-perturbative information via parton distribution
functions, choosing proper values for factorization/renormalization scales,
and much much more). The modern scattering amplitudes program is, in a word,
all about finding potential _simplicity_ hiding within this highly complex
process in order to lessen the calculational burden [5]111 In fact, it is now
quite common in this field to dispense with Lagrangians entirely, as they
carry so much physically irrelevant (i.e. off-shell) information! . More
narrowly, for the purpose of this thesis we may ponder the question: _do FIs
also exhibit any simple mathematical structure?_

###### IBPs and DEQs.

The first sign of "simplicity" was arguably the discovery of the _linear
relations_ among FIs called _integration-by-parts_ identities (IBPs) [6, 7,
8]. To mention a concrete application of these relations, take the recent
3-loop QCD amplitude computed in [9]. The original expression for the
amplitude contains of order $\mathcal{O}(10^{5})$ FIs to be evaluated.
However, only $\mathcal{O}(300)$ of these integrals are actually linearly
independent! In other words, the original $\mathcal{O}(10^{5})$ integrals can
all be written as linear combinations of $\mathcal{O}(300)$ _master integrals_
(MIs), where the coefficients in front of the MIs are rational functions in
the kinematic variables and the dimension of spacetime. When such linear
relations can be generated systematically and efficiently, it constitutes a
significant _simplification_ of a given calculation. There are now many public
codes dedicated to this purpose being used in virtually every current
amplitude calculation beyond leading order [10, 11, 12, 13, 14, 15]. Indeed,
we shall make frequent use of them in this thesis.

Suppose the many FIs in a given calculation have been reduced to a small set
of MIs. The question is then how to evaluate the MIs themselves. IBPs once
again lend a helping hand, as they can be used to derive _differential
equations_ for MIs [16, 17, 18, 19]. The second example of "simplicity" that
we highlight was discovered by Henn [20]: there is a freedom in choosing which
FIs that are promoted to MIs, and if this choice is made judiciously, then
their DEQs simplify considerably. Concretely, the simplification stems from
the decoupling, or factorization, of a parameter in the DEQs called
$\epsilon$, which appears because it is related to the dimension $\mathrm{D}$
of spacetime (as mentioned above, $\mathrm{D}$ is brought into the game due to
IBPs). Such a basis is dubbed to be _canonical_. When the DEQs are brought
into this canonical form, they are _much_ easier to solve.

######  $\mathcal{D}$-modules.

These two notions, IBPs and DEQs, thus appear to be central to the
(analytical, or symbolic) evaluation of FIs. To learn more about the
mathematical structure of these integrals, we argue that it is interesting to
reformulate these two notions in novel mathematical frameworks222 A prime
example hereof is the recent reformulation of IBPs within the framework of
_twisted cohomology_ [21]. , as new perspectives often tend to bring new
insights. One such mathematical framework, which will occupy a large portion
of this thesis, is that of $\mathcal{D}$-_modules_. $\mathcal{D}$-modules
study DEQs from an _algebraic_ point of view, as a supplement to the usual
tools from complex analysis. Its mathematical development goes back to the
Japanese school of the 70s and 80s spearheaded by M. Sato and M. Kashiwara.
(See [22] for a historical account and the textbooks [23, 24] for technical
details. Early studies on the relationship between FIs and
$\mathcal{D}$-modules include [25, 26] in the 1970s, and later [27, 28, 29] in
the 2010s.)

One of the goals of this thesis is thus to re-interpret IBPs and DEQs in this
language. To foreshadow the contents of future chapters, we shall interpret
IBPs in terms of so-called _Macaulay matrices_ , and solutions to DEQs via
_"restrictions" of $\mathcal{D}$-modules_.

###### Numerical integration.

Alongside the aforementioned developments in the analytical evaluation of FIs,
there is also a tradition for the _numerical_ evaluation thereof (see [30,
Section 3.3] for a historical account). However, because of limited computer
resources in the previous century, this tradition only goes back roughly two
decades (some of the earliest developments include [31, 32]). But why go down
this route? Although analytic control over FIs is arguably the most favorable
situation, sometimes the symbolic expressions are just too complicated even
with current tools! Nevertheless, blind numerical evaluation is rarely the
most efficient way to proceed. It is better to perform some kind of symbolic
_preprocessing_ , so as to reveal an expression that is simple to numerically
evaluate. For instance, one could still use IBPs to derive DEQs for FIs, but
now with every parameter except one fixed to a number in order to simplify the
IBPs. The DEQs are then written w.r.t. this single, unfixed parameter, and are
solved through numerical series expansions rather than the analytic method of
canonical bases. This approach, called _auxiliary mass flow_ , has proven to
be very fruitful in recent years [33, 34].

###### Tropical geometry.

With the goal of numerical evaluation in mind, we shall carry out such a
symbolic preprocessing by stepping into the world of _tropical geometry_. The
mathematical discipline of _algebraic geometry_ is perhaps already familiar to
the reader. The latter studies the "smooth" geometries associated to the
solutions of polynomial equations. Tropical geometry, instead, studies the
solutions to equations involving "piecewise-linear" polynomials. This
"piecewise" nature then translates into more "chunky" geometric figures, more
precisely _polytopes_ , which are comparatively easier to handle. (See the
textbook [35] for a pedagogical exposition on tropical geometry.) The
connection to FIs comes from representing them as integrals over certain
polynomials raised to rational powers. By studying the simpler, tropical
nature of these integrand polynomials, one can learn about the UV/IR divergent
behaviour of FIs [36], as well as the connection between QFT and string theory
amplitudes [37, 38].

For our purposes, the central notion will be that of a _tropical
approximation_ to a polynomial, envisioned by Panzer [39], which essentially
just replaces $"\\!+\\!"$ with $"\mathrm{max}"$. We obtain tropical versions
of FIs by applying this operation to the integrands. Despite their
_simplicity_ (or perhaps because of it!), tropical versions of FIs appear to
retain an uncanny amount of information about the original integral [39,
Section 5.1]. We shall employ this particular symbolic preprocessing, i.e. the
tropicalization of the Feynman integrand, to simplify the numerical evaluation
of FIs. This idea was originally put forward by Borinsky [40].

#### The content of this thesis

The above discussion hopefully motivates the main topic of this thesis, namely
the evaluation of FIs. We have divided the presentation into three parts.

###### Part I.

The first part of the thesis focuses on the evaluation of FIs through the
method of DEQs.

In Chapter (2), starting from a single FI, we show how to associate to it a
whole family of integrals $I_{\nu_{1}\ldots\nu_{n}}$, where each
$\nu_{i}\in\mathbb{Z}$. There exist linear relations among integrals with
different values of $(\nu_{1},\ldots,\nu_{n})$, namely the IBPs mentioned
above. We shall describe how to derive these relations, how to identify a set
of MIs $\vec{I}$, and lastly how set up a system of first-order DEQs

$\displaystyle\frac{\partial}{\partial
z_{i}}\vec{I}(z)=P_{i}(z)\cdot\vec{I}(z)\quad\text{for}\quad i=1,\ldots,N\,.$
(1.2)

$z=(z_{1},\ldots,z_{N})$ is a collection of kinematic variables, and
$P_{i}(z)$ are rational matrices. This is called a _Pfaffian system_ of DEQs.
Direct integration of FIs, often being too hard, is thereby sidestepped in
favor of the comparatively easier task of solving Pfaffian systems.

In Chapter (3), we use the method of IBPs and Pfaffian systems to calculate a
collection of 2-loop 3-point FIs with internally massive propagators. These
integrals are then fed into the form factors for a scattering amplitude that
is of relevance to dark matter detection.

###### Part II.

The second part of this thesis is much more mathematically formal (though
hopefully still accessible to a physicist audience), as it dives into the
theory of $\mathcal{D}$-_modules_.

In Chapter (4), we begin by describing those concepts in the theory of
$\mathcal{D}$-modules that are needed for developments in future chapters.
This includes the notion of a _holonomic_ $\mathcal{D}$-module, which carries
the structure of a finite-dimensional vector space. This mimics the vector
space structure of families of FIs due to IBPs.

In Chapter (5), we present a novel algorithm whose output is a Pfaffian system
for a holonomic $\mathcal{D}$-module. It takes as input a collection of
differential operators that annihilate a given function as well a basis of
operators for the associated vector space. The key concept of the algorithm is
that of a _Macaulay matrix_ , which encodes relations among differential
operators. This matrix is frequently used in the setting of commutative
polynomial rings to avoid Gröbner basis computations, as they tend to scale
poorly with the complexity of the problem. Indeed, one feature of the
algorithm presented here is that it too avoids the need for Gröbner bases, now
in a non-commutative ring of differential operators.

In Chapter (6), we draw our attention to a specific holonomic
$\mathcal{D}$-module called the _GKZ system_ , named after Gel’fand, Kapranov
and Zelevinsky [41, 42]. It is built from a collection of differential
operators that can immediately be written down given the integral
representation of an _Euler integral_. The benefit from employing the GKZ
system is two-fold: 1) a complete set of annihilating operators is known, and
2) an explicit vector space basis can be computed very fast.

In Chapter (7), the connection between Euler integrals and FIs is established
by generalizing the parameters present in the _Lee-Pomeransky representation_
[43] of a FI. The GKZ system can consequently be used to study this
_generalized_ FI. In particular, we use the GKZ system as input for the
Macaulay matrix algorithm in order to derive Pfaffian systems for a couple of
simple FIs.

In Chapter (8), we address the discrepancy between Euler integrals and FIs: an
Euler integral often has too many parameters compared to the number of
kinematic variables in a FI. To land on the $\mathcal{D}$-module for the FI
proper, one must take a kind of "limit" of the GKZ system, called a
_restriction_ in the language of $\mathcal{D}$-modules. We develop two
different methods for computing restrictions whose output is the Pfaffian
system for the restricted $\mathcal{D}$-module. The first method is based on
manipulating Pfaffian systems through series expansions and gauge
transformations. The second method simply takes the limits directly at the
level of the Macaulay matrix. The first method turns out to be applicable to
the computation of logarithmic series solutions to Pfaffian systems. The
second method turns out to allow for restrictions onto hyper _surfaces_ in
addition to just hyper _planes_ , which, to the best of our knowledge, has so
far not been achieved within the theory of $\mathcal{D}$-modules.

In Chapter (9), we showcase both restriction algorithms through several FI
examples. The examples concerning logarithmic solutions to Pfaffian systems
unfortunately only contain partial results, and so will have to be completed
in future work.

###### Part III.

The third and final part of this thesis shifts focus from analytical to
numerical methods for evaluating FIs. The mathematical framework is that of
_tropical geometry_.

In Chapter (10), we begin by describing the notion of a _tropical
approximation_ to a polynomial. For instance, given

$\displaystyle p(x)=2x_{1}x_{2}+3x_{1}^{2}-8x_{2}^{2}\,,$ (1.3)

then its tropical approximation is

$\displaystyle
p^{\mathrm{tr}}(x)=\mathrm{max}\\{x_{1}x_{2},\,x_{1}^{2},\,x_{2}^{2}\\}\,.$
(1.4)

We thereafter present a striking observation made by Borinsky [40]: for an
integral of a polynomial raised to some power $\int p(x)^{a}\mathrm{d}^{n}x$
(over a certain integration contour and modulo some technical conditions), one
can build a _probability measure_ from $p^{\mathrm{tr}}(x)$. By sampling
points from this measure, it is thereby possible to perform numerical Monte
Carlo integration. This can be used to evaluate FIs in parametric
representation. We combine this _tropical integration_ algorithm with the
$\epsilon$-expansion of FIs in the dimensional regularization scheme.
Feynman’s $i\varepsilon$-prescription is incorporated via a certain contour
deformation $x\to xe^{if(x)}$, which maps the Schwinger integration variables
$x$ from real line into a path in the complex plane.

In Chapter (11), the program Feyntrop is introduced. It is an implementation
of the tropical Monte Carlo integration scheme, and includes both the
$\epsilon$-expansion and the $i\varepsilon$-prescription. The program can
therefore numerically evaluate dimensionally regulated FIs for physical
kinematics in Minkowski space (rather than just Euclidean kinematics). After a
tutorial of the program, we then apply it to several non-trivial examples
between 1-5 loops and 0-5 legs.

In Chapter (12), we present the conclusions reached throughout this thesis and
speculate on developments for the future.

##  Part I Multi-loop Technology

### Chapter 2Feynman Integrals

This chapter is a summary of basic concepts pertaining to Feynman integrals.
It is far from an exhaustive review (for such a text, consult e.g. [44] and
[45]).

####  Section 2.1Definition

The main object of study in this thesis is the _Feynman integral_ (FI)
associated to an $L$-loop Feynman diagram $G$:

$\displaystyle I_{G}=\lim_{\varepsilon\to
0^{+}}\int_{(\mathbb{R}^{1,\mathrm{D}-1})^{L}}\,\mathrm{d}^{L\cdot\mathrm{D}}\ell\,\prod_{e=1}^{E_{\text{int}}}\frac{1}{\left[-q_{e}(p,\ell)^{2}+m_{e}^{2}-i\varepsilon\right]^{\nu_{e}}}\,.$
(2.1)

We often also write $I_{\nu_{1}\nu_{2}\ldots}=I(\nu_{1},\nu_{2},\ldots)$. Let
us go through the notation used above.

##### 2.1.1 Feynman diagram

A (scalar) _Feynman diagram_ $G=(E,V)$ is defined by a set of edges
$E=\\{e_{1},\ldots,e_{|E|}\\}$ and vertices $V=\\{v_{1},\ldots,v_{|V|}\\}.$
Given the number of edges and vertices, the number of loops is given by
$L=|E|-|V|+1\,.$ To each edge $e$, we associate a $\mathrm{D}$-dimensional
momentum vector flowing along the orientation of $e$. Momentum conservation is
imposed at each vertex $v$, meaning that all edge momenta flowing into $v$
equal those flowing out of $v$ (in analogy with Kirchhoff’s law for currents
in electrical circuits).

##### 2.1.2 The integrand

Let us partition the set of edges into _external_ and _internal_ ones, i.e.
$E=E_{\text{ext}}\,\sqcup\,E_{\text{int}}$. The external edges have one end
which is not connected to any other edge, while the internal edges are
connected to at least one other edge on each end.

Consider an edge $e$. If $e\in E_{\text{ext}}$ is external, then we label its
associated momentum vector by

$\displaystyle
p_{e}=\bigl{(}p_{e}^{0},\ldots,p_{e}^{\mathrm{D}-1}\bigr{)}\in\mathbb{R}^{1,\mathrm{D}-1}\,,$
(2.2)

where $\mathbb{R}^{1,\mathrm{D}-1}$ denotes $\mathrm{D}$-dimensional Minkowski
space equipped with the inner product

$\displaystyle a\cdot
b=a^{0}b^{0}-a^{1}b^{1}-\ldots-a^{\mathrm{D}-1}b^{\mathrm{D}-1}\,.$ (2.3)

On the other hand, if $e\in E_{\text{int}}$ is internal, then we associate to
it the vector $q_{e}\in\mathbb{R}^{1,\mathrm{D}-1}$ written as a sum of
external momenta $p_{i}$ and _loop momenta_ $\ell_{j}$:

$\displaystyle q_{e}(p,\ell)=\sum_{i}\pm\,p_{i}\hskip 5.69046pt+\hskip
5.69046pt\sum_{j}\pm\,\ell_{j}\,.$ (2.4)

The signs and summation ranges for $(i,j)$ depend on the chosen _momentum
routing_ for the graph $G$. The $L$ loop momenta
$\\{\ell_{1},\ldots,\ell_{L}\\}$ can be viewed as the leftover degrees of
freedom after imposing momentum conservation at each vertex. In (2.1), these
loop momenta are integrated against the measure

$\displaystyle\mathrm{d}^{L\cdot\mathrm{D}}\ell:=\mathrm{d}\ell_{1}^{(0)}\wedge\cdots\wedge\mathrm{d}\ell_{1}^{(\mathrm{D}-1)}\wedge\cdots\wedge\mathrm{d}\ell_{L}^{(0)}\wedge\cdots\wedge\mathrm{d}\ell_{L}^{(\mathrm{D}-1)}$
(2.5)

over the whole of Minkowski space $\mathbb{R}^{1,\mathrm{D}-1}$.

To any edge $e$, we additionally assign a scalar mass
$m_{e}\in\mathbb{R}_{\geq 0}.$ For an external momentum, we say that it is
_on-shell_ if it satisfies Einstein’s energy momentum relation:
$p_{e}^{2}=m_{e}^{2}$. An internal momentum is generally not on-shell, that is
$q_{e}^{2}\neq m_{e}^{2}$ (apart from special points within the integration
range).

The integrand of (2.1) is constructed by taking a product of _Feynman
propagators_

$\displaystyle\frac{1}{\left[-q_{e}(p,\ell)^{2}+m_{e}^{2}-i\varepsilon\right]^{\nu_{e}}}$
(2.6)

for each internal edge. The integers $\nu_{e}\in\mathbb{Z}$ are called _edge
weights_ or simply _propagator powers_. The infinitesimal imaginary part
$i\varepsilon$ is introduced to set a particular integration contour. We shall
have more to say about this in Section (10.2.1), but until then we can safely
leave the $i\varepsilon$ as well as the limit $\varepsilon\to 0^{+}$ in front
of (2.1) as implicit. In gauge theories such as QED and QCD, Feynman
propagators become dressed with numerators having some theory-dependent tensor
structure. This additional structure is often factored out in front of scalar
FIs (cf. Section (3.2)).

Example 2.1. Consider the Feynman diagram

$\displaystyle
G\qquad=\parbox[c]{0.0pt}{\includegraphics{Figs/introduction_example.pdf}}\,.$
(2.7)

It has $L=2$ loops by inspection. As a sanity check, we get the same value by
writing $L=|E_{\text{int}}|-|V_{\text{int}}|+1=5-4+1=2,$ where we just used
the internal edges and vertices for counting (the external lines add three
edges but technically also three vertices, so these cancel in the formula for
$L$).

We have chosen an arbitrary momentum routing by drawing arrows on each edge.
The external momenta are denoted by $\\{p_{1},p_{2},p_{3}\\}$, and the
internal ones are $\\{q_{1},\ldots,q_{5}\\}$. Imposing momentum conservation
at each vertex, we have

$\displaystyle\begin{array}[]{lll}&\text{vertex
$1$:}&p_{1}+q_{1}+q_{5}=q_{4}\\\ &\text{vertex $2$:}&p_{2}+q_{2}=q_{1}\\\
&\text{vertex $3$:}&p_{3}+q_{4}=q_{3}\\\ &\text{vertex
$4$:}&q_{3}=q_{2}+q_{5}\,.\end{array}$ (2.12)

There are not enough relations here to completely fix all of the momenta
$q_{e}$. Given that $L=2$, there will be two leftover degrees of freedom
$\ell_{1}$ and $\ell_{2}$. Arbitrarily choosing to call $q_{1}=\ell_{1}$ and
$q_{5}=\ell_{2}$, then a solution to the linear system above is

$\displaystyle\begin{array}[]{lll}q_{1}&=&\ell_{1}\\\
q_{2}&=&\ell_{1}-p_{2}\\\ q_{3}&=&\ell_{1}+\ell_{2}-p_{2}\\\
q_{4}&=&\ell_{1}+\ell_{2}-p_{2}-p_{3}\\\ q_{5}&=&\ell_{2}\,.\end{array}$
(2.18)

Suppose the internal masses are given by
$\\{m_{1},m_{2},m_{3},m_{4},m_{5}\\}=\\{m,0,0,0,m\\},$ and that the edge
weights are $\\{\nu_{1},\nu_{2},\nu_{3},\nu_{4},\nu_{5}\\}=\\{1,2,1,1,3\\}.$

The FI associated to this data is then

$\displaystyle
I_{G}=\int_{\left(\mathbb{R}^{(1,\mathrm{D}-1)}\right)^{2}}\frac{\mathrm{d}^{2\cdot\mathrm{D}}\ell}{\left[-q_{1}^{2}+m^{2}\right]\left[-q_{2}^{2}\right]^{2}\left[-q_{3}^{2}\right]\left[-q_{4}^{2}\right]\left[-q_{5}^{2}+m^{2}\right]^{3}}\,,$
(2.19)

where the $q_{e}$ are given by (2.18). $\blacksquare$

##### 2.1.3 General features

Having defined our object of study, let us now describe some general features
of FIs. The first thing to mention is that they are often _divergent_ when the
loop momenta are small and/or large, i.e. when $\ell\to 0$ and/or
$|\ell|\to\infty$. This should not deter us, since we can introduce a
regularization scheme to make these integrals well-defined. The most common
scheme, and the one we adopt throughout this thesis, is that of _dimensional
regularization_ (DR) [46]. Although (2.1) will generically diverge in, say,
$\mathrm{D}=4$ spacetime dimensions, we can promote $\mathrm{D}$ to an
indeterminate parameter, thereby considering $I_{G}=I_{G}(\mathrm{D})$ as a
function of $\mathrm{D}$. Divergences then show up as poles in
$(\mathrm{D}-4)$, leading us to interpret the FI a Laurent series in
$(\mathrm{D}-4)$. It is customary to write

$\displaystyle\mathrm{D}=\mathrm{D}_{0}-2\epsilon$ (2.20)

with $\mathrm{D}_{0}\in\mathbb{Z}_{>0}$ and $|\epsilon|\ll 1$, such that

$\displaystyle
I_{G}(\mathrm{D})=\sum_{i=i_{\text{min}}}^{\infty}I_{G}^{(i)}\epsilon^{i}\,,\quad
i_{\text{min}}\in\mathbb{Z}\,.$ (2.21)

The poles in $\epsilon$ happily cancel during the computations of physical
observables. There are two mechanisms for these cancellations: UV
renormalization and IR subtraction, but we shall not dive into these deep
topics here. Note that the mass dimension of the integral measure
$\left[\mathrm{d}^{\mathrm{D}\cdot L}\ell\right]=\mathrm{D}\cdot L$ changes in
DR. To counterbalance this, one introduces another arbitrary parameter $\mu$
with mass dimension $[\mu]=1$, leading to a new measure $(\mu^{2})^{\epsilon
L}\,\mathrm{d}^{\mathrm{D}\cdot L}\ell.$ We will often be sloppy with factors
of $\mu$ throughout this thesis, since they are easy to reinstate.

The end goal of a FI computation is to obtain the coefficient functions
$I_{G}^{(i)}$. What do they look like, qualitatively speaking? Inspecting
(2.1), observe that $I_{G}^{(i)}$ will be a functions of masses $m_{e}$ and
$(|E_{\text{ext}}|-1)/2$ scalar products $p_{e}\cdot p_{d}$. The factor of
$1/2$ is due to $p_{e}\cdot p_{d}=p_{d}\cdot p_{e}$, and the factor
$(|E_{\text{ext}}|-1)$ comes from momentum conservation

$\displaystyle p_{|\mathrm{Ext}|}=-\sum_{e=1}^{|E_{\text{ext}}|-1}p_{e}\,,$
(2.22)

with the convention that all external momenta are in-coming. In practical
calculations, it is useful to trade the dependence on scalar products for a
dependence on _Mandelstam variables_ $s_{ij}=(p_{i}+p_{j})^{2}.$

So the $I_{G}^{(i)}$ are functions of masses and external momenta, but what
kind of functions? The answer to this question is not yet settled in full
generality, though much is known at low loop-order [47]. Part of the answer
depends on what properties we would _like_ a space of functions to have. Take,
for instance, the famous $\Gamma$-function

$\displaystyle\Gamma(z)=\int_{0}^{\infty}x^{z-1}e^{-x}\mathrm{d}x\,,\quad\mathrm{Re}(z)>0\,.$
(2.23)

Because $\Gamma$ is _defined_ by its integral representation, we don’t think
of that integral as being "leftover work" to be done, as opposed to e.g.
$\int_{0}^{z}\,x^{2}\,\mathrm{d}x.$ The integral representation of $\Gamma$ is
perfectly fine because

  1. 1.

It allows for the derivation of many identities, such as the famous reflection
formula $\Gamma(1-z)\Gamma(z)=\frac{\pi}{\sin(\pi z)},$ which can simplify
calculations involving large expressions. For instance, there could be two
terms like $4\Gamma(1-z)\Gamma(z)-4\frac{\pi}{\sin(\pi z)}$ in some result,
but it would be wasteful to numerically evaluate this expression since it
gives $0$.

  2. 2.

It opens the door to other representations which allow for accurate numerical
evaluations, such as the Stirling series $\Gamma(1+z)=\sqrt{2\pi
z}\left(z/e\right)^{z}\exp\left(\sum_{i=1}^{\infty}\frac{B_{2i}}{2i(2i-1)z^{2i-1}}\right)$
written in terms of Bernoulli numbers $B_{i}$.

  3. 3.

It is a starting point for determining the analytic continuation of
$\Gamma(z)$ to larger regions of the complex plane.

Based on this analogy, we claim that "calculating" the $i$th
$\epsilon$-coefficient function $I_{G}^{(i)}$ is equivalent to rewriting it in
terms of functions with "nice" numerical and analytic properties. A natural
choice for such "nice" functions turns out to be the _Chen iterated integrals_
which we define in Section (2.4).

####  Section 2.2Integration-by-parts identities

In state-of-the-art computations at 2- or 3-loop orders, one typically
encounters $\mathcal{O}(10^{4})$ to $\mathcal{O}(10^{6})$ different FIs.
Fortunately, only a select few of these integrals, typically of order
$\mathcal{O}(10^{2})$, need to be calculated. This huge decrease in workload
is a consequence of linear relations among FIs called _integration-by-parts
identities_ (IBPs) [6, 7, 8]. These identities express the many integrals
appearing in a QFT computation in terms of a small basis of _master integrals_
(MIs). The MIs can then be computed via the method of differential equations
[16, 17, 18, 19], as outlined in Section (2.3).

IBPs are born from the fact that a total derivative vanishes in DR [44,
Section 6.1]:

$\displaystyle 0=\int\mathrm{d}^{\mathrm{D}\cdot
L}\ell\,\frac{\partial}{\partial\ell^{\mu}}\bigl{[}v^{\mu}\cdot
f(\ell)\bigr{]}\,,$ (2.24)

where $v^{\mu}$ is an arbitrary vector and $f(\ell)$ is some product of
propagators. Let us go through an example to see why (2.24) leads to linear
relations among FIs.

Example 2.1. Consider the following 1-loop bubble integral family111 The
solid, thick line has mass $m$, the solid, thin line is massless, and the
dashed lines have squared momentum $p^{2}$.

$\displaystyle I(\nu_{1},\nu_{2})\ =\
\parbox[c]{0.0pt}{\includegraphics[scale={0.15}]{Figs/MM_bubble.pdf}}\ =\
\int\frac{\mathrm{d}^{\mathrm{D}}\ell}{\left[-\ell^{2}+m^{2}\right]^{\nu_{1}}\left[-(\ell+p)^{2}\right]^{\nu_{2}}}\,.$
(2.25)

Say we choose $v^{\mu}=p^{\mu}$ in (2.24). Setting $D_{1}=-\ell^{2}+m^{2}$ and
$D_{2}=-(\ell+p)^{2},$ then the chain and product rules yield

$\displaystyle\frac{\partial}{\partial\ell^{\mu}}\,\frac{p^{\mu}}{D_{1}^{\nu_{1}}D_{2}^{\nu_{2}}}$
$\displaystyle=\frac{p^{\mu}}{D_{2}^{\nu_{2}}}\,\frac{\partial}{\partial\ell^{\mu}}\,\frac{1}{D_{1}^{\nu_{1}}}\hskip
5.69046pt+\hskip
5.69046pt\frac{p^{\mu}}{D_{1}^{\nu_{1}}}\,\frac{\partial}{\partial\ell^{\mu}}\,\frac{1}{D_{2}^{\nu_{2}}}$
(2.26)
$\displaystyle=\frac{p^{\mu}}{D_{2}^{\nu_{2}}}\,\frac{2\nu_{1}\ell_{\mu}}{D_{1}^{\nu_{1}+1}}\hskip
9.38945pt\hskip 5.69046pt+\hskip
5.69046pt\frac{p^{\mu}}{D_{1}^{\nu_{1}}}\,\frac{2\nu_{2}(\ell_{\mu}+p_{\mu})}{D_{2}^{\nu_{2}+1}}\,.$
(2.27)

The combination $2p^{\mu}\ell_{\mu}=2p\cdot\ell$ appears in both terms. It can
be rewritten in terms of the two propagators $D_{1}$ and $D_{2}$ as follows:

$\displaystyle\begin{array}[]{ll}&D_{2}={\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-\ell^{2}}-p^{2}-{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}2p\cdot\ell}={\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}D_{1}-m^{2}}-p^{2}-{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}2p\cdot\ell}\implies\\\\[3.0pt]
&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}2p\cdot\ell}=D_{1}-D_{2}-m^{2}-p^{2}\,.\end{array}$
(2.30)

Inserting this into (2.27) gives

$\displaystyle\nu_{1}\frac{D_{1}-D_{2}-m^{2}-p^{2}}{D_{2}^{\nu_{2}}\,D_{1}^{\nu_{1}+1}}\hskip
5.69046pt+\hskip
5.69046pt\nu_{2}\frac{D_{1}-D_{2}-m^{2}-p^{2}+2p^{2}}{D_{1}^{\nu_{1}}\,D_{2}^{\nu_{2}+1}}\,.$
(2.31)

Simplifying, reinstating the integral sign, and equating to zero, we find the
family of IBPs

$\displaystyle 0$ $\displaystyle=(\nu_{2}-\nu_{1})I(\nu_{1},\nu_{2})$
$\displaystyle\hskip 5.69046pt+\hskip
5.69046pt\nu_{1}\bigl{[}\bigl{(}m^{2}+p^{2}\bigr{)}I(\nu_{1}+1,\nu_{2})+I(\nu_{1}+1,\nu_{2}-1)\bigr{]}$
(2.32) $\displaystyle\hskip 5.69046pt+\hskip
5.69046pt\nu_{2}\bigl{[}\bigl{(}m^{2}-p^{2}\bigr{)}I(\nu_{1},\nu_{2}+1)-I(\nu_{1}-1,\nu_{2}+1)\bigr{]}\,.$

Choosing instead $v^{\mu}=\ell^{\mu},$ then a similar calculation yields

$\displaystyle 0$
$\displaystyle=(\mathrm{D}-2\nu_{1}-\nu_{2})I(\nu_{1},\nu_{2})$
$\displaystyle\hskip 5.69046pt+\hskip
5.69046pt2m^{2}\nu_{1}I(\nu_{1}+1,\nu_{2})$ (2.33) $\displaystyle\hskip
5.69046pt+\hskip
5.69046pt\nu_{2}\bigl{[}(m^{2}-p^{2})I(\nu_{1},\nu_{2}+1)-I(\nu_{1}-1,\nu_{2}+1)\bigr{]}\,.$

The latter derivation requires the additional identity

$\displaystyle\frac{\partial}{\partial\ell^{\mu}}\ell^{\mu}=(\partial_{0},-\partial_{i})^{\mu}\,\eta_{\mu\nu}\,(\ell_{0},\ell_{i})^{\nu}=\mathrm{D}.$
(2.34)

Let’s now showcase an _IBP reduction_. Setting $\nu_{1}=\nu_{1}=1$ in (2.32),
we have

$\displaystyle 0=(m^{2}+p^{2})I(2,1)+I(2,0)+(m^{2}-p^{2})I(1,2)\,.$ (2.35)

Here we set $I(0,2)=0$. This is an example of a _scaleless_ integral, i.e. it
does not depend on any variables with mass dimension. It can be shown that
scaleless integrals vanish in DR [44, Page 40]222 $I(0,2)$ naively does appear
to depend on the scale $p$. But we can shift $\ell\to\ell-p$ without altering
the integral since we are integrating over the whole space
$\mathbb{R}^{1,\mathrm{D}-1}$, thereby removing the $p$-dependence. .

Setting also $\nu_{1}=\nu_{2}=1$ in (2.33), we get

$\displaystyle 0=(\mathrm{D}-3)I(1,1)+2m^{2}I(2,1)+(m^{2}-p^{2})I(1,2)\,.$
(2.36)

The last equation we will need is (2.33) for $\nu_{1}=1,\,\nu_{2}=0$:

$\displaystyle 0=(\mathrm{D}-2)I(1,0)+2m^{2}I(2,0)\,.$ (2.37)

Combining the three equations (2.35), (2.36) and (2.37), a few lines of
algebra grant

$\displaystyle
I(2,1)=\frac{\mathrm{D}-2}{2m^{2}(p^{2}-m^{2})}I(1,0)+\frac{\mathrm{D}-3}{p^{2}-m^{2}}I(1,1)\,.$
(2.38)

This is an example of an IBP reduction: we have expressed an integral $I(2,1)$
with "large" denominator powers in terms of two simpler integrals, $I(1,0)$
and $I(1,1)$. In fact, the IBPs (2.32) and (2.33) allow us to express _any_
integral $I(\nu_{1},\nu_{2})$ in terms of $I(1,0)$ and $I(1,1)$. These two
integrals are thereby dubbed the _master integrals_ (MIs) for the integral
family $I(\nu_{1},\nu_{2})$. These MIs are of course not unique. The MI
$I(1,1)$ could for instance be swapped for $I(2,1)$ via (2.38). $\blacksquare$

IBP reduction has been automated in many public codes such has kira [10], Fire
[11], Reduze [12], NeatIBP [13], Blade [14], and LiteRed [15].

##### 2.2.1 Irreducible scalar products

There is an important detail to deriving IBPs which did not feature in Example
(2.2), namely that of _irreducible scalar products_ (ISPs). Recall (2.30),
where we expressed the scalar product $2p\cdot\ell$ in terms of the
propagators $D_{1}$ and $D_{2}$. This is always possible at 1-loop order
because there is a one-to-one correspondence between scalar products and
propagators. At 2-loop order or higher, this is no longer the case, as there
will generally be more scalar products than propagators. These "extra" scalar
products are called _irreducible_.

Example 2.2. Take for instance a 2-loop integral written in terms $3$
propagators:

$\displaystyle
I(\nu_{1},\nu_{2},\nu_{2})=\int\frac{\mathrm{d}^{2\cdot\mathrm{D}}\ell}{D_{1}^{\nu_{1}}\,D_{2}^{\nu_{2}}\,D_{3}^{\nu_{3}}}\quad\text{where}\quad$
(2.39) $\displaystyle D_{1}=-\ell_{1}^{2}\,,\quad D_{2}=-\ell_{2}^{2}\,,\quad
D_{3}=-(\ell_{1}+\ell_{2}+p)^{2}\,,$

with $p$ being the only external momentum. It is possible to form $5>3$ scalar
products

$\displaystyle\ell_{1}^{2}\,,\quad\ell_{1}^{2}\,,\quad\ell_{1}\cdot\ell_{2}\,,\quad\ell_{1}\cdot
p\,,\quad\ell_{2}\cdot p\,.$ (2.40)

Two out of the last three scalar products are irreducible; say we choose
$\ell_{1}\cdot\ell_{2}\text{ and }\ell_{1}\cdot p.$ The trick is simply to
enlarge the original integral family to include two new integrand factors
$D_{4}=\ell_{1}\cdot\ell_{2}$ and $D_{5}=\ell_{1}\cdot p$, leading to

$\displaystyle
I(\nu_{1},\ldots,\nu_{5})=\int\frac{\mathrm{d}^{2\cdot\mathrm{D}}\ell}{D_{1}^{\nu_{1}}\cdots
D_{5}^{\nu_{5}}}\,.$ (2.41)

In intermediate steps of an IBP reduction, $D_{4}$ and $D_{5}$ will generally
appear in numerators. This is the reason for letting the edge weights
$\nu_{e}\in\mathbb{Z}$ be integral rather than just non-negative.
$\blacksquare$

####  Section 2.3The method of differential equations

So IBPs allow us to write a given FI in terms of MIs. The next question is how
to compute the MIs themselves. Let $I(\nu)=I(\nu_{1},\nu_{2},\ldots)$ denote a
family of FIs. This family depends on some kinematic variables, call them
$z=(z_{1},z_{2},\ldots)$, such as masses $m_{e}^{2}$ and Mandelstam variables
$s_{ij}$. Suppose a vector of MIs $\vec{I}$ has been revealed after IBP
reduction. The _method of differential equations_ for computing the MIs
consists of the following steps:

  1. 1.

Take derivatives $\partial_{i}\vec{I},$ w.r.t. the kinematic variables
$z_{i}$. These derivatives lead to shifts in the edge weights $\nu_{e}$.

  2. 2.

IBP reduce the integrals with shifted $\nu_{e}$ down to MIs $\vec{I}$. The
reduction can be factored as $P_{i}(z)\cdot\vec{I}$ for some matrices $P_{i}$
whose entries are rational functions of $z$ and $\mathrm{D}$. The resulting
system of DEQs is

$\displaystyle\partial_{i}\vec{I}=P_{i}(z)\cdot\vec{I}\,.$ (2.42)

This is called a _Pfaffian system_ , and will be a major topic of interest for
this thesis. The matrices $P_{i}$ are called _Pfaffian matrices_.

  3. 3.

Solve the Pfaffian system in terms of a _path-ordered exponential_.

Let us illustrate the first three steps by an example. The fourth step will be
discussed in the next section.

Example 2.1. This is a continuation of Example (2.2), wherein we found two MIs

$\displaystyle\vec{I}=\left[\begin{array}[]{c}I(1,0)\\\
I(1,1)\end{array}\right]$ (2.45)

for a 1-loop integral family.

There are two kinematic variables $(z_{1},z_{2})=(m^{2},p^{2})$ w.r.t. which
we can take derivatives. The derivative w.r.t. $z_{1}$ works out to

$\displaystyle\partial_{1}\vec{I}$
$\displaystyle=\left[\begin{array}[]{c}-I(2,0)\\\ -I(2,1)\end{array}\right]$
(2.48) $\displaystyle\overset{\rm
IBP}{=}\left[\begin{array}[]{c}\frac{\mathrm{D}-2}{2z_{1}}\,I(1,0)\\\\[5.0pt]
\frac{\mathrm{D}-2}{2z_{1}(z_{1}-z_{2})}\,I(1,0)+\frac{\mathrm{D}-3}{z_{1}-z_{2}}\,I(1,1)\end{array}\right]$
(2.51) $\displaystyle=\left[\begin{array}[]{c
c}\frac{\mathrm{D}-2}{2z_{1}}&0\\\
\frac{\mathrm{D}-2}{2z_{1}(z_{1}-z_{2})}&\frac{\mathrm{D}-3}{z_{1}-z_{2}}\end{array}\right]\cdot\left[\begin{array}[]{c}I(1,0)\\\
I(1,1)\end{array}\right]$ (2.56) $\displaystyle=P_{1}\cdot\vec{I}\,.$ (2.57)

A similar computation can be carried out for $z_{2}$. The result is

$\displaystyle\partial_{2}\vec{I}=\left[\begin{array}[]{cc}0&0\\\
\frac{\mathrm{D}-2}{2z_{2}(z_{1}-z_{2})}&\frac{(2-\mathrm{D})z_{1}+(4-\mathrm{D})z_{2}}{2z_{2}(z_{1}-z_{2})}\end{array}\right]\cdot\vec{I}=P_{2}\cdot\vec{I}\,.$
(2.60)

Pfaffian matrices such as $P_{1}$ and $P_{2}$ satisfy interesting identities.
The _Euler homogeneity condition_

$\displaystyle
z_{1}P_{1}+z_{2}P_{2}=\left[\begin{array}[]{cc}\frac{\mathrm{D}-2}{2}&0\\\
0&\frac{\mathrm{D}-4}{2}\end{array}\right]$ (2.63)

reflects the fact that the MIs are certain homogeneous functions of $z$.
Furthermore, the _integrability condition_

$\displaystyle\partial_{1}P_{2}-\partial_{2}P_{1}=[P_{1},P_{2}]$ (2.64)

follows from the fact that the MIs are $C^{2}$ (twice continuously
differentiable). $\blacksquare$

##### 2.3.1 Canonical bases

Obtaining the Pfaffian system for a basis of MIs is only part of the work, the
next step is to calculate a solution. In practice, we are only interested in
computing solutions up a specific order in the $\epsilon$-expansion (2.21). As
discovered by Henn [20], it is particularly simple to find this expansion if
we first write the Pfaffian system in _canonical form_.

The starting point is the Pfaffian system
$\partial_{i}\vec{I}=P_{i}(z|\epsilon)\cdot\vec{I}$ for $N$ kinematic
variables $z=(z_{1},\ldots,z_{N})$. Note that we highlight both the $z$\- and
$\epsilon$-dependence of the matrix $P_{i}$. It is more clean to write this
system in terms of the total derivative
$\mathrm{d}_{z}=\mathrm{d}z_{1}\frac{\partial}{\partial
z_{1}}+\ldots+\mathrm{d}z_{N}\frac{\partial}{\partial z_{N}}$ as

$\displaystyle\mathrm{d}_{z}\vec{I}=P(z|\epsilon)\cdot\vec{I}\quad\text{where}\quad
P=P_{1}(z|\epsilon)\,\mathrm{d}z_{1}+\ldots+P_{N}(z|\epsilon)\,\mathrm{d}z_{N}\,.$
(2.65)

Suppose we would like to change to a new basis of MIs $\vec{J}$. Due to IBP
identities, there must necessarily exist an invertible matrix
$G=G(z|\epsilon)$ such that $\vec{I}=G\cdot\vec{J}.$ Inserting this expression
for $\vec{I}$ into (2.65), it follows from the chain rule that

$\displaystyle\mathrm{d}_{z}\vec{J}=Q(z|\epsilon)\cdot\vec{J}\quad\text{where}\quad
Q_{i}(z|\epsilon)=G^{-1}\cdot\big{[}P_{i}\cdot G-\partial_{i}G\big{]}\,.$
(2.66)

We say that the new Pfaffian matrices $Q_{i}$ are obtained from a _gauge
transformation_ of $P_{i}$ by $G$. Now, assume that the basis $\vec{J}$ has
the following special properties:

  1. 1.

The $\epsilon$-dependence _factorizes_ in its associated Pfaffian matrices,
that is

$\displaystyle Q(z|\epsilon)=\epsilon\,Q(z)\,.$ (2.67)

  2. 2.

The poles coming from the denominators of the Pfaffian matrices all have order
1. In particular, we can write

$\displaystyle\epsilon\,Q(z)=\epsilon\sum_{\eta\in
A}Q_{i}\,\mathrm{d}\log\eta(z)\,,$ (2.68)

where

$\displaystyle\mathrm{d}\log\eta(z)=\frac{1}{\eta(z)}\frac{\partial\eta(z)}{\partial
z_{1}}\mathrm{d}z_{1}+\ldots+\frac{1}{\eta(z)}\frac{\partial\eta(z)}{\partial
z_{N}}\mathrm{d}z_{N}\,.$ (2.69)

Every matrix entry $(Q_{i})_{ab}\in\mathbb{Q}$ is here a rational number, and
the (possibly algebraic or even transcendental) functions $\eta(z)$ are called
_letters_. The full set of letters, denoted here by $A$, is called the
_alphabet_.

By virtue of these two properties, $\vec{J}$ is called a _canonical basis_.
Many calculations in the literature have shown that when condition 1 is
satisfied, then condition 2 follows automatically - though, to the best of our
knowledge, this has yet to be proven rigorously. In practice, the effort in
searching for canonical bases is hence put into $\epsilon$-factorization.

The $\epsilon$-expansion of the canonical solution vector $\vec{J}$ is now
particularly simple:

$\displaystyle\vec{J}=\mathcal{P}\exp\left(\epsilon\int_{\Gamma}Q(z)\right)\cdot\vec{J}_{0}(\epsilon)\,,$
(2.70)

where the _path-ordered exponential_ is given by

$\displaystyle\mathcal{P}\exp\left(\epsilon\int_{\Gamma}Q(z)\right)=\mathbf{1}+\epsilon\int_{0}^{t}Q(\Gamma(t_{1}))+\epsilon^{2}\int_{0}^{t}Q(\Gamma(t_{1}))\int_{0}^{t_{1}}Q(\Gamma(t_{2}))+\mathcal{O}(\epsilon^{3})$
(2.71)

for some chosen integration contour $\Gamma$ parametrized by $t\in\mathbb{R}$.
The vector of boundary constants also has an $\epsilon$-expansion:

$\displaystyle\vec{J}_{0}(\epsilon)=\vec{J}_{0}\hskip
0.1424pt^{(0)}+\epsilon\vec{J}_{0}\hskip
0.1424pt^{(1)}+\epsilon^{2}\vec{J}_{0}\hskip
0.1424pt^{(2)}+\mathcal{O}(\epsilon^{3})\,,$ (2.72)

with each vector $\vec{J}_{0}\hskip 0.1424pt^{(n)}$ consisting solely of real
numbers. These boundary vectors are often determined by considering special
limits of the MIs; concrete examples are given in Section (3.3).

The solution above is written somewhat abstractly in terms of a matrix-valued
one-form $Q(z)$ (given in (2.68)) and some unspecified contour $\Gamma$. In
Section (2.4), we shall rewrite the iterated integrals (2.71) in form that is
more suitable for computer implementations.

##### 2.3.2 Magnus expansion

With the discussion above in mind, the key question is thus the following:
given a Pfaffian system $\partial_{i}\vec{I}=P_{i}\cdot\vec{I},$ is it
possible to algorithmically construct a gauge transformation $G$ such that the
vector $\vec{J}=G^{-1}\cdot\vec{I}$ satisfies a new system in canonical form?
Such an algorithm should preferably terminate in reasonable time, and be
applicable to a large category of FIs. This question has been, and is, subject
to intense study (see e.g. [44, Chapter 7] and the references therein). In
this section, we present a method for that calculating canonical forms that is
at least effective for the class of functions called _generalized
polylogarithms_ (GPLs), to be defined in Section (2.4). (There are now also
many methods that output canonical forms for FIs even when they cannot be
evaluated in terms of GPLs [48, 49, 50, 51, 52, 53, 54].)

The approach described here is based on the _Magnus expansion_ for solving
systems of DEQs [55]. 60 years after Magnus’ work, Argeri et al. [56]
discovered that his expansion can be used to obtain canonical forms. The
starting point is a Pfaffian system in _pre-canonical form_ :

$\displaystyle\mathrm{d}_{z}\vec{I}=\big{[}P^{(0)}(z)+\epsilon
P^{(1)}(z)\big{]}\cdot\vec{I}\,.$ (2.73)

I.e. a system linear in $\epsilon$ (the algorithm can be generalized to a
polynomial dependence on $\epsilon$ too [57, Section 5.2.4]). By choosing
judicial edge weights and simple $\epsilon$-dependent prefactors for the MIs
$\vec{I}$, it is often not too difficult to find a pre-canonical form.

The clever observation is now this: suppose that there exists a matrix
$G=G(z)$, depending only on $z$, subject to the matrix DEQ

$\displaystyle\mathrm{d}_{z}G=P^{(0)}(z)\cdot G\,.$ (2.74)

If we change basis to $\vec{I}=G\cdot\vec{J}$, then the gauge transformation
(2.66) simplifies to

$\displaystyle Q(z|\epsilon)$
$\displaystyle=G^{-1}\cdot\big{[}\big{(}P^{(0)}+\epsilon P^{(1)}\big{)}\cdot
G-\mathrm{d}_{z}G\big{]}$ (2.75)
$\displaystyle=G^{-1}\cdot\big{[}\big{(}P^{(0)}+\epsilon P^{(1)}\big{)}\cdot
G-P^{(0)}\cdot G\big{]}$ (2.76) $\displaystyle=\epsilon\,G^{-1}\cdot
P^{(1)}\cdot G\,,$ (2.77)

which is indeed $\epsilon$-factorized.

The question is thus how to find the solution $G$ to the matrix DEQ (2.74). It
turns out that the solution takes the form of a matrix exponential:

$\displaystyle G=e^{\Omega\left[P^{(0)}\right]}\cdot G(z=0)\,.$ (2.78)

For our purposes, we can simply take the boundary condition to be the identity
matrix, $G(z=0)=\mathbf{1}$. The $\Omega$ matrix in the exponent enjoys a
_Magnus expansion_

$\displaystyle\Omega\big{[}P^{(0)}\big{]}$
$\displaystyle=\sum_{n=1}^{\infty}\Omega_{n}\big{[}P^{(0)}\big{]}\,,$ (2.79)

with each term computed by integrating iterated commutators:

$\displaystyle\Omega_{1}\big{[}P^{(0)}\big{]}$
$\displaystyle=\int_{\bullet}^{z}\mathrm{d}z_{1}\,P^{(0)}(z_{1})$
$\displaystyle\Omega_{2}\big{[}P^{(0)}\big{]}$
$\displaystyle=\frac{1}{2}\int_{\bullet}^{z}\mathrm{d}z_{1}\int_{\bullet}^{z_{1}}\mathrm{d}z_{2}\,\left[P^{(0)}(z_{1}),\,P^{(0)}(z_{2})\right]$
(2.80) $\displaystyle\Omega_{3}\big{[}P^{(0)}\big{]}$
$\displaystyle=\frac{1}{6}\int_{\bullet}^{z}\mathrm{d}z_{1}\int_{\bullet}^{z_{1}}\mathrm{d}z_{2}\int_{\bullet}^{z_{2}}\mathrm{d}z_{3}\,\Big{\\{}\left[P^{(0)}(z_{1}),\,\left[P^{(0)}(z_{2}),\,P^{(0)}(z_{3})\right]\right]+$
$\displaystyle\hskip
145.10922pt\left[P^{(0)}(z_{3}),\,\left[P^{(0)}(z_{2}),\,P^{(0)}(z_{1})\right]\right]\Big{\\}}\,.$

Clean formulas for the $n$th term of this expansion can be found in [58]. Here
$P^{(0)}$ was assumed to depend on just a single variable $z$. For
multivariate problems, one can simply apply these formulas one variable at a
time [59, Page 108].

It turns out that the Magnus expansion truncates for FIs that can be evaluated
in terms of GPLs, i.e. $\Omega_{n_{\mathrm{max}}}=0$ for some
$n_{\mathrm{max}}\in\mathbb{Z}_{>0}$. In such cases, the matrix $P^{(0)}$ is
apparently sparse enough to become nilpotent, meaning that the iterated
commutators eventually evaluate to zero.

Example 2.2. Here we apply the Magnus expansion to the bubble integral from
Examples (2.2) and (2.3). Consider a basis with squared propagators:

$\displaystyle\vec{I}=\left[\begin{array}[]{c}I(2,0)\\\
I(2,1)\end{array}\right]\,.$ (2.83)

Using IBPs as in Example (2.3), it is possible to obtain two Pfaffian matrices
$P_{1}$ and $P_{2}$ in the variables $z_{1}=m^{2}$ and $z_{2}=p^{2}$. Working
in units such that $z_{1}=1$, then $P_{1}$ can be discarded. Call $z_{2}=z$.
We thence get a pre-canonical system

$\displaystyle\frac{\partial}{\partial
z}\vec{I}=\Big{[}P^{(0)}+\epsilon\,P^{(1)}\Big{]}\cdot\vec{I}\,,$ (2.84)

where the constant and linear matrices w.r.t. $\epsilon$ are

$\displaystyle P^{(0)}=\left[\begin{array}[]{cc}0&0\\\
0&-\frac{1}{z}\end{array}\right]\quad\text{and}\quad
P^{(1)}=\left[\begin{array}[]{cc}0&0\\\
\frac{1}{(1-z)z}&\frac{1+z}{(1-z)z}\end{array}\right]\,.$ (2.89)

The Magnus expansion truncates at the first term, giving

$\displaystyle\Omega\big{[}P^{(0)}\big{]}=\Omega_{1}\big{[}P^{(0)}\big{]}=\int_{\bullet}^{z}\left[\begin{array}[]{cc}0&0\\\
0&-\frac{1}{z^{\prime}}\end{array}\right]\mathrm{d}z^{\prime}=\left[\begin{array}[]{cc}0&0\\\
0&-\log(z)\end{array}\right]\,,$ (2.94)

wherefore

$\displaystyle
G=e^{\Omega\left[P^{(0)}\right]}=\left[\begin{array}[]{cc}1&0\\\
0&\frac{1}{z}\end{array}\right]\,.$ (2.97)

Gauge transforming with $G$ via (2.77) yields a new Pfaffian matrix

$\displaystyle Q=\epsilon\left[\begin{array}[]{cc}0&0\\\
\frac{1}{1-z}&\frac{2}{1-z}+\frac{1}{z}\end{array}\right]\,,$ (2.100)

which is indeed canonical. $\blacksquare$

Example 2.3. Every physicist’s favorite equation is the DEQ for the harmonic
oscillator:

$\displaystyle\frac{\mathrm{d}^{2}f(t)}{\mathrm{d}t^{2}}+\omega^{2}f(t)=0\,.$
(2.101)

The angular frequency $\omega$ is taken to be constant. Let us solve this
equation using the Magnus expansion. The initial step is to rewrite the DEQ as
a first-order matrix equation. Setting $f_{1}(t)=f(t)$ and
$f_{2}(t)=\frac{\mathrm{d}f(t)}{\mathrm{d}t}$, then the DEQ is equivalent to

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\left[\begin{array}[]{c}f_{1}(t)\\\
f_{2}(t)\end{array}\right]=\left[\begin{array}[]{cc}0&1\\\
-\omega^{2}&0\end{array}\right]\cdot\left[\begin{array}[]{c}f_{1}(t)\\\
f_{2}(t)\end{array}\right]\,.$ (2.108)

In addition to solving _matrix_ DEQs such as (2.74), the Magnus expansion also
grants solution _vectors_ in the form

$\displaystyle\left[\begin{array}[]{c}f_{1}(t)\\\
f_{2}(t)\end{array}\right]=e^{\Omega(t)}\cdot\left[\begin{array}[]{c}f_{1}(0)\\\
f_{2}(0)\end{array}\right]\,.$ (2.113)

The Magnus expansion truncates at the first term because $\omega$ is constant.
Thus,

$\displaystyle e^{\Omega(t)}$
$\displaystyle=\exp\left(\int_{0}^{t}\left[\begin{array}[]{cc}0&1\\\
-\omega^{2}&0\end{array}\right]\mathrm{d}t^{\prime}\right)$ (2.116)
$\displaystyle=\left[\begin{array}[]{cc}\cosh\left(\sqrt{-\omega^{2}}t\right)&\frac{\sinh\left(\sqrt{-\omega^{2}}t\right)}{\sqrt{-\omega^{2}}}\\\
\sqrt{-\omega^{2}}\sinh\left(\sqrt{-\omega^{2}}t\right)&\cosh\left(\sqrt{-\omega^{2}}t\right)\end{array}\right]$
(2.119) $\displaystyle=\left[\begin{array}[]{cc}\cos(\omega
t)&\frac{\sin(\omega t)}{\omega}\\\ -\omega\sin(\omega t)&\cos(\omega
t)\end{array}\right]\,.$ (2.122)

We obtain oscillating functions as expected. Further, this solution satisfies
$f_{2}(t)=\dot{f}_{1}(t)$. $\blacksquare$

####  Section 2.4Generalized polylogarithms

The path-ordered exponential (2.71) naturally leads to a class of functions
called _Chen iterated integrals_ [60]. Let us first say what we mean by an
integration contour for these integrals. Suppose that there are $N$ kinematic
variables $z=(z_{1},\ldots,z_{N})$ which parametrize an $N$-dimensional
(smooth) complex manifold $Z$. We write write a path in this space as

$\displaystyle\Gamma:[a,b]\to Z\,,$ (2.123)

where $\Gamma(a)=z_{a}\in Z$ and $\Gamma(b)=z_{b}\in Z$ are the start and end
points. One can think of $t\in[a,b]\subset\mathbb{R}$ as being a "time"
coordinate.

Let $\\{\omega_{1},\ldots,\omega_{M}\\}$ denote a collection of differential
1-forms on $Z$. For example, they could be the $\mathrm{d}\\!\log$ forms from
(2.68). These 1-forms are supposed to be integrated over the contour $\Gamma$.
So instead of parametrizing the $\omega_{i}$ in terms of $z$, we would like to
parametrize them using $t\in[a,b]$. This is called a _pullback_ to the
interval $[a,b]$, written as

$\displaystyle\Gamma^{*}\omega_{i}=f_{i}(t)\mathrm{d}t\,.$ (2.124)

This gives a concrete expression for the integral of 1-form along $\Gamma$:

$\displaystyle\int_{\Gamma}\omega_{i}=\int_{a}^{b}f_{i}(t)\mathrm{d}t\,.$
(2.125)

The question is how to compute the function $f_{i}(t)$. Writing

$\displaystyle\omega_{i}=\sum_{j=1}^{N}\omega_{ij}(z)\mathrm{d}z_{j}\quad\text{and}\quad\Gamma(t)=\left[\begin{array}[]{c}\Gamma_{1}(t)\\\
\vdots\\\ \Gamma_{N}(t)\end{array}\right]\,,$ (2.129)

then the formula for $f_{i}(t)$ follows from the chain rule:

$\displaystyle\Gamma^{*}\omega_{i}=\sum_{j=1}^{N}\omega_{ij}\big{(}\Gamma(t)\big{)}\frac{\mathrm{d}\Gamma_{j}(t)}{\mathrm{d}t}\mathrm{d}t\quad\implies\quad
f_{i}(t)=\sum_{j=1}^{N}\omega_{ij}\big{(}\Gamma(t)\big{)}\frac{\mathrm{d}\Gamma_{j}(t)}{\mathrm{d}t}\,.$
(2.130)

Example 2.1. To demystify the construction above, let us compute

$\displaystyle\int_{\Gamma}\omega$ (2.131)

with

$\displaystyle\omega=\frac{2}{z_{2}-2}\mathrm{d}z_{1}+\frac{3}{z_{1}+4}\mathrm{d}z_{2}\,,\quad\text{and}\quad\Gamma:[0,1]\to\mathbb{C}^{2}\quad\text{given
by}\quad\Gamma(t)=\left[\begin{array}[]{c}t^{2}\\\ 1-t\end{array}\right]\,.$
(2.134)

Inserting this data into (2.130) gives

$\displaystyle f(t)$
$\displaystyle=\omega_{1}\big{(}\Gamma(t)\big{)}\times\frac{\mathrm{d}\Gamma_{1}(t)}{\mathrm{d}t}\hskip
5.69046pt+\hskip
5.69046pt\omega_{2}\big{(}\Gamma(t)\big{)}\times\frac{\mathrm{d}\Gamma_{2}(t)}{\mathrm{d}t}$
(2.135) $\displaystyle=\frac{2}{(1-t)-2}\times 2t\hskip 5.69046pt+\hskip
5.69046pt\frac{3}{t^{2}+4}\times(-1)$ (2.136)
$\displaystyle=-\frac{4t^{3}+19t+3}{(1+t)(4+t^{2})}\,.$ (2.137)

Using any CAS, one can then easily integrate

$\displaystyle\int_{\Gamma}\omega\ =\ \int_{0}^{1}f(t)\mathrm{d}t\ \approx\
-1.923\,.$ (2.138)

$\blacksquare$

Given a $t\in[a,b]$, the _Chen iterated integral_ can now be defined as

$\displaystyle I_{\Gamma}(\omega_{1},\ldots,\omega_{M}|t)$
$\displaystyle=\int_{\Gamma}\omega_{1}\wedge\cdots\wedge\omega_{M}$ (2.139)
$\displaystyle=\int_{a}^{t}f_{1}(t_{1})\mathrm{d}t\int_{a}^{t_{1}}f_{2}(t_{2})\mathrm{d}t_{2}\cdots\int_{a}^{t_{M-1}}f_{M}(t_{M})\mathrm{d}t_{M}\,.$
(2.140)

An equivalent, recursive definition is

$\displaystyle I_{\Gamma}(\,|t)$ $\displaystyle=1$ (2.141) $\displaystyle
I_{\Gamma}(\omega_{1},\ldots,\omega_{M}|t)$
$\displaystyle=\int_{a}^{t}f_{1}(t_{1})I_{\Gamma}(\omega_{2},\ldots,\omega_{M}|t_{1})\mathrm{d}t_{1}\,.$
(2.142)

Chen proved that $I_{\Gamma}$ is independent of $\Gamma$, provided that
$\Gamma$ does not hit any singularities nor cross any branch points of the
integrand for $a<t<b$. If the integrand has singularities at the endpoints
$t=a$ or $t=b$, these are assumed to be integrable singularities (caveat:
endpoint singularities can also be regularized using so-called shuffle
identities [61]).

Consider now a special class of 1-forms with a simple pole at $t=z_{i}$:

$\displaystyle\omega_{\mathrm{GPL}}(z_{i})=\frac{\mathrm{d}t}{t-z_{i}}=\mathrm{d}\log(t-z_{i})\,.$
(2.143)

The complex number $z_{i}$ is treated as a parameter. Let $\Gamma$ be a path
along the non-negative real axis such that $\Gamma(a=0)=0$ and
$\Gamma(t)=y\in\mathbb{R}_{\geq 0}$. The special class of Chen iterated
integrals called _generalized multiple polylogarithms_ (GPLs) [62] is defined
by

$\displaystyle
G(z_{1},\ldots,z_{M}|y)=I_{\Gamma}\big{(}\omega_{\mathrm{GPL}}(z_{1}),\ldots,\omega_{\mathrm{GPL}}(z_{M})|y\big{)}\,.$
(2.144)

Explicitly,

$\displaystyle G(\underbrace{0,\ldots,0}_{M\text{ times }}|y)$
$\displaystyle=\frac{1}{M!}\log^{M}(y)$ (2.145) $\displaystyle
G(z_{1},\ldots,z_{M}|y)$
$\displaystyle=\int_{0}^{y}\frac{\mathrm{d}y_{1}}{y_{1}-z_{1}}G(z_{2},\ldots,z_{M}|y_{1})\,.$
(2.146)

The $z$-parameters are often called _weights_. The case $z=0$ above naively
constitutes a divergence of the integral representation, so the first formula
above is more of a stipulation than a derivation.

There is a huge literature on these functions, see [61] [44, Chapters 6,8] and
references therein. There are also plenty of powerful codes for symbolically
manipulating and numerically evaluating these functions: the Maple package
[63], Hyperint [64], HPLog [65], HPL [66], PolyLogTools [67], Ginac [68],
Chaplin [69], handyG [70], and FastGPL [71].

Now harken back to the path-ordered exponential (2.71). If the functions
$\eta(z)$ from (2.68) are _rational_ functions of the kinematic variables $z$,
then the solution for the FI in question can be written in terms of GPLs. This
is indeed a favorable form for the solution because so much is known about the
identities among GPLs, how to analytically continue them, and how to
numerically evaluate them.

One may encounter letters $\eta$ that depend on $z$ in terms of algebraic
roots. In those cases, one can sometimes find a change of a variables $w=f(z)$
that turns the roots into rational functions. The package RationalizeRoots
[72] automates the procedure for finding such an $f$.

### Chapter 3Two-loop Feynman Integrals for Dark Matter Detection

We now apply the multi-loop technology from the previous chapter to compute
form factors for a _leptophillic dark matter_ (DM) model. Form factors are
certain building blocks of scattering amplitudes. For phenomenological
implications of these results, we refer to the paper [73] on which this
chapter is based.

####  Section 3.1The model

Consider a DM particle $\chi$ that scatters elastically off a nucleon $A$:
$\chi A\to\chi A.$ This is called a _direct DM search_ , and is probed by
remarkable experiments having large tanks filled with e.g. oxygen, helium or
xenon nucleons.

Assuming that DM does not couple _directly_ to SM particles, this scattering
can only occur via a _mediator_ particle $\phi$ that _does_ couple to the SM.
Here we consider a _leptophillic_ DM model, wherein $\chi$ couples to $\phi$,
which in turn only couples to SM leptons; the interaction between leptons and
the quarks inside the nucleon $A$ is finally mediated via photons. The
mediator $\phi$ is taken to be a singlet scalar (called $\phi_{S}$) or a
pseudo-scalar (called $\phi_{P}$). The Lagrangian interaction terms are thus
of the form

$\displaystyle-\mathcal{L}_{S}\ \supseteq\ g_{S}\,\phi_{S}\,\bar{\ell}\,\ell\
+\ g_{\chi}\,\phi_{S}\,\Gamma_{\chi}\
\quad\text{or}\quad-\mathcal{L}_{P}\,\supseteq\
i\,g_{P}\,\phi_{P}\,\bar{\ell}\,\gamma_{5}\,\ell\ +\
g_{\chi}\,\phi_{P}\,\Gamma_{\chi}\,,$ (3.1)

where $\Gamma_{\chi}=\chi^{\dagger}\chi$ for scalar DM and
$\Gamma_{\chi}=\\{\bar{\chi}\,\chi,\,\bar{\chi}\,\gamma_{5}\,\chi\\}$ for
fermionic DM, $\\{g_{S},g_{P},g_{\chi}\\}$ are coupling constants, and
$\\{\ell,\bar{\ell}\\}$ denote a lepton and an anti-lepton.

Since this model assumes no tree-level couplings to quarks, the $\chi A$
scattering cross section is purely loop-induced. It turns out that the leading
contributions come from the diagrams111 These are only the so-called one-body
diagrams. The virtual photons can also scatter off two nucleons, yielding
additional one-loop two-body diagrams. See [73] for more details.  in Figure
(3.1).

Figure 3.1: Leading contributions to the scattering amplitude for the
leptophillic DM model.

This is a $t$-channel process, with an incoming mediator $\phi(q)$ (dashed
black lines) transferring momentum to a nucleon $A(p)$ (leftmost  red lines),
leading to a final state nucleon $A(p^{\prime})$ (rightmost  red lines).
Momentum conservation dictates that $q=p^{\prime}-p$. The  blue triangle loops
represent a lepton which couples to (the quarks inside) the nucleon via
photons (wiggly black lines).

Let $m_{A}$ denote the nucleon mass, meaning that it can either be the proton
or neutron mass. The kinematic invariants for this process are then defined by

$\displaystyle t=q^{2}=(p^{\prime}-p)^{2},\quad m_{A}^{2}=p^{2}=p^{\prime
2}\,.$ (3.2)

####  Section 3.2The scattering amplitude

According to the Feynman rules dictated by the Lagrangian terms (3.1), the
desired scattering amplitude takes the form

$\displaystyle\mathcal{A}_{S,P}(t,m_{\ell},m_{A})=i\,g_{S,P}\,Q_{A}^{2}\sum_{\ell=e,\,\mu,\,\tau}Q_{\ell}^{2}\Big{[}\bar{u}_{A}(p^{\prime})\cdot\Gamma_{S,P}(t,m_{\ell},m_{A})\cdot
u_{A}(p)\Big{]}\,.$ (3.3)

Let us parse this expression from left to right. The amplitude $\mathcal{A}$
has an index $S$ or $P$, depending on whether the mediator $\phi_{S,P}$ is
taken to be scalar or pseudo-scalar. $\mathcal{A}$ has kinematic dependence on
the momentum transfer variable $t$, the lepton masses
$m_{\ell}=\\{m_{e},m_{\mu},m_{\tau}\\}$ and the nucleon mass $m_{A}$. The
prefactor in front of the sum contains the coupling constants $g_{S,P}$ from
(3.1) and the electric charge of the nucleon, where $Q_{A}=1$ for the proton
and $Q_{A}=0$ for the neutron222 The expression (3.3) therefore vanishes when
$A$ is the neutron. This is an approximation which holds only in the non-
relativistic limit, which is the regime being probed by DM experiments. See
[73] for more details. . The lepton electric charge inside the sum is
$Q_{\ell}=1$. The sum over leptons contains an operator
$\Gamma_{S,P}(t,m_{\ell},m_{A})$ sandwiched between spinors
$\bar{u}(p^{\prime})$ and $u(p)$ for the final and initial state nucleons
respectively. This operator stems from the sum of the two diagrams above. Its
analytic expression is

$\displaystyle\Gamma_{S,P}(t,m_{\ell},m_{A})=-32\pi^{2}\alpha_{\mathrm{em}}\int\frac{\mathrm{d}^{\mathrm{D}}\ell_{1}}{(2\pi)^{\mathrm{D}}}\frac{\mathrm{d}^{\mathrm{D}}\ell_{2}}{(2\pi)^{\mathrm{D}}}\frac{g_{\mu\rho}\
\gamma^{\rho}\ (\not{\ell}_{2}-m_{A})\ g_{\nu\sigma}\ \gamma^{\sigma}\
\mathrm{Tr}_{S,P}^{\mu\nu}}{D_{1}D_{2}D_{3}D_{4}D_{5}D_{6}}\,,$ (3.4)

where the electromagnetic coupling is $\alpha_{\mathrm{em}}=e^{2}/(4\pi),$ and
the Dirac trace over spinor indices is given by

$\displaystyle\mathrm{Tr}_{S,P}^{\mu\nu}=\mathrm{Tr}\left[(\not{\ell}_{1}+\not{q}+m_{\ell})\cdot\Lambda_{S,P}\cdot(\not{\ell}_{1}+m_{\ell})\cdot\gamma^{\mu}\cdot(\not{\ell}_{1}+\not{\ell}_{2}+\not{p}^{\prime}+m_{\ell})\cdot\gamma^{\nu}\right]\,,$
(3.5)

with $\Lambda_{S}=\mathbf{1}$ and $\Lambda_{P}=\gamma_{5}$. The denominators
of the momentum space integral (3.4) are

$\displaystyle D_{1}$ $\displaystyle=\ell_{1}^{2}-m_{\ell}^{2},\quad$
$\displaystyle D_{2}$ $\displaystyle=(\ell_{1}+q)-m_{\ell}^{2},\quad$
$\displaystyle D_{3}$ $\displaystyle=\ell_{2}^{2}-m_{A}^{2}$ (3.6)
$\displaystyle D_{4}$ $\displaystyle=(\ell_{2}+p)^{2},\quad$ $\displaystyle
D_{5}$ $\displaystyle=(\ell_{2}+p^{\prime})^{2},\quad$ $\displaystyle D_{6}$
$\displaystyle=(\ell_{1}+k_{2}+p^{\prime})^{2}-m_{\ell}^{2}\,.$

It is standard to write such an amplitude containing tensor structure as

$\displaystyle\mathcal{A}_{S,P}(t,m_{\ell},m_{A})\propto\sum_{\ell=e,\mu,\tau}\mathcal{F}_{S,P}$
(3.7)

where the _form factors_ $\mathcal{F}_{S,P}$ are related to $\Gamma_{S,P}$ via
projection:

$\displaystyle\mathcal{F}_{S,P}(t,m_{\ell},m_{A})=\frac{1}{2(p^{\prime}\pm
p)^{2}}\mathrm{Tr}\left[\Lambda_{S,P}\cdot(\not{p}\pm
m_{A})\cdot\Gamma_{S,P}\cdot(\not{p}^{\prime}\pm m_{A})\right]\,.$ (3.8)

The signs are $(+1)$ and $(-1)$ for the scalar ($S$) and pseudo-scalar ($P$)
cases respectively. These form factors can in turn be expanded as sums over
scalar FIs:

$\displaystyle\mathcal{F}_{S,P}=\sum_{\vec{\nu}}c^{\vec{\nu}}_{S,P}(t,m_{\ell},m_{A}|\epsilon)\times
I_{\vec{\nu}}(t,m_{\ell},m_{A}|\epsilon),\quad\vec{\nu}\in\mathbb{N}^{7}\,.$
(3.9)

The integral family $I_{\vec{\nu}}$ is given by

$\displaystyle
I_{\nu_{1},\ldots,\nu_{7}}(t,m_{\ell},m_{A}|\epsilon)=\int\frac{\mathrm{d}^{\mathrm{D}}\ell_{1}}{(2\pi)^{\mathrm{D}}}\frac{\mathrm{d}^{\mathrm{D}}\ell_{2}}{(2\pi)^{\mathrm{D}}}\frac{D_{7}^{\nu_{7}}}{D_{1}^{\nu_{1}}\cdots
D_{6}^{\nu_{6}}}\,,$ (3.10)

with denominators defined in (3.6) and an ISP chosen to be

$\displaystyle D_{7}=(\ell_{1}-p)^{2}\,.$ (3.11)

The coefficients $c^{\vec{\nu}}_{S,P}$ are rational functions of the kinematic
variables and the DR parameter $\epsilon$, and stem from carrying out the
Dirac algebra333 The prescription for treating $\gamma_{5}$ in $\mathrm{D}\neq
4$ dimensions is taken from [74]. .

It turns out that $\mathcal{F}_{S}$ is a sum of $34$ integrals, and
$\mathcal{F}_{P}$ a sum of $5$ integrals. The $34$ integrals are not all
linearly independent due to IBP relations. The goal for the rest of this
chapter is to compute a linearly independent set of MIs, which in turn leads
to expressions for the form factors $\mathcal{F}_{S,P}$ via (3.9).

####  Section 3.3Computing the Feynman integrals

Although the diagrams in Figure (3.1) are finite in $\mathrm{D}=4$ dimensions,
we treat them in DR in order to benefit from IBPs and DEQs. The limit
$\epsilon\to 0$ in (3.9) then recovers the desired form factor results.

For the purpose of verification, we compute MIs for both the two-mass
($m_{A}\neq m_{\ell}$) and equal-mass ($m_{A}=m_{\ell}$) cases. The _ab
initio_ equal-mass calculation can then be compared with the limit $m_{A}\to
m_{\ell}$ taken on the two-mass form factors.

##### 3.3.1 Two-mass case

We begin by computing MIs for the family of FIs associated to Figure (3.2).

Figure 3.2: Two-mass Feynman diagram.

The integral family was defined in the formulas (3.6), (3.10) and (3.11). As
before, the kinematics are $p^{2}=p^{\prime 2}=m_{A}^{2}$ and
$t=(p-p^{\prime})^{2}$. MIs for this topology have already been computed in
[75] in the context of $H\to\bar{b}\,b$ decay; here we re-compute them
independently.

IBP reduction and differentiation of FIs is performed with the codes Reduze2
[12] and LiteRed [15]. Using the Magnus expansion from Section (2.3.2), the
following set of canonical MIs are identified:

$\displaystyle\begin{array}[]{lllllll}&I_{1}&=&\epsilon^{2}J_{1}&I_{2}&=&\epsilon^{2}J_{2}\\\
&I_{3}&=&\epsilon^{2}\lambda_{\ell}J_{3}&I_{4}&=&-\epsilon^{2}tJ_{4}\\\
&I_{5}&=&\epsilon^{2}\left[\frac{\lambda_{\ell}-t}{2}J_{4}+\lambda_{\ell}J_{4}\right]&I_{6}&=&-\epsilon^{2}tJ_{6}\\\
&I_{7}&=&\frac{\epsilon^{2}m_{A}^{2}\rho_{A}(t+\lambda_{\ell})}{\rho_{\ell}(t+\lambda_{A})}[J_{7}+2J_{8}]&I_{8}&=&\epsilon^{2}m_{A}^{2}J_{8}\\\
&I_{9}&=&\epsilon^{2}\lambda_{\ell}J_{9}&I_{10}&=&-\epsilon^{2}t\lambda_{\ell}J_{10}\\\
&I_{11}&=&\epsilon^{3}\lambda_{A}J_{11}&I_{12}&=&\frac{\epsilon^{2}\lambda_{\ell}}{4t}\left[(t-\lambda_{A})(J_{4}+2J_{5})-4m_{A}^{2}\lambda_{A}J_{12}\right]\\\
&I_{13}&=&\epsilon^{3}\lambda_{A}J_{13}&I_{14}&=&\epsilon^{3}(2\epsilon-1)tJ_{14}\\\
&I_{15}&=&\epsilon^{3}\lambda_{\ell}\lambda_{A}J_{15}&I_{16}&=&\epsilon^{3}\lambda_{A}J_{16}\\\
&I_{17}&=&\epsilon^{3}\lambda_{A}J_{17}&I_{18}&=&\eqref{two_mass_MI_I18}\\\
&I_{19}&=&\eqref{two_mass_MI_I19}&I_{20}&=&-\epsilon^{4}t\lambda_{A}J_{20}\,,\end{array}$
(3.22)

where

$\displaystyle I_{18}$
$\displaystyle=\frac{\epsilon^{2}}{t}\left[\lambda_{\ell}(\lambda_{A}-t)J_{9}+\epsilon\lambda_{A}(t-\lambda_{\ell})J_{17}+(2\epsilon-1)\lambda_{\ell}\lambda_{A}J_{18}\right]$
(3.23)

and

$\displaystyle I_{19}$
$\displaystyle=\frac{\epsilon^{2}}{2t}\Big{[}t(\lambda_{\ell}-t)J_{3}-2tm_{\ell}^{2}(J_{7}+2J_{8})+\big{(}4tm_{\ell}^{2}+\lambda_{\ell}(\lambda_{A}-t)\big{)}J_{9}+\frac{4t^{2}m_{A}^{2}}{\lambda_{A}+t}J_{16}$
(3.24)
$\displaystyle\quad+\epsilon\big{(}\lambda_{A}(t-\lambda_{\ell})4tm_{\ell}^{2}\big{)}J_{17}+(2\epsilon-1)(4tm_{\ell}^{2}+\lambda_{\ell}\lambda_{A}-t^{2})J_{18}+2t^{2}(m_{\ell}^{2}-m_{A}^{2})J_{19}\Big{]}\,.$

The integrals $J_{i}$ are graphically represented in Figure (3.3).

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

(m)

(n)

(o)

(p)

(q)

(r)

(s)

(t)

Figure 3.3:  The MIs $J_{i}$ appearing in the canonical basis (3.22). Squared
propagators are represented by dots.

The canonical integrals contain square root prefactors $\lambda_{i}$ and
$\rho_{i}$ defined by

$\displaystyle\lambda_{i}=\sqrt{-t}\sqrt{4m_{i}^{2}-t}\,,\quad\rho_{i}=\sqrt{\frac{2m_{i}^{2}-t-\lambda_{i}}{m_{i}^{2}}}\,,\quad
i=\ell,A\,.$ (3.25)

These roots are simultaneously rationalized by the change of variables

$\displaystyle t=-m_{\ell}^{2}\frac{(1-x)^{2}}{x^{2}}\,,\quad
m_{A}^{2}=m_{\ell}^{2}\frac{(1-x^{2})^{2}y^{2}}{(1-y^{2})^{2}x^{2}}\,,$ (3.26)

with a corresponding inverse map

$\displaystyle
x=\frac{\sqrt{4-\sigma_{\ell}}-\sqrt{-\sigma_{\ell}}}{2}\,,\quad
y=\frac{\sqrt{4-\sigma_{A}}-\sqrt{-\sigma_{A}}}{2}\,,\quad\sigma_{i}=\frac{t}{m_{i}^{2}}\,.$
(3.27)

The Euclidean kinematic regime444 It is standard to solve Pfaffian systems for
Euclidean kinematics and then analytically continue the resulting GPLs to
physical kinematic regimes.

$\displaystyle\big{\\{}t<0\big{\\}}\quad\bigcap\quad\big{\\{}0<m_{A}^{2}<m_{\ell}^{2}\quad\text{or}\quad
0<m_{\ell}^{2}<m_{A}^{2}\big{\\}}$ (3.28)

is parametrized by

$\displaystyle\big{\\{}0<x<1\big{\\}}\quad\bigcap\quad\big{\\{}0<y<x\quad\text{or}\quad
0<x<y\big{\\}}$ (3.29)

in the $(x,y)$ coordinates.

The basis (3.22) satisfies a canonical Pfaffian system

$\displaystyle\mathrm{d}\vec{I}(x,y|\epsilon)=\epsilon\,P(x,y)\cdot\vec{I}(x,y|\epsilon)\,,\quad
P(x,y)=\sum_{i=1}^{12}P_{i}\ \mathrm{d}\log\eta_{i}(x,y)\,,\quad
P_{i}\in\mathbb{Q}^{20\times 20}$ (3.30)

with letters

$\displaystyle\eta_{1}$ $\displaystyle=x\,,\quad$ $\displaystyle\eta_{2}$
$\displaystyle=1+x\,,\quad$ $\displaystyle\eta_{3}$ $\displaystyle=1-x\,,$
$\displaystyle\eta_{4}$ $\displaystyle=1+x^{2}\,,\quad$
$\displaystyle\eta_{5}$ $\displaystyle=y\,,\quad$ $\displaystyle\eta_{6}$
$\displaystyle=1+y\,,$ (3.31) $\displaystyle\eta_{7}$
$\displaystyle=1-y\,,\quad$ $\displaystyle\eta_{8}$
$\displaystyle=1+y^{2}\,,\quad$ $\displaystyle\eta_{9}$ $\displaystyle=x+y\,,$
$\displaystyle\eta_{10}$ $\displaystyle=x-y\,,\quad$ $\displaystyle\eta_{11}$
$\displaystyle=1+xy\,,\quad$ $\displaystyle\eta_{12}$ $\displaystyle=1-xy\,.$

The general solution to this Pfaffian system is given by the path-ordered
exponential (2.71). The result takes the form
$\vec{I}(x,y|\epsilon)=\left[\sum_{i=0}^{4}M_{i}(x,y)\,\epsilon^{i}+O(\epsilon^{5})\right]\cdot\vec{I}_{0}(\epsilon)\,,$
where the matrices $M_{i}(x,y)$ contain GPLs built from the letters
$\eta_{i}$, and $\vec{I}_{0}(\epsilon)$ is a vector of boundary constants. The
latter is fixed as follows:

  * •

By normalizing the integration measure as

$\displaystyle\left(\frac{m_{\ell}^{2}}{\mu^{2}}\right)^{\epsilon}\int\frac{\mathrm{d}^{\mathrm{D}}\ell_{i}}{i\pi^{\mathrm{D}/2}\Gamma(1+\epsilon)}\,,\quad
i=1,2\,,$ (3.32)

then direct integration (via Feynman parametrization) of the tadpoles
$\\{I_{1},\,I_{2}\\}$ as well as the factorized integral $I_{6}$ gives

$\displaystyle I_{1}$ $\displaystyle=1$ $\displaystyle I_{2}$
$\displaystyle=\left(\frac{(1-y^{2})x^{2}}{(1-x^{2})y^{2}}\right)^{\epsilon}$
(3.33) $\displaystyle I_{6}$
$\displaystyle=\left(\frac{x^{2}}{(1-x^{2})^{2}}\right)^{\epsilon}\left[1-\zeta_{2}\epsilon^{2}-2\zeta_{3}\epsilon^{3}-\frac{9}{4}\zeta_{4}\epsilon^{4}+O(\epsilon^{5})\right]\,,$

where $\zeta_{n}$ is the $n$th Riemann zeta value. These external inputs are
fed into the DEQ to give relations among the remaining boundary constants.

  * •

Boundary constants for the set of integrals
$\\{I_{3},I_{4},I_{5},I_{9},I_{10},I_{11},I_{12},I_{14},I_{16},I_{17},I_{18},\allowbreak
I_{19}\\}$ are fixed by requiring regularity at the pseudo-threshold limit
$t\to 0$. In particular for
$\\{I_{3},I_{4},I_{5},I_{9},I_{10},I_{14},I_{16},I_{17}\\},$ the prefactors in
(3.22) make them vanish in this limit555 The vanishing of $I_{14}$ as $t\to 0$
follows from the results of [76]. . Further,
$\\{I_{11},I_{12},I_{18},I_{19}\\}$ are seen to vanish as $t\to 0$ by
analyzing their DEQs. Finally, the factorized double-bubble integral $I_{10}$
vanishes in the limit because the _massive_ bubble factor (drawn in  blue in
Figure (3.3)) goes to $0$ at a faster rate than the divergence of the
_massless_ bubble factor (drawn in black).

  * •

Boundary constants for $\\{I_{7},I_{8}\\}$ are determined by regularity at
$m_{A}^{2}\to 0$. Due to the prefactors in (3.22), both integrals vanish in
this limit.

  * •

Regularity at the limit $t\to 4m_{A}^{2}$ fixes boundary constants for the
final three integrals $\\{I_{13},I_{15},I_{20}\\}.$ Thanks to the prefactors
of (3.22), they all vanish in this limit.

  * •

In some of the cases listed above, the equations for the boundary constants
only led to solutions at very high precision, as opposed to analytic
expressions in terms of transcendental constants. Fortunately, using the
implementation of the PSLQ algorithm [77] in PolyLogTools [67], we can
reconstruct the analytic constants by fitting to an ansatz.

##### 3.3.2 Equal-mass case

The equal-mass topology is shown in Figure (3.4). Its associated FI family is
given by (3.10) upon setting $m_{A}^{2}=m_{\ell}^{2}=m^{2}$ in (3.6) and
(3.11).

Figure 3.4: Equal-mass Feynman diagram.

IBP reduction reveals 15 MIs. The Magnus expansion then yields the canonical
basis

$\displaystyle\begin{array}[]{lllllll}&I_{1}&=&\epsilon^{2}J_{1}&I_{2}&=&\epsilon^{2}\lambda_{m}J_{2}\\\
&I_{3}&=&-\epsilon^{2}tJ_{3}&I_{4}&=&\epsilon^{2}\left[\frac{\lambda_{m}+t}{2}J_{3}+\lambda_{m}J_{4}\right]\\\
&I_{5}&=&-\epsilon^{2}tJ_{5}&I_{6}&=&\epsilon^{2}m^{2}J_{6}\\\
&I_{7}&=&-\epsilon^{2}t\lambda_{m}J_{7}&I_{8}&=&\epsilon^{3}\lambda_{m}J_{8}\\\
&I_{9}&=&\epsilon^{2}\left[\frac{4m^{2}-\lambda_{m}-t}{4}(J_{3}+2J_{4})+\
m^{2}(4m^{2}-t)J_{9}\right]&I_{10}&=&\epsilon^{3}\lambda_{m}J_{10}\\\
&I_{11}&=&\epsilon^{3}(1-2\epsilon)tJ_{11}&I_{12}&=&\epsilon^{3}t(t-4m^{2})J_{12}\\\
&I_{13}&=&\epsilon^{3}\lambda_{m}J_{13}&I_{14}&=&\eqref{dm_equal_mass_MI_I14}\\\
&I_{15}&=&-\epsilon^{4}\lambda_{m}tJ_{15}\,,\end{array}$ (3.42)

with

$\displaystyle I_{14}=\epsilon^{2}\left[(4m^{2}-\lambda_{m}-t)(J_{2}-\epsilon
J_{13})+(2\epsilon-1)(4m^{2}-t)J_{14}\right]\,.$ (3.43)

Diagrams for the integrals $J_{i}$ can be seen in Figure (3.5).

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

(m)

(n)

(o)

Figure 3.5:  The MIs $J_{i}$ appearing in the canonical basis (3.42). Squared
propagators are represented by dots.

The canonical integrals contain an algebraic function $\lambda_{m}$ given by

$\displaystyle\lambda_{m}=\sqrt{-t}\sqrt{4m^{2}-t}\,.$ (3.44)

Both of these roots are simultaneously rationalized by the coordinate change

$\displaystyle t=-m^{2}\frac{(1-w)^{2}}{w}\,,$ (3.45)

whose inverse map is given by

$\displaystyle
w=\frac{\sqrt{4m^{2}-t}-\sqrt{-t}}{\sqrt{4m^{2}-t}+\sqrt{-t}}\,.$ (3.46)

The Euclidean kinematic regime

$\displaystyle\big{\\{}t<0\big{\\}}\quad\bigcap\quad\big{\\{}m^{2}>0\big{\\}}$
(3.47)

is then given by

$\displaystyle 0<w<1\,.$ (3.48)

The canonical Pfaffian system for the basis (3.42) takes the form

$\displaystyle\mathrm{d}\vec{I}(w,\epsilon)=\epsilon\,P(w)\cdot\vec{I}(w,\epsilon)\,,\quad
P(w)=\sum_{i=1}^{3}P_{i}\,\mathrm{d}\log\eta_{i}(w)\,,\quad
P_{i}\in\mathbb{Q}^{15\times 15}\,,$ (3.49)

where the three letters are

$\displaystyle\eta_{1}=w\,,\quad\eta_{2}=1+w\,,\quad\eta_{3}=1-w\,.$ (3.50)

This means that the Pfaffian system has a solution in terms _harmonic
polylogarithms_ (HPLs) denoted by $H(\vec{z}|w)$, i.e. GPLs with weights
$z_{i}\in\\{-1,0,1\\}$ [78].

Boundary constants are determined as follows:

  * •

The integrals $\\{I_{1},I_{5},I_{6}\\}$ can be directly integrated, and are
provided as external input. The first integral is normalized to unity by the
measure, and the second integral is swiftly computed with Feynman parameters,
giving

$\displaystyle I_{1}$ $\displaystyle=1$ (3.51) $\displaystyle I_{5}$
$\displaystyle=\left(\frac{w}{(1-w)^{2}}\right)^{\epsilon}\left[1-\zeta_{2}\epsilon^{2}-2\zeta_{3}\epsilon^{3}-\frac{9}{4}\zeta_{4}\epsilon^{4}+O(\epsilon^{5})\right]\,.$
(3.52)

The equal-mass sunrise integral $I_{6}$ is computed by following the strategy
of [79]. Since the square of its external momentum equals its internal masses,
$p^{2}=m^{2}$, this integral evaluates to a number rather than a function,
wherefore it cannot immediately be computed from a Pfaffian system. The trick
is to change two of the propagator masses to $M^{2}$ rather than $m^{2}$, and
then solve a Pfaffian system in the auxiliary variable $x=m^{2}/M^{2}$.
Boundary constants are fixed by regularity at $x=0$. The desired integral is
obtained by setting $x=1$, giving the result

$\displaystyle I_{6}=$
$\displaystyle-\frac{\zeta_{2}}{2}\epsilon^{2}+\left[3\zeta_{2}\log(2)-\frac{7}{4}\zeta_{3}\right]\epsilon^{3}$
(3.53)
$\displaystyle+\left[\frac{31}{4}\zeta_{4}-6\zeta_{2}\log^{2}(2)-\frac{\log^{4}(2)}{2}-12\mathrm{Li}_{4}(\tfrac{1}{2})\right]\epsilon^{4}+O(\epsilon^{5})\,.$

  * •

Boundary values for the integrals
$\\{I_{2},I_{3},I_{4},I_{7},I_{8},I_{9},I_{11},I_{13},I_{14}\\}$ are fixed by
regularity at the pseudo-threshold $t\to 0$. In particular, due to the
prefactors in (3.42), the integrals
$\\{I_{2},I_{3},I_{4},I_{8},I_{11},I_{13}\\}$ go to $0$ in this limit. By
analyzing the Pfaffian system as $t\to 0$, we obtain the relations
$I_{9}|_{t\to 0}=3I_{6}$ and $I_{14}|_{t\to 0}=6I_{6},$ which fixes boundary
values for $I_{9}$ and $I_{14}$ because $I_{6}$ is known by the formula
(3.53). The integral $I_{7}$ vanishes because the _massive_ canonical bubble
factor goes to zero faster than the divergence of _massless_ bubble factor.

  * •

For $I_{10}$, we use regularity at the pseudo-threshold $t\to 4m^{2}$. Thanks
to the prefactor in (3.42), the integral vanishes in this limit.

  * •

Boundary constants for the integral $I_{12}$ follow from comparing the left-
and right-hand sides of $I_{12}=I_{2}\times I_{10},$ since $I_{2}$ and
$I_{10}$ are known. This relation can be gleaned from Figure (3.42) by
recalling that the double-tadpole $I_{1}$ is normalized to unity due to the
integration measure.

  * •

We are able to determine boundary constants for the top sector integral
$I_{15}$ by comparing our result with [80].

The final results for the canonical MIs are quite compact. For instance, the
top sector integral reads

$\displaystyle I_{15}=\epsilon^{4}$
$\displaystyle\Big{[}2\zeta_{3}H(0|w)+4\zeta_{2}H(0,0|w)+8H(0,0,-1,0|w)+H(0,0,0,0|w)$
(3.54)
$\displaystyle+6H(0,0,0,1|w)+4H(0,1,0,0|w)+3\zeta_{4}\Big{]}+O(\epsilon^{5})\,.$

##### 3.3.3 Soft limit

The soft limit $t\to 0$ of the form factors in Section (3.4.3) is of
phenomenological interest. For this reason, we now compute the two-mass MIs in
this limit. The original two-mass topology simplifies in this limit, as
depicted in Figure (3.6).

Figure 3.6: Diagrammatic representation of the soft limit.

The reason for this simplification is as follows. Recall that $t=q^{2}$, where
$q^{\mu}$ is the momentum of the mediator particle. Setting $q^{\mu}=0$ in the
integral family (3.10) implies that $D_{2}\to D_{1}\text{ and }D_{5}\to
D_{4},$ wherefore we only need to consider IBPs and DEQs for the simplified
integral family

$\displaystyle
I_{\nu_{1}\nu_{3}\nu_{4}\nu_{6}\nu_{7}}=\int\frac{\mathrm{d}^{\mathrm{D}}\ell_{1}}{(2\pi)^{\mathrm{D}}}\frac{\mathrm{d}^{\mathrm{D}}\ell_{2}}{(2\pi)^{\mathrm{D}}}\frac{D_{7}^{\nu_{7}}}{D_{1}^{\nu_{1}}D_{3}^{\nu_{3}}D_{4}^{\nu_{4}}D_{6}^{\nu_{6}}}$
(3.55)

having four denominators

$\displaystyle\begin{array}[]{lllllll}&D_{1}&=&\ell_{1}^{2}-m_{\ell}^{2}\,,&D_{3}&=&\ell_{2}^{2}-m_{A}^{2}\,,\\\
&D_{4}&=&(\ell_{2}+p)^{2}\,,&D_{6}&=&(\ell_{1}+\ell_{2}+p)^{2}-m_{\ell}^{2}\,,\end{array}$
(3.58)

and an ISP $D_{7}=(\ell_{1}-p)^{2}$. Note that $p^{\prime}$ has been set equal
to $p$ due to momentum conservation.

IBP reduction for the family (3.55) reveals 4 MIs. By the Magnus exponential,
we obtain the canonical basis

$\displaystyle\begin{array}[]{lllllll}&I_{1}&=&\epsilon^{2}J_{1}\,,&I_{2}&=&\epsilon^{2}J_{2}\,,\\\
&I_{3}&=&\epsilon^{2}m_{\ell}m_{A}\big{[}J_{3}+2J_{4}\big{]}\,,&I_{4}&=&\epsilon^{2}m_{A}^{2}J_{4}\,,\end{array}$
(3.61)

in terms of the integrals $\\{J_{1},J_{2},J_{3},J_{4}\\}$ shown in Figure
(3.7).

(a)

(b)

(c)

(d)

Figure 3.7:  The integrals $J_{1}$ appearing in the canonical basis (3.61).
Squared propagators are depicted by dots.

Defining the mass-ratio variable

$\displaystyle z=\frac{m_{A}}{m_{\ell}}\,,$ (3.62)

the canonical Pfaffian system takes the form

$\displaystyle\mathrm{d}\vec{I}(z|\epsilon)=\epsilon\,P(z)\cdot\vec{I}(z|\epsilon)\,,\quad
P(z)=\sum_{i=1}^{3}P_{i}\,\mathrm{d}\log\eta_{i}(z)\,,\quad
P_{i}\in\mathbb{Q}^{4\times 4}\,,$ (3.63)

where the $\eta_{i}$ are HPL letters

$\displaystyle\eta_{1}=z\,,\quad\eta_{2}=1+z\,,\quad\eta_{3}=1-z\,.$ (3.64)

Boundary constants for the solution vector are fixed as follows:

  * •

The two double-tadpoles $\\{I_{1},I_{2}\\}$ can be directly integrated to

$\displaystyle I_{1}=1\quad\text{and}\quad I_{2}=z^{-2\epsilon}\,.$ (3.65)

  * •

Boundary values for the two sunrises $\\{I_{3},I_{4}\\}$ are fixed by
regularity at $m_{A}\to 0$ as in [79]. In particular, due to the prefactors of
(3.61), both integrals vanish in this limit.

####  Section 3.4Form factor results

We now present form factor results for the following cases:

  * •

Two masses ($m_{A}\neq m_{\ell}$): scalar and pseudo-scalar.

  * •

Equal-mass ($m_{A}=m_{\ell}$): scalar and pseudo-scalar.

  * •

Soft limit ($t\to 0$): scalar and pseudo-scalar.

Let us emphasize again that we only include results for one-body interactions.
The analytic form factors in this section are all rescaled as
$\mathcal{F}_{S,P}\to\frac{\pi^{2}}{\alpha_{\mathrm{em}}}\mathcal{F}_{S,P}$
for the sake of readability.

##### 3.4.1 Two-mass case

The analytic expressions for the two-mass form factors
$\mathcal{F}_{S,P}(t,m_{\ell},m_{A})$ would span several pages, so we omit
them here. Instead, we showcase the results as numerical plots in Figures
(3.8(a)) and (3.8(b)).

(a) Two-mass pseudoscalar form factor as a function of momentum transfer $t$.
The nucleon mass $m_{A}$ is set equal to the proton mass. The dashed lines
represent individual contributions from the electron ($e$), the muon ($\mu$),
and the tau lepton ($\tau$). The solid  red line shows the sum of all
contributions. The grey shaded band indicates the range of $t$ that is
expected to be relevant for DM particles arriving from our local galaxy.

(b) Two-mass scalar form factor as a function of momentum transfer $t$. See
the caption of Figure (3.8(a)) for additional details.

##### 3.4.2 Equal-mass case

Define the shorthand notation $G(\vec{z}|w):=G_{\vec{z}}(w).$

###### Scalar form factor

The analytic expression for the equal-mass form factor with a scalar mediator
particle is

$\displaystyle\mathcal{F}_{S}(w)={-}\frac{(20\,w^{2}{+}4\,w)}{\,(w{-}1)\,(w{+}1)^{2}}\,G_{0}(w)\;{+}\;\frac{w}{\,(w{-}1)\,(w{+}1)}\,\big{[}24\,G_{{-}1,0}(w){+}4\,G_{1,0}(w)\big{]}$
$\displaystyle{+}\frac{w}{\,(w{-}1)^{2}\,(w{+}1)^{3}}[{-}12\,w^{3}{-}38\,w^{2}{+}8\,w{+}2]\,G_{0,0}(w)\;{-}\;\frac{w}{\pi^{2}\,(w{-}1)^{2}}\,\big{[}2\,G_{0,0,0}(w){+}4\,\zeta_{3}\big{]}$
$\displaystyle{+}\frac{w\,\zeta_{3}}{\,(w{-}1)^{3}\,(w{+}1)}\,[5\,w^{2}{-}6\,w{+}5]\,G_{0}(w)\;{+}\;\frac{w}{\pi^{2}\,(w{+}1)^{2}}\,\left[12\,G_{0,0,1}(w){-}\frac{4\,\pi^{2}}{3}\,G_{0}(w){+}16\,G_{1}(w)\right]$
$\displaystyle{+}\frac{w}{\,(w{-}1)^{2}\,(w{+}1)^{2}}\,\Big{[}({-}16\,w^{2}{+}64\,w{-}16)\,G_{0,{-}1,0}(w){+}(8\,w^{2}{-}48\,w{+}8)\,G_{1,0,0}(w){+}16\,w\,G_{0,1,0}(w)\Big{]}$
$\displaystyle{+}\frac{w}{\,(w{-}1)\,(w{+}1)^{3}}\Bigg{[}({-}14\,w^{2}{+}4\,w{-}14)\,G_{0,1}(w){+}\left(\frac{w^{2}}{3}{+}2\,w{+}\frac{1}{3}\right)\,G_{1,0}(w)\,\pi^{2}$
(3.66)
$\displaystyle{+}(2\,w^{2}{+}12\,w{+}2)\,\big{(}G_{1,0,0,0}(w)\,{-}\,G_{0,0,1,0}(w)\big{)}{+}\left(2\,w^{2}{-}\frac{20\,w}{3}{+}2\right)\pi^{2}\Bigg{]}$
$\displaystyle{+}\frac{w}{\,(w{-}1)^{3}\,(w{+}1)^{3}}\Bigg{[}\left(\frac{5\,w^{4}}{6}{-}2\,w^{3}{+}\frac{23\,w^{2}}{3}{-}2\,w{+}\frac{5}{6}\right)\,G_{0,0}(w)\,\pi^{2}$
$\displaystyle{+}\left(8\,w^{4}{-}32\,w^{3}{+}112\,w^{2}{-}32\,w{+}8\right)\,G_{0,0,{-}1,0}(w){+}\left(\frac{3\,w^{4}}{2}{-}2\,w^{3}{+}9\,w^{2}{-}2\,w{+}\frac{3}{2}\right)\,G_{0,0,0,0}(w)$
$\displaystyle{+}\left({-}9\,w^{4}{+}12\,w^{3}{-}54\,w^{2}{+}12\,w{-}9\right)\,G_{0,0,0,1}(w){+}\left({-}2\,w^{4}{+}24\,w^{3}{-}76\,w^{2}{+}24\,w{-}2\right)\,G_{0,1,0,0}(w)$
$\displaystyle{+}\left(\frac{w^{4}}{18}{-}\frac{2\,w^{3}}{45}{+}\frac{11\,w^{2}}{45}{-}\frac{2\,w}{45}{+}\frac{1}{18}\right)\,\pi^{4}\Bigg{]}.$

This expression can be turned into an $s$-channel form factor by crossing
symmetry. Then it is related to a form factor for Higgs decay into two heavy
quarks [81].

###### Pseudo-scalar form factor

The equal-mass form factor for the pseudoscalar reads

$\displaystyle\mathcal{F}_{P}(w)={-}\frac{w^{2}}{\,(w{-}1)^{3}\,(w{+}1)}\,\left[\frac{2\pi^{2}}{3}\,G_{0}(w){+}4\,G_{0,0,0}(w)\right]$
$\displaystyle{+}\frac{w}{\,(w{-}1)^{2}}\,\left[G_{0,0}(w){-}4\,G_{0,0,1}(w){-}4\,G_{0,1,0}(w){+}8\,G_{1,0,0}(w){+}12\,\zeta_{3}{+}\frac{\pi^{2}}{3}\right]$
(3.67)
$\displaystyle{+}\frac{w}{\,(w{-}1)\,(w{+}1)}\,\left[{-}3\,\zeta_{3}\,G_{0}(w){-}\frac{\pi^{2}}{6}\,G_{0,0}(w){-}\frac{\pi^{2}}{3}\,G_{1,0}(w){-}\frac{1}{2}\,G_{0,0,0,0}(w)\right.$
$\displaystyle\hskip
89.626pt\left.{+}3\,G_{0,0,0,1}(w){+}2\,G_{0,0,1,0}(w){-}2\,G_{0,1,0,0}(w){-}2\,G_{1,0,0,0}(w){-}\frac{\pi^{4}}{45}\right].$

Using crossing symmetry, this form factor is related to the $s$-channel QCD
form factor dubbed $A_{R}(s,m,m)$ in [74].

##### 3.4.3 Soft limit

The soft limit $t\to 0$ is important for phenomenological applications, as the
DM particles are expected to scatter elastically. It is therefore useful to
have simple expressions for the form factors in this limit that can be easily
numerically evaluated.

###### Scalar form factor

The scalar form factor is finite at $t=0$, so the result
$\mathcal{F}_{S}(t=0,m_{\ell},m_{A})$ will only depend on the mass ratio
$z=m_{A}/m_{\ell}$. Writing the form factor as an expansion in terms of FIs as
in (3.9), we then IBP reduce those integrals to the soft limit MIs from
Section (3.3.3). The smooth $4$-dimensional limit $\epsilon\to 0$ finally
yields the result

$\displaystyle\mathcal{F}_{S}(t=0,m_{\ell},m_{A})$
$\displaystyle=-\frac{2}{z}\left[1-\frac{\log(z)}{2}+f_{S}(z)+f_{S}(-z)\right]$
(3.68) $\displaystyle\text{with}\quad f_{S}(z)$
$\displaystyle=\frac{4+3z+z^{3}}{4z^{2}}\big{[}\log|z|\log(1+z)+\mathrm{Li}_{2}(-z)\big{]}\,.$

###### Pseudo-scalar form factor

The pseudo-scalar form factor is logarithmically divergent as $t\to 0$, so we
compute an asymptotic expansion in terms of the variables $(x,y)$ from (3.26).
Keeping $m_{A,\ell}$ fixed, the soft limit then corresponds to
$(x,y)\to(1,1)$.

Our strategy for computing the asymptotic expansion is as follows:

  1. 1.

Staring with the full expression for the form factor
$\mathcal{F}_{P}(t,m_{\ell},m_{A})$ in terms of GPLs, we use PolyLogTools to
rewrite every $G(\vec{z}\,|\,v)$ such that $\vec{z}\in\mathbb{Z}^{M}$ and
$v=x$ or $v=y$. The weight vectors $\vec{z}$ will thereby not contain any
$x$’s or $y$’s, which makes it easier to expand the GPLs.

  2. 2.

PolyLogTools has a built-in function for asymptotically expanding GPLs.
However, it only supports an expansion around $v\sim 0$, not $v\sim 1$.
Introducing a dummy variable $v^{\prime}$ by $v=1-v^{\prime}$, we hence seek
to rewrite every $G(\vec{z}\,|\,v)$ in terms of $G(\vec{z}\,|\,v^{\prime})$.
This is done by employing the well-known recursion

$\displaystyle G(z_{1},\ldots,z_{M}|v)$
$\displaystyle=G(z_{1},\ldots,z_{M}|1-v^{\prime})$ (3.69)
$\displaystyle=G(z_{1},\ldots,z_{M}|1)+\int_{0}^{v^{\prime}}\frac{\mathrm{d}t}{t-(1-z_{1})}G(z_{2},\ldots,z_{M}|1-t)\,,$

which produces GPLs having argument $v^{\prime}$ due to the upper integration
bound in the second term. Experimentally, we find that this recursion admits a
closed-form solution

$\displaystyle G(z_{1},\ldots,z_{M}|1-v^{\prime})=\sum_{i=0}^{M}$
$\displaystyle G(1-z_{1},\ldots,1-z_{i}|v^{\prime})$ (3.70)
$\displaystyle\times$ $\displaystyle
G(z_{i+1},\ldots,z_{M}|1)\,\big{|}_{\text{shuffle reg.}}\,,$

where we emphasized the need for _shuffle regularization_ to deal with
divergent GPLs [61, Section 5.3]. We use the convention
$G(1-z_{1},\ldots,1-z_{i}|v^{\prime})=1$ when $i=0$.

As an aside, although it is not a deep result, we have not been able to locate
this closed-form solution elsewhere in the literature.

  3. 3.

PolyLogTools can now expand the terms $G(1-z_{1},\ldots,1-z_{i}|v^{\prime})$
from (3.70) around $v^{\prime}\sim 0$ up the desired order
$\mathcal{O}\big{(}(1-x)^{n},(1-y)^{m}\big{)}.$ The latter notation means that
we first expand in $(1-x)$ up to $n$th order, and then in $(1-y)$ up to $m$th
order. The reversed expansion order is written as
$\mathcal{O}\big{(}(1-y)^{n},(1-x)^{m}\big{)}.$ The maximal powers $(n,m)$ are
chosen so as to give sufficiently accurate numerical results in comparison
with the full form factor $\mathcal{F}_{P}(t,m_{\ell},m_{A})$ when the latter
is evaluated for small $t$.

A double asymptotic expansion in $(x,y)$ will depend on the orderings
$x<y\text{ or }y<x,$ so we treat these two cases separately. The orderings
correspond to the mass hierarchies $m_{\ell}<m_{A}\text{ and }m_{A}<m_{\ell}$
respectively. By the steps outlined above, we obtain the following results for
the pseudoscalar form factor in the soft limit:

  * •

Mass hierarchy $m_{\ell}<m_{A}\,(x<y)$:

$\displaystyle\mathcal{F}_{P}(t\sim 0,\,m_{\ell}<m_{A})=-$ $\displaystyle
2\zeta_{2}+\frac{(1-x)}{(1-y)}\big{[}-2\log(1-x)-2\log(2)+3\big{]}$
$\displaystyle+$
$\displaystyle\frac{(1-x)^{2}}{(1-y)}\left[-\log(1-x)-\log(2)+\frac{1}{2}\right]$
$\displaystyle+$
$\displaystyle(1-x)^{2}\left[\frac{\log(1-x)}{2}+\frac{\log(2)}{2}+\frac{37\zeta_{2}}{60}-\frac{1}{4}\right]$
(3.71) $\displaystyle+$
$\displaystyle\frac{6(1-x)}{\zeta_{2}}\left[\log(1-x)+\log(2)+\frac{3\zeta_{2}}{2}-\frac{3}{2}\right]$
$\displaystyle+$ $\displaystyle\mathcal{O}\big{(}(1-y),\,(1-x)^{3}\big{)}\,.$

  * •

Mass hierarchy $m_{A}<m_{\ell}\,(y<x)$:

$\displaystyle\mathcal{F}_{P}(t\sim 0,\,m_{A}<m_{\ell})=$
$\displaystyle\frac{(1-x)}{(1-y)}\left[-\frac{3\log(1-x)}{2}-\frac{(1-y)}{2}-2\log(2)+\frac{17}{8}\right]$
(3.72) $\displaystyle+$
$\displaystyle(1-x)\left[\frac{3\log(1-x)}{4}+\frac{\log(1-y)}{4}+\log(2)+\frac{3\zeta_{2}}{2}-\frac{21}{16}\right]$
$\displaystyle+$ $\displaystyle\mathcal{O}\big{(}(1-x)^{2},\,(1-y)\big{)}\,.$

These formulas differ by roughly $1\%$ in comparison with the full analytic
results evaluated at $t=-10^{-7}\text{GeV}^{2}$.

##  Part II $\mathcal{D}$-modules

### Chapter 4Introduction to $\mathcal{D}$-modules

The preceding chapters of this thesis dealt with tried and tested methods for
evaluating FIs. We now venture into somewhat uncharted territory from the
perspective of the physics literature.

Based on the papers [82, 83], the next six chapters of this thesis study FIs
within the framework of $\mathcal{D}$-_modules_. A $\mathcal{D}$-module can,
informally, be thought of as a vector space of linear differential operators
that annihilate a given function. As a supplement to the tools coming from
complex _analysis_ , $\mathcal{D}$-modules provide _algebraic_ machinery to
uncover properties of PDEs. For instance, one can learn the dimension of the
solution space, construct bases of differential operators, obtain relations
between basis elements, and more. We will only scratch the surface of this
vast topic, so we refer to the classic textbooks [23] and [24] for exhaustive
discussions on $\mathcal{D}$-modules.

Our overarching goal in Part II of this thesis is to build a bridge between
these two independently well-studied topics, namely FIs and
$\mathcal{D}$-modules. But, as it so often happens in blue skies research,
practical applications will spring out as well.

####  Section 4.1What is a $\mathcal{D}$-module?

How might one represent a function? Take $\sin(x)$ as an example. It can be
represented in several ways, each one having different advantages. For
instance, it can be analytically represented as a power series, or
geometrically via the unit circle. In the philosophy of $\mathcal{D}$-modules,
$\sin(x)$ is represented by the set of linear PDEs which annihilate the
function; in this case simply $\\{\partial^{2}+1\\}.$

We note that this philosophy would not work for non-linear PDEs, since a
differential operator $L$ cannot be "pulled out" in the form $L\bullet f=0$ if
there is a term like $f(x)^{2}$.

Example 4.1. A striking example of this philosophy comes from calculating
roots of an $N$th order polynomial

$\displaystyle
z_{0}+z_{1}x+z_{2}x^{2}+\cdots+z_{N-1}x^{N-1}+z_{N}x^{N}=(x-r_{1})\cdots(x-r_{N})\,.$
(4.1)

Is there an equation for the roots $r_{k}=r_{k}(z_{1},\ldots,z_{N})$ in terms
of the indeterminate polynomial coefficients $z_{k}$? There is, for instance,
a well-known formula for $r_{1,2}(z_{0},z_{1},z_{2})$ in the quadratic case
$N=2$.

Alas, Galois theory teaches us that no such formula exists in terms of
radicals when $N\geq 5$. But if we settle for an infinite series
representation of $r_{k}$, then such a formula _does_ exist!

Fixing $k$, it turns out that $r_{k}$ is annihilated, indeed _represented_ ,
by the following set of partial differential operators:

$\displaystyle\frac{\partial^{2}}{\partial z_{i}\partial
z_{j}}-\frac{\partial^{2}}{\partial z_{k}\partial z_{l}}\quad$
$\displaystyle\text{whenever}\quad i+j=k+l\,,$
$\displaystyle\sum_{i=0}^{N}i\,z_{i}\frac{\partial}{\partial z_{i}}+1\quad$
$\displaystyle\hskip 11.38092pt\text{and}\hskip
11.38092pt\quad\sum_{i=0}^{N}z_{i}\frac{\partial}{\partial z_{i}}\,.$ (4.2)

Many features of polynomial roots can then be learned from studying these
differential operators. For one, it is possible to construct power series
solutions of the form $r_{k}=\sum_{i_{0}\ldots
i_{N}}^{\infty}f_{k}(i_{0},\ldots,i_{N})z_{0}^{i_{0}}\cdots z_{N}^{i_{N}},$
where the $f_{k}$’s are various combinations of factorials [84].
$\blacksquare$

So, what is a $\mathcal{D}$-module? It is a module over a ring of
$\mathcal{D}$ifferential operators. Let us spell this out more formally.

Let $R$ be a ring with a unit $"1"$, and $M$ an abelian group w.r.t. addition
$"+"$. Given $r\in R$ and $m\in M$, we assume the existence of a map

$\displaystyle\psi:R\times M\to M$ $\displaystyle\psi(r,m)=r\bullet m\,,$
(4.3)

where $r\bullet m$ denotes the _action_ of $r$ on $m$. For instance, if
$R=\mathbb{C}$ and $M=\mathbb{C}^{n}$, then $\psi$ could represent the scalar
multiplication of a complex number $r$ on a complex vector $m$.

######  Definition 4.1.

An abelian group $M$ becomes a _left_ 111 _Right_ $R$-modules are defined by
the same axioms with $r$’s acting from the right, i.e. $m\bullet r$.
$R$-_module_ when the following three axioms are fulfilled222 We omit $"+"$
subscripts of kind $(r_{1}+_{R}r_{2})\bullet m=r_{1}\bullet m+_{M}r_{2}\bullet
m$ for notational cleanliness. :

  1. 1.

For any $r\in R$ and $m_{i}\in M$, we have $r\bullet(m_{1}+m_{2})=r\bullet
m_{1}+r\bullet m_{2}.$

  2. 2.

For any $r_{i}\in R$ and $m\in M$, we have $(r_{1}r_{2})\bullet
m=r_{1}\bullet(r_{2}\bullet m)$ and $(r_{1}+r_{2})\bullet m=r_{1}\bullet
m+r_{2}\bullet m.$

  3. 3.

For any $m\in M$, we have $1\bullet m=m.$

If $R$ is a commutative field such as $\mathbb{C}$, then a left $R$-module is
nothing but an $R$-vector space.

The $N$th _Weyl algebra_

$\displaystyle\mathcal{D}_{N}=\mathbb{C}[z_{1},\ldots,z_{N}]\langle\partial_{1},\ldots,\partial_{N}\rangle$
(4.4)

is generated by symbols $z_{i}$ and $\partial_{i}$ subject to the commutation
relations

$\displaystyle[z_{i},z_{j}]=0\,,\quad[\partial_{i},\partial_{j}]=0\,,\quad[\partial_{i},z_{j}]=\delta_{ij}\,.$
(4.5)

We shall often suppress the subscript in $\mathcal{D}_{N}$ when the number of
variables is understood. The square brackets $[\ldots]$ in (4.4) mean that the
$z_{i}$ are thought of as coefficients of the generators $\partial_{i}$. Note
that the last commutator of (4.5) encodes Leibniz’ rule of differentiation
when we insert a dummy function $f(z)$: $\frac{\partial}{\partial
z}[zf(z)]-z\frac{\partial}{\partial z}f(z)=f(z).$ We therefore interpret
$\mathcal{D}$ as the collection of all differential operators having
polynomial coefficients. For instance,
$(z_{1}^{2}+z_{3}+1)\partial_{1}^{2}\in\mathcal{D}_{3}.$ We emphasize that
$\mathcal{D}$ has the structure of a ring.

We also have the $N$th _rational_ Weyl algebra

$\displaystyle\mathcal{R}_{N}=\mathbb{C}(z_{1},\ldots,z_{N})\langle\partial_{1},\ldots,\partial_{N}\rangle\,,$
(4.6)

wherein the coefficients are promoted to rational functions. We shall often
shorten $\mathcal{R}_{N}$ to $\mathcal{R}$. In this case, the commutators are
extended to

$\displaystyle[z_{i},z_{j}]=0\,,\quad[\partial_{i},\partial_{j}]=0\,,\quad[\partial_{i},r(z)]=\frac{\partial
r(z)}{\partial z_{i}}\,,$ (4.7)

for a rational function $r(z)\in\mathbb{C}(z_{1},\ldots,z_{N}).$ The rational
Weyl algebra is also endowed with a ring structure.

######  Definition 4.2.

A $\mathcal{D}$-module is a left $R$-module (see Definition (4.1)) with the
specific choice of ring $R=\mathcal{D}$. An $\mathcal{R}$-module is defined
similarly for $R=\mathcal{R}$.

We shall be careful in distinguishing between ring _multiplication_ and ring
_action_ $"\bullet"$, because there is a difference between e.g.

$\displaystyle\partial_{i}z_{i}=z_{i}\partial_{i}+1\quad\text{and}\quad\partial_{i}\bullet
z_{i}=1\,.$ (4.8)

Definition (4.2) can be instantiated in many different ways because it only
specifies the ring $R$. _Particular_ $\mathcal{D}$\- and $\mathcal{R}$-modules
follow from further specifying the abelian group $M$ and the action map
$\psi$. Let us list some examples.

Example 4.2. The Weyl algebras $\mathcal{D}$ and $\mathcal{R}$ are themselves
$\mathcal{D}$\- and $\mathcal{R}$-modules respectively. The abelian group $M$
is formed from the addition of differential operators. The algebras acts on
themselves according to (4.5) in the case of $\mathcal{D}$ and (4.7) in the
case of $\mathcal{R}$, which specifies $\psi$. $\blacksquare$

Example 4.3. Many function spaces are equipped with the structure of a
$\mathcal{D}$-module. Take $M$ to be some suitable space of holomorphic
functions on a domain of $\mathbb{C}^{n}$. Assume that $M$ is closed under
differentiation, meaning that if $f(z)\in M$ then also $\frac{\partial
f}{\partial z_{i}}\in M$. The addition of such functions indeed forms an
abelian group. The action map $\psi:R\times M\to M$ is defined by

$\displaystyle\psi(\partial_{i},f)=\partial_{i}\bullet f:=\frac{\partial
f}{\partial z_{i}}\,,\quad\psi(z_{i},f)=z_{i}\bullet f:=z_{i}\cdot f\,,$ (4.9)

where $"\cdot"$ in the last equation stands for ordinary scalar
multiplication. $\blacksquare$

Example 4.4. The most important example of a $\mathcal{D}$-module (or
$\mathcal{R}$-module) for the study of PDEs is arguably that of the _quotient_
ring $\mathcal{D}/\mathcal{I}$ (or $\mathcal{R}/\mathcal{I}$), where
$\mathcal{I}$ is a so-called _left ideal_ to be defined shortly. The abelian
group $M$ is given by the addition of operators, and the action is that of the
Weyl algebra.

To motivate the notion of a left ideal $\mathcal{I}$, consider a collection of
differential operators $\\{L_{1},L_{2},\ldots\\}$ that annihilate some
function $f$: $L_{1}\bullet f=L_{2}\bullet f=\ldots=0.$ Note that we still get
zero if we act with any other operator $L\in\mathcal{D}$ from the left:
$LL_{1}\bullet f=LL_{2}\bullet f=\ldots=0.$ The definition333 This is more of
a ”physicist’s definition”. Formally, a left ideal $\mathcal{I}$ of a ring $R$
satisfies 1) it is a subgroup under addition, and 2) for $r\in R$ and
$L\in\mathcal{I}$ we have $r\,L\in\mathcal{I}$.  of a left ideal $\mathcal{I}$
in the Weyl algebra encapsulates this fact:

$\displaystyle\mathcal{I}=\langle
L_{1},L_{2},\ldots\rangle=\\{LL_{1},\,LL_{2},\ldots\ |\ \forall\
L\in\mathcal{D}\\}\,.$ (4.10)

The notation $\langle L_{1},L_{2},\ldots\rangle$ means that the $L$’s
_generate_ $\mathcal{I}$.

The quotient ring $\mathcal{D}/\mathcal{I}$ is then given by all operators
$L\in\mathcal{D}$ modulo the equivalence relation $\sim$ determined by

$\displaystyle L\sim K\iff L=K+K_{\mathcal{I}}\quad\text{where}\quad
K\in\mathcal{D}\,,\ K_{\mathcal{I}}\in\mathcal{I}\,.$ (4.11)

Take, for instance, the left ideal generated by
$\mathcal{I}=\langle\partial_{1}\partial_{2}\rangle\,.$ Define two operators
by

$\displaystyle L\ $ $\displaystyle=\ \partial_{2}\quad\text{and}\quad$
$\displaystyle K\ $ $\displaystyle=\ \partial_{1}z_{1}\partial_{2}\
\overset{\eqref{Weyl_algebra_commutators}}{=}\
(1+z_{1}\partial_{1})\partial_{2}\ =\
\partial_{2}+z_{1}\partial_{1}\partial_{2}\,.$ (4.12)

$L$ and $K$ represent _the same_ element inside $\mathcal{D}_{2}/\mathcal{I}$
because they only differ by a term
$z_{1}\partial_{1}\partial_{2}\in\mathcal{I}$. We say that $L$ and $K$ lie in
the same _equivalence class_ , written as $[L]=[K].$

We shall see that the algebraic structure of the $\mathcal{D}$-module
$\mathcal{D}/\mathcal{I}$ carries a wealth of information about the functions
that are annihilated by the ideal $\mathcal{I}$. This $\mathcal{D}$-module
hence embodies the intuition, argued for in the beginning of this section,
that a function can be represented by the differential operators that
annihilate it. $\blacksquare$

####  Section 4.2Holonomic $\mathcal{D}$-modules

A $\mathcal{D}$-module $\mathcal{D}/\mathcal{I}$ has particularly nice
features when it is _holonomic_. Holonomic $\mathcal{D}$-modules describe FIs,
trigonometric functions, various probability distributions, hypergeometric
functions, Bessel functions, and GPLs. It is necessary to introduce a few
concepts before we can define the notion of holonomicity (we refer to the
seminal text [85] for many more details).

  * •

Using the Weyl algebra commutation relations, any operator $L\in\mathcal{D}$
can be uniquely expressed in the _normally ordered form_

$\displaystyle
L\quad=\quad\sum_{(p,q)\,\in\,\mathrm{supp}\left(L\right)}c_{p,q}\,z^{p}\,\partial^{q}\,,\quad
c_{p,q}\in\mathbb{C}\,.$ (4.13)

The summation range $\mathrm{supp}\left(L\right)$ denotes a finite set of non-
zero multi-indices $(p,q)=(p_{1},\ldots,p_{N},q_{1},\ldots,q_{N}),$ and we
introduced the multi-index exponent notation

$\displaystyle z^{p}=z_{1}^{p_{1}}\cdots
z_{N}^{p_{N}}\quad\text{and}\quad\partial^{q}=\partial_{1}^{q_{1}}\cdots\partial_{N}^{q_{N}}\,,$
(4.14)

which will be used frequently throughout this text.

  * •

For a given a vector $(u,v)\in\mathbb{R}^{2N}$, the _$(u,v)$ -weight_ of
$z^{p}\,\partial^{q}$ is defined by $u\cdot p+v\cdot q.$

  * •

The differential operator $\mathrm{in}_{(u,v)}(L)$ is constructed by only
keeping the largest $(u,v)$-weight terms in $L$.

  * •

The _initial form_ is defined by $\mathrm{in}_{(u,v)}(L)$ for the specific
choice $u=(0,\ldots,0)=\mathbf{0}$ and $v=(1,\ldots,1)=\mathbf{1}.$ I.e. one
assigns weight $0$ to all $z$’s, and weight $1$ to all $\partial$’s. The
standard notation for the initial form is

$\displaystyle\mathrm{in}_{(\mathbf{0},\mathbf{1})}(L)=\sum_{\begin{subarray}{c}(p,q)\,\in\,\mathrm{supp}\left(L\right),\\\
|q|\text{ is max}\end{subarray}}c_{p,q}\,z^{p}\,\xi^{q}\,,$ (4.15)

where $|q|=q_{1}+\ldots+q_{N}.$ The $\partial$’s have here been replaced by
commuting $\xi$ variables, but this replacement may only be done _after_
writing $L$ in normally ordered form. As an example, take the normally ordered
operator

$\displaystyle
L=z_{1}^{2}+2\partial_{1}^{2}\partial_{2}+z_{2}^{3}\partial_{1}\partial_{2}+z_{1}z_{2}\partial_{1}^{3}.$
(4.16)

From left to right, the terms have weights $(0,3,2,3).$ So we only keep the
2nd and 4th terms, giving

$\displaystyle\mathrm{in}_{(\mathbf{0},\mathbf{1})}(L)=2\xi_{1}^{2}\xi_{2}+z_{1}z_{2}\xi_{1}^{3}\,.$
(4.17)

  * •

Given a left ideal $\mathcal{I}$ in $\mathcal{D}_{N}$, we set

$\displaystyle\mathrm{in}_{(\mathbf{0},\mathbf{1})}(\mathcal{I})=\mathbb{C}\cdot\\{\mathrm{in}_{(\mathbf{0},\mathbf{1})}(L)\
|\ L\in\mathcal{I}\\}\,,$ (4.18)

i.e. the ideal generated by the highest
$(\mathbf{0},\mathbf{1})\text{-weight}$ terms in $\mathcal{I}$. It is a
_polynomial_ ideal since it is written in terms of commuting variables
$(z,\xi)$.

  * •

The _characteristic variety_ of an ideal is defined as

$\displaystyle\mathrm{char}\\!\left(\mathcal{I}\right)=\left\\{(z,\xi)\ |\
p(z,\xi)=0\ \text{ for all }\
p\,\in\,\mathrm{in}_{(\mathbf{0},\mathbf{1})}(\mathcal{I})\right\\}\ \subset\
\mathbb{C}^{2N}\,.$ (4.19)

Namely, it the zero-set of all the polynomials in the ideal
$\mathrm{in}_{(\mathbf{0},\mathbf{1})}(\mathcal{I}).$ The computation of
$\mathrm{char}\\!\left(\mathcal{I}\right)$ typically requires knowing a
_Gröbner basis_ for $\mathcal{I}$ (to be defined in the next section).

With all of these notions in place, we can finally state

######  Definition 4.3.

An ideal $\mathcal{I}$ is _holonomic_ when its characteristic variety
$\mathrm{char}\\!\left(\mathcal{I}\right)$ has dimension444 According to a
result called Bernstein’s inequality, all connected components of
$\mathrm{char}\\!\left(\mathcal{I}\right)$ have dimension at least $N$. An
ideal is hence holonomic when its characteristic variety has the smallest
possible dimension (the maximal dimension would be $2N$).  $N$. A
$\mathcal{D}$-module $\mathcal{D}/\mathcal{I}$ is holonomic if and only if
$\mathcal{I}$ is holonomic555 There is actually an independent definition of
holonomicity for $\mathcal{D}$-modules, depending on so-called Bernstein
filtrations. The fact that holonomic $\mathcal{D}$-modules coincide with
holonomic ideals is then a non-trivial theorem. See [86, Chapter 6] for
details. . The _holonomic rank_ is defined as the following vector space
dimension over $\mathbb{C}(z)$:

$\displaystyle\mathrm{rank}\\!\left(\mathcal{I}\right)$
$\displaystyle=\dim_{\mathbb{C}(z)}\big{(}\mathbb{C}(z)[\xi]\,/\,\mathbb{C}(z)[\xi]\cdot\mathrm{char}\\!\left(\mathcal{I}\right)\big{)}\,.$

The above definition of $\mathrm{rank}\\!\left(\mathcal{I}\right)$ relies on
the following important fact: a holonomic $\mathcal{D}$-module carries the
structure of a _finite-dimensional_ vector space [86, 85]. It follows that
there exists a _basis_ of differential operators for the vector space
$\mathcal{D}/\mathcal{I}$ having exactly
$\mathrm{rank}\\!\left(\mathcal{I}\right)$ many elements. Quite remarkably,
$\mathrm{rank}\\!\left(\mathcal{I}\right)$ also counts the number of
independent holomorphic solutions to the system of PDEs defined by
$\mathcal{I}$. When $\mathcal{I}$ is the annihilating ideal for a FI, then
$\mathrm{rank}\\!\left(\mathcal{I}\right)$ is indeed the number of MIs!

##### 4.2.1 Regular holonomic

Throughout the rest of this thesis, when we say an ideal $\mathcal{I}$ is
holonomic we shall always mean _regular holonomic_. The extra adjective
"regular" signifies, loosely speaking, that the solutions to the PDEs dictated
by $\mathcal{I}$ only have mild singularities.

Suppose that $z=0$ is a singularity of a solution $f(z)$ in the univariate
case. It is a _regular singularity_ if there exists a $k\in\mathbb{Z}_{>0}$
such that $|z|^{k}f(z)\to 0$ as $z\to 0$ at a fixed angle $\mathrm{arg}(z)$ in
the complex plane.

The definition of regular holonomic is much more complicated in the
multivariate case. We refer to [85, Section 2.4] [23, Chapter 6] for details.
A working definition for our purposes is to say that the solutions to
$\mathcal{I}$ admit logarithmic series expansions near the singular
hypersurfaces. Equivalently, the Pfaffian system born from $\mathcal{I}$, to
be discussed in Section (4.4), will have at most simple poles (higher-order
poles can appear, but they are "fake" in the sense they can be turned into
simple poles via a gauge transformation).

####  Section 4.3Gröbner bases

How does one compute a basis for a holonomic $\mathcal{D}$\- or
$\mathcal{R}$-module? The answer, as with almost any computational question in
commutative and non-commutative algebra, is via _Gröbner bases_. A couple of
extra concepts are needed before the definition can be stated.

  * •

Given two monomials $\partial^{p}$ and $\partial^{q}$, a _term order_ $\prec$
discerns which of them is the largest one. We refer to [86] for the technical
definition, and content ourselves with giving an example, namely the _graded
reverse lexicographic order_ 666 The speed with which Gröbner bases are
calculated can heavily depend on the chosen term order (of which there are
many). The graded reverse lexicographic order tends to be the fastest. .
Setting $|p|=p_{1}+\cdots+p_{N}$, then $\partial^{p}\prec\partial^{q}$ iff

$\displaystyle|p|<|q|\quad\text{or}\quad|p|=|q|\ \text{and the \emph{last}
non-zero entry of}\ q-p\in\mathbb{Z}^{N}\ \text{is \emph{negative}.}$

So the monomial with the largest degree wins, unless there is a tie, in which
case the one with the smallest variable index wins. In the case of $N=3$
variables and degree less than 3, we have

$\displaystyle
1\prec\partial_{3}\prec\partial_{2}\prec\partial_{1}\prec\partial_{3}^{2}\prec\partial_{2}\partial_{3}\prec\partial_{1}\partial_{3}\prec\partial_{2}^{2}\prec\partial_{1}\partial_{2}\prec\partial_{1}^{2}\,.$
(4.20)

  * •

Consider an operator
$L=\sum_{q\in\mathrm{supp}\left(L\right)}c_{q}(z)\,\partial^{q}\in\mathcal{R}$
in normally ordered form. Fixing a term order, the _initial monomial_ 777 The
similarity between this notation and $\mathrm{in}_{(u,v)}$ from (4.15) is not
coincidental. There is indeed a strong connection between comparing monomials
according to weights $(u,v)\in\mathbb{R}^{2N}$ and according to term orders
$\prec$ [85, Chapter 1].  $\mathrm{in}_{\prec}(L)$ is defined as the largest
monomial $\partial^{q}$ (stripped of its rational function
$c_{q}(z)\in\mathbb{C}(z)$) w.r.t. $\prec$.

  * •

The _initial ideal_ of $\mathcal{I}$ is the ideal generated by its initial
monomials:

$\displaystyle\mathrm{in}_{\prec}(\mathcal{I})=\mathbb{C}\cdot\big{\\{}\mathrm{in}_{\prec}(L)\
|\ L\in\mathcal{I}\big{\\}}\ \subset\
\mathbb{C}[\partial_{1},\ldots,\partial_{N}]\,.$ (4.21)

######  Definition 4.4.

A finite subset $G=\\{g_{1},\ldots,g_{M}\\}\in\mathcal{R}$ is called a
_Gröbner basis_ for the ideal $\mathcal{I}$ w.r.t. the term order $\prec$ iff

  1. 1.

$G$ generates $\mathcal{I}$, i.e. $\langle
g_{1},\ldots,g_{M}\rangle=\mathcal{I}$.

  2. 2.

$\mathrm{in}_{\prec}(\mathcal{I})=\langle\mathrm{in}_{\prec}(g_{1}),\ldots,\mathrm{in}_{\prec}(g_{M})\rangle.$

A Gröbner basis containing the minimal number of generators $g_{i}$ is called
_reduced_ , and is unique for a given term order. A monomial $\partial^{q}$ is
called a _standard monomial_ w.r.t. $G$ when it does _not_ belong to
$\mathrm{in}_{\prec}(\mathcal{I})$.

Absent any other methods, the standard monomials only become known at the end
of a Gröbner basis computation.

We can now answer the question posed in the very beginning of this section:
there are exactly $\mathrm{rank}\\!\left(\mathcal{I}\right)$ many standard
monomials, and they form an operator basis for the $\mathcal{D}$-module
$\mathcal{D}/\mathcal{I}$!

Gröbner basis computations are typically performed using Buchberger’s
algorithm888 Bruno Buchberger was in fact himself the inventor of Gröbner
bases, but he named them after his doctoral supervisor Wolfgang Gröbner as a
sign of respect.  [86]. The intensity of such calculations grows significantly
with the number of variables. For this reason, the main theme of the next
chapter will be to develop a more efficient algorithm that avoids Gröbner
bases.

####  Section 4.4Pfaffian systems

Let $\mathcal{I}$ be a holonomic ideal born from a system of PDEs, and suppose
$f=f(z)$ is a solution to this system. Assume a basis of standard monomials is
known for the $\mathcal{R}$-module $\mathcal{R}/\mathcal{I}$. It is here
denoted by

$\displaystyle\mathrm{Std}=\\{\partial^{q_{1}},\partial^{q_{2}},\ldots,\partial^{q_{R}}\\}\,,$
(4.22)

where $R:=\mathrm{rank}\\!\left(\mathcal{I}\right)$ and each
$q_{i}\in\mathbb{Z}^{N}_{\geq 0}$. Without loss of generality, we may set
$\mathrm{Std}_{1}=1$. A basis of _functions_ for the PDE solution space is
then given by the action of each standard monomial on $f$:

$\displaystyle e=\left[\begin{array}[]{c}\mathrm{Std}_{1}\bullet f=f\\\
\mathrm{Std}_{2}\bullet f\\\ \vdots\\\ \mathrm{Std}_{R}\bullet
f\end{array}\right]\,.$ (4.27)

We are interested in finding a set of first-order PDEs satisfied by $e$. To
this end, consider the Weyl algebra element $\partial_{i}\mathrm{Std}_{j}$.
Since the standard monomials form a basis for $\mathcal{R}/\mathcal{I}$, it
must be possible to write

$\displaystyle\partial_{i}\mathrm{Std}_{j}=\sum_{k=1}^{R}p_{ij}^{k}(z)\,\mathrm{Std}_{k}\quad\text{mod}\
\mathcal{I}\,,\quad p_{ij}^{k}(z)\in\mathbb{C}(z)\,.$ (4.28)

The notation $"\text{mod}\ \mathcal{I}"$ means that this equation holds modulo
the ideal $\mathcal{I}$, in the sense of (4.11). Most of the time we leave the
$"\text{mod}\ \mathcal{I}"$ as implicit.

Combining (4.27) with (4.28), we then have

$\displaystyle\partial_{i}\bullet
e=\left[\begin{array}[]{c}(\partial_{i}\mathrm{Std}_{1})\bullet f\\\
(\partial_{i}\mathrm{Std}_{2})\bullet f\\\ \vdots\\\
(\partial_{i}\mathrm{Std}_{R})\bullet f\end{array}\right]=P_{i}(z)\cdot
e\quad\text{mod}\ \mathcal{I}\,,\quad i=1,\ldots,N\,,$ (4.33)

where $P_{i}(z)$ is an $R\times R$ matrix whose $(j,k)$th entry equals the
rational function $p_{ij}^{k}(z)$. The set of first-order PDEs (4.33) is
dubbed a _Pfaffian system_. By manipulating the standard monomials, it can be
shown [86] that the _Pfaffian matrices_ $P_{i}$ satisfy the _integrability
condition_

$\displaystyle\partial_{i}\bullet P_{j}+P_{i}\cdot P_{j}=\partial_{j}\bullet
P_{i}+P_{j}\cdot P_{i}\,.$ (4.34)

The next chapter will be dedicated to deriving Pfaffian systems from the point
of view of $\mathcal{D}$-modules. Such a calculation requires two steps:

  1. 1.

Find a set of standard monomials $\mathrm{Std}$.

  2. 2.

Derive the decomposition (4.28).

Both steps are possible when a Gröbner basis is known (the second step is then
performed via an algorithm analogous to that of polynomial division). However,
as mentioned earlier, we shall seek to avoid Gröbner bases as their
computation often requires huge resources. As an aside, if a Pfaffian system
is known by other means, then one can immediately construct a Gröbner basis
too [82, Appendix B].

We end this section with the following observation: if $f$ were a single FI,
then (4.33) is exactly the kind of Pfaffian systems that we studied in Part I
of this thesis! The decomposition of $\partial_{i}\mathrm{Std}_{j}$ in (4.28)
is then analogous to an IBP reduction.

### Chapter 5Macaulay Matrices

Based on [82], this chapter presents a novel algorithm for the computation of
Pfaffian systems. It avoids the use of Gröbner bases in favor of linear
algebra methods, making it quite efficient.

Throughout this chapter we assume that

a set of standard monomials
$\mathrm{Std}=\\{\partial^{q_{i}}\\}_{i=1,\ldots,R}$ is given, (5.1)

where the cardinality $|\mathrm{Std}|=R$ equals the holonomic rank of a given
left ideal $\mathcal{I}$. While this is a strong assumption, it will actually
be fulfilled for the so-called GKZ system that we study in Chapter (6).

####  Section 5.1From Pfaffian to Macaulay matrix

Let $\mathcal{I}=\langle d_{1},\ldots,d_{D}\rangle$ be a holonomic ideal
generated by a collection of $D$ operators $d_{i}\in\mathcal{R}$. In this
section we suppose to know the Pfaffian matrices $P_{i}$ associated to
$\mathcal{I}$. It holds that

$\displaystyle\partial_{i}\mathrm{Std}_{j}=\sum_{k=1}^{R}(P_{i})_{jk}\,\mathrm{Std}_{k}\quad\text{in}\quad\mathcal{R}/\mathcal{I}\,.$
(5.2)

This expression can be lifted from the quotient ring $\mathcal{R}/\mathcal{I}$
to the whole ring $\mathcal{R}$ by writing