aainstitutetext: Physik-Institut, Universität Zürich, Winterthurerstrasse 190,
CH-8057 Zürich, Switzerlandbbinstitutetext: Dipartimento di Fisica e
Astronomia, Universitá di Bologna e INFN, Sezione di Bologna, via Irnerio 46,
I-40126 Bologna, Italy

# Reduction to master integrals via intersection numbers and polynomial
expansions

Gaia Fontana b Tiziano Peraro<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Intersection numbers are rational scalar products among functions that admit
suitable integral representations, such as Feynman integrals. Using these
scalar products, the decomposition of Feynman integrals into a basis of
linearly independent master integrals is reduced to a projection. We present a
new method for computing intersection numbers that only uses rational
operations and does not require any integral transformation or change of
basis. We achieve this by systematically employing the polynomial series
expansion, namely the expansion of functions in powers of a polynomial. We
also introduce a new prescription for choosing dual integrals, de facto
removing the explicit dependence on additional analytic regulators in the
computation of intersection numbers. We describe a proof-of-concept
implementation of the algorithm over finite fields and its application to the
decomposition of Feynman integrals at one and two loops.

††preprint: ZU-TH 19/23

##  1 Introduction

Scattering amplitudes at the loop level are linear combinations of Feynman
integrals. These integrals, in dimensional regularization, obey linear
relations – most notably integration by parts identities Tkachov:1981wb ;
Chetyrkin:1981qh , Lorentz invariance identities Gehrmann:1999as  and symmetry
relations. Using these we can reduce a large number of integrals into a much
smaller set of linearly independent ones called _master integrals_. This
reduction is an essential step of most modern multi-loop predictions.

The most successful and widely used integral decomposition method is the
Laporta algorithm Laporta:2000dsw , which consists in generating and solving a
large and sparse system of equations obeyed by loop integrals. This approach
produced extraordinary results and has been a key ingredient of the bulk of
multi-loop predictions in the last two decades, also thanks to efficient
public implementations vonManteuffel:2012np ; Smirnov:2008iw ;
Maierhofer:2017gsa . More recently, this approach has been combined with
finite fields and functional reconstruction techniques vonManteuffel:2014ixa ;
Peraro:2016wsq , where complex algebraic manipulations are replaced with
efficient numerical evaluations over machine-size integers and analytic
results are reconstructed out of several numerical evaluations. This method –
also implemented in both private and public Klappert:2019emp ; Smirnov:2019qkx
; Peraro:2019svx ; Klappert:2020nbg  programs – significantly pushed the state
of the art of multi-loop integral decomposition and, more generally, multi-
loop theoretical predictions.

Despite its tremendous success, the reduction to master integrals is still a
major bottleneck of higher-order theoretical predictions and it very often
determines whether a multi-loop calculation is doable or not. Moreover, this
approach doesn’t make manifest the vector-space structure obeyed by Feynman
integrals. Indeed, since all the Feynman integrals in a family (i.e. having a
common denominator structure) can be decomposed into a finite Smirnov:2010hn
basis of master integrals, they can be viewed as the elements of a finite-
dimensional vector space over the field of rational functions of the
invariants that describe a process. In traditional reduction approaches, this
property only indirectly comes out of a large number of linear identities and
it is not directly exploited. The search for decomposition methods that more
directly exploit this structure is definitely of interest from a purely
theoretical point of view, as it is going to improve our understanding of
Feynman integrals. From a more pragmatic and phenomenological viewpoint, this
improved understanding may inspire the development of more efficient
decomposition methods, opening up a number of important theoretical and
phenomenological studies that are currently beyond our technical capability.

An approach that has been gaining a lot of interest is the study of Feynman
integrals within the mathematical framework of _intersection theory_
Mizera:2017rqa ; Mastrolia:2018uzb ; Frellesvig:2019kgj ; Frellesvig:2019uqt ;
Mizera:2019gea ; Mizera:2019vvs ; Frellesvig:2020qot ; Weinzierl:2020xyy ;
Chestnov:2022alh ; Chestnov:2022xsy . The latter offers a method for defining
rational scalar products, called _intersection numbers_ , in vector spaces of
functions with a suitable integral representation – among which are Feynman
integrals. Using a scalar product, in turn, it is straightforward to project
any element of a vector space into its components with respect to a vector
basis. Hence, the problem of decomposing any Feynman integral into a basis of
master integrals is reduced to the problem of computing scalar products.

This approach, while promising, is currently not competitive with more
traditional decomposition methods. One of the main drawbacks of current
approaches to the calculation of intersection numbers is the appearance of
non-rational expressions in intermediate stages of the calculation. Indeed,
intersection numbers are typically computed as sums of residues of local
solutions of differential equations around the poles of some rational
functions. These poles can generally appear at non-rational locations. Because
intersection numbers are rational, all non-rational contributions will
eventually combine into a rational expression after summing over all residues.
The handling of non-rational expressions is however a major algebraic
bottleneck. Moreover, their appearance prevents us from using the available
finite-fields based techniques and programs in this context, since these apply
to algorithms that are a sequence of rational operations. Finally, since both
the inputs and the outputs of intersection numbers are rational, it would be
mathematically and theoretically more satisfactory to avoid non-rational
expressions in all steps of the calculation. A possible solution to this
problem was proposed in Weinzierl:2020xyy , where a rational algorithm for
computing intersection numbers was formulated in terms of a sequence of
changes of bases, integral transformations and the application of the global
residue theorem. The required changes of bases and integral transformations
are however highly non-trivial and can themselves become a bottleneck,
especially for more complex integrals.

In this paper, we present a new purely rational method for computing
intersection numbers which does not involve any integral transformation or
change of basis. The explicit appearance of irrational poles is sidestepped
via the systematic usage of the _polynomial series expansion_ , or
$p(z)$_-adic expansion_ for brevity, namely the series expansion of rational
functions in powers of a prime polynomial (it is conceptually similar to
$p$-adic numbers, which are instead series expansions of rational numbers in
powers of a prime integer). A $p(z)$-adic expansion is effectively equivalent
to a local expansion around all the (potentially irrational) zeroes of a
polynomial $p(z)$ at the same time. However, it can simply be computed in
terms of a polynomial division algorithm – which is obviously rational – and
does not require to know the location of the roots of $p(z)$. We thus locally
solve the aforementioned differential equations in terms of $p(z)$-adic
(rather than Laurent) expansions of their solution. For the final sum over
residues we then use a simple generalization of the global residue theorem
that applies to functions with a known $p(z)$-adic expansion.

We implement the new method over finite fields using the Mathematica package
FiniteFlow and test it on several one- and two-loop examples, finding
agreement with reductions obtained with the Laporta algorithm.

Another drawback of intersection theory, when applied to loop integrals, is
the need of introducing additional analytic regulators that significantly
increase the complexity of the calculation and obfuscate the typical block
triangular structure of integral reductions. A mathematical formalism for
dealing with this was proposed in Caron-Huot:2021xqj ; Caron-Huot:2021iev . In
this paper we propose an alternative approach based on a suitable choice of
master integrals for the dual space of Feynman integrals and on performing
operations on the coefficients of the leading terms in an expansion where the
additional regulators disappear. This “effectively” removes any analytic
dependence on the regulators in the calculation and provides a number of major
simplifications, besides recovering the aforementioned block triangular
structure.

The paper is structured as follows. In section 2 we give a pedagogical review
of intersection theory and its application to Feynman integrals, setting the
notation and the main concepts used afterwards. Section 3 contains the main
result of this paper, namely the description of a new rational algorithm for
computing intersection numbers. In section 4 we describe our new approach for
dealing with analytic regulators. An implementation of our new method over
finite fields is described in section 5. Finally, in section 6 we validate our
results with some examples and in section 7 we draw our conclusions and
discuss possible future developments. We give more details about the
$p(z)$-adic expansion and the global residue theorem in the appendixes.

##  2 Intersection numbers of Feynman integrals

In this section we set the notation and review some basic concepts of
intersection theory and its application to the reduction of Feynman integrals
into a basis of master integrals. This overview is meant to be accessible to
readers that are not familiar with intersection theory and we do not attempt
to be mathematically rigorous or thorough. For a more comprehensive treatment
of the subject and some of its applications, we refer readers to
Mizera:2017rqa ; Mastrolia:2018uzb ; Frellesvig:2019kgj ; Frellesvig:2019uqt ;
Mizera:2019gea ; Mizera:2019vvs ; Frellesvig:2020qot ; Caron-Huot:2021xqj ;
Caron-Huot:2021iev ; Weinzierl:2020xyy ; Weinzierl:2022eaz ; Chestnov:2022alh
; Chestnov:2022xsy ; Cacciatori:2021nli ; Cacciatori:2022mbi ; Chen:2020uyk ;
Chen:2022lzr ; amslaurea27132  and references therein.

We are interested in analyzing classes of integrals, such as Feynman
integrals, that can be viewed as elements of a vector space. In particular, we
consider $n$-fold integrals in the variables $\mathbf{z}=(z_{1},\ldots z_{n})$
of the form111A more rigorous definition considers the elements
$|\varphi_{R}\rangle$ and $\langle\varphi_{L}|$ as equivalence classes of
$n$-forms such that forms in the same class would yield integrands that differ
by total derivatives, thus implying eq.s (5). In this paper, with an abuse of
notation, we identify these with the integrals themselves, since we are only
interested in finding linear relations between integrals that can be written
as integration by parts identities, i.e. in the form of eq. (5).

$|\varphi_{R}\rangle=\int\differential z_{1}\cdots\differential
z_{n}\,\frac{1}{u(\mathbf{z})}\,\varphi_{R}(\mathbf{z})$ (1)

and _dual integrals_

$\langle\varphi_{L}|=\int\differential z_{1}\cdots\differential
z_{n}\,u(\mathbf{z})\,\varphi_{L}(\mathbf{z}).$ (2)

We also refer to these as _right_ and _left_ integrals respectively. In the
previous equations $\varphi_{R}$ and $\varphi_{L}$ are _rational functions_ in
$\mathbf{z}$ while $u(\mathbf{z})$ is a multivalued function that regulates
the singularities of $\varphi_{R}$ and $\varphi_{L}$. In particular, in our
applications, $u$ takes the form

$u(\mathbf{z})=\prod_{j}B_{j}(\mathbf{z})^{\gamma_{j}},$ (3)

where $B_{j}(\mathbf{z})$ are polynomials and $\gamma_{j}$ are _generic
exponents_ , i.e. exponents that are not identically equal to any integer. For
Feynman integrals, $\gamma_{j}$ are functions of symbolic regulators, such as
the dimensional regulator. The explicit form of $u$ (or $B_{j}$) and the
rational functions $\varphi_{R,L}$ depends on the specific problem one is
interested in – for our identifications see subsection 2.2. We study classes
of integrals with the same $u(\mathbf{z})$ and different rational functions
$\varphi_{L,R}$, namely Feynman integrals within the same family.

We also assume that these (regulated) integrands vanish at the boundaries of
integration. Hence, we can write _integration by parts identities (IBPs)_ as

$\displaystyle\sum_{j=1}^{n}\int\differential z_{1}\cdots\differential
z_{n}\,\partial_{z_{j}}\Big{(}\frac{1}{u}\,\xi^{(R)}_{j}\Big{)}={}0$
$\displaystyle\sum_{j=1}^{n}\int\differential z_{1}\cdots\differential
z_{n}\,\partial_{z_{j}}\Big{(}u\,\xi^{(L)}_{j}\Big{)}={}0$ (4)

or equivalently,

$\displaystyle\sum_{j=1}^{n}\left|\Big{(}\partial_{z_{j}}-(\partial_{z_{j}}u)/u\Big{)}\xi^{(R)}_{j}\right\rangle=0$
$\displaystyle\sum_{j=1}^{n}\left\langle\Big{(}\partial_{z_{j}}+(\partial_{z_{j}}u)/u\Big{)}\xi^{(L)}_{j}\right|=0$
(5)

for a list of $n$ rational functions $\xi^{(R)}_{j}$ and $\xi^{(L)}_{j}$. By
generating and solving a large number of these identities, one can decompose
any integral $|\varphi_{R}\rangle$ into a _basis_ of linearly independent
integrals $\left\\{|e_{j}^{(R)}\rangle\right\\}_{j=1}^{\nu}$, commonly called
_master integrals_ , where $\nu$ is the (finite) dimension of the vector space
of integrals. One can similarly reduce dual integrals $\langle\varphi_{L}|$
into a dual basis $\left\\{\langle e_{j}^{(L)}|\right\\}_{j=1}^{\nu}$. In this
work we are however interested in a more direct method for achieving such
decompositions.

Our main goal is the computation of _rational scalar products_

$\langle\varphi_{L}|\varphi_{R}\rangle$ (6)

between left and right integrals. These are referred to as _intersection
numbers_. Intersection numbers are obviously linear in both $\varphi_{L}$ and
$\varphi_{R}$ and consistent with the identities in eq.s (5).

By calculating these scalar products and combining them, we can project left
and right integrals into master integrals in the same way a vector belonging
to a vector space can be decomposed into a basis via projections. Therefore we
express the decomposition of any integral $|\varphi_{R}\rangle$ and dual
integral $\langle\varphi_{L}|$ in their bases as

$\displaystyle|\varphi_{R}\rangle={}$
$\displaystyle\sum_{i=1}^{\nu}c^{(R)}_{i}|e_{i}^{(R)}\rangle$
$\displaystyle\langle\varphi_{L}|={}$
$\displaystyle\sum_{i=1}^{\nu}c^{(L)}_{i}\langle e_{i}^{(L)}|$ (7)

with

$\displaystyle c^{(R)}_{i}={}$
$\displaystyle\sum_{j=1}^{\nu}\quantity(\mathbf{C}^{-1})_{ij}\langle
e_{j}^{(L)}|\varphi_{R}\rangle$ $\displaystyle c^{(L)}_{i}={}$
$\displaystyle\sum_{j=1}^{\nu}\langle\varphi_{L}|e_{j}^{(R)}\rangle\quantity(\mathbf{C}^{-1})_{ji},$
(8)

where we introduced the metric $\mathbf{C}_{ij}$

$\mathbf{C}_{ij}=\langle e_{i}^{(L)}|e_{j}^{(R)}\rangle,$ (9)

that contains the intersection numbers between all the master integrals.

