# Feynman parameter integration through differential equations

Martijn Hidding  Johann Usovitsch

###### Abstract

We present a new method for numerically computing generic multi-loop Feynman
integrals. The method relies on an iterative application of Feynman’s trick
for combining two propagators. Each application of Feynman’s trick introduces
a simplified Feynman integral topology which depends on a Feynman parameter
that should be integrated over. For each integral family, we set up a system
of differential equations which we solve in terms of a piecewise collection of
generalized series expansions in the Feynman parameter. These generalized
series expansions can be efficiently integrated term by term, and segment by
segment. This approach leads to a fully algorithmic method for computing
Feynman integrals from differential equations, which does not require the
manual determination of boundary conditions. Furthermore, the most complicated
topology that appears in the method often has less master integrals than the
original one. We illustrate the strength of our method with a five-point two-
loop integral family.

##  1 Introduction

The computation of multi-loop Feynman integrals is a crucial component in
generating predictions for particle processes at high-energy colliders such as
the LHC. Many modern techniques for computing Feynman integrals rely on the
differential equation method [1, 2]. This method involves setting up a system
of differential equations for the master integrals of a Feynman integral
family, and solving the differential equations either analytically or
numerically. Analytic solutions can often be found in terms of special classes
of iterated integrals such as multiple polylogarithms [3, 4] or
generalizations thereof (see e.g. [5, 6]), but are hard to obtain in general.
In addition, there are many Feynman integrals which do not evaluate to well-
studied classes of iterated integrals. Therefore, renewed interest has been
expressed in solving differential equations numerically without reference to
an intermediate function space.

Many works in the literature have studied methods for solving differential
equations for Feynman integrals numerically (see e.g. [7, 8, 9, 10, 11]). In
this paper, we employ the strategy outlined in [12]. In this approach,
differential equations are repeatedly solved in terms of generalized power
series expansions (which may contain powers of logarithms) along connected
line segments in phase-space. The Mathematica package DiffExp [13] provides a
public implementation of the method that works on user-provided differential
equations. (Another recent public implementation is the SeaSyde [14]
Mathematica package.) The approach yields a precise and efficient way of
computing Feynman integrals at any point in phase-space, given an initial set
of boundary conditions as input. The computation of boundary conditions can be
performed analytically in special limits with the aid of expansion by regions
[15, 16, 17], but generally requires some manual work.

In [18] the Mathematica package AMFlow was presented which is based on the
method of auxiliary mass flow [8, 19]. The package uses expansion by regions
to automatically determine boundary conditions in a special limit of infinite
complex mass, and a customized series solution solver to transport the
boundary conditions to physical points in phase-space. The result is an
efficient and fully automated package for computing Feynman integrals
numerically at high precision. In [20] we explored a different automated
strategy which relies on sector decomposition [21, 22, 23] to compute boundary
conditions in the Euclidean region for a special basis of (quasi-)finite
master integrals [24]. We then obtained results in physical regions by
transporting the boundary conditions with DiffExp. This strategy is fully
automated as well, but it is not able to obtain results at a precision as high
as AMFlow.

The main result of the present paper is a new automated approach for computing
Feynman integrals numerically, which relies on repeatedly combining two
propagators using Feynman’s trick:

$\displaystyle\frac{1}{D_{i}^{\nu_{i}}D_{j}^{\nu_{j}}}=\frac{\Gamma(\nu_{i}+\nu_{j})}{\Gamma(\nu_{i})\Gamma(\nu_{j})}\int_{0}^{1}dx\,\frac{x^{\nu_{i}-1}(1-x)^{\nu_{j}-1}}{(D_{i}x+D_{j}(1-x))^{\nu_{i}+\nu_{j}}}\,,$
(1)

where $\nu_{i}$ and $\nu_{j}$ are assumed to be positive integers. We will
show that the integration with respect to the Feynman parameter $x$ can be
performed by solving an associated system of differential equations in terms
of generalized series expansions, and by integrating the series expansions
term by term. This method has the computational advantage that the first
application of Feynman’s trick typically gives a simplified integral family
with less master integrals than the original family [25]. This is generally
not the case for the auxiliary mass flow method, which introduces (complex)
masses on the propagators. The reduction of the number of master integrals in
our approach simplifies the integration by parts (IBP) reduction, which can
present a bottleneck in complex calculations.

We remark that the use of (1) was also considered in [26] and [25] in the
context of simplifying Feynman integral computations. In [26], the trick was
used to study differential equations of elliptic Feynman integrals by
deferring the elliptic type integration to the last integration parameter. In
[25], the formula was shown to simplify an integral family when the combined
propagators $D_{i}$ and $D_{j}$ have the same internal loop momentum, and the
method was called the ‘internal reduction method’. Furthermore, it was already
observed in [25] that this procedure reduces the number of master integrals.

This paper expands on previous ideas in a number of novel ways. Firstly, we
discuss how (1) can be employed recursively, and how the necessary
integrations can be performed by solving differential equations with series
expansion methods. Furthermore, we present a formula for resolving poles in
dimensional regularization which is needed to regulate the integrand, and we
discuss how to apply it to an integrand which is in a generalized series
representation. We also discuss some first steps into obtaining
$i\delta$-prescriptions for the Feynman parameters, in order to obtain results
outside the Euclidean region.

The remainder of the paper is organized as follows. In section 2 we set up
conventions, and we review the Feynman parametrization and Feynman’s trick for
combining propagators. In section 3 we describe the heart of our method, which
is to perform the integrals in (13) by making use of differential equations
and series solution methods. We describe in section 3.2 how to regulate the
Feynman parameter integrals in our approach when the master integrals are not
finite in $\varepsilon$, and in section 3.3 we discuss first steps into
crossing threshold singularities in our approach, in order to obtain results
directly in physical regions. Lastly, we discuss in section 4 two non-trivial
pedagogical two-loop examples which serve as a proof of concept for our
method. Finally, we summarize our methods and results in section 5.

##  2 Preliminaries

Here we review a few general definitions of Feynman integrals. For a thorough
review, see e.g. [27]. The main object of interest is a family of Feynman
integrals given by

$\displaystyle
I_{\nu_{1},\ldots,\nu_{n}}\left(\\{s_{i},m_{i}\\},d\right)=\int\left(\prod_{j=1}^{l}\frac{d^{d}k_{j}}{i\pi^{\frac{d}{2}}}\right)\,\prod_{j=1}^{n}\frac{1}{D_{j}^{\nu_{j}}}\,.$
(2)

Here $\\{s_{i},m_{i}\\}$ schematically denotes the set of external scales and
internal masses on which the integral family depends. We work in dimensional
regularization $d=d_{\mathrm{int}}-2\varepsilon$, where $d_{\mathrm{int}}$ is
a positive (even) integer. In the remainder of this paper, we will typically
leave the dependence on the external scales, masses and dimension implicit.
The propagator exponents $\nu_{i}$ are assumed to be integers, and we will
denote their sum by $\nu=\nu_{1}+\ldots+\nu_{n}$. Typically some of the
exponents are taken to be negative, which gives numerator factors. We use the
convention $D_{j}=-q_{j}^{2}+m_{j}^{2}-i\delta$, where $q_{j}$ denotes the
momentum of the $j$-th propagator. The momentum $q_{j}$ can be written as a
linear combination of $l$ loop momenta $k_{i}$ and $E$ external momenta
$p_{j}$, such that

$\displaystyle q_{j}=\sum_{n=1}^{l}a_{jn}k_{n}+\sum_{n=1}^{E}b_{jn}p_{n}\,,$
(3)

for some integers $a_{jn}$ and $b_{jn}$. The Feynman $i\delta$-prescription
tells us how to deform around poles in the integration region, and can be
viewed as adding an infinitesimally small complex mass in the lower half-plane
to each propagator.

It is well known that IBP relations exist between Feynman integrals whose
propagator exponents are related by integer shifts [28, 29]. Furthermore,
Feynman integrals in different dimensions may be related through dimensional
raising and lowering relations [30, 31]. In order to perform IBP reductions in
closed-form, the propagators $D_{j}$ should form a basis for the span of dot
products of the form $k_{i}\cdot k_{j}$, and $k_{i}\cdot p_{j}$. A minimal set
of independent integrals may be found in terms of which all other integrals
can be expressed, which is called a basis of master integrals. Finding the
basis of master integrals can be done through IBP-reduction programs such as
LiteRed [32], Reduze [33], FIRE [34] and KIRA [35]. The computational
complexity of the IBP reductions is typically closely connected to the number
of master integrals that the integral family has.

###  2.1 Feynman parametrization

It is common to write Feynman integrals in the so-called Feynman
parametrization. For brevity, we will assume in the following that all
propagator exponents $\nu_{1},\ldots,\nu_{n}$ are positive, and that we do not
necessarily have a basis of propagators and numerators suitable for IBP
reduction. For the case with numerators, we refer the reader to [32]. First,
we consider the following formula, which was first observed by Feynman, and
which we will refer to as Feynman’s trick:

$\displaystyle\frac{1}{D_{1}^{\nu_{1}}\ldots
D_{n}^{\nu_{n}}}=\frac{\Gamma(\nu)}{\Gamma(\nu_{1})\ldots\Gamma(\nu_{n})}\int_{0}^{1}d^{n}\vec{x}\,\frac{x_{1}^{\nu_{1}-1}\ldots
x_{n}^{\nu_{n}-1}\delta\left(1-\sum_{j=1}^{n}x_{j}\right)}{(x_{1}D_{1}+\ldots+x_{n}D_{n})^{\nu}}\,,$
(4)

where $\nu=\nu_{1}+\ldots+\nu_{n}$. Upon plugging the formula into (2), we may
integrate out the loop momenta, which gives the Feynman parametrization:

$\displaystyle I_{\nu_{1},\ldots,\nu_{n}}=\frac{\Gamma(\nu-
ld/2)}{\prod_{j=1}^{n}\Gamma\left(\nu_{j}\right)}\int\left(\prod_{j=1}^{n}dx_{j}\,x_{j}^{\nu_{j}-1}\right)\frac{\mathcal{U}^{\nu-(l+1)d/2}}{\mathcal{F}^{\nu-
ld/2}}\delta\left(1-\sum_{j=1}^{n}x_{j}\right)\,.$ (5)

The graph polynomial $\mathcal{U}$ is of degree $l$ in the Feynman parameters
$x_{i}$, and $\mathcal{F}$ is of degree $l+1$. These are also called the first
and second Symanzik polynomials. The Symanzik polynomials can be constructed
from the graph $G$ of the Feynman integral. It holds that

$\displaystyle\mathcal{U}$ $\displaystyle=\sum_{T\in T(G)}\prod_{e_{i}\notin
T}x_{i}\,,\quad\quad\mathcal{F}_{0}=\sum_{\left(T_{1},T_{2}\right)\in
T_{2}(G)}\left(\prod_{e_{i}\notin\left(T_{1}\cup
T_{2}\right)}x_{i}\right)s_{\left(T_{1},T_{2}\right)},\quad$
$\displaystyle\mathcal{F}$
$\displaystyle=-\mathcal{F}_{0}+\mathcal{U}\left(\sum
x_{i}m_{i}^{2}\right)-i\delta\,.$ (6)

The sum goes over the spanning trees $T(G)$ of $G$, and the spanning 2-forests
$T_{2}(G)$ of $G$. We have added a $-i\delta$ to $\mathcal{F}$, which takes
care of the Feynman prescription.

More generally, we may wish to consider Feynman integrals with generalized
propagators, which are quadratic functions in the internal and external
momenta. In this case a graph representation may not exist for the Feynman
integral, but the Symanzik polynomials can still be defined through
determinants. We refer the reader to [36, 27] for explicit formulas.

###  2.2 Iterating Feynman’s trick

Let us derive an alternative form of the Feynman parametrization by repeatedly
applying (4) to combine two propagators. For example, we may choose to combine
every time the leftmost two propagators, which we will do next. More
generally, we may choose to combine propagators in some other sequence, which
we will do for our examples. Let us consider the following set of generalized
propagators

$\displaystyle D_{12}$ $\displaystyle=x_{1}D_{1}+(1-x_{1})D_{2}\,,$
$\displaystyle D_{123}$ $\displaystyle=x_{2}D_{12}+(1-x_{2})D_{3}\,,$
$\displaystyle\ldots$ $\displaystyle D_{1\ldots n}$
$\displaystyle=x_{n-1}D_{1\ldots(n-1)}+(1-x_{n-1})D_{n}\,.$ (7)

Next, we define a set of integral families $I_{\vec{\nu}}^{(\kappa)}$, such
that for $\kappa=0$ we have the original family, and for each successive step
$\kappa$ we combine two propagators. Like before, we will assume in the
following that all propagator exponents are positive, and we leave out the
numerator factors which are needed to obtain a basis for IBP reductions. We
define the integral families by

$\displaystyle I_{\nu_{1},\ldots,\nu_{n-\kappa}}^{(\kappa)}$
$\displaystyle=\int\left(\prod_{j=1}^{l}\frac{d^{d}k_{j}}{i\pi^{\frac{d}{2}}}\right)D_{1\ldots(\kappa+1)}^{-\nu_{1}}\prod_{j=\kappa+2}^{n}D_{j}^{-\nu_{j-\kappa}}\quad\text{
for }0\leq\kappa<n-2,$ $\displaystyle I_{\nu}^{(n-1)}$
$\displaystyle=\int\left(\prod_{j=1}^{l}\frac{d^{d}k_{j}}{i\pi^{\frac{d}{2}}}\right)D_{1\ldots
n}^{-\nu}\,.$ (8)

We have chosen our conventions such that the integrals
$I_{\vec{\nu}}^{(\kappa)}$ depend on the (unintegrated) Feynman parameters
$x_{1},\ldots,x_{\kappa}$. We remark that an integral
$I_{\vec{\nu}}^{(\kappa)}$ can be represented as a Feynman diagram if we
combine propagators with the same internal momentum [25]. For example, suppose
that

$\displaystyle D_{1}=-(k+p)^{2}+m_{1}^{2}\,,\quad\quad
D_{2}=-(k+q)^{2}+m_{2}^{2}\,.$ (9)

Then we have

$\displaystyle xD_{1}+(1-x)D_{2}=-(k+P)^{2}+M^{2}\,,$ (10)

where $P=xp+(1-x)q$ and $M^{2}=xm_{1}^{2}+(1-x)m_{2}^{2}-x(1-x)(p-q)^{2}$.
Note that $I_{\nu}^{(n-1)}$ can be loosely seen as a generalized tadpole
diagram, which evaluates to:

$\displaystyle I_{\nu}^{(n-1)}$
$\displaystyle=\int\left(\prod_{j=1}^{l}\frac{d^{d}k_{j}}{i\pi^{\frac{d}{2}}}\right)D_{1\ldots
n}^{-\nu}=\frac{\Gamma(\nu-
ld/2)}{\Gamma\left(\nu\right)}\frac{\tilde{\mathcal{U}}^{\nu-(l+1)d/2}}{\tilde{\mathcal{F}}^{\nu-
ld/2}}\,,$ (11)

where $\tilde{\mathcal{U}}$ and $\tilde{\mathcal{F}}$ are obtained from
$\mathcal{U}$ and $\mathcal{F}$ by letting:

$\displaystyle x_{1}$ $\displaystyle\rightarrow
x_{1}^{\prime}=\prod_{i=1}^{n-1}x_{i}\,,$ $\displaystyle x_{j}$
$\displaystyle\rightarrow
x_{j}^{\prime}=(1-x_{j-1})\prod_{i=j}^{n-1}x_{i}\quad\text{ for
}j=2,\ldots,n-1\,,$ $\displaystyle x_{n}$ $\displaystyle\rightarrow
x_{n}^{\prime}=(1-x_{n-1})\,.$ (12)

In short, this follows from observing that $D_{1\ldots
n}=x_{1}^{\prime}D_{1}+\ldots+x_{n}^{\prime}D_{n}$. More generally, if we
choose to repeatedly combine propagators in some different way, we can look at
the generalized tadpole propagator to understand how the Feynman parameters
are rescaled from the standard Feynman parametrization.

Using (4), we may write

$\displaystyle
I_{\nu_{1},\ldots,\nu_{n-(\kappa-1)}}^{(\kappa-1)}=\frac{\Gamma(\nu_{1}+\nu_{2})}{\Gamma(\nu_{1})\Gamma(\nu_{2})}\int_{0}^{1}dx_{\kappa}\,x_{\kappa}^{\nu_{1}-1}(1-x_{\kappa})^{\nu_{2}-1}I_{\nu_{1}+\nu_{2},\nu_{3},\ldots,\nu_{n-\kappa}}^{(\kappa)}\,,$
(13)