In the following, we will thus focus on methods for computing intersection
numbers and their application to the decomposition of Feynman integrals to
master integrals.

###  2.1 Computing intersection numbers

We now review a known method for computing multivariate intersection numbers
Frellesvig:2019kgj ; Frellesvig:2019uqt , which is recursive in the
integration variables. We start from the univariate case, which is the base
case of the recursion, and then describe the recursive step, which computes
multivariate intersection numbers of $n$-fold integrals in terms of
intersection numbers of $(n-1)$-fold integrals.

While intersection numbers are rational, the method reviewed in this section
generally involves non-rational intermediate expressions, namely poles of
rational functions. In section 3 we will propose a new method that avoids non-
rational expressions in all steps of the calculation.

#### Univariate case

For the univariate case, we consider intersection numbers between $1$-fold
integrals of the form

$\displaystyle|\varphi_{R}\rangle={}$ $\displaystyle\int\differential
z\,\frac{1}{u(z)}\,\varphi_{R}(z)$ $\displaystyle\langle\varphi_{L}|={}$
$\displaystyle\int\differential z\,u(z)\,\varphi_{L}(z).$ (10)

Intersection numbers are calculated as a sum of residues

$\langle\varphi_{L}|\varphi_{R}\rangle=\sum_{p\in\mathcal{P}_{\omega}}\text{Res}_{z=p}\quantity(\psi\,\varphi_{R}),$
(11)

where the sum is extended over all $p\in\mathcal{P}_{\omega}$

$\mathcal{P}_{\omega}=\left\\{\,z\,|\,z\,\text{is a pole
of}\,\omega\,\right\\}\bigcup\left\\{\infty\right\\}$ (12)

and $\psi$ is the local solution around the poles $p$ of the differential
equation

$(\partial_{z}+\omega)\psi=\varphi_{L},$ (13)

where

$\omega\equiv\frac{\partial_{z}u}{u}.$ (14)

The local solution can be computed as a Laurent expansion around $z=p$

$\psi=\sum_{i=\text{min}}^{\text{max}}c_{i}\quantity(z-p)^{i}+O((z-p)^{\text{max}+1}),$
(15)

which inserted in eq. (13) yields a linear system of equations for the unknown
coefficients $c_{i}$. The Laurent expansion can be truncated at a finite order
that is sufficient for the purpose of taking the residue in eq. (11).

#### Multivariate case

We now consider $n$-fold integrals of the form in eq.s (1) and (2). We rewrite
these as

$\displaystyle|\varphi_{R}\rangle={}$ $\displaystyle\int\differential
z_{n}\,|\varphi_{R}\rangle_{n-1}$ $\displaystyle\langle\varphi_{L}|={}$
$\displaystyle\int\differential z_{n}\,\langle\varphi_{L}|_{n-1},$ (16)

where $|\varphi_{R}\rangle_{n-1}$ and $\langle\varphi_{L}|_{n-1}$ are
$(n-1)$-fold integrals which have a parametric dependence on $z_{n}$. Since
this is a recursive method, we assume to be able to compute intersection
numbers of $(n-1)$-fold integrals, which we write as

$\langle\varphi_{L}|\varphi_{R}\rangle_{n-1}$ (17)

and are _rational functions_ of $z_{n}$. These, in turn, can be used to
project any $(n-1)$-fold right and left integral into a basis of master
integrals $\left\\{|e_{j}^{(R)}\rangle_{n-1}\right\\}_{j=1}^{\nu_{(n-1)}}$ and
$\left\\{\langle e_{j}^{(L)}|_{n-1}\right\\}_{j=1}^{\nu_{(n-1)}}$
respectively,

$\displaystyle|\varphi_{R}\rangle_{n-1}={}$
$\displaystyle\sum_{j=1}^{\nu_{(n-1)}}\varphi_{R,j}\,|e_{j}^{(R)}\rangle_{n-1}$
$\displaystyle\langle\varphi_{L}|_{n-1}={}$
$\displaystyle\sum_{j=1}^{\nu_{(n-1)}}\,\varphi_{L,j}\langle
e_{j}^{(L)}|_{n-1},$ (18)

where the coefficients $\varphi_{R,j}$ and $\varphi_{L,j}$ are also rational
functions of $z_{n}$.

Intersection numbers of $n$-fold integrals can be computed as

$\langle\varphi_{L}|\varphi_{R}\rangle=\sum_{p\in\mathcal{P}_{\mathbf{\Omega}}}\text{Res}_{z_{n}=p}\quantity(\sum_{j}\psi_{j}\,\langle
e_{j}^{(L)}|\varphi_{R}\rangle_{n-1}),$ (19)

where $\psi_{j}$ are solutions of the system of differential equations

$\partial_{z_{n}}\psi_{j}+\sum_{k=1}^{\nu_{(n-1)}}\psi_{k}\,\mathbf{\Omega}_{kj}=\varphi_{L,j}\quad(j=1,\dots,\nu_{(n-1)}),$
(20)

with $\varphi_{L,j}$ defined as in eq. (18), while the matrix
$\mathbf{\Omega}_{jk}$ is the differential equation matrix of $(n-1)$-fold
left master integrals with respect to $z_{n}$,

$\partial_{z_{n}}\langle
e_{j}^{(L)}|_{n-1}=\sum_{k=1}^{\nu_{(n-1)}}\mathbf{\Omega}_{jk}\,\langle
e_{k}^{(L)}|_{n-1}.$ (21)

In practice, one can take the derivative of the $(n-1)$-fold master integrals
with respect to $z_{n}$ under the integral sign and decompose the resulting
expression into the same master integrals to obtain the matrix
$\mathbf{\Omega}_{jk}$. The sum in eq. (19) runs over the poles of the matrix
elements $\mathbf{\Omega}_{jk}$

$\mathcal{P}_{\mathbf{\Omega}}=\\{\,z\,|\,z\,\text{is a pole
of}\,\mathbf{\Omega}\,\\}\bigcup\\{\infty\\}.$ (22)

On a side note, we observe that eq. (19) can be rewritten as

$\langle\varphi_{L}|\varphi_{R}\rangle=\sum_{p\in\mathcal{P}_{\mathbf{\Omega}}}\text{Res}_{z_{n}=p}\langle\psi|\varphi_{R}\rangle_{n-1},$
(23)

where $\langle\psi|_{n-1}$ solves

$\partial_{z_{n}}\langle\psi|_{n-1}=\langle\varphi_{L}|_{n-1},$ (24)

from which the more practically useful eq. (20) is obtained by decomposing
both $\langle\varphi_{L}|_{n-1}$ and $\langle\psi|_{n-1}$ into $(n-1)$-fold
left master integrals, defining $\psi_{j}$ as

$\langle\psi|_{n-1}=\sum_{j=1}^{\nu_{(n-1)}}\,\psi_{j}\langle
e_{j}^{(L)}|_{n-1}.$ (25)

###  2.2 Application to Feynman integrals

Multiloop Feynman integrals in dimensional regularization are usually written,
as a first step, in momentum representation. After a tensor decomposition, a
generic $L$-loop integral of a given _integral family_ with $E+1$ external
momenta $p_{1},\ldots,p_{E+1}$ (of which only $E$ are independent because of
momentum conservation) takes the form

$I[\alpha_{1},\dots,\alpha_{n}]=\int\prod_{j=1}^{L}\differential^{d}k_{j}\frac{1}{z_{1}^{\alpha_{1}}\dots
z_{n}^{\alpha_{n}}},$ (26)

where $z_{1},\ldots,z_{n}$ are a complete set of denominators of loop
propagators and auxiliary denominators, such that any scalar product of the
form $k_{i}\cdot k_{j}$ or $k_{i}\cdot p_{j}$ is a linear combination of the
$z_{i}$. It follows that $n=L(L+1)/2+LE$. The exponents $\alpha_{j}$ are
_integers_ , which can be positive, negative or zero (negative exponents thus
correspond to numerators).

Because of the vanishing of total derivatives in dimensional regularization,
we can write IBP identities Tkachov:1981wb ; Chetyrkin:1981qh

$\int\prod_{j=1}^{L}\differential^{d}k_{j}\partialderivative{k_{i}^{\mu}}\quantity(\frac{v^{\mu}}{z_{1}^{\alpha_{1}}\dots
z_{n}^{\alpha_{n}}})=0,$ (27)

where $v^{\mu}=k_{i}^{\mu}$ or $v^{\mu}=p_{i}^{\mu}$. We can exploit these
identities, possibly combined with Lorentz invariance identities
Gehrmann:1999as  and symmetry relations, to perform the decomposition of a
given set of Feynman integrals to master integrals. The most widely used
method, i.e. the Laporta algorithm Laporta:2000dsw , consists in writing a
long list of linear identities of the form of eq. (27) (and similar for
Lorentz invariance identities and symmetries) for different choices of
$v^{\mu}$ and $\alpha_{j}$ and combining them into a large and sparse system
of linear equations satisfied by the integrals. The solution of such system
yields the reduction to master integrals. As mentioned, solving this system is
one of the main bottlenecks of many higher-order predictions. Moreover, the
vector-space structure obeyed by Feynman integrals is not manifest in this
approach.

In order to work in the framework of intersection theory, we need to rewrite
eq. (26) into a representation that mimics the one in eq. (1). Throughout this
work, we use the Baikov representation Baikov:1996iu , which consists in
changing the integration variables from loop momenta $k_{j}^{\mu}$ to
propagators $z_{i}$. Other representations are also usable in this context,
see e.g. Lee:2013hzt ; Caron-Huot:2021xqj ; Caron-Huot:2021iev . In the Baikov
representation, Feynman integrals read

$I[\alpha_{1},\dots,\alpha_{n}]=K\int\differential z_{1}\dots\differential
z_{n}\,B(z_{1},\ldots,z_{n})^{\gamma}\,\frac{1}{z_{1}^{\alpha_{1}}\dots
z_{n}^{\alpha_{n}}}$ (28)

with $\gamma=\frac{d-L-E-1}{2}$. The Baikov polynomial $B$ is the Gram
determinant of the loop momenta and the independent external momenta,
expressed as a function of the $z_{i}$. $K$ is a constant that does not depend
on the integration variables nor the exponents. Feynman integrals in the
Baikov representation also obey IBP identities, that are the vanishing of
total derivatives with respect to the $z_{j}$.

Comparing eq. (28) to the integrals defined in section 2, we identify the
multivalued function $u$ and the rational functions $\varphi_{R,L}$ as

$u(\mathbf{z})=B^{-\gamma}\,\prod_{j=1}^{n}z_{j}^{\rho_{j}}\,\qquad\varphi_{R,L}(\mathbf{z})=\frac{1}{z_{1}^{\alpha_{1}}\dots
z_{n}^{\alpha_{n}}},$ (29)

respectively. As in Frellesvig:2019uqt , we inserted factors
$z_{j}^{\rho_{j}}$ in $u(\mathbf{z})$ because $\varphi_{R,L}$ can have
singularities in $z_{j}=0$ that are generally not regulated by the Baikov
polynomial alone, while having them regulated is a requirement for applying
the framework of intersection theory, as we stated in section 2. Indeed, one
can easily check that the differential equations we need to locally solve to
compute intersection numbers have no solution around unregulated
singularities. It is understood that we are interested in the limit
$\rho_{j}\to 0$ of a decomposition to master integrals. We will refer to these
additional variables $\rho_{j}$ as _analytic regulators_. Analytic regulators
can be dropped for auxiliary denominators $z_{i}$ which only appear with a
negative exponent in $\varphi_{R,L}$.

The need of analytic regulators is an additional complication of the
application of intersection theory to Feynman integrals, since it adds
additional parameters into the problem and obfuscates the typical block
triangular structure of integral reduction identities. We will show in section
4, however, that if we are interested in the decomposition of Feynman
integrals, we can effectively remove the explicit dependence on $\rho_{j}$
from the calculation of intersection numbers, by making a suitable choice of
basis for dual integrals and systematically work in the $\rho_{j}\to 0$ limit.
An alternative mathematical framework for dual Feynman integrals which also
removes the need of regulators was proposed in Caron-Huot:2021xqj ; Caron-
Huot:2021iev .

In this paper we identify Feynman integrals in the Baikov representation as
right integrals $|\varphi_{R}\rangle$ and their duals as left integrals
$\langle\varphi_{L}|$. We therefore apply the framework of intersection theory
for computing intersection numbers between Feynman integrals and their duals.
This approach allows to directly decompose an integral into the basis of
chosen master integrals via the projections in (7), without resorting to the
construction of an IBP system. For a given integral, this approach implies the
necessity to calculate intersection numbers of $n$-fold integrals, with $n$
being the number of independent denominators.

##  3 Intersection numbers via polynomial series expansions

In this section, we present the main result of this paper, namely a new
rational method for computing intersection numbers via the systematic use of
_polynomial series expansions_ , or $p(z)$-adic expansions for brevity.

As mentioned, one of the main drawbacks of the method described in section 2
is the appearance of non-rational poles in the calculation of intersection
numbers. While intersection numbers are rational, irrational intermediate
expressions represent a bottleneck in complicated algebraic manipulations and
make the usage of efficient techniques based on finite fields impractical.

Because we always deal with rational integrands, all their (rational and non-
rational) poles are roots of rational polynomial factors of their
denominators. In the following, we will thus consider each of these polynomial
factors, which are _prime_ 222We can, in principle, relax the requirement of
having prime polynomial factors to the one of having a list of pairwise co-
prime polynomial factors. over $\mathbb{Q}$, rather than individual poles, as
the building blocks of our calculation.

###  3.1 $p(z)$-adic series expansions

We elaborate on the expansion of any rational function in terms of a
polynomial $p(z)$, a necessary step to calculate intersection numbers with our
newly-proposed rational algorithm.

Consider a univariate polynomial $p(z)$ of degree $\deg p$ in the variable $z$
that is prime over the field of rational numbers. We define the $p(z)$-adic
series expansion of a rational function $f(z)$ as

$f(z)=\sum_{i=\textrm{min}}^{\textrm{max}}c_{i}(z)\,p^{i}(z)+\mathcal{O}\left(p(z)^{\textrm{max}+1}\right),$
(30)

where the coefficients $c_{i}(z)$ are polynomials in $z$ over the rational
field, whose degree is lower than the one of $p(z)$