where $\nu_{1}$ and $\nu_{2}$ are assumed to be positive. If we iterate the
recursion, we obtain

$\displaystyle I_{\nu_{1},\ldots,\nu_{n}}$
$\displaystyle=\frac{\Gamma(\nu)}{\Gamma\left(\nu_{1}\right)\ldots\Gamma\left(\nu_{n}\right)}\left(\prod_{j=1}^{n-1}\int_{0}^{1}dx_{j}\,x_{j}^{\mu_{j}-1}\left(1-x_{j}\right)^{\nu_{j+1}-1}\right)I_{\nu}^{(n-1)}\,,$
$\displaystyle=\frac{\Gamma(\nu-
ld/2)}{\Gamma\left(\nu_{1}\right)\ldots\Gamma\left(\nu_{n}\right)}\left(\prod_{j=1}^{n-1}\int_{0}^{1}dx_{j}\,x_{j}^{\mu_{j}-1}\left(1-x_{j}\right)^{\nu_{j+1}-1}\right)\frac{\tilde{\mathcal{U}}^{\nu-(l+1)d/2}}{\tilde{\mathcal{F}}^{\nu-
ld/2}}\,,$ (14)

where $\mu_{j}=\nu_{1}+\ldots+\nu_{j}$, and $\nu=\mu_{n}$. This result could
have also been obtained directly from (5) by performing the change of
variables in (2.2) (including the Jacobian determinant). This makes it clear
that (13) implements a recursion identity that leads to (an alternate form of)
the Feynman parametrization. Note that we have the following special cases
when $\nu_{1}$ and/or $\nu_{2}$ is zero

$\displaystyle I_{0,0,\nu_{3},\ldots,\nu_{n-(\kappa-1)}}^{(\kappa-1)}$
$\displaystyle=I_{0,\nu_{3},\ldots,\nu_{n-\kappa}}^{(\kappa)}\,,$
$\displaystyle I_{\nu_{1},0,\nu_{3},\ldots,\nu_{n-(\kappa-1)}}^{(\kappa-1)}$
$\displaystyle=\lim_{x_{\kappa}\rightarrow
1}I_{\nu_{1},\nu_{3},\ldots,\nu_{n-\kappa}}^{(\kappa)}\,,$ $\displaystyle
I_{0,\nu_{2},\nu_{3},\ldots,\nu_{n-(\kappa-1)}}^{(\kappa-1)}$
$\displaystyle=\lim_{x_{\kappa}\rightarrow
0}I_{\nu_{2},\nu_{3},\ldots,\nu_{n-\kappa}}^{(\kappa)}\,.$ (15)

We elaborate on these limits at the end of section 3.2.

##  3 Direct integration through differential equations

###  3.1 Main method

Now that we have obtained the recursion identities in (13) and (2.2), we will
discuss how each step in the recursion can be performed by solving a system of
differential equations. First, we add a set of numerator factors to each
family $I_{\vec{\nu}}^{(\kappa)}$, in order to have a complete linearly
independent set of propagators for the IBP reduction. Let us denote the set of
master integrals by $\vec{I}^{(\kappa)}$. We may then obtain a closed-form
linear system of differential equations of the form [1, 37]

$\displaystyle\partial_{x_{\kappa}}\vec{I}^{(\kappa)}=M_{x_{\kappa}}\vec{I}^{(\kappa)}\,.$
(16)

Note that the master integrals $\vec{I}^{(\kappa)}$ and the corresponding
differential matrix $M_{x_{\kappa}}$, depend on Feynman parameters
$x_{1},\ldots,x_{\kappa}$, and on the external scales and internal masses of
the original integral family. Let us assume we have boundary conditions for
the master integrals $\vec{I}^{(\kappa)}$ at a given point in phase-space, and
for some numeric values of the Feynman parameters $x_{1},\ldots,x_{\kappa}$.
In the remainder of the paper, we will typically choose $x_{j}=11/23$ for all
$j$. The precise choice is not important, as long as the external scales,
masses, and $x_{j}$ do not lie on a singularity of the differential equations.
Note that boundary conditions in step $\kappa=n-1$ can be obtained by
evaluating (11) in a chosen numerical point for the external scales, masses,
and Feynman parameters.

We may use DiffExp [13] to obtain a solution for the master integrals
$\vec{I}^{(\kappa)}$ in the interval $0<x_{\kappa}<1$, while keeping all other
variables at a fixed numerical value. The solution is obtained as a piece-wise
collection of generalized series expansions up to a given order, on a covering
set of line segments. We seek to use this representation to perform the
integrals in (13). In principle, this is straightforward since we may
integrate the series expansions term-by-term and segment by segment in an
efficient manner. We may also use the series expansion representation to
evaluate the limits in (2.2). However, there are some obstacles which we
discuss next.

The first obstacle is that the integrals which appear on the right-hand side
of (13) and (2.2) may not necessarily correspond to the choice of master
integrals $\vec{I}^{(\kappa)}$. Therefore, we should use IBP relations to
write all integrals in terms of the given master integrals. Since we have
representations for the masters in terms of collections of generalized series
expansions, we have to convert the IBP coefficients to the same representation
as well, which can be done in an automatic way.

A second obstacle comes from the (dimensional) regularization of the
integrals. In particular, integrals of the form (13) do not converge in
general when $\varepsilon$ is close to zero. Instead, it may be necessary to
obtain the result by analytically continuing from larger values of
$\varepsilon$. Practically, this can be achieved through a formula for
resolving the poles (18). Similarly, the limits in (2.2) have to be taken in a
way that is consistent with dimensional regularization. We explain in more
detail how to deal with these issues in section 3.2.

The last difficulty is the need for $i\delta$-prescriptions to perform the
transport past (threshold) singularities in the bulk of the interval
$0<x_{\kappa}<1$. If we choose our kinematics in the Euclidean region, there
are no such singularities by definition. If we choose kinematics in physical
regions, there are singularities and it is necessary to make choices which are
in correspondence with the Feynman prescription. We will discuss this point
further in section 3.3.

Taking care of the previously mentioned subtleties, we may use (13) and (2.2)
to obtain new boundary conditions for the master integrals
$\vec{I}^{(\kappa-1)}$. By iterating this procedure, we eventually obtain
numerical results for the master integrals of the original integral family
$\vec{I}=\vec{I}^{(0)}$.

###  3.2 Regularization

Let us consider the regularization of (13). In general, Feynman parameter
integrals may not converge when $\varepsilon$ is close to zero, due to the
presence of non-integrable singularities at the integration boundaries
$x_{\kappa}=0$ and $x_{\kappa}=1$. There may also be singularities in the bulk
of the integration region $0<x_{\kappa}<1$, which can be integrated over by
deforming the contour using suitable $i\delta$-prescriptions. We will discuss
this case in section 3.3.

Singularities at the boundary can be dealt with using sector decomposition.
What makes our computation different from normal is that our integrand is
represented by a collection of series expansions along a partitioning of the
interval $x_{\kappa}\in[0,1]$. Let us first distinguish two ways to deal with
the dimensional regulator $\varepsilon$ in the series expansion approach. The
first option is to work order-by-order in $\varepsilon$. DiffExp was
originally written with this approach in mind. However, we found it difficult
to apply regularization formulas while the series expansions are in this
representation. Instead, we may find a (quasi-)finite basis of master
integrals for the differential equations [24], which can be done using for
example the functionality in Reduze 2 [33]. Non-integrable singularities at
the boundary disappear when choosing (quasi-)finite master integrals for each
family $I_{\vec{\nu}}^{(\kappa)}$. A downside to this approach is that finite
basis integrals may contain high degree numerators and dots, or may be shifted
in dimension, and therefore may require more complicated IBP reductions.

The second option for the transport with DiffExp is to use a fixed numeric
value for $\varepsilon$. By repeating the computation multiple times for
different values of $\varepsilon$, it is possible to reconstruct the
dependence on $\varepsilon$ up to a given order [18]. Although DiffExp has not
been designed with this approach in mind, the package also works on a
differential matrix with a numeric value plugged in for
$\varepsilon$.111However, because DiffExp uses Mathematica’s SeriesData
format, it is not able to work with very small values of $\varepsilon$. In
particular, the SeriesData representation stores many redundant zero
coefficients when a power series is multiplied by an overall fractional power
$x^{b\varepsilon}$. For the examples in section 4, we chose the values in
(27). The $\varepsilon$-sampling approach is furthermore suitable for analytic
regularization of the poles [38], which we discuss next. Let us consider an
integral of the form

$\displaystyle\int_{0}^{c}dx\,x^{a+b\varepsilon}g(x,\varepsilon)\,.$ (17)

Here we assume that $a$ and $b$ are rational numbers, that $b\neq 0$, and that
the upper bound is chosen such that there are no non-integrable singularities
in the interval $0<x\leq c$, except possibly at $x=0$. Furthermore, we assume
that $g(x,\varepsilon)$ is a Taylor series in $x$ with a non-zero finite part,
and whose coefficients may depend on $\varepsilon$. When $a\leq-1$ and
$\varepsilon\rightarrow 0$, there is a non-integrable singularity at $x=0$. We
may resolve it by (possibly repeatedly) applying the following integral
identity:

$\displaystyle\int_{0}^{c}dx\,x^{a+b\varepsilon}g(x)=\int_{0}^{c}dx\,\frac{x^{a+b\varepsilon+1}}{(1+a+b\varepsilon)}\left(\frac{(2+a+b\varepsilon)}{c}g(x)-\left(1-\frac{x}{c}\right)g^{\prime}(x)\right)\,,$
(18)

which is similar to the formula for extracting poles in the sector
decomposition approach (see e.g. [39]). It is clear that this formula
increases the overall scaling exponent in the limit $x\rightarrow 0$, such
that repeated application of the formula resolves non-integrable
singularities.

We discuss next how to apply this formula to integrate results from DiffExp.
The solutions from DiffExp are given as a collection of line segments. In
particular, there is a segment around $x=0$, and a segment around
$x^{\prime}=1-x=0$. For illustrative purposes, we will consider the segment
around $x=0$, but the discussion is the same for the segment at the upper
integration boundary with local line parameter $x^{\prime}$. We decompose the
expansions in the following way:

$\displaystyle
g(x)=g_{0}(x,\epsilon)+x^{a_{1}+b_{1}\varepsilon}g_{1}(x,\varepsilon)+\ldots+x^{a_{k}+b_{k}\varepsilon}g_{k}(x,\varepsilon)\,,$
(19)

where $g_{0}(x,\epsilon)$ is a Taylor series in $x$, and where each
contribution $x^{a_{i}+b_{i}\varepsilon}g_{i}(x,\varepsilon)$ with
$i=1,\ldots,k,$ leads to an integration of the form of (17). We may assume
that a numerical value for $\varepsilon$ is already plugged into each
$g_{i}(x,\varepsilon)$. The exponents $a_{i}$ and $b_{i}$ may be read off from
the indicial equations of the differential equations.

The fact that a decomposition of the form of (19) is possible, and that there
are no contributions with $a_{i}\leq-1$ and $b_{i}=0$, was established
empirically by computing a number of examples. Furthermore, we might expect
logarithms in the expansions, but we did not observe these in our examples.
(However, logarithms of the type $\log(x)$ may be obtained by expanding in
$\epsilon$.) Lastly, we found in our examples that the exponents $a_{i}$ and
$b_{i}$ were always integers.

After performing the decomposition, we may apply (18) at integrand level to
each factor in (19) until we reach a form where every $a_{i}\geq 0$.
Afterwards, we can perform the integration as usual. Lastly, we remark on how
to perform the limits in (2.2) for integrals which do not depend on one of the
combined propagators. To evaluate the limits, we take the segment centered at
$x=0$ or at $x^{\prime}=1-x=0$, and we filter out the finite coefficient of
the Taylor series $g_{0}(x,\epsilon)$ in (19). This means we put any
contributions of the form $x^{a_{i}+b_{i}\varepsilon}g_{i}(x,\varepsilon)$
with $b_{i}\neq 0$ to zero (even when $a_{i}<0$).

###  3.3 Threshold singularities

When crossing threshold singularities (with DiffExp) it is necessary to
provide $i\delta$-prescriptions (contour deformations) which are in agreement
with the Feynman prescription. Since our method involves solving differential
equations of Feynman parameters, instead of external scales and masses, it is
not always clear how to proceed with a suitable choice of $i\delta$. One
option is to first compute a point in the Euclidean region (when a Euclidean
region exists for the given integral family.) One may then obtain results in
the physical region by transporting from the Euclidean region using the system
of differential equations of the (original) integral family.

Alternatively, we may seek to compute results directly in the physical region
with our method. Let us first recall that if a graph representation exists for
an integral family, we may identify possible thresholds by drawing unitarity
cuts. In practice, we can use a simple automatic strategy which was also
employed in [20]. We take the set of external legs $p_{j}$ and compute all
possible subsets (including the empty set.) Next, we do the same for the set
of all internal masses $m_{j}$. Lastly, we generate a list of all
$i\delta$-prescriptions of the form

$\displaystyle s-M^{2}+i\delta\,,$ (20)

where $s$ is a sum of momenta, and $M$ is a sum of internal masses. We then
provide DiffExp with the resulting set of $i\delta$-prescriptions. The list is
overcomplete (and includes pseudo-thresholds), but any redundant prescriptions
will be of no effect. For topologies which contain combined propagators, a
graphical representation often still exist and we can employ the above
strategy on the (generalized) external momenta and masses (cf. (10)). In this
case, $s$ and $M$ will depend on the unintegrated Feynman parameter(s).

If we choose to combine only propagators which have the same loop momenta in
the first steps of the iteration, we have a diagrammatic representation for
the integral families in these steps. For example, for the two-loop examples
in section 4 we eventually obtain a (generalized) sunrise topology (see Figure
7). In the final steps of the iteration, we have to combine propagators with
different loop momenta, and here we do not currently have a general
understanding of which $i\delta$-prescriptions to give to DiffExp. We note
that we do not always cross physical thresholds for every integration. For
example, in section 4.1 we compute a non-planar double pentagon family in a
kinematic point that lies outside the Euclidean region. However, we found that
we do not cross any physical thresholds during the transports of
$x_{4},x_{5},x_{6}$ and $x_{7}$. We elaborate on this point in section 4.1.

More generally, we may need to cross physical thresholds of integral families
that do not have a diagrammatic depiction. On the one hand, we could choose to
compute boundary conditions for the last diagrammatic topology in the
iteration using another method such as auxiliary mass flow, or by deriving
boundary conditions manually using expansion by regions. On the other hand, we
believe that a suitable contour prescription exist in principle for each
Feynman parameter. For example, we managed to obtain consistent results in the
physical region of the sunrise integral family with our method, after trying a
few different choices of $i\delta$-prescriptions. To get consistent results,
we also included a small $-i\delta$ in the second Symanzik polynomial
appearing in the boundary conditions of the generalized tadpole topology (cf.
(6), (11)).

We hope that a more rigorous analysis in the future will shed further light on
the general choice of $i\delta$-prescriptions for the Feynman parameters,
especially when computing integral families that do not have a diagrammatic
representation.

##  4 Pedagogical examples

In this section, we define two integral families which we will use for our
examples. The computation of the second integral family will be discussed in
detail in section 4.1. For the examples, we have opted to flip the sign of the
momentum and mass in the propagators compared to section 2. The first example
family is defined by

$\displaystyle
I^{\text{topo7}}_{\nu_{1},\nu_{2},\nu_{3},\nu_{4},\nu_{5},\nu_{6},\nu_{7},\nu_{8},\nu_{9}}=\int\frac{d^{d}k_{1}}{i\pi^{\frac{d}{2}}}\frac{d^{d}k_{2}}{i\pi^{\frac{d}{2}}}\,\frac{D_{8}^{-\nu_{8}}D_{9}^{-\nu_{9}}}{D_{1}^{\nu_{1}}D_{2}^{\nu_{2}}D_{3}^{\nu_{3}}D_{4}^{\nu_{4}}D_{5}^{\nu_{5}}D_{6}^{\nu_{6}}D_{7}^{\nu_{7}}}\,,$
(21)

where $d=4-2\varepsilon$, and where the propagators are given by