$c_{i}(z)=\sum_{j=0}^{\deg p-1}c_{ij}z^{j},\quad c_{ij}\in\mathbb{Q}.$ (31)

The notation $\mathcal{O}((p(z)^{k})$ stands for terms proportional to the
$k$-th power of the polynomial $p(z)$. If $z$ is such that $p(z)=\delta$,
where $\delta$ is some small quantity, then we obviously have
$\mathcal{O}(p(z)^{k})=\mathcal{O}(\delta^{k})$. Hence, the $p(z)$-adic
expansion allows us to consider the limit $p(z)\to 0$, without resorting to
non-rational operations nor knowing the explicit location of the (irrational)
roots of $p(z)$. In particular we have

$f(z)|_{p(z)=\delta}=\mathcal{O}(\delta^{\textrm{min}}),$ (32)

where $i=\textrm{min}$ is the starting order of the expansion. Similarly to
Laurent expansions, the latter can be deduced from polynomial factors of the
form $p(z)^{k}$ in $f(z)$, but unlike Laurent expansions around irrational
points, here they can be detected using only operations over the rational
field.

When calculating intersection numbers, we are often required to expand around
the irrational poles of $\omega$ or $\mathbf{\Omega}$. These irrational poles
are the roots of some polynomials $p(z)$ irreducible over $\mathbb{Q}$.
Therefore a $p(z)$-adic expansion is effectively equivalent to considering all
expansions of a function around all the roots of the polynomial at once,
ultimately avoiding their explicit appearance.

The expansion in (30) can be obtained by repeated polynomial divisions modulo
$p(z)$ with remainder, as illustrated in appendix A. As a simpler alternative,
we observe that we can implement the same expansion via a shortcut that avoids
repeated polynomial divisions. We introduce an auxiliary parameter $\delta$
and define, for a generic function $f(z)$, its remainder333As explained more
in detail in appendix A, this also involves computing the multiplicative
inverse of the denominator of $f(z)$ modulo $p(z)-\delta$. Some computer
algebra systems have builtin functions for polynomial reminders that are also
applicable to rational functions, such as the PolynomialRemainder function in
Mathematica. modulo $p(z)-\delta$

$\lfloor f(z)\rfloor_{p(z)-\delta}\equiv
f(z)\,\text{mod}\,\big{(}p(z)-\delta\big{)},$ (33)

that is a polynomial of degree $\deg p-1$ in $z$ and has a rational dependence
on the parameter $\delta$. We highlight that this step amounts to the
substitution $p(z)=\delta$ and is exact, i.e. we obtain back the original
function by substituting $\delta\to p(z)$ in the right-hand-side. After doing
this, we series expand for small $\delta$ obtaining

$\lfloor f(z)\rfloor_{p(z)-\delta}\evaluated{}_{\delta\rightarrow
0}=\sum_{i=\textrm{min}}^{\text{max}}\sum_{j=0}^{\deg
p-1}c_{ij}z^{j}\delta^{i}+O(\delta^{\text{max}+1}),$ (34)

where the coefficients $c_{ij}$ are identified with those appearing in the
$p(z)$-adic expansion as in eq.s (30) and (31). In this way we obtain the
$p(z)$-adic expansion of a rational function up to the order
$O(p(z)^{\text{max}+1})$.

Polynomial series expansions can be combined with the _univariate global
residue theorem_ (see appendix B) to compute the sum of the residues of a
rational function around the poles of a rational polynomial, with purely
rational operations. A polynomial $p(z)$ can always be factorized over the
complex field as

$p(z)=l_{c}\,\prod_{k=1}^{\textrm{deg }p}(z-y_{k}),$ (35)

where $l_{c}$ is the _leading coefficient_ , i.e. the coefficient of
$z^{\textrm{deg }p}$ in $p(z)$, and $y_{k}$ are the roots of $p(z)$. For a
generic function $f(z)$ we define

$\text{Res}_{p(z)}\quantity(f(z))\equiv\sum_{k=1}^{\textrm{deg
}p}\text{Res}_{z=y_{k}}\quantity(f(z)).$ (36)

i.e. the sum of the residues of $f(z)$ at the roots of $p(z)$. If $f(z)$
admits the $p(z)$-adic expansion in eq. (30), then the sum of residues in the
previous equation can be computed without resorting to irrational operations
and without explicitly knowing the location of the roots $y_{k}$, using the
following generalization of the global residue theorem

$\text{Res}_{p(z)}\quantity(f(z))=\frac{c_{-1,\textrm{deg }p-1}}{l_{c}},$ (37)

where the coefficients $c_{ij}$ are defined as in eq. (31).

We now apply $p(z)$-adic expansions to the computation of univariate and
multivariate intersection numbers.

###  3.2 Univariate intersection numbers

For the computation of intersection numbers of univariate integrals in the
framework of $p(z)$-adic expansion, we rewrite eq. (11) as

$\langle\varphi_{L}|\varphi_{R}\rangle=\sum_{p(z)\in\mathcal{P}_{\omega}[z]}\langle\varphi_{L}|\varphi_{R}\rangle_{p(z)},$
(38)

where the sum runs over all factors $p(z)$ of the denominators of $\omega$

$\mathcal{P}_{\omega}[z]=\left\\{\textrm{factors of the denominator of
$\omega$}\right\\}\bigcup\left\\{\infty\right\\}.$ (39)

The contribution $\langle\varphi_{L}|\varphi_{R}\rangle_{\infty}$ is computed
exactly as the contribution at $p=\infty$ to eq. (11). Each other addend of
the form $\langle\varphi_{L}|\varphi_{R}\rangle_{p(z)}$ is instead the sum of
all contributions to (11) coming from the roots of $p(z)$.

For each factor $p(z)$, we make an ansatz for the local solution $\psi$ of the
differential equation (13) close to all the roots of $p(z)$. This takes the
form of a $p(z)$-adic expansion

$\psi=\sum_{i=\text{min}}^{\text{max}}\sum_{j=0}^{\deg
p-1}c_{ij}z^{j}p(z)^{i}+O(p(z)^{\text{max}+1}).$ (40)

We stress that, even when $p(z)$ has irrational roots, these do not explicitly
appear in eq. (40) which, once plugged into eq. (13) and $p(z)$-adic expanded
again, yields a linear system of equations with rational coefficients that is
then solved for the unknowns $c_{ij}$.

The solution for $\psi$ is then multiplied by $\varphi_{R}$ and the product is
again $p(z)$-adic expanded with respect to $p(z)$. We obtain

$\psi\,\varphi_{R}=\sum_{i}\sum_{j=0}^{\deg
p-1}\tilde{c}_{ij}z^{j}p(z)^{i}+O(p(z)^{0})$ (41)

and apply the univariate global residue theorem as in eq. (37)

$\langle\varphi_{L}|\varphi_{R}\rangle_{p(z)}=\frac{\tilde{c}_{-1,\deg
p-1}}{l_{c}},$ (42)

where $l_{c}$ is the leading coefficient of $p(z)$. As we stated, this is the
sum of the contributions to the intersection number between
$\langle\varphi_{L}|\varphi_{R}\rangle$ given by all the roots of $p(z)$.

###  3.3 Multivariate intersection numbers

The generalization to the multivariate case is straightforward. In this
subsection we let $z\equiv z_{n}$. We rewrite eq. (19) as

$\langle\varphi_{L}|\varphi_{R}\rangle=\sum_{p(z)\in\mathcal{P}_{\mathbf{\Omega}}[z]}\langle\varphi_{L}|\varphi_{R}\rangle_{p(z)},$
(43)

where the sum is over the factors $p(z)$ of the denominators of
$\mathbf{\Omega}$.

$\mathcal{P}_{\mathbf{\Omega}}[z]=\left\\{\textrm{factors of the denominator
of $\mathbf{\Omega}_{ij}$}\right\\}\bigcup\left\\{\infty\right\\}.$ (44)

As before, the contribution $\langle\varphi_{L}|\varphi_{R}\rangle_{\infty}$
is computed exactly as the contribution at $p=\infty$ to eq. (19). Each other
addend of the form $\langle\varphi_{L}|\varphi_{R}\rangle_{p(z)}$ is instead
the sum of all contributions to (19) coming from the roots of $p(z)$.

We make an ansatz for the local solutions $\psi_{i}$ to the differential
equation (20) around all the roots of $p(z)$ that, as in the univariate case,
take the form of $p(z)$-adic expansions

$\psi_{i}=\sum_{j=\text{min}}^{\text{max}}\sum_{k=0}^{\deg
p-1}c_{ijk}z^{k}p(z)^{j}+O(p(z)^{\text{max}+1})\quad(i=1,\dots,\nu_{n-1}).$
(45)

Eq. (20) thus yields a linear system that we solve for the unknowns $c_{ijk}$.

We then take the scalar product of $\psi_{j}$ and $\langle
e_{j}^{(L)}|\varphi_{R}\rangle_{n-1}$ and $p(z)$-adic expand with respect to
$p(z)$. We obtain

$\sum_{i=1}^{\nu_{(n-1)}}\psi_{i}\,\langle
e_{i}^{(L)}|\varphi_{R}\rangle_{n-1}=\sum_{j}\sum_{k=0}^{\deg
p-1}\tilde{c}_{jk}z^{k}p(z)^{j}+O(p(z)^{0}).$ (46)

The contribution of the roots of $p(z)$ to the intersection number
$\langle\varphi_{L}|\varphi_{R}\rangle$ is again given by

$\langle\varphi_{L}|\varphi_{R}\rangle_{p(z)}=\frac{\tilde{c}_{-1,\deg
p-1}}{l_{c}}.$ (47)

##  4 Dual integrals and analytic regulators

As explained in subsection 2.2, one of the drawbacks of the application of
intersection theory to Feynman integrals is the need of introducing _analytic
regulators_ $\rho_{j}$ for each integration variable $z_{i}$ that may appear
as a denominator for $\varphi_{R}$ or $\varphi_{L}$. However, since we are
interested in computing intersection numbers for the purpose of reducing
Feynman integrals $\varphi_{R}$ to master integrals, we can exploit the
freedom of choice for the dual basis to simplify our problem. In this section
we discuss a simple strategy that “effectively” removes the explicit
dependence on analytic regulators from the calculation, without substantially
changing the algorithm for computing intersection numbers. An alternative
mathematical formalism for choosing dual integrals and dealing with this issue
was presented in Caron-Huot:2021xqj ; Caron-Huot:2021iev . A comparison
between the two approaches, as well as an in-depth study of all implications
and possible limitations of our new approach, is outside the purposes of this
work and is left to future investigations.

###  4.1 Choice of dual integrals

Our approach simply consists in choosing dual integrals of the form

$\varphi_{L}(\mathbf{z})=\rho_{1}^{\Theta(\alpha_{1}-\frac{1}{2})}\cdots\rho_{n}^{\Theta(\alpha_{n}-\frac{1}{2})}\frac{1}{z_{1}^{\alpha_{1}}\dots
z_{n}^{\alpha_{n}}}$ (48)

and _systematically_ work in the $\rho_{j}\to 0$ limit. In other words, we
multiply each factor $z_{j}^{-\alpha_{j}}$ by the analytic regulator
$\rho_{j}$ if (and only if) $\alpha_{j}>0$, i.e. when $z_{j}$ appears in the
denominator of $\varphi_{L}$. Then, each step in the algorithm is worked out
by keeping only the leading contributions in the $\rho_{j}\to 0$ limit.

One can easily see that the dependence on $\rho_{j}$ of intersection numbers
only appears at the $j$-th level of the recursion. On the other hand,
intersection numbers of $j$-fold integrals are systematically computed in the
$\rho_{j}\to 0$ limit, which is _finite_ when we use left integrands of the
form (48). Hence, we only have to consider one regulator at a time, namely
$\rho_{j}$ (and only $\rho_{j}$) at the $j$-th step of the recursion.

Consider a generic integration variable $z=z_{j}$ and the corresponding
regulator $\rho=\rho_{j}$. When solving eq. (13) for $\psi$, for each order of
the Laurent or $p(z)$-adic series expansion in $z$ of the equation, we only
consider the leading term in a $\rho\to 0$ expansion (similar statements hold
for $\psi_{j}$ solving eq. (20)). Since we systematically work on the leading
coefficients of a $\rho\to 0$ expansion and obtain _finite_ intersection
numbers after each step of the recursion, the calculation is effectively as
simple as one without any dependence on $\rho$ – or even simpler in most
cases.

For a more detailed discussion, we distinguish the two cases of exponents
$\alpha\leq 0$ and exponents $\alpha>0$ inside $\varphi_{L}\sim z^{-\alpha}$.
If $\alpha\leq 0$, then $\varphi_{L}$ has no singularity in $z=0$ hence we do
not expect to need an analytic regulator at all and one may set $\rho=0$
before solving for $\psi$. If instead $\alpha>0$, since poles in $z=0$ are
generally _not_ regulated by $\omega$, eq. (13) has no solution unless we add
an analytic regulator. The solution for $\psi$, in turn, develops a pole in
$\rho=0$, which cancels against the prefactor $\rho$ we add according to the
prescription in (48), yielding a finite solution for $\psi$ in the $\rho\to 0$
limit. We observe that, therefore, for $\alpha>0$ the solution for $\psi$ is
only non-zero (in the $\rho\to 0$ limit) around poles at $z=0$ and possibly at
$z=\infty$, hence we can skip the calculation at all other poles or
denominator factors.

This strategy has a number of advantages.

  * •

The calculation is _effectively_ independent of $\rho_{j}$ since we are
directly working on the leading coefficients of expansion around $\rho_{j}\to
0$. The full dependence on $\rho_{j}$ is never considered and, when we
implement the algorithm over finite fields, we do not sample over or
reconstruct the $\rho_{j}$ dependence.

  * •

Intermediate analytic expressions (if reconstructed) are dramatically simpler.

  * •

The matrix of intersection numbers $\langle\varphi_{L}|\varphi_{R}\rangle$ for
a list of $\varphi_{L}$ and $\varphi_{R}$, as well as the matrices
$\mathbf{C}_{ij}$ and $\mathbf{\Omega}_{ij}$ have a block triangular
structure, similar to what we expect from integral decompositions in
traditional approaches, where blocks correspond to so-called _sectors_. In
particular, $\langle\varphi_{L}|\varphi_{R}\rangle\neq 0$ only if the set of
indexes $k$ such that $\alpha_{k}>0$ for the integrand $\varphi_{L}$ is either
the same set or a subset of the one for $\varphi_{R}$ – in the language of
Feynman integrals, we say that $\varphi_{L}$ belongs to a sector that is a
subsector of $\varphi_{R}$. The same structure appears in the matrices
$\mathbf{C}_{ij}$ and $\mathbf{\Omega}_{ij}$ (examples are given in section
6). We note that this structure is transposed to the one that is observed in
the reduction or differential equations for (right) Feynman integrals
$|\varphi_{R}\rangle$.

  * •

A large set of intersection numbers and contributions of poles to intersection
numbers vanishes, for the reasons illustrated above.

  * •

At intermediate steps of the recursive algorithm, sometimes there are _fewer
master integrals_ with respect to the case where the full $\rho$ dependence is
considered.

###  4.2 Subtleties of analytic regulators

A few subtleties need to be discussed. The system of equations for the
coefficients of the expansion of $\psi$ in $z$ can sometimes be
underdetermined, regardless of the sign of the exponent $\alpha$ of $z$. In
this case, we add more constraints coming from higher-order terms in the
$\rho\to 0$ expansion of the equations. This can also happen when $\alpha\leq
0$ – in this case a regulator is temporarily added – as well as for
contributions to intersection numbers from denominator factors $p(z)\neq
0,\infty$ where we can regulate the system with the prescription

$\omega\to\omega+\rho_{p}\frac{\partial_{z}\,p(z)}{p(z)}\qquad\textrm{or}\qquad\mathbf{\Omega}_{ij}\to\mathbf{\Omega}_{ij}+\rho_{p}\frac{\partial_{z}\,p(z)}{p(z)}\,\delta_{ij}.$
(49)

Even in such cases, we still solve the equation for the coefficients of $\psi$
in the $\rho_{p}\to 0$ limit, keeping higher-order terms in such limit only as
necessary. In most cases the additional constraints simply amount to setting
to zero the previously undetermined coefficients, hence the next-to-leading
equations in the $\rho\to 0$ (or $\rho_{p}\to 0$) expansion can almost always
be dropped after a numerical check.

Another subtlety is that, even if a variable $z$ is not a factor of the Baikov
polynomial, it is still possible that $\mathbf{\Omega}$ may have $z$ as a
denominator factor at the corresponding step of the recursion. This, in turn,
might regulate the singularity at $z\to 0$ for at least some of the integrals
$\varphi_{L}$, without the need of the corresponding analytic regulator
$\rho$. In such cases $\psi$ does not develop the singularity $1/\rho$ that
cancels the $\rho$ prefactor in eq. (48) and therefore we obtain $\psi=0$ for
all poles – which for the purpose of computing intersection numbers is
equivalent to setting $\langle\varphi_{L}|=0$. While in principle we may
remove the $\rho$ prefactor from eq. (48) when this happens, in all our tests
this was not needed. Indeed, even in such cases we have always found enough
integrals of the form in (48) to obtain a complete basis of dual integrals in
the $\rho\to 0$ limit.

As a final remark, our approach is only meant to deal with singularities at
$z=0$ or underdetermined systems for expansions around arbitrary factors
$p(z)$. In cases where a system with no solution appears for a $p(z)\neq
z,\infty$, both this approach and the traditional one that considers the full
dependence of intersection numbers on the analytic regulators in (29) are
problematic and deserve further investigation. In principle the prescription
in (49) may fix the singular point, but the solution for $\psi$ generally
develops a singularity in $\rho_{p}=0$ which only cancels after the
intersection numbers are combined into a decomposition as in eq. (8). It is
possible that a generalization of the strategy proposed here – namely
multiplying the offending integrands by $\rho_{p}$ and systematically work in
the $\rho_{p}\to 0$ limit after applying eq. (49) – may also simplify the
calculation in these cases.

##  5 Implementation over finite-fields

One of the main motivations for developing a rational algorithm for
intersection numbers, besides avoiding algebraic difficulties in handling non-
rational roots, is the possibility of implementing it over finite fields.
Finite fields are numerical fields with a finite number of elements. In the
context of scattering amplitudes, the most common choice is the field of
integers modulo a prime number, namely
$\mathbb{Z}_{p}=\left\\{0,\ldots,p-1\right\\}$ with $p$ a prime. In
$\mathbb{Z}_{p}$ we can perform efficient numerical evaluations using machine-
size integers (for primes $p<2^{64}$) which are also exact, i.e. not affected
by numerical inaccuracies. We can implement over $\mathbb{Z}_{p}$ any
algorithm that is rational, i.e. a sequence of rational arithmetic operations.
Out of these numerical evaluations, we can thus use rational reconstruction
techniques to infer full analytic results as rational functions of the free
parameters of the problem (see e.g. vonManteuffel:2014ixa ; Peraro:2016wsq ;
Klappert:2019emp ; Peraro:2019svx  for more details).

In this section we describe a proof-of-concept implementation of our new
algorithm over finite fields that uses the Mathematica package FiniteFlow
Peraro:2019svx . This package offers the possibility of building new
algorithms by combining core algorithms in computational graphs (dataflow
graphs) from a high-level interface and reconstruct analytic results out of
numerical evaluations.

We give an overview of how our implementation of multivariate intersection
numbers works. At first, for simplicity, we will assume to know in advance a
basis of master integrals and dual master integrals for all levels of the
recursive algorithm. We will then address this point more thoroughly in
subsection 5.4.

###  5.1 Setting up the recursive calculation

Before doing any explicit calculation, we make a list of intersection numbers
we need to compute at each step of the recursive algorithm. Our starting point
is a list of $n$-variate intersection numbers we wish to compute. This will
typically include all the intersection numbers needed for the decomposition of
a selection of integrals $|\varphi_{R}\rangle$ via eq.s (7) and (8), namely
$\langle e^{(L)}_{j}|\varphi_{R}\rangle$ and the ones appearing in the metric
$\mathbf{C}_{ij}=\langle e^{(L)}_{i}|e^{(R)}_{j}\rangle$, for a basis
$\left\\{|e_{j}^{(R)}\rangle\right\\}_{j=1}^{\nu}$ and a dual basis
$\left\\{\langle e_{j}^{(L)}|\right\\}_{j=1}^{\nu}$ – but we consider a
generic list of $n$-variate intersection numbers
$\langle\varphi_{L}|\varphi_{R}\rangle$ for the sake of generality. From
these, by reviewing the algorithm in subsection 2.1, we can easily make a list
of $(n-1)$-variate intersection numbers which are required as input for the
recursive algorithm. More precisely, if
$\left\\{|e_{j}^{(R)}\rangle_{n-1}\right\\}_{j=1}^{\nu_{(n-1)}}$ and
$\left\\{\langle e_{j}^{(L)}|_{n-1}\right\\}_{j=1}^{\nu_{(n-1)}}$ are the
chosen bases for right and left $(n-1)$-fold integrals respectively, the
calculation of the $n$-variate intersection numbers in the list requires the
following $(n-1)$-variate intersection numbers

  * •

the matrix elements of the $(n-1)$-variate metric $\langle
e^{(L)}_{i}|e^{(R)}_{j}\rangle_{n-1}$

  * •

$\langle\varphi_{L}|e^{(R)}_{j}\rangle_{n-1}$ for the decomposition of
$\langle\varphi_{L}|$ needed in the differential equation (20)

  * •

$\langle d^{(L)}_{i}|e^{(R)}_{j}\rangle_{n-1}$ with $\langle
d^{(L)}_{i}|_{n-1}\equiv\partial_{z_{n}}\langle e^{(L)}_{i}|_{n-1}$ needed for
computing the matrix elements $\mathbf{\Omega}_{ij}$ defined in eq. (21)

  * •

$\langle e^{(L)}_{j}|\varphi_{R}\rangle_{n-1}$ appearing in eq. (19).

This makes up a new list of required $(n-1)$-variate intersection numbers. We
proceed recursively, by applying the same analysis to this new list, obtaining
all the required $(n-2)$-variate intersection numbers, until we reach the
univariate case.

Using this method, we easily generate a list of required intersection numbers
for each step of the recursive algorithm, before any explicit calculation is
done. We can thus proceed with the implementation of the algorithm described
in section 3, starting from univariate intersection numbers and moving up
until the $n$-variate case.

###  5.2 Laurent and $p(z)$-adic expansions

In this subsection, we describe how we implement Laurent and $p(z)$-adic
expansions over finite fields, which are extensively used in all steps of the
recursion. While FiniteFlow already has algorithms for reconstructing Laurent
expansions, we find it is more convenient to adopt a different method that
uses the sparse linear solvers implemented in the package.

The starting point is a rational function $f(z)$ for which we want to compute
the series expansion, which we write as

$f(z)=\frac{n(z)}{d(z)}=\frac{\sum_{j}n_{j}\,z^{j}}{\sum_{j}d_{j}z^{j}}$ (50)

and for which we have reconstructed the full dependence on $z$, while the
coefficients $n_{j}$ and $d_{j}$ are generally just known numerically, i.e.
they can be evaluated via a rational algorithm.

A Laurent expansion for $f(z)$ around $z=0$ can be easily obtained by making
an ansatz for it

$f(z)=\sum_{j=\min}^{\max}c_{j}\,z^{j}+\mathcal{O}(z^{\max+1}),$ (51)

where $\min$ is easily determined from the lowest-degree non-vanishing terms
in the denominator and numerator of $f$, while $\max$ is the order at which
the expansion is needed. We thus plug eq. 51 into the relation

$f(z)\,d(z)-n(z)=0,$ (52)

truncated at the suitable order determined by $\max$, and impose that the
coefficient of all terms proportional to $z^{k}$ in the previous equation must
vanish. This gives a linear system of equations for the coefficients $c_{i}$
in eq. 51. It is easy to see that this system is triangular, hence equivalent
to a list of substitutions. Expansions around $z=\infty$ can be obtained
similarly, by considering instead the function

$g(\delta)=-\frac{1}{\delta^{2}}\,f\left(\frac{1}{\delta}\right)$ (53)

where the prefactor $-1/\delta^{2}$ is the Jacobian of the transformation.
Hence we expand $g(\delta)$ for $\delta\to 0$, using the same method
illustrated above. Note that the coefficients in the numerator and the
denominator of $g$ coincide with the ones for $f$, except that they multiply
different monomials, hence there’s no need of shifting variables in this case.

The calculation of $p(z)$-adic expansions is instead done in two steps. Let
$p(z)$ be the polynomial with respect to which we want to expand our functions

$p(z)=\sum_{k=0}^{\deg p}a_{k}\,z^{k}.$ (54)

In the first step we compute

$z^{k}\mod\big{(}p(z)-\delta\big{)}$ (55)

for all needed values of the exponent $k$. This is done by writing a linear
system of equations of the form

$\sum_{k=0}^{\deg p}a_{k}\,z^{k+\alpha}\,\delta^{\beta}-\delta^{\beta+1}=0,$
(56)

for integers $\alpha$ and $\beta$, where the unknowns are monomials
$z^{k}\delta^{j}$. We solve these linear relations for the monomials,
expressing monomials with higher powers in $z$ as linear combinations of
monomials with lower powers in $z$ (and higher powers in $\delta$). This
yields the $p(z)$-adic expansion of every monomial we need, that takes the
form

$z^{k}=\sum_{j}\sum_{l=0}^{\deg p-1}b_{kjl}z^{l}\,\delta^{j}.$ (57)

We stress that this expansion is exact, i.e. it will produce an exact identity
after substituting $\delta\to p(z)$, although in practice it can also be
truncated to the needed order in $\delta$. We also notice that the system is
triangular, hence equivalent to a sequence of linear substitutions. In the
second step, we proceed similarly to the expansion around $z=0$, by making an
ansatz for the $p(z)$-adic expansion of $f$

$f(z)=\sum_{j=\min}^{\max}\sum_{k=0}^{\deg
p-1}c_{jk}z^{k}\,\delta^{j}+\mathcal{O}(\delta^{\max+1}),$ (58)

where we have set $p(z)=\delta$ for convenience. The minimum order can be
inferred by the order of the expansion of the denominator $d(z)$, which in
turn can be easily derived by substituting eq. (57) into its expression. We
thus substitute this expansion into eq. (52), perform the substitutions in
(57) and impose that all the coefficients of terms $z^{k}\,\delta^{j}$ of the
(appropriately truncated) expansion of (57) vanish. This yields a linear
system of equations for the coefficients $c_{jk}$. Notice that, in the special
case where $p(z)$ is linear, eq. (57) is equivalent to a shift of variables
and the second step becomes identical to the expansion around zero discussed
above.

We observe that we typically need to expand several rational functions having
the same denominator. As a shortcut, we expand a generic numerator where the
coefficients $n_{j}$ are symbolic. In practice, we treat $n_{j}$ as additional
unknowns in our linear system, so that the solution for the coefficients of
the expansion will be a linear combination of the coefficients $n_{j}$.
Moreover, for expansions around $z=0$ and $z=\infty$, we remove factors of $z$
in the denominator by allowing the numerator $n(z)$ to contain terms with
negative powers of $z$. This allows us to combine the expansion of a larger
set of functions together.

###  5.3 Polynomial factors and multivariate intersection numbers

We now describe the most important step of the implementation, namely the
recursive step in the calculation of multivariate intersection numbers. Both
the univariate and the multivariate algorithm illustrated in section 3 consist
in the computation of Laurent or $p(z)$-adic expansions which are then plugged
into eq. (13) to set up linear systems for the coefficients of $\psi$, which
are ultimately used to take the sum of residues contributing to the
intersection numbers. All of these are clearly combinations of rational
operations. The non-trivial part, that we discuss in the following, is how to
generate the input for each step of the recursion without reconstructing
(potentially large) intermediate results of the previous steps.

For the univariate case, the input is just the analytic expression for
$u(\mathbf{z})$, or equivalently the polynomials $B_{j}(\mathbf{z})$ and the
generic exponents $\gamma_{j}$.

For the multivariate case, we start from a numerical implementation of the
required list of all $(n-1)$-variate intersection numbers. These, via
straightforward algebraic operations, are combined to obtain the following
quantities

  * •

the matrix elements $\mathbf{\Omega}_{ij}$ defined in eq. (21)

  * •

the coefficients $\varphi_{L,j}$ of the decomposition into left master
integrals of $\langle\varphi_{L}|_{n-1}$ needed in the differential equation
(20)

  * •

the intersection numbers $\langle e^{(L)}_{j}|\varphi_{R}\rangle_{n-1}$
appearing in eq. (19).

In the following we call the union of these three lists of quantities
$\mathcal{X}_{n}$, which we are thus able to evaluate numerically (over finite
fields) as a function of $z_{n}$ and the other parameters it depends on,
including other integration variables $z_{k}$ with $k>n$. It is easy to see
that $\mathcal{X}_{n}$ is the input required for the computation of
$n$-variate intersection numbers via the recursive algorithm we illustrated in
this paper. More precisely, we need to know the dependence on $z_{n}$ of each
function $f\in\mathcal{X}_{n}$. Moreover, we would like to have the list of
factors $p(z_{n})$ appearing in eq. (44). In general, we want to avoid to
compute the polynomial factors $p(z_{n})$ of the denominators of
$\mathcal{X}_{n}$ every time we evaluate intersection numbers. Knowing the
full analytic dependence of the factors $p(z_{n})$ (not just on $z_{n}$ but on
every other parameter) avoids the need to perform polynomial factorization at
each numerical evaluation and also allows to significantly simplify the
reconstruction in $z_{n}$. Moreover, as we will now see, the reconstruction of
the full analytic dependence of $p(z_{n})$ is generally very simple and can be
performed in a relatively small number of numerical evaluations.

In order to set up the recursion in our numerical implementation, we follow
these three steps

  1. 1.

We first reconstruct every $f\in\mathcal{X}_{n}$ in the variable $z_{n}$
modulo a prime number $p$, while every other parameter is set to a random
numerical value. We thus factorize the denominators of the reconstructed
functions, modulo $p$, to obtain a semi-numerical list of factors $p(z_{n})$.

  2. 2.

We proceed to reconstruct the full analytic dependence of the denominator
factors $p(z_{n})$. We first identify a “simple” subset $\mathcal{S}_{n}$ of
$\mathcal{X}_{n}$ such that i) the union of the $z_{n}$-dependent denominator
factors in $\mathcal{S}_{n}$ coincides with the one in $\mathcal{X}_{n}$ and
ii) the degrees of numerators and denominators of the functions in
$\mathcal{S}_{n}$ are as low as possible. The subset $\mathcal{S}_{n}$ is
easily identified from the $z_{n}$ reconstruction of point 1. We thus proceed
to reconstruct the $z_{n}$-dependent denominator factors in $\mathcal{S}_{n}$
as follows. For each $f\in\mathcal{S}_{n}$, we make an ansatz of its
dependence on $z_{n}$ of the form