$\displaystyle D_{1}$ $\displaystyle=(k_{1}+p_{3})^{2}-m_{2}^{2}\,,$
$\displaystyle D_{4}$ $\displaystyle=k_{2}^{2}\,,$ $\displaystyle D_{7}$
$\displaystyle=(k_{1}-k_{2})^{2}\,,$ $\displaystyle D_{2}$
$\displaystyle=k_{1}^{2}\,,$ $\displaystyle D_{5}\,$
$\displaystyle=(k_{2}-p_{2})^{2}\,,$ $\displaystyle D_{8}$
$\displaystyle=(k_{1}-p_{2})^{2}\,,$ $\displaystyle D_{3}$
$\displaystyle=(k_{1}-p_{1}-p_{2})^{2}\,,$ $\displaystyle D_{6}$
$\displaystyle=(k_{2}-p_{1}-p_{2})^{2}\,,$ $\displaystyle D_{9}$
$\displaystyle=(k_{2}-p_{3})^{2}\,,$ (22)

where $D_{8}$ and $D_{9}$ are assumed to be numerators ($\nu_{8},\nu_{9}\leq
0$). The top sector is depicted in Figure 1. The kinematics is given by

$\displaystyle p_{1}^{2}=p_{2}^{2}=p_{3}^{2}=0\,,\quad\quad p_{1}\cdot
p_{2}=s/2\,,\quad\quad p_{1}\cdot p_{3}=(t-s+m_{1}^{2})/2\,,\quad\quad
p_{2}\cdot p_{3}=-t/2\,,$ (23)

where the momenta satisfy the conservation identity
$p_{1}+p_{2}+p_{3}+p_{4}=0$.

Figure 1:  A double box topology named topo7 which in this paper is dressed
with one massive external leg and one massive propagator.

Our second example is of higher complexity, and defined by

$\displaystyle
I^{\text{5p}}_{\nu_{1},\nu_{2},\nu_{3},\nu_{4},\nu_{5},\nu_{6},\nu_{7},\nu_{8},\nu_{9},\nu_{10},\nu_{11}}=\int\frac{d^{d}k_{1}}{i\pi^{\frac{d}{2}}}\frac{d^{d}k_{2}}{i\pi^{\frac{d}{2}}}\,\frac{D_{9}^{-\nu_{9}}D_{10}^{-\nu_{10}}D_{11}^{-\nu_{11}}}{D_{1}^{\nu_{1}}D_{2}^{\nu_{2}}D_{3}^{\nu_{3}}D_{4}^{\nu_{4}}D_{5}^{\nu_{5}}D_{6}^{\nu_{6}}D_{7}^{\nu_{7}}D_{8}^{\nu_{8}}}\,,$
(24)

where $d=4-2\varepsilon$, and where the propagators are

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

We assume that $D_{9}$, $D_{10}$ and $D_{11}$ are numerators
($\nu_{9},\nu_{10},\nu_{11}\leq 0$). The kinematics is given by:

$\displaystyle p_{1}^{2}=p_{2}^{2}=p_{3}^{2}=p_{4}^{2}=0\,,\quad\quad
p_{1}\cdot p_{2}=s_{12}/2\,,\quad\quad p_{1}\cdot p_{3}=s_{13}/2\,,\quad\quad
p_{1}\cdot p_{4}=s_{14}/2\,,$ $\displaystyle p_{2}\cdot
p_{3}=s_{23}/2\,,\quad\quad p_{2}\cdot
p_{4}=-(s_{12}+s_{13}+s_{14}+s_{23}+s_{34}-s_{55})/2\,,\quad\quad p_{3}\cdot
p_{4}=s_{34}/2\,,$ (26)

where the momenta are conserved according to
$p_{1}+p_{2}+p_{3}+p_{4}+p_{5}=0$. The family is depicted in Figure 2.

Figure 2:  A double pentagon topology named 5p which in this paper is dressed
with one massive external leg.

###  4.1 Non-planar double pentagon family

In the following section we describe the application of our method to the non-
planar double pentagon family $I^{\text{5p}}_{\vec{\nu}}$, with a focus on the
computational aspects of the calculation. We obtained our results using the
$\varepsilon$-sampling strategy, which we borrowed from [18]. More
specifically, we considered a set of small numerical values
$\\{\varepsilon_{i}\\}$ inserted for the dimensional regulator $\varepsilon$,
and we performed the computation every time for a fixed $\varepsilon_{i}$. At
the end, we reconstructed our results as a power series in $\varepsilon$ up to
order $\varepsilon^{4}$.

We aimed to reach a precision of 40 digits for the finite part of all master
integrals of the uncombined integral topology defined in (24). To do so, we
proceeded with some small trial and error. We found it sufficient to choose
numerical values

$\displaystyle\varepsilon_{i}=(-1)^{i}/(50+i)\quad\text{where
$i=1,\dots,40$\,,}$ (27)

and to perform the evaluation of all integrals at a precision of 60 digits.
Our samples $\varepsilon_{i}$ are relatively large in magnitude, because
DiffExp can currently not handle very small samples. Note that when
reconstructing the dependence on $\varepsilon$ up to a given order, there is a
precision loss involved. We reconstructed our results a few times by leaving
out a sample from the full set. By comparing the different results, we found a
variation that was within the desired precision of 40 digits for the finite
part. A more systematic analysis for how to choose the precision of the
evaluations, the values for $\varepsilon_{i}$, and the number of samples, is
described in [18].

Table 1: We iterate Feynman’s trick seven times to reduce the non-planar double pentagon topology to a trivial form. Each step $\kappa$ combines two propagators in the input column. The new combined propagator is given in the output column. The number of master integrals of the integral family with the combined propagator(s) is given in the last column. A large reduction of master integrals is seen for each iteration. $\kappa$ | input | output | Number of master integrals
---|---|---|---
0 | - | uncombined | 142
1 | $\\{D_{1},D_{2}\\}$ | $D_{12}=D_{1}x_{1}+(1-x_{1})D_{2}$ | 69
2 | $\\{D_{4},D_{5}\\}$ | $D_{45}=D_{4}x_{2}+(1-x_{2})D_{5}$ | 32
3 | $\\{D_{7},D_{8}\\}$ | $D_{78}=D_{7}x_{3}+(1-x_{3})D_{8}$ | 16
4 | $\\{D_{12},D_{3}\\}$ | $D_{123}=D_{12}x_{4}+(1-x_{4})D_{3}$ | 8
5 | $\\{D_{45},D_{6}\\}$ | $D_{456}=D_{45}x_{5}+(1-x_{5})D_{6}$ | 4
6 | $\\{D_{123},D_{456}\\}$ | $D_{123456}=D_{123}x_{6}+(1-x_{6})D_{456}$ | 2
7 | $\\{D_{123456},D_{78}\\}$ | $D_{12345678}=D_{123456}x_{7}+(1-x_{7})D_{78}$ | 1

Our computation required seven iterations, where we repeatedly combine two
propagators. We provide the order in which the propagators are combined in
Table 1. The ancillary file 5-point-FT-iteration.wl contains the definitions
of the topologies with combined propagators. There exists a diagrammatic
representation for the integral families in steps $\kappa=0,\dots,5$, which
are depicted in figures 2, 3, 4, 5, 6, 7. In step $\kappa=5$, we have a
generalized sunrise family, which is computed with two further iterations
where the topologies do not have a diagrammatic representation anymore. The
integral family $I_{\nu}^{\text{5p},(\kappa=7)}$ in step $\kappa=7$ can be
loosely viewed as a generalized tadpole with a single master integral, which
evaluates to a combination of (rescaled) Symanzik polynomials and gamma
functions (cf. (11)). We obtain the first boundary conditions by setting all
external scales and Feynman parameters to numerical values. In particular, we
consider the point $s_{14}=3$, $s_{13}=-11/17$, $s_{23}=-13/17$,
$s_{12}=-7/17$, $s_{34}=-7/13$, $s_{55}=-1$, and we set all Feynman parameters
to the value $x_{j}=11/23$. We then use DiffExp to transport in $x_{7}$, and
to obtain a generalized series representation in the interval $0<x_{7}<1$. By
integrating the series representation in accordance with Feynman’s trick (1),
we obtain boundary conditions for the integral family in step $\kappa=6$. We
repeat the procedure of transporting and integrating, until we reach the last
Feynman parameter integration in step $\kappa=1$.