$f(z_{n})=\frac{n(z_{n})}{d(z_{n})}=\frac{\sum_{j}n_{j}\,z_{n}^{j}}{\sum_{j}d_{j}\,z_{n}^{j}}$
(59)

which, again, is easily built from the result of point $1$. We stress that the
coefficients $n_{j}$ and $d_{j}$ will have a rational dependence on other
parameters, that is implicit in the previous equation. We also fix the
arbitrary normalization of numerator and denominator by setting the lowest
degree coefficient of the denominator to be equal to one, i.e.
$d_{\textrm{min}(j)}=1$. This choice also has the effect of moving any factor
of $f$ that does _not_ depend on $z_{n}$ into the numerator coefficients
$n_{j}$, hence keeping the coefficients $d_{j}$ very simple. We can thus write
the previous equation as

$\sum_{j}n_{j}\,z_{n}^{j}-\sum_{j}d_{j}\,z_{n}^{j}\,f(z_{n})=0,$ (60)

which evaluated for several values of of $z_{n}$ yields a linear system that
we solve for $n_{j}$ and $d_{j}$. From numerical solutions of this linear
systems, obtained with different values of the parameters they depend on, we
reconstruct the full analytic dependence of the denominator coefficients
$d_{j}$ only. Given their relative simplicity, this usually requires a small
number of evaluations and the reconstruction is thus quite efficient. By
factorizing the full analytic form of these denominators in a computer algebra
system, we obtain the full list of factors defined in eq. (44).

  3. 3.