After performing the last integration, we obtain numerical results for all
master integrals of the original family ($\kappa=0$) defined in (24). These
results may also be used as boundary conditions for the differential equations
of the original integral family, which allows one to reach other points in
phase space without going through the iteration again. We remark that we did
not cross any physical thresholds for the first transports in the Feynman
parameters $x_{7},\ldots,x_{4}$. This can be motivated by the fact that the
second Symanzik polynomial, obtained from the generalized tadpole in step
$\kappa=7$, is positive in our kinematic point with $x_{i}=11/23$ for
$i=1,2,3,$ and with $0<x_{j}<1$ for $j=4,\ldots,7$. For the transports of
$x_{3}$ and $x_{1}$ we found that we have to cross physical thresholds, while
for the transport in $x_{2}$ we did not observe any physical thresholds. We
chose $i\delta$-prescriptions for $x_{3}$ and $x_{1}$ in the manner described
in section 3.3. The result for one of the 142 most complicated master
integrals is given below at the numerical point $s_{14}=3$, $s_{13}=-11/17$,
$s_{23}=-13/17$, $s_{12}=-7/17$, $s_{34}=-7/13$, $s_{55}=-1$:

$\displaystyle I^{\text{5p}}_{13111111000}$
$\displaystyle=\frac{1}{\varepsilon^{-4}}\big{(}-80991.44634941832815855134956686330134244459\big{)}$
$\displaystyle+$
$\displaystyle\frac{1}{\varepsilon^{-3}}\big{(}-1176854.140501650857516200908950071824160111-$
$\displaystyle\quad\quad 303701.8453350029342400125918254935316349429i\big{)}$
$\displaystyle+$
$\displaystyle\frac{1}{\varepsilon^{-2}}\big{(}-13432835.8477692962185637394931604891797674-$
$\displaystyle\quad\quad 4251651.64965980166114774272201533676580580i\big{)}$
$\displaystyle+$
$\displaystyle\frac{1}{\varepsilon^{-1}}\big{(}-111346171.63704503288070435527859004232921-$
$\displaystyle\quad\quad 32927342.395688330300021665788556801968176i\big{)}$
$\displaystyle+$
$\displaystyle\big{(}-763045644.5561305442093867867513427731742-$
$\displaystyle\quad\quad 183231121.4048774146788661490531205282119i\big{)}$
$\displaystyle+$
$\displaystyle\varepsilon\big{(}-4428755434.16119754697555927652734791719-$
$\displaystyle\quad\quad 816059490.912195429388068459166197648719i\big{)}$
$\displaystyle+$
$\displaystyle\varepsilon^{2}\big{(}-23085640630.259889520777994526537639199-$
$\displaystyle\quad\quad 3082908606.7551294811504215473642629605i\big{)}$
$\displaystyle+$
$\displaystyle\varepsilon^{3}\big{(}-110164352209.7092412652451256610943938-$
$\displaystyle\quad\quad 10252510409.42185691550687766152353640i\big{)}$
$\displaystyle+$
$\displaystyle\varepsilon^{4}\big{(}-497649560130.015209279192098631531920-$
$\displaystyle\quad\quad 30796992268.3516086870566559550754104i\big{)}.$ (28)

We have numerically cross-checked the master integrals with AMFlow.

###  4.2 Computational complexity

The main advantage of our method is that the integral families with combined
propagators typically have significantly fewer master integrals for each step
in the iteration. We also remark that there is a significant overlap between
the master integrals of different steps in the iteration. This feature could
be utilized in a future specialized integration code. We motivate with four
different examples in Table 2 that the computational complexity of our method
is in many cases reduced compared to auxiliary mass flow [18]. In particular,
if we assume that the same algorithms are used to perform the IBP reductions
and the transportation of boundary conditions from the differential equations,
we expect that the computational complexity is mostly determined by the
maximal number of master integrals appearing in the computation of the Feynman
integrals. Our example family 5p serves as an ideal candidate to showcase a
significant reduction in the computational complexity compared to the
auxiliary mass flow method. We found in this case that the CPU hours for
performing the IBP-reductions are reduced by a factor of 66 between our method
and auxiliary mass flow.

Table 2: The maximal number of master integrals appearing in the computation of a few integral families. Integral family | No deformation | Combined propagators | AMFlow
---|---|---|---
topo7 | 31 | 19 | 31
topo7 with $m_{1}=0$, $m_{2}=0$ | 8 | 12 | 21
5p | 142 | 69 | 191
5p with $s_{55}=0$ | 108 | 69 | 174

There are a few possible caveats to the above analysis. Firstly, in our
current implementation it takes relatively long to compute the series
expansions with DiffExp in the intervals $0<x_{j}<1$. Furthermore, the
integration and regularization of the series expansions (cf. (18) and (13))
takes some time in our current implementation. For example, the final
transport with DiffExp and the subsequent integration took 340 CPU hours in
total for the 5p integral family. The whole bottom-up iteration for the 5p
integral family took 700 CPU hours (which includes the computation of all 40
numerical samples in the dimensional regulator.) We expect that incorporating
ideas from the differential solver used in AMFlow will improve the performance
of DiffExp significantly.

##  5 Conclusion

In this work we demonstrated a novel method for computing Feynman integrals.
The method is similar to the direct integration method, in which the Feynman
parametrization is integrated one Feynman parameter at a time. However, in our
approach each integration is done numerically by solving a system of
differential equations of a simplified Feynman integral family. The solutions
are obtained in terms of generalized series expansions, which are subsequently
integrated term-by-term. We furthermore showed that non-integrable
singularities at the integration boundaries can be dealt with using a
regularization formula. This leads to a computationally efficient and precise
method for performing Feynman parameter integrals numerically. Lastly, we have
made a first exploration of how to obtain results outside the Euclidean region
with our method. (Although a full understanding of the necessary contour
deformations of the Feynman parameters is left for the future.)

The approach discussed in this paper gives a fully algorithmic way to
numerically compute master integrals belonging to a generic Feynman integral
family. Furthermore, we expect that our method is computationally more
efficient than auxiliary mass flow, as it typically requires the computation
of integral families with fewer master integrals (see Table 2). This
significantly reduces the complexity of the IBP reductions which are required.
As a proof of concept, we have discussed in detail the computation of a non-
planar double pentagon integral family using our method. We leave the
publication of a general public code that implements our method to a future
publication.

## Acknowledgements

JU and MH thank Xiao Liu for fruitful discussions during the 2 weeks workshop
Precision calculations for future e+e– colliders: targets and tools at CERN.
JU and MH thank Oliver Schlotterer and Stefan Weinzierl for short comments on
the manuscript. JU thanks Janusz Gluza and his group for providing computing
resources with the support in part by the Polish National Science Center (NCN)
under grant 2017/25/B/ST2/01987. MH is supported by the European Research
Council under ERC-STG-804286 UNISCAMP.

## Appendix A: Feynman trick iterations of the double pentagon family

Figure 3:  Step $\kappa=1$: Combining propagators for the first time, reducing
the complexity of the initial topology 5p.  Figure 4: Step $\kappa=2$:
Combining propagators for the second time, reducing the complexity further.
Figure 5: Step $\kappa=3$: Combining propagators for the third time, reducing
the complexity further.  Figure 6: Step $\kappa=4$: Combining propagators for
the fourth time, reducing the complexity further.  Figure 7: Step $\kappa=5$:
Combining propagators for the fifth time, reducing the complexity further, and
giving a sunrise topology. Further iterations do not have a diagrammatic
representation.