By combining the results obtained in steps $1$ and $2$ we can easily map each
denominator factor (or each suitable product of factors) appearing in the
semi-numerical reconstruction of point $1$ into its fully analytic counterpart
reconstructed in step $2$. This allows us to make an ansatz for each
$f\in\mathcal{X}_{n}$ of the form

$f(z_{n})=\sum_{\alpha_{1},\ldots,\alpha_{m},k}c_{\alpha_{1},\ldots,\alpha_{m},k}\,\frac{z_{n}^{k}}{\prod_{j=1}^{m}p_{j}(z_{n})^{\alpha_{j}}},$
(61)

where the unknown coefficients $c_{\alpha_{1},\ldots,\alpha_{m},k}$ are
independent of $z_{n}$. By evaluating $\mathcal{X}_{n}$ for several values of
$z_{n}$ we obtain a linear system of equations that we solve for the unknown
coefficients. This list of numerically-evaluated coefficients, together with
the functional form in $z_{n}$ given by the ansatz, are used as input for the
calculation of $n$-variate intersection numbers, which proceeds following
subsection 3.3.

###  5.4 Master integrals

In each step of the recursive algorithm, we need to select a basis and a dual
basis of master integrals, which for simplicity was assumed to be known in
advance in the discussion above. In ref. Lee:2013hzt  it is shown that the
number of master integrals can be found as the number of solutions of a system
of rational equations. As observed in Frellesvig:2019uqt  the same number can
also be computed as the dimension of a polynomial ideal, which in turn can be
found via rational operations. A guess for the list of master integrals can be
made from the list of independent monomials in such ideal, see e.g. ref.s
Weinzierl:2020xyy ; amslaurea27132  for more details.

As an alternative, we can also find bases of master integrals with a more
pragmatic approach that is based on the computation of intersection numbers.
We may start from a list of integrals which form overcomplete bases

$\mathcal{\tilde{E}}_{R}=\left\\{|\tilde{e}^{(R)}_{j}\rangle\right\\},\qquad\mathcal{\tilde{E}}_{L}=\left\\{\langle\tilde{e}^{(L)}_{j}|\right\\}$
(62)

and compute an “enlarged metric” $\mathbf{\tilde{C}}$

$\mathbf{\tilde{C}}_{ij}=\langle\tilde{e}^{(L)}_{i}|\tilde{e}^{(R)}_{j}\rangle$
(63)

which is obviously not invertible since the basis is overcomplete. We then
column reduce the matrix (or better, a numerical evaluation of it over finite
fields) to find a minimal list of independent columns, which corresponds to a
basis of independent dual integrals

$\mathcal{E}_{L}=\left\\{\langle
e^{(L)}_{j}|\right\\}_{j}=\left\\{\langle\tilde{e}^{(L)}_{j}|\right\\}_{j\in\textrm{indep.\
columns}}\subset\mathcal{\tilde{E}}_{L}.$ (64)

The independent columns can thus be row-reduced and the linearly independent
rows correspond to a basis of master integrals

$\mathcal{E}_{R}=\left\\{|e^{(R)}_{j}\rangle\right\\}_{j}=\left\\{|\tilde{e}^{(R)}_{j}\rangle\right\\}_{j\in\textrm{indep.\
rows}}\subset\mathcal{\tilde{E}}_{R}.$ (65)

In each step of the recursive algorithm, this strategy is used to find a basis
$\mathcal{E}_{R}$ and a dual basis $\mathcal{E}_{L}$ of master integrals.

##  6 Examples

We present some one- and two-loop examples of reduction to master integrals
that use the new rational algorithm presented in this paper. All the examples
we showcase also use our new prescription of section 4 for dealing with
analytic regulators. The two methods are however independent of each other and
our new rational algorithm has also been extensively tested on several
examples retaining the full dependence on analytic regulators (see also
amslaurea27132 ).

In the following subsections we report some of the decompositions we obtained
and the bases of master integrals we used, including the ones for intermediate
layers of integration. We also include, in some cases, the reconstructed
metric of intersection numbers between master integrals. We stress that one
never needs to analytically reconstruct the metric when performing a reduction
and we report it in some examples for the sole purpose of showing their
simplicity and block triangular structure, a feature that is absent when
keeping the full dependence on analytic regulators.

In the following, we use the notation in eq. (26) to identify Feynman
integrals $|\varphi_{R}\rangle$. As in the previous sections, the integration
variables are $z_{1},\ldots,z_{n}$ and in the recursive algorithm we choose
the ordering of integration where $z_{1}$ is the innermost integration
variable and $z_{n}$ the outermost one. Even if we do not explicitly write it,
to avoid cluttering the notation, we understand that all dual integrals
$\langle\varphi_{L}|$ are multiplied by analytic regulators according to the
prescription in section 4. In other words, the replacement

$\frac{1}{z_{1}^{\alpha_{1}}\dots
z_{n}^{\alpha_{n}}}\to\rho_{1}^{\Theta(\alpha_{1}-\frac{1}{2})}\cdots\rho_{n}^{\Theta(\alpha_{n}-\frac{1}{2})}\frac{1}{z_{1}^{\alpha_{1}}\dots
z_{n}^{\alpha_{n}}}$ (66)

is understood for all the integrands of _left integrals_ (i.e.  _dual
integrals_). For all families we produced reductions for a wide selection of
integrals with various combinations of positive and negative powers of
denominators. All of these have been successfully checked against the
decomposition obtained using the traditional Laporta algorithm. A small
selection among the simplest ones is quoted for illustration purposes.

###  6.1 One-loop families

At one-loop, we consider the integral families of the massless box (Fig. 1)
and the one-loop box with two off-shell external legs – also known as two-mass
hard box (Fig. 2). We compute the reduction of several integrals belonging to
these two families.

#### Massless box

Figure 1: One-loop massless box.

We consider the family of one-loop massless box with four external legs
$p_{1},\ldots,p_{4}$ satisfying

$p_{i}^{2}{}=0\quad(i=1,\dots,4),\\\ $ (67)

and momentum conservation

$p_{1}+p_{2}+p_{3}+p_{4}=0.$ (68)

The integrals in this family depend on two Mandelstam invariants

$\displaystyle s=$ $\displaystyle(p_{1}+p_{2})^{2},$ $\displaystyle t=$
$\displaystyle(p_{1}+p_{3})^{2}.$ (69)

The loop propagators $z_{i}$ are

$\displaystyle z_{1}$ $\displaystyle=k_{1}^{2},$ $\displaystyle z_{2}$
$\displaystyle=(k_{1}-p_{1})^{2},$ $\displaystyle z_{3}$
$\displaystyle=(k_{1}-p_{1}-p_{2})^{2},$ $\displaystyle z_{4}$
$\displaystyle=(k_{1}-p_{1}-p_{2}-p_{3})^{2}$ (70)

and yield the Baikov polynomial

$\displaystyle B=\frac{1}{16}\Big{(}$ $\displaystyle
s^{4}+2s^{3}t-2s^{3}z_{1}+2s^{3}z_{2}-2s^{3}z_{3}+2s^{3}z_{4}+s^{2}t^{2}-4s^{2}tz_{1}+2s^{2}tz_{2}$
$\displaystyle-4s^{2}tz_{3}+2s^{2}tz_{4}+s^{2}z_{1}^{2}+s^{2}z_{2}^{2}+s^{2}z_{3}^{2}+s^{2}z_{4}^{2}-2s^{2}z_{1}z_{2}+2s^{2}z_{1}z_{3}$
$\displaystyle-2s^{2}z_{2}z_{3}-2s^{2}z_{1}z_{4}+2s^{2}z_{2}z_{4}-2s^{2}z_{3}z_{4}-2st^{2}z_{1}-2st^{2}z_{3}+2stz_{1}^{2}$
$\displaystyle+2stz_{3}^{2}-2stz_{1}z_{2}-2stz_{2}z_{3}-2stz_{1}z_{4}+4stz_{2}z_{4}-2stz_{3}z_{4}+t^{2}z_{1}^{2}$
$\displaystyle+t^{2}z_{3}^{2}-2t^{2}z_{1}z_{3}\Big{)},$ (71)

We start with the one-fold integral in $z_{1}$. The input for univariate
intersection numbers is the Baikov polynomial and its exponent
$\gamma=1/2(d-5)$. The bases for both left and right integrals are given by

$\left\\{e^{(L)}_{j}\right\\}=\left\\{e^{(R)}_{j}\right\\}=\left\\{\frac{1}{z_{1}},1\right\\},$
(72)

where, as stated above, for left integrals the replacement in eq. (66) is
understood. The metric obtained by calculating the univariate intersection
numbers between the first-layer master integrals is

$\mathbf{C}=\left(\begin{array}[]{cc}1&0\\\
-\frac{s^{2}+s(t+z_{2}-z_{3}+z_{4})+tz_{3}}{(d-6)(s+t)}&-\frac{4(d-5)st(sz_{3}+tz_{3}+(z_{2}-z_{3})(z_{3}-z_{4}))}{(d-6)(d-4)(s+t)^{2}}\\\
\end{array}\right).$ (73)

Using one-fold intersection numbers, we are able to reduce the one-fold
integrals in $z_{1}$ to master integrals and compute the differential equation
matrix w.r.t. $z_{2}$, defined as in eq. (21), which reads

$\mathbf{\Omega}=\left(\begin{array}[]{cc}-\frac{(d-5)s\left(s^{2}+s(t+z_{2}-z_{3}+z_{4})-t(z_{3}-2z_{4})\right)}{p_{2}(\mathbf{z})}&0\\\
\frac{z_{4}(s(s+t-z_{2})+tz_{3})+z_{3}(s(s+t+z_{2})-z_{3}(s+t))+sz_{4}^{2}}{2(s+t)p_{1}(\mathbf{z})}&-\frac{(d-6)(z_{3}-z_{4})}{2p_{1}(\mathbf{z})}\\\
\end{array}\right)$ (74)

where the two $z_{2}$-dependent polynomial factors in the denominators are

$\displaystyle p_{1}(\mathbf{z})=$ $\displaystyle
z_{3}(s+t-z_{3}+z_{4})+z_{2}(z_{3}-z_{4}),$ $\displaystyle p_{2}(\mathbf{z})=$
$\displaystyle s^{2}z_{2}^{2}+2sz_{2}(s^{2}+st-
sz_{3}+sz_{4}-tz_{3}+2tz_{4})+(s^{2}+st-sz_{3}+sz_{4}-tz_{3})^{2}.$ (75)