## References

  * [1] A. V. Kotikov, “Differential equations method: New technique for massive Feynman diagrams calculation,” Phys. Lett. B 254 (1991) 158–164.
  * [2] J. M. Henn, “Multiloop integrals in dimensional regularization made simple,” Phys. Rev. Lett. 110 (2013) 251601, arXiv:1304.1806 [hep-th].
  * [3] A. B. Goncharov, “Multiple polylogarithms, cyclotomy and modular complexes,” Math. Res. Lett. 5 (1998) 497–516, arXiv:1105.2076 [math.AG].
  * [4] A. B. Goncharov, “Multiple polylogarithms and mixed Tate motives,” arXiv:math/0103059.
  * [5] J. Broedel, C. Duhr, F. Dulat, and L. Tancredi, “Elliptic polylogarithms and iterated integrals on elliptic curves II: an application to the sunrise integral,” Phys. Rev. D 97 no. 11, (2018) 116009, arXiv:1712.07095 [hep-ph].
  * [6] L. Adams and S. Weinzierl, “The $\varepsilon$-form of the differential equations for Feynman integrals in the elliptic case,” Phys. Lett. B 781 (2018) 270–278, arXiv:1802.05020 [hep-ph].
  * [7] S. Pozzorini and E. Remiddi, “Precise numerical evaluation of the two loop sunrise graph master integrals in the equal mass case,” Comput. Phys. Commun. 175 (2006) 381–387, arXiv:hep-ph/0505041.
  * [8] X. Liu, Y.-Q. Ma, and C.-Y. Wang, “A Systematic and Efficient Method to Compute Multi-loop Master Integrals,” Phys. Lett. B 779 (2018) 353–357, arXiv:1711.09572 [hep-ph].
  * [9] R. N. Lee, A. V. Smirnov, and V. A. Smirnov, “Solving differential equations for Feynman integrals by expansions near singular points,” JHEP 03 (2018) 008, arXiv:1709.07525 [hep-ph].
  * [10] M. K. Mandal and X. Zhao, “Evaluating multi-loop Feynman integrals numerically through differential equations,” JHEP 03 (2019) 190, arXiv:1812.03060 [hep-ph].
  * [11] M. L. Czakon and M. Niggetiedt, “Exact quark-mass dependence of the Higgs-gluon form factor at three loops in QCD,” JHEP 05 (2020) 149, arXiv:2001.03008 [hep-ph].
  * [12] F. Moriello, “Generalised power series expansions for the elliptic planar families of Higgs + jet production at two loops,” JHEP 01 (2020) 150, arXiv:1907.13234 [hep-ph].
  * [13] M. Hidding, “DiffExp, a Mathematica package for computing Feynman integrals in terms of one-dimensional series expansions,” Comput. Phys. Commun. 269 (2021) 108125, arXiv:2006.05510 [hep-ph].
  * [14] T. Armadillo, R. Bonciani, S. Devoto, N. Rana, and A. Vicini, “Evaluation of Feynman integrals with arbitrary complex masses via series expansions, ,” arXiv:2205.03345 [hep-ph].
  * [15] M. Beneke and V. A. Smirnov, “Asymptotic expansion of Feynman integrals near threshold,” Nucl. Phys. B 522 (1998) 321–344, arXiv:hep-ph/9711391.
  * [16] V. A. Smirnov, “Problems of the strategy of regions,” Phys. Lett. B 465 (1999) 226–234, arXiv:hep-ph/9907471.
  * [17] B. Jantzen, “Foundation and generalization of the expansion by regions,” JHEP 12 (2011) 076, arXiv:1111.2589 [hep-ph].
  * [18] X. Liu and Y.-Q. Ma, “AMFlow: a Mathematica package for Feynman integrals computation via Auxiliary Mass Flow,” arXiv:2201.11669 [hep-ph].
  * [19] X. Liu and Y.-Q. Ma, “Multiloop corrections for collider processes using auxiliary mass flow,” Phys. Rev. D 105 no. 5, (2022) L051503, arXiv:2107.01864 [hep-ph].
  * [20] I. Dubovyk, A. Freitas, J. Gluza, K. Grzanka, M. Hidding, and J. Usovitsch, “Evaluation of multi-loop multi-scale Feynman integrals for precision physics,” arXiv:2201.02576 [hep-ph].
  * [21] T. Binoth and G. Heinrich, “An automatized algorithm to compute infrared divergent multi-loop integrals,” Nucl. Phys. B585 (2000) 741–759, hep-ph/0004013v.2.
  * [22] S. Borowka, G. Heinrich, S. Jahn, S. P. Jones, M. Kerner, and J. Schlenk, “A GPU compatible quasi-Monte Carlo integrator interfaced to pySecDec,” Comput. Phys. Commun. 240 (2019) 120–137, arXiv:1811.11720 [physics.comp-ph].
  * [23] A. V. Smirnov, N. D. Shapurov, and L. I. Vysotsky, “FIESTA5: Numerical high-performance Feynman integral evaluation,” Comput. Phys. Commun. 277 (2022) 108386, arXiv:2110.11660 [hep-ph].
  * [24] A. von Manteuffel, E. Panzer, and R. M. Schabinger, “A quasi-finite basis for multi-loop Feynman integrals,” JHEP 02 (2015) 120, arXiv:1411.7392 [hep-ph].
  * [25] C. G. Papadopoulos and C. Wever, “Internal Reduction method for computing Feynman Integrals,” JHEP 02 (2020) 112, arXiv:1910.06275 [hep-ph].
  * [26] M. Hidding and F. Moriello, “All orders structure and efficient computation of linearly reducible elliptic Feynman integrals,” JHEP 01 (2019) 169, arXiv:1712.04441 [hep-ph].
  * [27] S. Weinzierl, “Feynman Integrals,” arXiv:2201.03593 [hep-th].
  * [28] K. Chetyrkin and F. Tkachov, “Integration by parts: The algorithm to calculate $\beta$ functions in 4 loops,” Nucl. Phys. B192 (1981) 159–204.
  * [29] S. Laporta, “High precision calculation of multiloop Feynman integrals by difference equations,” Int. J. Mod. Phys. A 15 (2000) 5087–5159, arXiv:hep-ph/0102033.
  * [30] O. V. Tarasov, “Connection between Feynman integrals having different values of the space-time dimension,” Phys. Rev. D 54 (1996) 6479–6490, arXiv:hep-th/9606018.
  * [31] R. N. Lee, “Space-time dimensionality D as complex variable: Calculating loop integrals using dimensional recurrence relation and analytical properties with respect to D,” Nucl. Phys. B 830 (2010) 474–492, arXiv:0911.0252 [hep-ph].
  * [32] R. N. Lee, “LiteRed 1.4: a powerful tool for reduction of multiloop integrals,” J. Phys. Conf. Ser. 523 (2014) 012059, arXiv:1310.1145 [hep-ph].
  * [33] A. von Manteuffel and C. Studerus, “Reduze 2 - Distributed Feynman Integral Reduction,” arXiv:1201.4330 [hep-ph].
  * [34] A. V. Smirnov and F. S. Chuharev, “FIRE6: Feynman Integral REduction with Modular Arithmetic,” Comput. Phys. Commun. 247 (2020) 106877, arXiv:1901.07808 [hep-ph].
  * [35] P. Maierhöfer, J. Usovitsch, and P. Uwer, “Kira - a Feynman integral reduction program,” Comput. Phys. Commun. 230 (2018) 99–112, arXiv:1705.05610 [hep-ph].
  * [36] C. Bogner and S. Weinzierl, “Feynman graph polynomials,” Int. J. Mod. Phys. A 25 (2010) 2585–2618, arXiv:1002.3458 [hep-ph].
  * [37] E. Remiddi, “Differential equations for Feynman graph amplitudes,” Nuovo Cim. A 110 (1997) 1435–1452, arXiv:hep-th/9711188.
  * [38] E. Panzer, “On hyperlogarithms and Feynman integrals with divergences and many scales,” JHEP 03 (2014) 071, arXiv:1401.4361 [hep-th].
  * [39] G. Heinrich, “Sector Decomposition,” Int. J. Mod. Phys. A 23 (2008) 1457–1486, arXiv:0803.4177 [hep-ph].


# STATUS OF QUBIC, THE Q&U BOLOMETRIC INTERFEROMETER FOR COSMOLOGY

L. MOUSSET on behalf of the QUBIC collaboration  IRAP, Université de Toulouse,
CNRS, CNES, UPS, Toulouse, France

###### Abstract

The Q&U Bolometric Interferometer for Cosmology (QUBIC) is a novel kind of
polarimeter optimized for the measurement of the $B$-mode polarization of the
Cosmic Microwave Background (CMB), which is one of the major challenges of
observational cosmology. The signal is expected to be of the order of a few
tens of nK, prone to instrumental systematic effects and polluted by various
astrophysical foregrounds which can only be controlled through multichroic
observations. QUBIC is designed to address these observational issues with a
novel approach that combines the advantages of interferometry in terms of
control of instrumental systematics with those of bolometric detectors in
terms of wide-band, background-limited sensitivity.

##  1 Introduction

The quest for $B$-mode polarization of the Cosmic Microwave Background (CMB)
is among the major challenges of observational cosmology. Cosmic inflation
predicts primordial scalar perturbations of the metric (density fluctuations),
but also tensor perturbations, equivalent to primordial gravitational waves.
These tensor modes should be imprinted in the CMB polarization fluctuations
with a very specific signature: odd-parity patterns (curl term), known as
$B$-modes [6]. The amplitude of the tensor modes, relative to the scalar
modes, is parametrized by the so called tensor-to-scalar ratio $r$. The
expected signal is weak, requiring high sensitivity detectors. In addition,
astrophysical foregrounds produce non-primordial $B$-mode polarization, such
as thermal emission from dust grains in the Galaxy.