We notice that both the metric and $\mathbf{\Omega}_{ij}$ present a triangular
structure. This is a general feature of our approach for dealing with analytic
regulators. We also observe that $p_{2}$ is quadratic in $z_{2}$, as well as
in the other variables $z_{i}$ it depends on. Hence, even in this simple
example, $\mathbf{\Omega}_{ij}$ has non-rational poles.

After calculating the relevant univariate intersection numbers, we move to the
second integration variable, $z_{2}$, and we compute the needed bivariate
intersection numbers. The bases for both left and right integrals for the
second layer are given by

$\left\\{e^{(L)}_{j}\right\\}=\left\\{e^{(R)}_{j}\right\\}=\left\\{\frac{1}{z_{1}z_{2}},\frac{1}{z_{1}},\frac{1}{z_{2}}\right\\}.$
(76)

The metric of bivariate intersection numbers reads

$\mathbf{C}=\left(\begin{array}[]{ccc}1&0&0\\\
\frac{s^{2}+s(t-z_{3}+z_{4})-t(z_{3}-2z_{4})}{(d-6)s}&-\frac{4(d-5)tz_{4}(s+t)(s-z_{3}+z_{4})}{(d-6)(d-4)s^{2}}&0\\\
-\frac{s^{2}+s(t-z_{3}+z_{4})+tz_{3}}{(d-6)(s+t)}&0&-\frac{4(d-5)stz_{3}(s+t-z_{3}+z_{4})}{(d-6)(d-4)(s+t)^{2}}\\\
\end{array}.\right)$ (77)

For the third integration variable $z_{3}$, the bases for both left and right
integrals are

$\left\\{e^{(L)}_{j}\right\\}=\left\\{e^{(R)}_{j}\right\\}=\left\\{\frac{1}{z_{1}z_{2}z_{3}},\frac{1}{z_{1}z_{3}},\frac{1}{z_{2}}\right\\}$
(78)

and the metric is

$\mathbf{C}=\left(\begin{array}[]{ccc}1&0&0\\\
\frac{s^{2}+s(t+z_{4})+2tz_{4}}{(d-6)s}&-\frac{4(d-5)tz_{4}(s+t)(s+z_{4})}{(d-6)(d-4)s^{2}}&0\\\
\frac{s(s+t+z_{4})^{2}}{(d-7)(d-6)(s+t)}&0&\frac{st(s+t+z_{4})^{4}}{4(d-7)(d-3)(s+t)^{2}}\\\
\end{array}.\right)$ (79)

We finally move on to the fourth and last integration variable, $z_{4}$. The
left and right bases for this integral family are

$\left\\{e^{(L)}_{j}\right\\}=\left\\{e^{(R)}_{j}\right\\}=\left\\{\frac{1}{z_{1}z_{2}z_{3}z_{4}},\frac{1}{z_{1}z_{3}},\frac{1}{z_{2}z_{4}}\right\\},$
(80)

corresponding to two bubbles and one box. This last step produces the
following metric

$\mathbf{C}=\left(\begin{array}[]{ccc}1&0&0\\\
\frac{s(s+t)}{(d-7)(d-6)}&-\frac{s^{2}t(s+t)}{4(d-7)(d-3)}&0\\\
\frac{s(s+t)}{(d-7)(d-6)}&0&\frac{st(s+t)^{2}}{4(d-7)(d-3)}\\\
\end{array}\right).$ (81)

An example of a simple reduction identity is

$\displaystyle I[1,2,1,2]=$
$\displaystyle-\frac{4(d-8)(d-5)(d-3)}{(d-6)s^{2}(s+t)^{2}}I[1,0,1,0]$
$\displaystyle-\frac{8(d-5)(d-3)}{(d-6)s(s+t)^{3}}I[0,1,0,1]$
$\displaystyle+\frac{(d-5)((d-6)s+2t)}{s(s+t)^{2}}I[1,1,1,1].$ (82)

For testing purposes, we also repeated this calculation (and several others)
without using the method of analytic regulators described in section 4 and
retaining the full dependence on the regulators, which we set to be equal to
each other, i.e. $\rho_{j}=\rho$ for all $j=1,\ldots,4$. In this case the
metric $\mathbf{C}_{ij}$ and the matrix $\mathbf{\Omega}_{ij}$ are dense and
their entries significantly more complicated. Moreover, at the third layer of
integration we have four master integrals instead of just three. These
features are _not_ particular to these example and display some of the general
advantages of the strategy illustrated in section 4.

#### Two-mass hard box

Figure 2: One-loop two-mass hard box.

We consider the family of the one-loop massless box with adjacent two massive
external legs $p_{3}$ and $p_{4}$, with $p_{3}^{2}\neq p_{4}^{2}$. This is
known in the literature as _two-mass hard box_. The external momenta
$p_{1},\dots,p_{4}$ satisfy

$\displaystyle p_{1}^{2}{}$ $\displaystyle=0,\,p_{2}^{2}{}=0,$ $\displaystyle
p_{3}^{2}{}$ $\displaystyle\neq 0,\,p_{4}^{2}{}\neq 0,$ $\displaystyle
p_{3}^{2}{}$ $\displaystyle\neq p_{4}^{2},$ (83)

and momentum conservation

$p_{1}+p_{2}+p_{3}+p_{4}=0.$ (84)

The integrals belonging to this family depend on the external masses
$p_{3}^{2}$ and $p_{4}^{2}$ and the following additional independent
invariants

$\displaystyle s=$ $\displaystyle(p_{1}+p_{2})^{2},$ (85) $\displaystyle t=$
$\displaystyle(p_{1}+p_{3})^{2}.$ (86)

The denominators $z_{i}$ are defined as in eq. (70) (but note that now
$p_{3}^{2}\neq 0$ and $p_{4}^{2}\neq 0$). The Baikov polynomial is

$\displaystyle B=\frac{1}{16}$
$\displaystyle\Bigl{(}p_{3}^{4}(s+z_{2}-z_{3})^{2}+2p_{3}^{2}\bigl{(}p_{4}^{2}(s^{2}-s(z_{1}+z_{3})+(z_{1}-z_{2})(z_{2}-z_{3}))-s^{3}$
$\displaystyle-s^{2}(t-z_{1}+2z_{2}-2z_{3}+z_{4})+st(z_{1}-z_{2}+2z_{3})+s(z_{2}-z_{3})(z_{1}-z_{2}+z_{3}-z_{4})$
$\displaystyle+t(z_{1}-z_{3})(z_{3}-z_{2})\bigr{)}+p_{4}^{4}(s-z_{1}+z_{2})^{2}-2p_{4}^{2}\bigl{(}s^{3}+s^{2}(t-2z_{1}+2z_{2}-z_{3}+z_{4})$
$\displaystyle+st(-2z_{1}+z_{2}-z_{3})+s(z_{1}-z_{2})(z_{1}-z_{2}+z_{3}-z_{4})+t(z_{1}-z_{2})(z_{1}-z_{3})\bigr{)}$
$\displaystyle+s^{4}+2s^{3}t-2s^{3}z_{1}+2s^{3}z_{2}-2s^{3}z_{3}+2s^{3}z_{4}+s^{2}t^{2}-4s^{2}tz_{1}+2s^{2}tz_{2}-4s^{2}tz_{3}$
$\displaystyle+2s^{2}tz_{4}+s^{2}z_{1}^{2}+s^{2}z_{2}^{2}+s^{2}z_{3}^{2}+s^{2}z_{4}^{2}-2s^{2}z_{1}z_{2}+2s^{2}z_{1}z_{3}-2s^{2}z_{2}z_{3}-2s^{2}z_{1}z_{4}$
$\displaystyle+2s^{2}z_{2}z_{4}-2s^{2}z_{3}z_{4}-2st^{2}z_{1}-2st^{2}z_{3}+2stz_{1}^{2}+2stz_{3}^{2}-2stz_{1}z_{2}-2stz_{2}z_{3}$
$\displaystyle-2stz_{1}z_{4}+4stz_{2}z_{4}-2stz_{3}z_{4}+t^{2}z_{1}^{2}+t^{2}z_{3}^{2}-2t^{2}z_{1}z_{3}\Bigr{)},$
(87)

There are four integration variables $z_{1},\ldots,z_{4}$. The bases of
$k$-fold master integrals for each layer $k$ of integration are given by

$\displaystyle 1-\text{fold}$
$\displaystyle\left\\{e^{(L)}_{j}\right\\}=\left\\{e^{(R)}_{j}\right\\}=\left\\{\frac{1}{z_{1}},1\right\\},$
$\displaystyle 2-\text{fold}$
$\displaystyle\left\\{e^{(L)}_{j}\right\\}=\left\\{e^{(R)}_{j}\right\\}=\left\\{\frac{1}{z_{1}z_{2}},\frac{1}{z_{2}},\frac{1}{z_{1}},1\right\\},$
$\displaystyle 3-\text{fold}$
$\displaystyle\left\\{e^{(L)}_{j}\right\\}=\left\\{e^{(R)}_{j}\right\\}=\left\\{\frac{1}{z_{1}z_{2}z_{3}},\frac{1}{z_{1}z_{3}},\frac{1}{z_{3}},\frac{1}{z_{2}},\frac{1}{z_{1}},1\right\\},$
$\displaystyle 4-\text{fold}$
$\displaystyle\left\\{e^{(L)}_{j}\right\\}=\left\\{e^{(R)}_{j}\right\\}=\left\\{\frac{1}{z_{1}z_{2}z_{3}z_{4}},\frac{1}{z_{1}z_{3}z_{4}},\frac{1}{z_{3}z_{4}},\frac{1}{z_{2}z_{4}},\frac{1}{z_{1}z_{4}},\frac{1}{z_{1}z_{3}}\right\\}.$
(88)

The master integrals for the last layer coincide with the one of the integral
family and consist of one box, one triangle and four bubbles. The metric
presents the following triangular structure

$\displaystyle\mathbf{C}=\left(\begin{array}[]{cccccc}C_{11}&0&0&0&0&0\\\
C_{21}&C_{22}&0&0&0&0\\\ C_{31}&C_{32}&C_{33}&0&0&0\\\
C_{41}&0&0&C_{44}&0&0\\\ C_{51}&C_{52}&0&0&C_{55}&0\\\
C_{61}&C_{62}&0&0&0&C_{66}\\\ \end{array}\right),$ (95)

where the $C_{ij}$ are rational functions of the kinematic invariants and of
the dimensional regulator. We report a simple representative reduction for
this family

$\displaystyle I[1,1,1,2]=$
$\displaystyle\frac{(d-4)\left(p_{3}^{4}-2p_{3}^{2}(p_{4}^{2}+s)+(p_{4}^{2}-s)^{2}\right)I[1,0,1,1]}{2p_{3}^{2}p_{4}^{2}s(p_{3}^{2}+p_{4}^{2}-s-t)}$
$\displaystyle+\frac{(d-3)(p_{3}^{2}-p_{4}^{2}-s)}{p_{3}^{2}p_{4}^{2}s(p_{3}^{2}+p_{4}^{2}-s-t)}I[0,0,1,1]$
$\displaystyle+\frac{(d-3)(p_{3}^{2}-p_{4}^{2}+s)}{p_{3}^{2}p_{4}^{2}s(-p_{3}^{2}-p_{4}^{2}+s+t)}I[1,0,0,1]$
$\displaystyle+\frac{(d-3)(p_{3}^{2}+p_{4}^{2}-s)}{p_{3}^{2}p_{4}^{2}s(-p_{3}^{2}-p_{4}^{2}+s+t)}I[1,0,1,0]$
$\displaystyle+\frac{(d-5)}{-p_{3}^{2}-p_{4}^{2}+s+t}I[1,1,1,1].$ (96)

###  6.2 Two-loop families

We consider the decomposition of integrals belonging to the families of two-
loop massless kite (Fig. 3), the two-loop massive sunrise with equal internal
masses (Fig. 4) and the two-loop massless pentabox on the maximal cut (Fig.
5).

#### Massless two-loop kite

Figure 3: Two loop massless kite with off-shell external legs.

We consider the family of two-loop massless kite (Fig. 3) with off-shell
external momentum $p$. The integrals depend on the invariant

$\displaystyle s=p^{2}.$ (97)

The propagators are

$\displaystyle z_{1}$ $\displaystyle=k_{1}^{2},$ $\displaystyle z_{2}$
$\displaystyle=(k_{1}-p_{1})^{2},$ $\displaystyle z_{3}=(k_{1}-k_{2})^{2},$
$\displaystyle z_{4}$ $\displaystyle=k_{2}^{2},$ $\displaystyle z_{5}$
$\displaystyle=(k_{2}-p_{1})^{2}$ (98)

and the Baikov polynomial is given by

$\displaystyle B=\frac{1}{4}\Big{(}$
$\displaystyle-s^{2}z_{3}-sz_{3}^{2}+sz_{1}z_{3}+sz_{2}z_{3}+sz_{2}z_{4}+sz_{3}z_{4}+sz_{1}z_{5}+sz_{3}z_{5}-sz_{4}z_{5}$
$\displaystyle-
sz_{1}z_{2}-z_{2}z_{4}^{2}-z_{1}z_{5}^{2}-z_{2}^{2}z_{4}+z_{1}z_{2}z_{4}-z_{1}z_{3}z_{4}+z_{2}z_{3}z_{4}-z_{1}^{2}z_{5}+z_{1}z_{2}z_{5}$
$\displaystyle+z_{1}z_{3}z_{5}-z_{2}z_{3}z_{5}+z_{1}z_{4}z_{5}+z_{2}z_{4}z_{5}\Big{)}$
(99)

We present the basis of $k$-fold master integrals for each layer $k$ of
integration, highlighting that in some cases the left and right bases may not
be the same

$\displaystyle 1-\text{fold}$
$\displaystyle\left\\{e^{(L)}_{j}\right\\}=\left\\{e^{(R)}_{j}\right\\}=\left\\{1,\frac{1}{z_{1}}\right\\},$
$\displaystyle 2-\text{fold}$
$\displaystyle\left\\{e^{(L)}_{j}\right\\}=\left\\{e^{(R)}_{j}\right\\}=\left\\{\frac{1}{z_{1}z_{2}},\frac{1}{z_{1}},\frac{1}{z_{2}},1\right\\},$
$\displaystyle 3-\text{fold}$
$\displaystyle\left\\{e^{(L)}_{j}\right\\}=\left\\{e^{(R)}_{j}\right\\}=\left\\{\frac{1}{z_{1}z_{2}z_{3}},\frac{1}{z_{2}z_{3}},\frac{1}{z_{1}z_{3}},\frac{1}{z_{1}z_{2}}\right\\},$
$\displaystyle 4-\text{fold}$
$\displaystyle\left\\{e^{(L)}_{j}\right\\}=\left\\{\frac{1}{z_{2}z_{3}z_{4}},\frac{1}{z_{1}z_{2}z_{4}},\frac{1}{z_{2}z_{3}},\frac{1}{z_{1}z_{3}},\frac{1}{z_{1}z_{2}},\frac{1}{z_{3}}\right\\}$
$\displaystyle\left\\{e^{(R)}_{j}\right\\}=\left\\{\frac{1}{z_{2}z_{3}z_{4}},\frac{1}{z_{1}z_{2}z_{4}},\frac{1}{z_{1}z_{2}z_{3}},\frac{1}{z_{2}z_{3}},\frac{1}{z_{1}z_{3}},\right\\},$
$\displaystyle 5-\text{fold}$
$\displaystyle\left\\{e^{(L)}_{j}\right\\}=\left\\{e^{(R)}_{j}\right\\}=\left\\{\frac{1}{z_{1}z_{2}z_{4}z_{5}},\frac{1}{z_{1}z_{3}z_{5}},\frac{1}{z_{2}z_{3}z_{4}}\right\\}.$
(100)

The metric for the last layer of master integrals reads

$\mathbf{C}=\left(\begin{array}[]{ccc}-\frac{(d-4)s^{2}}{4(d-5)(d-3)}&0&0\\\
0&\frac{4(d-4)^{2}s^{4}}{3(3d-16)(3d-14)(3d-10)(3d-8)}&0\\\
0&0&\frac{4(d-4)^{2}s^{4}}{3(3d-16)(3d-14)(3d-10)(3d-8)}\\\
\end{array}\right)$ (101)

and presents a diagonal structure due to the fact that no master integral
belongs to a subsector of the others. We quote a simple example of a reduction
for this family

$\displaystyle I[1,1,2,1,1]={}$
$\displaystyle-\frac{3(d-5)(d-2)(3d-10)(3d-8)}{(d-6)^{2}(d-4)s^{3}}I[0,1,1,1,0]$
$\displaystyle-\frac{3(d-5)(d-2)(3d-10)(3d-8)}{(d-6)^{2}(d-4)s^{3}}I[1,0,1,0,1]$
$\displaystyle+\frac{4(d-3)}{(d-6)s^{2}}I[1,1,0,1,1].$ (102)

#### Massive two-loop sunrise

Figure 4: Two loop massive sunrise with equal internal masses and off-shell
external momenta.

We consider the family of two-loop massive sunrise (Fig. 4) with equal-mass
internal lines and off-shell external momentum $p$. The integrals depend on
the internal mass $m$ and on the external invariant

$\displaystyle s=p^{2}.$ (103)

The propagators are

$\displaystyle z_{1}$ $\displaystyle=(k_{1}+p_{1})^{2},$ $\displaystyle z_{2}$
$\displaystyle=k_{1}^{2}-m^{2},$ $\displaystyle z_{3}$
$\displaystyle=(k_{1}+k_{2}+p_{1})^{2}-m^{2},$ $\displaystyle z_{4}$
$\displaystyle=(k_{2}+p_{1})^{2},$ $\displaystyle z_{5}$
$\displaystyle=k_{2}^{2}-m^{2}$ (104)

and the Baikov polynomial is given by

$\displaystyle B=\frac{1}{4}\Big{(}$
$\displaystyle-m^{6}+2m^{4}s-2m^{4}z_{2}+m^{4}z_{3}-2m^{4}z_{5}-m^{2}s^{2}+2m^{2}sz_{2}+2m^{2}sz_{5}$
$\displaystyle-m^{2}z_{2}^{2}-m^{2}z_{5}^{2}+m^{2}z_{1}z_{2}-m^{2}z_{1}z_{3}+m^{2}z_{2}z_{3}+3m^{2}z_{1}z_{4}-m^{2}z_{3}z_{4}$
$\displaystyle-3m^{2}z_{2}z_{5}+m^{2}z_{3}z_{5}+m^{2}z_{4}z_{5}-s^{2}z_{3}-sz_{3}^{2}-sz_{1}z_{2}+sz_{1}z_{3}+sz_{2}z_{3}$
$\displaystyle+sz_{1}z_{4}+sz_{3}z_{4}+sz_{2}z_{5}+sz_{3}z_{5}-sz_{4}z_{5}-z_{1}z_{4}^{2}-z_{2}z_{5}^{2}-z_{1}^{2}z_{4}+z_{1}z_{2}z_{4}$
$\displaystyle+z_{1}z_{3}z_{4}-z_{2}z_{3}z_{4}-z_{2}^{2}z_{5}+z_{1}z_{2}z_{5}-z_{1}z_{3}z_{5}+z_{2}z_{3}z_{5}+z_{1}z_{4}z_{5}+z_{2}z_{4}z_{5}\Big{)}.$
(105)

We consider integrals of the form in (26) with the constrain

$\alpha_{1}\leq 0,\qquad\alpha_{4}\leq 0$ (106)

i.e. we identify $z_{1}$ and $z_{4}$ as auxiliary denominators. We report the
bases for each layer of integration, highlighting that in some cases the basis
for left and right integrals may not be the same.

$\displaystyle 1-\text{fold}$
$\displaystyle\left\\{e^{(L,1)}_{j}\right\\}=\left\\{e^{(R,1)}_{j}\right\\}=\left\\{1\right\\},$
(107) $\displaystyle 2-\text{fold}$
$\displaystyle\left\\{e^{(L,2)}_{j}\right\\}=\left\\{e^{(R,2)}_{j}\right\\}=\left\\{\frac{1}{z_{2}},1\right\\},$
(108) $\displaystyle 3-\text{fold}$
$\displaystyle\left\\{e^{(L,3)}_{j}\right\\}=\left\\{e^{(R,3)}_{j}\right\\}=\left\\{\frac{1}{z_{2}z_{3}},\frac{1}{z_{3}},\frac{1}{z_{2}}\right\\},$
(109) $\displaystyle 4-\text{fold}$
$\displaystyle\left\\{e^{(L,4)}_{j}\right\\}=\left\\{\frac{1}{z_{2}^{2}},\frac{1}{z_{2}z_{3}},\frac{1}{z_{3}},\frac{1}{z_{2}},z_{3},1\right\\}$
$\displaystyle\left\\{e^{(R,4)}_{j}\right\\}=\left\\{\frac{1}{z_{2}z_{3}^{2}},\frac{1}{z_{2}^{2}z_{3}},\frac{z_{4}}{z_{2}z_{3}},\frac{1}{z_{2}z_{3}},\frac{1}{z_{3}},\frac{1}{z_{2}}\right\\},$
(110) $\displaystyle 5-\text{fold}$
$\displaystyle\left\\{e^{(L,5)}_{j}\right\\}=\left\\{e^{(R,5)}_{j}\right\\}=\left\\{\frac{1}{z_{2}z_{3}z_{5}^{2}},\frac{1}{z_{2}^{2}z_{3}z_{5}},\frac{1}{z_{2}z_{3}z_{5}},\frac{1}{z_{3}z_{5}},\frac{1}{z_{2}z_{5}},\frac{1}{z_{2}z_{3}}\right\\}.$
(111)

The $5$-fold metric presents the following block triangular structure

$\displaystyle\mathbf{C}=\left(\begin{array}[]{cccccc}C_{11}&C_{12}&C_{13}&0&0&0\\\
C_{21}&C_{22}&C_{23}&0&0&0\\\ C_{31}&C_{32}&C_{33}&0&0&0\\\
C_{41}&C_{42}&C_{43}&C_{44}&0&0\\\ C_{51}&C_{52}&C_{53}&0&C_{55}&0\\\
C_{61}&C_{62}&C_{63}&0&0&C_{66}\\\ \end{array}\right)$ (118)

Here is the example of a reduction

$\displaystyle I[0,1,1,-1,1]=$
$\displaystyle\frac{1}{48(d-5)(d-3)m^{4}s^{2}}\Big{(}-\left((d-4)(d-2)m^{8}\right)+2(d-2)(3d-14)m^{6}s$
$\displaystyle+16(d-5)(d-3)m^{4}s^{2}-2(d-2)(3d-14)m^{2}s^{3}$
$\displaystyle+(d-4)(d-2)s^{4}\Big{)}I[0,0,1,0,1]$
$\displaystyle+\frac{1}{48(d-5)(d-3)m^{4}s^{2}}\Big{(}(d-4)(d-2)m^{8}-2(d-2)(3d-14)m^{6}s$
$\displaystyle+16(d-5)(d-3)m^{4}s^{2}+2(d-2)(3d-14)m^{2}s^{3}$
$\displaystyle-(d-4)(d-2)s^{4}\Big{)}I[0,1,0,0,1]$
$\displaystyle+\frac{4m^{2}(s-m^{2})I[0,1,1,0,2]}{3(d-2)}+\frac{4m^{2}(m^{2}-s)I[0,2,1,0,1]}{3(d-2)}$
$\displaystyle+\left(m^{2}+\frac{s}{3}\right)I[0,1,1,0,1]+\frac{1}{3}I[0,1,1,0,0].$
(119)

We observe that the basis of master integrals in this example is not minimal.
Indeed, the framework of intersection theory does not always take into account
symmetry relations among integrals, but only IBP relations. This is not an
issue, because these additional relations can be easily identified and added
afterwards. In this case it is easy to see that integrals of the form
$I[0,\alpha_{2},\alpha_{3},0,\alpha_{5}]$ are invariant under permutations of
the exponents $\alpha_{j}$. This reduces the number of master integrals from
six to three, since

$I[0,1,1,0,2]=I[0,2,1,0,1],\qquad I[0,1,1,0,0]=I[0,1,0,0,1]=I[0,0,1,0,1].$
(120)

#### Massless pentabox on the maximal cut

Figure 5: Massless pentabox, in the example we consider its maximal cut.

We consider the family of massless pentabox on the maximal cut. The external
legs $p_{1},\dots,p_{5}$ satisfy

$\displaystyle p_{i}^{2}{}=0,$ (121)

and momentum conservation

$p_{1}+p_{2}+p_{3}+p_{4}+p_{5}=0.$ (122)

The integrals belonging to this family depend on five independent invariants
which we choose as

$\left\\{s_{12},s_{23},s_{34},s_{45},s_{51}\right\\},$ (123)

where

$s_{ij}=(p_{i}+p_{j})^{2}.$ (124)

The pentabox family depends on $11$ generalized denominators

$\displaystyle z_{1}$ $\displaystyle=(k_{2}+p_{1}+p_{2})^{2},$ $\displaystyle
z_{2}$ $\displaystyle=(k_{2}+p_{1})^{2},$ $\displaystyle z_{3}$
$\displaystyle=(k_{1}+p_{5})^{2},$ $\displaystyle z_{4}$
$\displaystyle=k_{1}^{2},$ $\displaystyle z_{5}$
$\displaystyle=(k_{1}-p_{1})^{2},$ $\displaystyle z_{6}$
$\displaystyle=(k_{1}-p_{1}-p_{2})^{2},$ $\displaystyle z_{7}$
$\displaystyle=(k_{1}-p_{1}-p_{2}-p_{3})^{2},$ $\displaystyle z_{8}$
$\displaystyle=k_{2}^{2},$ $\displaystyle z_{9}$
$\displaystyle=(k_{2}+p_{1}+p_{2}+p_{3}+p_{4})^{2},$ $\displaystyle z_{10}$
$\displaystyle=(k_{2}+p_{1}+p_{2}+p_{3})^{2},$ $\displaystyle z_{11}$
$\displaystyle=(k_{1}+k_{2})^{2}.$ (125)

of which the first three are irreducible scalar products, i.e. we consider
integrals of the form (26) with

$\alpha_{1}\leq 0,\quad\alpha_{2}\leq 0,\quad\alpha_{3}\leq 0.$ (126)

We consider the reduction of integrals belonging to the maximal cut of this
family, which corresponds to putting on-shell all the denominators via the
substitution

$\displaystyle\frac{1}{z_{i}}\rightarrow\delta(z_{i})\quad(i=4,\dots,11)$
(127)

in the expressions for $\varphi_{R,L}$. This still yields the correct
coefficients of reductions multiplying the master integrals of the top sector.
Therefore the integrals we reduce are function only of the auxiliary
denominators $z_{1},z_{2},z_{3}$, meaning that we calculate intersection
numbers of three-fold integrals. The Baikov polynomial of pentabox on the
maximal-cut is given by

$\displaystyle B=\frac{1}{64}$
$\displaystyle\Big{(}-z_{3}(s_{12}^{2}z_{3}(s_{23}-z_{2})^{2}+2s_{12}(-s_{23}^{2}z_{1}z_{3}+s_{23}(s_{45}z_{2}z_{3}+s_{45}z_{1}(2z_{2}+z_{3})$
$\displaystyle+z_{1}z_{2}z_{3})+s_{45}(z_{1}-z_{2})z_{2}z_{3})+z_{3}(s_{23}z_{1}+s_{45}(z_{2}-z_{1}))^{2})$
$\displaystyle+2s_{15}s_{45}z_{1}(z_{3}(s_{12}(s_{23}+z_{2})-s_{23}z_{1}+s_{45}z_{1}-s_{45}z_{2})+s_{34}s_{45}z_{2})$
$\displaystyle+2s_{34}s_{45}z_{2}z_{3}(s_{12}(s_{23}-z_{2})+s_{23}z_{1}+s_{45}(z_{2}-z_{1}))-s_{15}^{2}s_{45}^{2}z_{1}^{2}-s_{34}^{2}s_{45}^{2}z_{2}^{2}\Big{)},$
(128)

and its exponent is $\gamma=1/2(d-7)$. The bases of master integrals for each
layer of integration are