The Q&U Bolometric Interferometer for Cosmology (QUBIC) was designed to
address the $B$-mode detection challenge [1]. Characterization and calibration
of a Technical Demonstrator (TD) started in 2018, at Astroparticle Physics &
Cosmology (APC) laboratory. In May 2021, the instrument has been sent to
Argentina. A second calibration phase is undergoing and the instrument will be
installed on the observation site in the next months. In Hamilton et al. [1],
we give forecasts for typical observations and measurements: with three years
of integration on the sky and assuming perfect foreground removal as well as
stable atmospheric conditions from our site in Argentina, our simulations show
that we can achieve a statistical sensitivity to the effective tensor-to-
scalar ratio (including primordial and foreground $B$-modes)
$\sigma(r)=0.015$.

##  2 Control of systematic effects with QUBIC

As a bolometric interferometer, QUBIC combines the advantages of
interferometry in terms of control of instrumental systematic effects with
those of bolometric detectors in terms of wide-band, background-limited
sensitivity. A picture of the instrument with a sketch of the optical design
is shown in Figure 1.

Figure 1: Left: Picture of the QUBIC instrument at APC during the integration
phase. Right: Optical sketch of the instrument.

###  2.1 Low cross-polarization

QUBIC bolometers are full power detectors. The measurement of the polarization
is only done by modulating the signal amplitude with a rotating Half-Wave
Plate (HWP) and a fixed polarizing grid. Every optical element has its own
systematic effects which could induce cross-polarization. This means that the
two polarization directions $(x,y)$ of the signal can be mixed when
interacting with the instrument and this is a major issue. The $x$-axis is
defined as the transmission axis of the polarizing grid. By putting the HWP
and the polarizing grid right after the window, the $x$ polarization is
selected as early as possible. In this way, any cross-polarization occurring
after the polarizing grid, for example generated by the horn-array, the
mirrors or any reflection in the instrument, has no impact.

The modulation of the polarization was tested during the calibration phase and
a very low cross-polarization was indeed detected. Those measurements are
presented in Torchinsky et al. [3] and D’Alessandro et al. [4] and the cross-
polarization contamination at 150 GHz is compatible with zero to within 0.6%.

###  2.2 Self-calibration technique

Interferometry offers the possibility to self-calibrate the instrument
systematic effects. This technique has been used for a long time in radio
astronomy [8]. Self-calibration relies on the concept of equivalent baselines,
one baseline $b$ being formed with two horns. The full horn-array is shown in
Figure 2 (left).

The self-calibration technique is based on the fact that, in case of an ideal
instrument without any systematic effect, equivalent baselines produce the
same interference pattern on the focal plane, for an observation at infinity,
in the Fraunhofer regime. Thus, by measuring the differences, one can
calibrate the systematics of the instrument. The demonstration of this
technique for QUBIC was done in Bigot-Sazy et al. [7].

Figure 2: Left: Picture of the TD horn-array ($8\times 8$ back-to-back horns),
taken from Cavaliere et al. [5]. Right: Superimposition of three monochromatic
beams measured in the laboratory at 130 GHz (red), 150 GHz (green) and 170 GHz
(blue), taken from Torchinsky et al. [3].

##  3 Spectral imaging capability

The instrument beam pattern, shown in Figure 2 (right), is given by the
geometric distribution of the horn-array. It contains multiple peaks whose
angular separation is linearly dependent on the wavelength. As a result, and
after a non-trivial map-making process, a bolometric interferometer such as
QUBIC can simultaneously produce sky maps at multiple frequency sub-bands with
data acquired over a single wide frequency band. The demonstration of this
technique, called spectral imaging, and the characterization of its
performance are presented in Mousset et al. [2].

###  3.1 Tested on end-to-end simulations

We demonstrate spectral imaging capabilities by trying to recover the
frequency dependence of the thermal galactic dust emission with simulated
observations. We simulate an observation in a sky patch of 15 degree radius.
The parameters of the pipeline are set in such a way that the simulated
instrument has a single focal plane operating either at 150 GHz or at 220 GHz
with a 25% bandwidth each. From these wide-band TOD, we are able to
reconstruct several numbers of sub-bands using spectral imaging.

The reconstructed intensity as a function of frequency is studied in a given
pixel. Figure 3 shows the intensity of the input sky convolved with the
instrument beam, and the reconstructed intensity for a given pixel,
considering 5 sub-bands in each wide band at 150 (red) and 220 (blue) GHz.

Figure 3: Intensity as a function of the frequency for $N_{\rm rec}=5$ sub-
bands in each wide band at 150 (red) and 220 (blue) GHz for a given pixel. The
grey regions correspond to the unobserved frequencies outside our physical
bands. Two sky pixels are shown as red stars, one in a patch centered at the
Galactic center and one in the patch that QUBIC plans to observe centered in
[0, -57 deg]. Red and blue dots: Input sky convolved with the instrument beam.
In both cases are shown in light color the 68% CL regions for a modified
black-body spectrum reconstructed with a MCMC from our simulated measurements
and sub-band covariance matrices. Maps are in $\mu$K CMB and $N_{\rm
side}=32$.

###  3.2 First trial on real data

Spectral imaging has been applied on real data for the first time during the
calibration campaign at the APC laboratory. The QUBIC instrument was placed on
an alt-azimuth mount in order to scan a calibration source tuned at 150 GHz
(with 144 Hz bandwidth) and placed in the far field. The corresponding
analysis is presented in Torchinsky et al. [3]. We performed a scan in azimuth
and elevation with the instrument, obtaining a TOD for each bolometer. We then
applied our spectral imaging map-making algorithm with five sub-bands to a
selection of 26 bolometers that do not exhibit saturation. The synthesized
beam for each bolometer is realistically modeled in our map-making through a
series of Gaussian whose amplitude, width and locations are fit from a
measured map of the synthesized beam for each bolometer (see Figure 20 from
Torchinsky et al. [3] for an example). We were able to reconstruct a map of
the point-like artificial calibration source as well as its location in
frequency space. In Figure 4, we show the reconstruction onto 5 sub-bands. The
expected point-source shape is clearly visible in the central frequency sub-
band containing the emission frequency of the source at 150 GHz, it is fainter
in adjacent bands, and not visible in the furthest bands. On the right, we
show the detected amplitude in the central pixel as a function of the
frequency. The measurement in red is compared to the expected value spectrum
in blue.

Figure 4: Left: Calibration data with the source at 150 GHz projected on the
sky using our map-making software to deconvolve from the multiple peaked
synthesized beam and split the physical band of the instrument into 5 sub-
bands. Right: Measurement of the flux of the source in reconstructed sub-
bands. The measurement (simple aperture photometry) in red is compared to the
expected value spectrum in blue.

##  4 Conclusion and perspectives

The QUBIC Technical Demonstrator was successfully tested at APC laboratory and
sent to Argentina to be installed on the observation site in the following
year. The QUBIC instrument relies on an innovative design which allows a very
low cross-polarization, a high control of systematic effects thanks to self-
calibration and a spectroscopic capability. Spectral imaging technique was
demonstrated on end-to-end simulations and applied to calibration data
successfully. Testing it on real astrophysical sources is the next determining
step.

## References

## References

  * [1] J.-Ch. Hamilton et al. (QUBIC Collaboration), arXiv:2011.02213, to appear on JCAP.
  * [2] L. Mousset et al (QUBIC Collaboration), arXiv:2010.15119, to appear on JCAP.
  * [3] S.A. Torchinsky et al (QUBIC Collaboration), arXiv:2008.10056, submitted to JCAP.
  * [4] G. D’Alessandro et al (QUBIC Collaboration), arXiv:2008.10667, submitted to JCAP.
  * [5] F. Cavaliere et al (QUBIC Collaboration), arXiv:2008.12721, submitted to JCAP.
  * [6] M. Zaldarriaga and U. Seljak, Phys. Rev. D 55, 1830 (1997)
  * [7] M.-A. Bigot-Sazy et al, A&A, vol. 550, A59, Feb. 2013, arXiv:1209.4905
  * [8] T. Cornwell and P. Wilkinson, Mon. Not. Roy. Astron. Soc., vol. 196, 1067, Sep. 1981.