$\displaystyle 1-\text{fold}$
$\displaystyle\left\\{e^{(L)}_{j}\right\\}=\left\\{e^{(R)}_{j}\right\\}=\left\\{1\right\\},$
$\displaystyle 2-\text{fold}$
$\displaystyle\left\\{e^{(L)}_{j}\right\\}=\left\\{e^{(R)}_{j}\right\\}=\left\\{1\right\\},$
$\displaystyle 3-\text{fold}$
$\displaystyle\left\\{e^{(L)}_{j}\right\\}=\left\\{e^{(R)}_{j}\right\\}=\left\\{z_{2},z_{3},1\right\\}.$
(129)

Since we are on the maximal cut, all the integrals we are considering belong
to the same sector and thus all entries of the metric $\mathbf{C}_{ij}$ are
non-vanishing. We reconstructed the reduction of all integrals up to degree 5
in the variables $z_{i}$ (which are those that contribute e.g. to QCD
amplitudes) and successfully checked them against the Laporta method.

##  7 Conclusions and Outlook

We presented a new method for computing intersection numbers via a purely
rational algorithm that does not require any change of basis or integral
transformation. This is achieved via the systematic use of the $p(z)$-adic
series expansion. The latter expands functions in powers of a polynomial,
allowing to study their behaviour close to potentially irrational singular
points, without having to perform any irrational operation nor knowing the
explicit location of these points.

This result represents significant progress in the application of intersection
theory to the decomposition of Feynman integrals. The new algorithm is elegant
and satisfactory from a mathematical and theoretical point of view, since
everything is rational in all steps of the calculation. It also sidesteps the
algebraic bottleneck of dealing with irrational expressions and opens up the
possibility of combining intersection theory with efficient computational
techniques, such as finite fields and functional reconstruction. As a proof of
concept, we implemented the new algorithm over finite fields, using the
Mathematica package FiniteFlow and tested it on several one- and two-loop
examples.

We also proposed a new strategy for dealing with analytic regulators, de facto
removing any analytic dependence of the calculation on them and showing
dramatic simplifications on both analytic expressions and the structure of the
results.

These results open up several possible lines of research aimed at making
intersection theory a viable approach for phenomenological applications.
Besides optimizing and simplifying our proof-of-concept implementation,
exploring different integral representations which reduce the number of
integration variables may lead to substantial improvements. Examples are the
loop-by-loop Baikov representation and the Lee-Pomeransky representation
Lee:2013hzt . Even more interestingly, a method for the direct calculation of
multivariate intersection numbers – thus stepping away from the recursive
univariate method – has recently been proposed Chestnov:2022xsy . A possible
generalization of our method that builds on Chestnov:2022xsy  and a
multivariate generalization of the $p(z)$-adic expansion, combined with the
multivariate global residue theorem, is worth investigating in the near
future.

Finally, we expect the method of $p(z)$-adic expansions may also find uses in
a broader spectrum of applications, where one needs to study the behaviour of
certain functions close to the zeroes of certain polynomials. For instance,
they might be employed in studies of the singular structure of amplitudes for
a more efficient reconstruction of them – as recently done using the similar
concept of $p$-adic numbers DeLaurentis:2022otd . The possibility of
systematically employing $p(z)$-adic expansions using publicly available
finite-field frameworks Peraro:2019svx  and the fact that these are not
affected by numerical instabilities, makes this technique a suitable candidate
for future theoretical and phenomenological investigations.

##  Appendix A $p(z)$-adic expansions via polynomial division

We illustrate how we can obtain the coefficients $c_{i}(z)$ of the $p(z)$-adic
expansion in (30) as reminders of polynomial divisions. Since $f(z)$ is a
rational function, we can write it as

$f(z)=\frac{n(z)}{p(z)^{k}\,d(z)},$ (130)

where $d(z)$ is co-prime (i.e. has no common factor) with $p(z)$. Since the
factor $1/p(z)^{k}$ can be trivially added back after the expansion, it is
sufficient to consider a function of the form

$f^{\star}(z)=\frac{n(z)}{d(z)},$ (131)

whose denominator $d(z)$ is co-prime with $p(z)$.

Let $\tilde{d}(z)$ be the multiplicative inverse of $d(z)$ with respect to
$p(z)$, namely the polynomial of degree lower than $\deg p$ such that

$d(z)\,\tilde{d}(z)=1\mod p(z).$ (132)

One can show that, if $d(z)$ and $p(z)$ are co-prime, the multiplicative
inverse of $d(z)$ exists and is unique. It can be computed using the extended
Euclidean algorithm, or simply by making an ansatz for it and fixing the
unknown coefficients using eq. (132).

The first term of the $p(z)$-adic expansion of $f^{\star}(z)$ is simply its
polynomial reminder modulo $p(z)$

$c_{0}(z)\equiv
f^{\star}(z)\,\text{mod}\,p(z)=n(z)\tilde{d}(z)\,\text{mod}\,p(z).$ (133)

We can thus rewrite the numerator $n(z)$ as

$n(z)=c_{0}(z)d(z)+q_{1}(z)p(z),$ (134)

for some polynomial $q_{1}(z)$, hence

$f^{\star}(z)=c_{0}(z)+\frac{q_{1}(z)}{d(z)}\,p(z).$ (135)

To find the next term in the expansion we repeat the procedure considering
$q_{1}(z)/d(z)$ as a new rational function. In general we have the recursive
formula

$c_{i}(z)=\frac{q_{i}(z)}{d(z)}\,\text{mod}\,p(z),$ (136)

where $q_{i}(z)$ for $i>0$ can be computed from $q_{i-1}(z)$ and $c_{i-1}(z)$
using

$q_{i-1}(z)=c_{i-1}(z)d(z)+q_{i}(z)p(z),$ (137)

and $q_{0}(z)=n(z)$. Following this procedure we can compute $c_{i}(z)$ up to
the order we need.

Readers familiar with the concept of $p$-adic numbers may recognize that the
construction of the $p(z)$-adic expansion of a rational function is analogous
to the construction of the $p$-adic expansion of a rational number, where $p$
is a prime number rather than a prime polynomial over the rational field and
the remainder of integer division is used instead of the polynomial remainder.

##  Appendix B The univariate global residue theorem

We give a brief overview of the univariate global residue theorem, whose
generalization is a key ingredient of the calculation of intersection numbers
with the rational algorithm, as described in subsection 3.1. Given a rational
function

$f(z)=\frac{n(z)}{d(z)},$ (138)

whose denominator is coprime with a polynomial $p(z)$, we can calculate the
polynomial reminder of $f(z)$ modulus $p(z)$ with eq. (133). It takes the
explicit form

$\tilde{f}(z)\equiv
f(z)\,\,\text{mod}\,\,p(z)=\sum_{j=1}^{\text{deg}\,p-1}f_{j}z^{j}.$ (139)

The univariate global residue theorem states that the global residue of
$f(z)/p(z)$, namely the sum of all the local residues of $f(z)/p(z)$ at the
zeroes of $p(z)$, is

$\text{Res}_{p(z)}\left(\frac{f(z)}{p(z)}\right)=\frac{f_{\text{deg}\,p-1}}{l_{c}},$
(140)

where $\text{Res}_{p(z)}$ is defined in eq. (36).

###### Acknowledgements.

We thank Vsevolod Chestnov, Federico Gasparotto and Pierpaolo Mastrolia for
many valuable discussions and comments on this work. We also grateful to
Thomas Gehrmann for feedback on the draft. This work was supported by the
European Research Council (ERC) under the European Union’s Horizon $2020$
research and innovation programme grant agreements $101019620$ (ERC Advanced
Grant TOPUP) and 101040760 (ERC Starting Grant FFHiggsTop).

## References

  * (1) F. V. Tkachov, A Theorem on Analytical Calculability of Four Loop Renormalization Group Functions, Phys. Lett. B 100 (1981) 65–68.
  * (2) K. G. Chetyrkin and F. V. Tkachov, Integration by Parts: The Algorithm to Calculate beta Functions in 4 Loops, Nucl. Phys. B 192 (1981) 159–204.
  * (3) T. Gehrmann and E. Remiddi, Differential equations for two loop four point functions, Nucl. Phys. B 580 (2000) 485–518, [hep-ph/9912329].
  * (4) S. Laporta, High precision calculation of multiloop Feynman integrals by difference equations, Int. J. Mod. Phys. A 15 (2000) 5087–5159, [hep-ph/0102033].
  * (5) A. von Manteuffel and C. Studerus, Reduze 2 - Distributed Feynman Integral Reduction, arXiv:1201.4330.
  * (6) A. V. Smirnov, Algorithm FIRE – Feynman Integral REduction, JHEP 10 (2008) 107, [arXiv:0807.3243].
  * (7) P. Maierhöfer, J. Usovitsch, and P. Uwer, Kira—A Feynman integral reduction program, Comput. Phys. Commun. 230 (2018) 99–112, [arXiv:1705.05610].
  * (8) A. von Manteuffel and R. M. Schabinger, A novel approach to integration by parts reduction, Phys. Lett. B 744 (2015) 101–104, [arXiv:1406.4513].
  * (9) T. Peraro, Scattering amplitudes over finite fields and multivariate functional reconstruction, JHEP 12 (2016) 030, [arXiv:1608.01902].
  * (10) J. Klappert and F. Lange, Reconstructing rational functions with FireFly, Comput. Phys. Commun. 247 (2020) 106951, [arXiv:1904.00009].
  * (11) A. V. Smirnov and F. S. Chuharev, FIRE6: Feynman Integral REduction with Modular Arithmetic, Comput. Phys. Commun. 247 (2020) 106877, [arXiv:1901.07808].
  * (12) T. Peraro, FiniteFlow: multivariate functional reconstruction using finite fields and dataflow graphs, JHEP 07 (2019) 031, [arXiv:1905.08019].
  * (13) J. Klappert, F. Lange, P. Maierhöfer, and J. Usovitsch, Integral reduction with Kira 2.0 and finite field methods, Comput. Phys. Commun. 266 (2021) 108024, [arXiv:2008.06494].
  * (14) A. V. Smirnov and A. V. Petukhov, The Number of Master Integrals is Finite, Lett. Math. Phys. 97 (2011) 37–44, [arXiv:1004.4199].
  * (15) S. Mizera, Scattering Amplitudes from Intersection Theory, Phys. Rev. Lett. 120 (2018), no. 14 141602, [arXiv:1711.00469].
  * (16) P. Mastrolia and S. Mizera, Feynman Integrals and Intersection Theory, JHEP 02 (2019) 139, [arXiv:1810.03818].
  * (17) H. Frellesvig, F. Gasparotto, S. Laporta, M. K. Mandal, P. Mastrolia, L. Mattiazzi, and S. Mizera, Decomposition of Feynman Integrals on the Maximal Cut by Intersection Numbers, JHEP 05 (2019) 153, [arXiv:1901.11510].
  * (18) H. Frellesvig, F. Gasparotto, M. K. Mandal, P. Mastrolia, L. Mattiazzi, and S. Mizera, Vector Space of Feynman Integrals and Multivariate Intersection Numbers, Phys. Rev. Lett. 123 (2019), no. 20 201602, [arXiv:1907.02000].
  * (19) S. Mizera, Aspects of Scattering Amplitudes and Moduli Space Localization.  PhD thesis, Princeton, Inst. Advanced Study, 2020.  arXiv:1906.02099.
  * (20) S. Mizera and A. Pokraka, From Infinity to Four Dimensions: Higher Residue Pairings and Feynman Integrals, JHEP 02 (2020) 159, [arXiv:1910.11852].
  * (21) H. Frellesvig, F. Gasparotto, S. Laporta, M. K. Mandal, P. Mastrolia, L. Mattiazzi, and S. Mizera, Decomposition of Feynman Integrals by Multivariate Intersection Numbers, JHEP 03 (2021) 027, [arXiv:2008.04823].
  * (22) S. Weinzierl, On the computation of intersection numbers for twisted cocycles, J. Math. Phys. 62 (2021), no. 7 072301, [arXiv:2002.01930].
  * (23) V. Chestnov, F. Gasparotto, M. K. Mandal, P. Mastrolia, S. J. Matsubara-Heo, H. J. Munch, and N. Takayama, Macaulay matrix for Feynman integrals: linear relations and intersection numbers, JHEP 09 (2022) 187, [arXiv:2204.12983].
  * (24) V. Chestnov, H. Frellesvig, F. Gasparotto, M. K. Mandal, and P. Mastrolia, Intersection Numbers from Higher-order Partial Differential Equations, arXiv:2209.01997.
  * (25) S. Caron-Huot and A. Pokraka, Duals of Feynman integrals. Part I. Differential equations, JHEP 12 (2021) 045, [arXiv:2104.06898].
  * (26) S. Caron-Huot and A. Pokraka, Duals of Feynman Integrals. Part II. Generalized unitarity, JHEP 04 (2022) 078, [arXiv:2112.00055].
  * (27) S. Weinzierl, Feynman Integrals.  1, 2022.
  * (28) S. L. Cacciatori, M. Conti, and S. Trevisan, Co-Homology of Differential Forms and Feynman Diagrams, Universe 7 (2021), no. 9 328, [arXiv:2107.14721].
  * (29) S. L. Cacciatori and P. Mastrolia, Intersection Numbers in Quantum Mechanics and Field Theory, arXiv:2211.03729.
  * (30) J. Chen, X. Jiang, X. Xu, and L. L. Yang, Constructing canonical Feynman integrals with intersection theory, Phys. Lett. B 814 (2021) 136085, [arXiv:2008.03045].
  * (31) J. Chen, X. Jiang, C. Ma, X. Xu, and L. L. Yang, Baikov representations, intersection theory, and canonical Feynman integrals, JHEP 07 (2022) 066, [arXiv:2202.08127].
  * (32) G. Fontana, Rational algorithms for the decomposition of feynman integrals via intersection theory, Master’s thesis, 2022.
  * (33) P. A. Baikov, Explicit solutions of the multiloop integral recurrence relations and its application, Nucl. Instrum. Meth. A 389 (1997) 347–349, [hep-ph/9611449].
  * (34) R. N. Lee and A. A. Pomeransky, Critical points and number of master integrals, JHEP 11 (2013) 165, [arXiv:1308.6676].
  * (35) G. De Laurentis and B. Page, Ansätze for scattering amplitudes from p-adic numbers and algebraic geometry, JHEP 12 (2022) 140, [arXiv:2203.04269].