11institutetext: Higgs Centre for Theoretical Physics, University of
Edinburgh,
James Clerk Maxwell Building, Peter Guthrie Tait Road, Edinburgh, EH9 3FD,
United Kingdom

# Feynman Integrals from Positivity Constraints

Mao Zeng<EMAIL_ADDRESS>
###### Abstract

We explore inequality constraints as a new tool for numerically evaluating
Feynman integrals. A convergent Feynman integral is non-negative if the
integrand is non-negative in either loop momentum space or Feynman parameter
space. Applying various identities, all such integrals can be reduced to
linear sums of a small set of master integrals, leading to infinitely many
linear constraints on the values of the master integrals. The constraints can
be solved as a semidefinite programming problem in mathematical optimization,
producing rigorous two-sided bounds for the integrals which are observed to
converge rapidly as more constraints are included, enabling high-precision
determination of the integrals. Positivity constraints can also be formulated
for the $\epsilon$ expansion terms in dimensional regularization and reveal
hidden consistency relations between terms at different orders in $\epsilon$.
We introduce the main methods using one-loop bubble integrals, then present a
nontrivial example of three-loop banana integrals with unequal masses, where
11 top-level master integrals are evaluated to high precision.

###### Contents

  1. 1 Introduction
  2. 2 One-loop bubble integrals
    1. 2.1 Review: integration-by-parts (IBP) and dimensional-shifting identities
    2. 2.2 Positivity constraints in loop momentum space
      1. 2.2.1 A crude first attempt
      2. 2.2.2 Formulation in terms of matrix eigenvalues
      3. 2.2.3 High-precision evaluation using semidefinite programming
    3. 2.3 Positivity constraints in Feynman parameter space
    4. 2.4 Constraints for expansions in dimensional regularization parameter
      1. 2.4.1 Generic constraints
      2. 2.4.2 Taylored constraints for specific Feynman integrals
      3. 2.4.3 Bubble integral up to second order in $\epsilon$
    5. 2.5 $\epsilon$ expansion from numerical differentiation w.r.t. spacetime dimension
  3. 3 Three-loop banana integrals with unequal masses
    1. 3.1 Definitions and conventions
    2. 3.2 Positivity constraints
    3. 3.3 Numerical results
  4. 4 Discussions
  5. A Analytic results for bubble integrals

##  1 Introduction

Evaluation of Feynman integrals is both a classic problem and an expanding
frontier of research in quantum field theory, with many areas of applications
such as collider physics, statistical mechanics, and gravitational physics. A
widely used method for numerical evaluations of Feynman integrals is Monte
Carlo integration with sector decomposition Binoth:2000ps ; Bogner:2007cr ;
Kaneko:2009qx ; Borowka:2015mxa ; Borowka:2017idc ; Smirnov:2013eza ;
Smirnov:2015mct . Another method, numerical series solutions of differential
equations Moriello:2019yhu ; Hidding:2020ytt ; Liu:2017jxz ; Liu:2021wks ;
Liu:2022chg ; Armadillo:2022ugh ; Hidding:2022ycg , enjoyed intense
developments in the last few years. Some other methods, which are either
inherently numerical or can be used numerically, include Mellin-Barnes
representations Usyukina:1992jd ; Usyukina:1993ch ; Smirnov:1999gc ;
Tausk:1999vh ; Czakon:2005rk ; Smirnov:2009up ; Gluza:2007rt ;
Belitsky:2022gba , difference equations from dimensional recursion
Laporta:2000dsw ; Lee:2012te ; Lee:2022art , Loop-tree duality Catani:2008xa ;
Runkel:2019yrs ; Capatti:2019ypt , and tropical Monte Carlo integration
Borinsky:2020rqs ; Borinsky:2023jdv . See also books Weinzierl:2022eaz ;
Smirnov:2012gma  which comprehensively review both numerical and analytic
methods. Despite these developments, numerical evaluation of Feynman integrals
still presents challenges in the ongoing quest for higher precisions in
perturbative calculations, and new explorations are warranted.

A fruitful recent idea in theoretical physics is the use of positivity
constraints, e.g. arising from the unitarity of a Hilbert space, to constrain
unknown parameters from first principles, sometimes reaching great accuracy.
Prominent examples include the conformal bootstrap Rattazzi:2008pe , the non-
perturbative S-matrix bootstrap (see review Kruczenski:2022lot  and references
within), and EFT positivity bounds Adams:2006sv ; Bellazzini:2020cot ; Caron-
Huot:2020cmc ; Tolley:2020gtv ; Arkani-Hamed:2020blm ; Caron-Huot:2021rmr ;
Caron-Huot:2022ugt . Some of the predictions in the S-matrix and EFT contexts
have been checked against explicit perturbative calculations involving Feynman
integrals Bern:2021ppb ; Chen:2022nym ; Bellazzini:2022wzv , so it is natural
to explore positivity properties for Feynman integrals themselves. What
directly inspired this paper, though, is recent applications of positivity
constraints to bootstrapping simple quantum mechanical systems Lin:2020mme ;
Han:2020bkb ; Berenstein:2021dyf ; Berenstein:2022unr  and lattice models
Anderson:2016rcw ; Kazakov:2022xuh ; Cho:2022lcj , which crucially use various
linear identities that are reminiscent of integration-by-parts identities for
Feynman integrals Chetyrkin:1981qh  as well as a numerical technique known as
semidefinite programming vandenberghe1996semidefinite  which will be adapted
to our calculations. Semidefinite programming was introduced to theoretical
physics by Ref. Poland:2011ey  in the conformal bootstrap and subsequently
applied to wider contexts.

In this work, we will formulate positivity constraints for _Euclidean Feynman
integrals_ , which can be either (1) integrals in Euclidean spacetime or can
be rewritten in Euclidean spacetime after a trivial Wick rotation, or (2)
integrals in Minkowskian spacetime but with kinematics in the so called
_Euclidean region_ , i.e. with center-of-mass energy of incoming momenta below
the particle production threshold. In case (1), the loop integrand in
Euclidean momentum space has non-negative propagator denominators, and the
integrand is non-negative as long as the numerator is non-negative. In case
(2), the integral is real due to the absence of Cutkosky cuts. After Feynman
parametrization, the parametric integrand involves non-negative graph
polynomials and stays non-negative when multiplied by a positive function of
the Feynman parameters. In both cases, a non-negative integrand implies a non-
negative integral as long as the integral is convergent, i.e. has no
ultraviolet or infrared divergences. The initial restriction to convergent
integrals is not a fundamental limitation, as divergent integrals can be
rewritten as linear sums of convergent integrals multiplied by divergent
coefficients vonManteuffel:2014qoa , and then it suffices to evaluate the
convergent integrals as a Taylor series in the dimensional regularization
parameter $\epsilon$.

We will see that positivity constraints, expressed in the language of
semidefinite programming, are strong enough to precisely determine the values
of Feynman integrals. Moreover, rigorous error bounds can be obtained. Our
machinery relies on linear relations between Feynman integrals, including
integration-by-parts identities Chetyrkin:1981qh  and dimension-shifting
identities Tarasov:1996br , to express all positivity constraints as linear
constraints on the values of a set of master integrals.

The outline of the paper is as follows. In Section 2, we introduce the main
methods using the simple example of massive one-loop bubble integrals.
Specifically, Subsection 2.1 gives a short review of two kinds of linear
identities for Feynman integrals, integration-by-parts identities and
dimension shifting identities. Subsection 2.2 discusses positivity constraints
in Euclidean momentum space, starting from ad hoc constraints that narrow down
the allowed value of the bubble master integral to $\sim 50\%$ accuracy at an
example kinematic point, then developing the machinery of semidefinite
programming to reach an accuracy of $10^{-14}$. We switch to Feynman parameter
space in Subsection 2.3, which allows calculating the integrals in a wider
region of kinematic parameters below the two-particle cut threshold, while
many steps of the calculations are unchanged from the momentum-space case.
Subsection 2.4 formulates positivity constraints for $\epsilon$ expansions of
Feynman integrals in dimensional regularization. Subsection 2.5 presents an
alternative method for calculating the $\epsilon$ expansion based on numerical
differentiation of results w.r.t. the spacetime dimension. Having laid out
most of the methods, in Section 3, we present a nontrivial application to
three-loop banana integrals with four unequal masses. Due to a large number of
undetermined master integrals, semidefinite programming becomes essential in
efficiently solving the positivity constraints. We calculate all master
integrals at an example kinematic point up to the second order in $\epsilon$
expansion in $d=2-2\epsilon$, and a detailed account of the numerical
accuracies is given. We end with some discussions in Section 4.

##  2 One-loop bubble integrals

We consider the following family of Feynman integrals in Minkowski spacetime
parametrized by two integers $a_{1}$ and $a_{2}$,

$I_{a_{1},a_{2}}^{d}\equiv\int\frac{d^{d}l\,e^{\gamma_{E}\epsilon}}{i\pi^{d/2}}\frac{1}{(-l^{2}+m^{2})^{a_{1}}[-(p+l)^{2}+m^{2}]^{a_{2}}},$
(1)

which correspond to the diagram in Fig. but with the propagator denominators
raised to general integer powers.

Figure 1: The one-loop bubble integral with external legs of virtuality
$p^{2}$ and two internal massive line with the same squared mass $m^{2}$.

The external mass is $\sqrt{p^{2}}$ and the internal mass for the two
propagators is $m$, with the two propagators raised to powers $a_{1}$ and
$a_{2}$. If either $a_{1}$ or $a_{2}$ is non-positive, Eq. (1) becomes a
massive tadpole integral possibly with a numerator. If $a_{1}$ and $a_{2}$ are
both non-positive, the integral vanishes in dimensional regularization. The
spacetime dimension $d$ is equal to $4-2\epsilon$, with $\epsilon$ being the
dimensional regularization parameter. Eq. (1) can be re-written, by Wick
rotation of the integration contour, as the following integrals in _Euclidean_
spacetime,

$I_{a_{1},a_{2}}^{d}\equiv\int\frac{d^{d}\bm{l}\,e^{\gamma_{E}\epsilon}}{\pi^{d/2}}\frac{1}{(\bm{l}^{2}+m^{2})^{a_{1}}[(\bm{p}+\bm{l})^{2}+m^{2}]^{a_{2}}},\,.$
(2)

where $\bm{p}^{2}\equiv-p^{2}$.

We will use positivity constraints to numerically evaluate bubble integrals
Eq. (1) for $p^{2}<4m^{2}$, i.e. below the two-particle production threshold.
We first consider the case $p^{2}<0$, i.e. $\bm{p}^{2}>0$. In this case,
$\bm{p}$ can be literally embedded in Euclidean spacetime as a vector with
real-valued components, and we will derive positivity constraints starting
from the loop momentum integral Eq. (2). We will subsequently present a
treatment applicable to all $p^{2}<4m^{2}$ by using Feynman parameter
representations of the integrals. Before actual calculations, we first briefly
review linear identities for the Feynman integrals involved, arising from
_integration by parts_ and _dimension shifting_. Efficiently solving IBP
identities for more complicated Feynman integrals is a major research problem,
and we refer readers to Refs. Laporta:2000dsw ; Laporta:1996mq ;
Anastasiou:2004vj ; vonManteuffel:2012np ; Smirnov:2008iw ; Smirnov:2019qkx ;
Lee:2012cn ; Lee:2013mka ; Maierhofer:2017gsa  for relevant computational
algorithms and software.

###  2.1 Review: integration-by-parts (IBP) and dimensional-shifting
identities

$I^{d}_{a_{1},a_{2}}$ at different values of $a_{1}$ and $a_{2}$, as defined
in Eq. (1), are related through integration-by-parts (IBP) identities
Chetyrkin:1981qh , as total derivatives integrate to zero without boundary
terms in dimensional regularization. To derive the identities, we will write
dot products as linear combinations of denominators and constants,

$l^{2}=-(-l^{2}+m^{2})+m^{2},\qquad 2p\cdot
l=-[-(p+l)^{2}+m^{2}]+[-l^{2}+m^{2}]-p^{2}\,.$ (3)

The actual identities are

$\displaystyle\int\frac{d^{d}l\,e^{\gamma_{E}\epsilon}}{i\pi^{d/2}}\frac{\partial}{\partial
l^{\mu}}\frac{p^{\mu}}{(-l^{2}+m^{2})^{a_{1}}[-(p+l)^{2}+m^{2}]^{a_{2}}}=0$
$\displaystyle\implies-
a_{1}p^{2}I^{d}_{a_{1}+1,a_{2}}+a_{2}p^{2}I^{d}_{a_{1},a_{2}+1}-a_{1}I^{d}_{a_{1}+1,a_{2}-1}+(a_{1}-a_{2})I^{d}_{a_{1},a_{2}}+a_{2}I^{d}_{a_{1}-1,a_{2}+1}=0\,,$
(4)

and

$\displaystyle\int\frac{d^{d}l\,e^{\gamma_{E}\epsilon}}{\pi^{d/2}}\frac{\partial}{\partial
l^{\mu}}\frac{l^{\mu}}{(-l^{2}+m^{2})^{a_{1}}[-(p+l)^{2}+m^{2}]^{a_{2}}}=0\implies$
$\displaystyle
2a_{1}m^{2}I^{d}_{a_{1}+1,a_{2}}+a_{2}(-p^{2}+2m^{2})I^{d}_{a_{1},a_{2}+1}+(d-2a_{1}-a_{2})I^{d}_{a_{1},a_{2}}-a_{2}I^{d}_{a_{1}-1,a_{2}+1}=0\,,$
(5)

We also need the diagram reflection symmetry relations from relabeling
$l\rightarrow-(p+l)$,

$I^{d}_{a_{1},a_{2}}=I^{d}_{a_{2},a_{1}}\,,$ (6)

and the boundary condition in dimensional regularization,

$I^{d}_{a_{1},a_{2}}=0,\ \text{when both }a_{1}\leq 0\ \text{and }a_{2}\leq
0\,.$ (7)

Solving eqs. (4)-(7), e.g. by iteratively eliminating $I^{d}_{a_{1},a_{2}}$
with the largest values of $|a_{1}|+|a_{2}|$, allow us to rewrite any
$I^{d}_{a_{1},a_{2}}$ with fixed spacetime dimension $d$ in terms of two
_master integrals_ ,

$I^{d}_{1,1},\quad I^{d}_{1,0},$ (8)

i.e. a bubble integral and a tadpole integral. We will then change to an
alternative basis $I_{1}$ and $I_{2}$, which is ultraviolet finite as $d$
approaches $4$ and in fact for any $d<6$,

$\displaystyle I^{d}_{2,1}$
$\displaystyle=\frac{d-3}{p^{2}-4m^{2}}I^{d}_{1,1}+\frac{d-2}{2m^{2}(p^{2}-4m^{2})}I^{d}_{1,0}\,,$
(9) $\displaystyle I^{d}_{3,0}$
$\displaystyle=\frac{(d-2)(d-4)}{8m^{4}}I^{d}_{1,0}\,.$

_For the rest of this paper, we will impose $d<6$ in the treatment of the
bubble integrals and pay special attention to values of $d$ near $4$._

The process of carrying out the above calculation and rewriting any integral
as linear sums of master integrals is referred to as _IBP reduction_. Here are
two examples showing how other integrals are expressed as linear combinations
of the basis Eq. (9),

$\displaystyle I^{d}_{2,2}$
$\displaystyle=\frac{1}{p^{2}(4m^{2}-p^{2})}\Big{(}[(6-d)p^{2}-4m^{2}]I^{d}_{2,1}+4m^{2}I^{d}_{3,0}\Big{)}$
(10) $\displaystyle I^{d}_{3,1}$
$\displaystyle=I^{d}_{1,3}=\frac{1}{2p^{2}(4m^{2}-p^{2})}\Big{(}[(4-d)p^{2}+4m^{2}]I^{d}_{2,1}+2(p^{2}-2m^{2})I^{d}_{3,0}\Big{)}\,.$
(11)

The tadpole integral $I^{d}_{3,0}$ is well known in textbooks. The bubble
integral with one of the propagator raised to a higher power, $I^{d}_{2,1}$,
is also known analytically, but we will use positivity constraints to
numerically evaluate it for the purpose of demonstrating our methods. We will
write

$\displaystyle I^{d}_{2,1}=\hat{I}^{d}_{2,1}\cdot I^{d}_{3,0}\,,$ (12)

and treat the normalized bubble integral

$\hat{I}^{d}_{2,1}=I^{d}_{2,1}/I^{d}_{3,0}$ (13)

as an unknown quantity to be bounded by positivity constraints. More
generally, we define a “hatted” notation for integrals normalized against the
finite tadpole integral,

$\hat{I}^{d}_{a_{1},a_{2}}=I^{d}_{a_{1},a_{2}}/I^{d}_{3,0},\,$ (14)

under which Eqs. (10) and (11) become

$\displaystyle\hat{I}^{d}_{2,2}$
$\displaystyle=\frac{1}{p^{2}(4m^{2}-p^{2})}\Big{(}[(6-d)p^{2}-4m^{2}]\hat{I}^{d}_{2,1}+4m^{2}\Big{)}$
(15) $\displaystyle\hat{I}^{d}_{3,1}$
$\displaystyle=I^{d}_{1,3}=\frac{1}{2p^{2}(4m^{2}-p^{2})}\Big{(}[(4-d)p^{2}+4m^{2}]\hat{I}^{d}_{2,1}+2(p^{2}-2m^{2})\Big{)}\,.$
(16)

We also need dimensional-shifting identities. Using the Schwinger
parametrization, Eq. (1) is rewritten as

$I^{d}_{a_{1},a_{2}}=\frac{e^{\gamma_{E}\epsilon}}{\Gamma(a_{1})\Gamma(a_{2})}\int_{0}^{\infty}dx_{1}\int_{0}^{\infty}dx_{2}\,x_{1}^{a_{1}-1}x_{2}^{a_{2}-1}(x_{1}+x_{2})^{-d/2}\exp(-i\mathcal{U}/\mathcal{F})\,,$
(17)

where $\mathcal{U}$ and $\mathcal{F}$ are graph polynomials depending on
$x_{1}$ and $x_{2}$,

$\mathcal{U}(x_{1},x_{2})=x_{1}+x_{2},\quad\mathcal{F}(x_{1},x_{2})=m^{2}(x_{1}+x_{2})^{2}-p^{2}x_{1}x_{2}-i0^{+}\,.$
(18)

Therefore, $(d-2)$ dimensional integrals can be written as

$\displaystyle I^{d-2}_{a_{1},a_{2}}$
$\displaystyle=\frac{e^{\gamma_{E}\epsilon}}{\Gamma(a_{1})\Gamma(a_{2})}\int_{0}^{\infty}dx_{1}\int_{0}^{\infty}dx_{2}\,x_{1}^{a_{1}-1}x_{2}^{a_{2}-1}(x_{1}+x_{2})^{-(d-2)/2}\exp(-iU/F)$
$\displaystyle=\frac{e^{\gamma_{E}\epsilon}}{\Gamma(a_{1})\Gamma(a_{2})}\int_{0}^{\infty}dx_{1}\int_{0}^{\infty}dx_{2}\,x_{1}^{a_{1}-1}x_{2}^{a_{2}-1}(x_{1}+x_{2})(x_{1}+x_{2})^{-d/2}\exp(-iU/F)$
$\displaystyle=\frac{e^{\gamma_{E}\epsilon}}{\Gamma(a_{1})\Gamma(a_{2})}\int_{0}^{\infty}dx_{1}\int_{0}^{\infty}dx_{2}\,\big{(}x_{1}^{a_{1}}x_{2}^{a_{2}-1}+x_{1}^{a_{1}-1}x_{2}^{a_{2}}\big{)}(x_{1}+x_{2})^{-d/2}\exp(-iU/F)$
$\displaystyle=a_{1}I^{d}_{a_{1}+1,a_{2}}+a_{2}I^{d}_{a_{1},a_{2}+1}\,,$ (19)

where the last line used Eq. (18) to map different monomials in $x_{1}$ and
$x_{2}$ to bubble integrals with different propagator powers. Applying Eq.
(19) to the finite bubble master integral $I^{d}_{2,1}$ on the LHS of Eq. (9)
in $(d-2)$ spacetime dimensions and then performing IBP reduction, we obtain

$I^{d-2}_{2,1}=\frac{2(d-5)}{p^{2}-4m^{2}}I^{d}_{2,1}-\frac{2}{p^{2}-4m^{2}}I^{d}_{3,0}\,.$
(20)

The dimension-shifting formula for the tadpole integral can be obtained
easily, e.g. by using closed-form results for tadpole integrals. The result is

$I^{d-2}_{3,0}=\frac{6-d}{2m^{2}}I^{d}_{3,0}\,.$ (21)

Eqs. (20) and (21) are the dimension-shifting identities that express the two
master integrals in Eq. (9) in lower spacetime dimensions to the same master
integrals in higher spacetime dimensions. By inverting Eqs. (20) and (21), we
obtain dimension shifting identities in the reverse direction, i.e. expressing
master integrals in higher spacetime dimensions in terms of master integrals
in lower spacetime dimensions,

$\displaystyle I^{d+2}_{2,1}$
$\displaystyle=\frac{p^{2}-4m^{2}}{2(d-3)}I^{d}_{2,1}-\frac{2m^{2}}{(d-3)(d-4)}I^{d}_{3,0}\,,$
(22) $\displaystyle I^{d+2}_{3,0}$
$\displaystyle=\frac{2m^{2}}{4-d}I^{d}_{3,0}\,.$ (23)

###  2.2 Positivity constraints in loop momentum space

Here we focus on the case $p^{2}<0$ and use the Wick-rotated expression for
the integrals, Eq. (2), with

$\bm{p}^{2}=M^{2}=-p^{2}\,.$ (24)

We will later present a Feynman parameter-space treatment that seems more
powerful and in particular works for all $p^{2}<4m^{2}$, but the loop momentum
space treatment here will help build up intuitions.

####  2.2.1 A crude first attempt

We first consider the following two convergent integrals with non-negative
integrands in loop momentum space,

$\displaystyle I^{d}_{2,2}$
$\displaystyle=\int\frac{d^{d}\bm{l}\,e^{\gamma_{E}\epsilon}}{\pi^{d/2}}\frac{1}{(\bm{l}^{2}+m^{2})^{2}[(\bm{p}+\bm{l})^{2}+m^{2}]^{2}}\,,$
(25) $\displaystyle I^{d}_{3,1}+I^{d}_{1,3}-2I^{d}_{2,2}$
$\displaystyle=\int\frac{d^{d}\bm{l}\,e^{\gamma_{E}\epsilon}}{\pi^{d/2}}\frac{1}{(\bm{l}^{2}+m^{2})[(\bm{p}+\bm{l})^{2}+m^{2}]}$
$\displaystyle\quad\times\left(\frac{1}{\bm{l}^{2}+m^{2}}-\frac{1}{(\bm{p}+\bm{l})^{2}+m^{2}}\right)^{2}\,.$

We have used the notation Eq. (1) and the subsequent Wick rotation Eq. (2).
Since the integrals are massive and have no infrared divergence, a simple
ultraviolet power-counting shows that the above integrals are convergent near
$4$ dimensions. These integrals therefore have finite non-negative values,

$I^{d}_{2,2}\geq 0,\qquad I^{d}_{3,1}+I^{d}_{1,3}-2I^{d}_{2,2}\geq 0\,.$ (26)

By IBP reduction as described in Section 2.1, the above inequalities translate
into constraints on the finite master integrals on the LHS of Eq. (9). The
needed IBP reduction results are shown in Eqs. (10) and (11), in which we will
rewrite $p^{2}=-M^{2}$. We use the parametrization Eq. (13) to factor out the
positive tadpole integral $I^{d}_{3,0}$, finally arriving at

$\displaystyle\frac{I^{d}_{3,0}}{M^{2}(M^{2}+4m^{2})}\big{[}\left((6-d)M^{2}+4m^{2}\right)\hat{I}^{d}_{2,1}-4m^{2}\big{]}\geq
0\,,$
$\displaystyle\frac{I^{d}_{3,0}}{M^{2}(M^{2}+4m^{2})}\big{[}-\left((8-d)M^{2}+12m^{2}\right)\hat{I}^{d}_{2,1}+(2M^{2}+12m^{2})\big{]}\geq
0\,,$ $\displaystyle\text{for any }d<6\,,$ (27)

i.e.

$\frac{4m^{2}}{(6-d)M^{2}+4m^{2}}\leq\hat{I}^{d}_{2,1}=I^{d}_{2,1}/I^{d}_{3,0}\leq\frac{2M^{2}+12m^{2}}{(8-d)M^{2}+12m^{2}}\,,\quad\text{for
any }d<6\,.$ (28)

It is easily shown that the LHS of the inequality above is always less than
the RHS when $d<6$, so the normalized bubble integral
$\hat{I}^{d}_{2,1}=I^{d}_{2,1}/I^{d}_{3,0}$ is bounded in a finite range. When
$d$ tends to $6$ from below, the LHS and RHS of the inequality both approach
$1$, leading to the prediction that $I^{d}_{2,1}/I^{d}_{3,0}\to 1$ as $d\to
6$; this is exactly as expected since both $I^{d}_{2,1}$ and $I^{d}_{3,0}$
have an ultraviolet pole $1/(d-6)$ with the same coefficient. Specializing to
the case $d=4$, Eq. (28) becomes

$\frac{2m^{2}}{M^{2}+2m^{2}}\leq\hat{I}^{d=4}_{2,1}\leq\frac{M^{2}+6m^{2}}{2M^{2}+6m^{2}}\,,\quad\text{for
}d=4\,.$ (29)

Let us arbitrarily choose an example numerical point,

$M^{2}=2,\quad m^{2}=1\,,$ (30)

At this point, Eq. (29) becomes

$0.5\leq\hat{I}^{d=4}_{2,1}\leq 0.8\,,\quad\text{for }M^{2}=2,m^{2}=1\,.$ (31)

This is consistent with the analytic result given in Appendix A with
$p^{2}=-M^{2}$,

$\displaystyle\hat{I}^{d=4}_{2,1}$
$\displaystyle=I^{d=4}_{2,1}/I^{d=4}_{3,0}=\frac{2m^{2}}{\beta
M^{2}}\log\frac{\beta+1}{\beta-1}\,,$ (32) $\displaystyle\beta$
$\displaystyle\equiv\sqrt{1+\frac{4m^{2}}{M^{2}}}\ ,$ (33)

which evaluates to

$\hat{I}^{d=4}_{2,1}\approx 0.7603459963\,,\quad\text{at }M^{2}=2,m^{2}=1\,.$
(34)

Since our crude positivity bound Eq. (31) and the analytic result Eq. (32)
only depend on the dimensionless ratio $M^{2}/m^{2}$, we plot them in Fig. 2.
It is not surprising that all three curves in the plot tend to $1$ as
$M^{2}/m^{2}\to 0$, since in this case we can set the external momenta to 0,
and the bubble and tadpole integrals in Eq. (9) then become identical. The
general observation is that our positivity bounds can become exact in special
limits of kinematics or the spacetime dimension.

Figure 2: Comparison between ad hoc positivity bounds Eq. (29) and the
analytic result for the bubble integral $\hat{I}^{d=4}_{2,1}$ normalized
according to Eq. (14).

####  2.2.2 Formulation in terms of matrix eigenvalues

In preparation for the introduction of semidefinite programming, we will first
recast positivity constraints in terms of eigenvalues of an appropriate
matrix. To simplify the notation of Eq. (1), let us define the following
shorthand notations for the denominators,

$\rho_{1}=-l^{2}+m^{2},\qquad\rho_{2}=-(p+l)^{2}+m^{2}\,.$ (35)

Let us consider a class of positive-semidefinite integrals of the bubble
family, parametrized by three real numbers $\alpha_{1},\alpha_{2},\alpha_{3}$,

$\int\frac{d^{d}\bm{l}\,e^{\gamma_{E}\epsilon}}{i\pi^{d/2}}\frac{1}{\rho_{1}^{2}\rho_{2}}\left(\alpha_{1}+\frac{\alpha_{2}}{\rho_{1}}+\frac{\alpha_{3}}{\rho_{2}}\right)^{2}=\sum\limits_{i,j}\alpha_{i}M_{ij}\alpha_{j}\\\
={\vec{\alpha}}^{T}\,\mathbb{M}\,\vec{\alpha}\,,$ (36)

where the last line switched to a notation involving a length-3 column vector

$\vec{\alpha}=\begin{pmatrix}\alpha_{1}\\\ \alpha_{2}\\\
\alpha_{3}\end{pmatrix}$ (37)

and a $3\times 3$ symmetric matrix $\mathbb{M}$, given by

$\mathbb{M}=\begin{pmatrix}I^{d}_{2,1}&I^{d}_{3,1}&I^{d}_{2,2}\\\
I^{d}_{3,1}&I^{d}_{4,1}&I^{d}_{3,2}\\\
I^{d}_{2,2}&I^{d}_{3,2}&I^{d}_{2,3}\end{pmatrix}\,,$ (38)

using the index notation Eq. (1). The expression Eq. (36) is non-negative for
any choice of $(\alpha_{1},\alpha_{2},\alpha_{3})$ because after Wick
rotation, $\rho_{1}$ and $\rho_{2}$ are non-negative and the squared
expression is also non-negative. Therefore, the symmetric matrix $\mathbb{M}$
must be positive-semidefinite, represented by the shorthand notation

$\mathbb{M}\succcurlyeq 0\,,\,.$ (39)

By IBP reduction, the matrix entries of $\mathbb{M}$ are integrals which can
be re-expressed as linear sums of the two finite master integrals on the LHS
of Eq. (9). As an example,

$\displaystyle\mathbb{M}_{23}$
$\displaystyle=\mathbb{M}_{32}=\int\frac{d^{d}\bm{l}\,e^{\gamma_{E}\epsilon}}{i\pi^{d/2}}\frac{1}{\rho_{1}^{2}\rho_{2}^{2}}=I^{d}_{2,2},$
(40)

which is then reduced to the finite master integrals according to Eq. (10).
Therefore $\mathbb{M}$ can be written as the sum of two individual master
integral contributions,

$\mathbb{M}=I_{3,0}^{d}\mathbb{M}_{1}+I_{2,1}^{d}\mathbb{M}_{2}=I_{3,0}^{d}\left(\mathbb{M}_{1}+\hat{I}^{d}_{2,1}\mathbb{M}_{2}\right)\,,$
(41)

using the notation $\hat{I}^{d}_{2,1}=I^{d}_{2,1}/I^{d}_{3,0}$ introduced in
Eq. (13). Since $I_{3,0}$ is itself positive, the positive-semidefiniteness of
$\mathbb{M}$ implies

$\widetilde{\mathbb{M}}\equiv\mathbb{M}/I^{d}_{3,0}=\mathbb{M}_{1}+\hat{I}^{d}_{2,1}\mathbb{M}_{2}\succcurlyeq
0\,,$ (42)

again employing the shorthand notation introduced in Eq. (39) to indicate
positive-semidefiniteness. Equivalently, all the eigenvalues of
$\mathbb{M}_{1}+\hat{I}^{d}_{2,1}\mathbb{M}_{2}$ must be non-negative. In Eq.
(41), the matrices $\mathbb{M}_{1}$ and $\mathbb{M}_{2}$ contain entries that
are rational functions of the spacetime $d$ and kinematic variables
$p^{2},m^{2}$, since IBP reduction always produces rational coefficients for
master integrals. Although the general $d$ dependence is not complicated, we
will present $\mathbb{M}_{1}$ and $\mathbb{M}_{2}$ in the $d=4$ case for
brevity of presentation,

$\displaystyle\mathbb{M}_{1}\big{|}_{d=4}$
$\displaystyle=\begin{pmatrix}0&\frac{2m^{2}+M^{2}}{M^{2}\left(4m^{2}+M^{2}\right)}&-\frac{4}{M^{2}\left(4m^{2}+M^{2}\right)}\\\
\frac{2m^{2}+M^{2}}{M^{2}\left(4m^{2}+M^{2}\right)}&\frac{\left(m^{2}+M^{2}\right)\left(6m^{2}+M^{2}\right)}{3m^{2}M^{2}\left(4m^{2}+M^{2}\right)^{2}}&-\frac{2m^{2}-M^{2}}{m^{2}M^{2}\left(4m^{2}+M^{2}\right)^{2}}\\\
-\frac{4m^{2}}{M^{2}\left(4m^{2}+M^{2}\right)}&\frac{M^{2}-2m^{2}}{M^{2}\left(4m^{2}+M^{2}\right)^{2}}&-\frac{2m^{2}-M^{2}}{m^{2}M^{2}\left(4m^{2}+M^{2}\right)^{2}}\\\
\end{pmatrix}\,,$ (43) $\displaystyle\mathbb{M}_{2}\big{|}_{d=4}$
$\displaystyle=\begin{pmatrix}1&-\frac{2m^{2}}{M^{2}\left(4m^{2}+M^{2}\right)}&\frac{2\left(2m^{2}+M^{2}\right)}{m^{2}M^{2}\left(4m^{2}+M^{2}\right)}\\\
-\frac{2m^{2}}{M^{2}\left(4m^{2}+M^{2}\right)}&-\frac{2m^{2}}{M^{2}\left(4m^{2}+M^{2}\right)^{2}}&\frac{2\left(m^{2}+M^{2}\right)}{m^{2}M^{2}\left(4m^{2}+M^{2}\right)^{2}}\\\
\frac{2\left(2m^{2}+M^{2}\right)}{M^{2}\left(4m^{2}+M^{2}\right)}&\frac{2\left(m^{2}+M^{2}\right)}{M^{2}\left(4m^{2}+M^{2}\right)^{2}}&\frac{2\left(m^{2}+M^{2}\right)}{m^{2}M^{2}\left(4m^{2}+M^{2}\right)^{2}}\\\
\end{pmatrix}\,.$ (44)

Now we treat $\hat{I}^{d}_{2,1}$ as an undetermined parameter to be
constrained by Eq. (42). Again we look at the example numerical point as in
Eq. (30),

$M^{2}=2,\quad m^{2}=1\,,$ (45)

and plot the three eigenvalues of the $3\times 3$ matrix
$\widetilde{\mathbb{M}}=\mathbb{M}_{1}+\hat{I}^{d}_{2,1}\mathbb{M}_{2}$ in
Fig. 3.

Figure 3: Three eigenvalues of $\widetilde{\mathbb{M}}$ defined in Eq. (42) as
a function of $\hat{I}^{d}_{2,1}$, at the kinematic point Eq. (30). The lowest
eigenvalue corresponds to the bottom orange curve and the remaining two
eigenvalues correspond to the upper black curves.

It can be seen in the figure that most of the parameter range shown is ruled
out due to the presence of a negative eigenvalue indicated by the lowest
orange curve. In Fig. 4, we zoom in to a smaller parameter region and only
plot the smallest eigenvalue, since the matrix is positive semidefinite as
long as the smallest eigenvalue is non-negative.

Figure 4: Magnified version of the vicinity of a small region of Fig. 3 in
which the lowest eigenvalue of $\widetilde{\mathbb{M}}$, shown in the curve,
is non-negative.

The allowed parameter region shown in the plot is

$0.630\leq\hat{I}^{d=4}_{2,1}\leq 0.847\,,$ (46)

which provides a more stringent lower bound than the previous result Eq. (31).
In fact, far better lower and upper bounds can be achieved in this approach by
using an ansatz larger than the one in Eq. (36). For example, let us replace
the squared term in Eq. (36) by an arbitrary degree-3 polynomial in
$1/\rho_{1}$ and $1/\rho_{2}$, parametrized by 10 undetermined free
coefficients. We obtain the constraint,

$0<\int\frac{d^{d}\bm{l}\,e^{\gamma_{E}\epsilon}}{i\pi^{d/2}}\frac{1}{\rho_{1}^{2}\rho_{2}}\left(\alpha_{1}+\frac{\alpha_{2}}{\rho_{1}}+\frac{\alpha_{3}}{\rho_{2}}+\frac{\alpha_{4}}{\rho_{1}^{2}}+\frac{\alpha_{5}}{\rho_{1}\rho_{2}}+\frac{\alpha_{6}}{\rho_{2}^{2}}+\frac{\alpha_{7}}{\rho_{1}^{3}}+\frac{\alpha_{8}}{\rho_{1}^{2}\rho_{2}}+\frac{\alpha_{9}}{\rho_{1}\rho_{2}^{2}}+\frac{\alpha_{10}}{\rho_{2}^{3}}\right)^{2}\,.$
(47)

Repeating the above analysis, we obtain a constraint for $\hat{I}^{d}_{2,1}$
similar to Eq. (42), except that the matrices involved have $10\times 10$
sizes. The ten eigenvalues as functions of $\hat{I}^{d}_{2,1}$ are plotted in
Fig. 5, as an “upgraded” version of Fig. 3 with a larger ansatz size. Some of
the ten eigenvalues are too close to each other on the plot to be seen
individually, but only the smallest eigenvalue (represented by the lowest
curve in orange color) matters as it determines whether we can satisfy the
constraint that all eigenvalues are positive.

Figure 5: Eigenvalues as a function $\hat{I}^{d}_{2,1}$, for the $10\times 10$
symmetric matrix that represent the quadratic dependence of the RHS of Eq.
(47) on the $\alpha_{i}$ parameters after factoring out $I^{d}_{3,0}$.

In Fig. (6), we zoom in to the small allowed range for the parameter
$\hat{I}^{d}_{2,1}$.

Figure 6: A magnified version of Fig. 5 around the small region of
$\hat{I}^{d}_{2,1}$ where all eigenvalues are positive, showing only the
smallest eigenvalue.

The allowed region, as shown in the plot, is

$0.7598\leq\hat{I}^{d=4}_{2,1}\leq 0.7610\,,$ (48)

which tightly constrains $\hat{I}^{d=4}_{2,1}$ around its true value
$\hat{I}^{d=4}_{2,1}\approx 0.7603$ from evaluating the known analytic result
at $d=4,M^{2}=-p^{2}=2$, with a relative error of around $10^{-3}$.

The above two-sided bounds are rigorous, but we will also explore a
prescription to assign a “central value”, or “best estimate”, of the value of
the integrals. The prescription described below, though not justified from
first principles, empirically achieve a closer agreement with true values of
the integrals than the rigorous bounds in the examples in this paper. The
prescription is simply finding the value of $\hat{I}^{d}_{2,1}$ which
maximizes the smallest eigenvalue of the matrix that is required to be
positive semidefinite, e.g. the matrix $\widetilde{M}$ of Eq. (41). For the
positivity constraint Eq. (47) with 10 free parameters, the prescription picks
the value of $\hat{I}^{d}_{2,1}$ corresponding to the maximum of the curve in
Fig. 6, which deviates from the exact result, again at the example point
$d=4,m=1,M^{2}=-p^{2}=2$, by only a relative error of about $10^{-6}$. In this
case, the prescription happens to produce a value that is very close to the
exact result, but typically we observe the prescription to give a “central
value” that is one to two orders of magnitude better than the accuracy
indicated by rigorous bounds.

####  2.2.3 High-precision evaluation using semidefinite programming

We have formulated positivity constraints in terms of eigenvalues of a matrix
which, in the toy example of the one-loop bubble integrals, depends linearly
on only one undetermined parameter, as shown in Eq. (42). For more complicated
Feynman integrals, there will be more than one master integrals to be
evaluated and all of them will be considered as undetermined parameters. So a
search in higher-dimensional space is needed to locate the region in which all
eigenvalues of the matrix are non-negative, and this can become
computationally challenging. Fortunately, very efficient algorithms exist to
solve _semidefinite programming problems_ in mathematical optimization
vandenberghe1996semidefinite . Loosely speaking, semidefinite programs are
generalizations of linear programs allowing not only linear constraints but
also positive semidefiniteness constraints on matrices that have linear
dependence on the optimization variables. Here we show how our problem of
constraining unknown master integrals can be stated as semidefinite
programming problems. Following the treatment of Section 2.2.2 above, finding
the minimum allowed value of $\hat{I}^{d}_{2,1}$ can be formulated as

minimize $\displaystyle\hat{I}^{d}_{2,1}\,,$ (49) subject to
$\displaystyle\mathbb{M}_{1}+\hat{I}^{d}_{2,1}\cdot\mathbb{M}_{2}\succcurlyeq
0\,,$

which is in the form of a semidefinite program. We used the $\succcurlyeq 0$
notation, already introduced in Eq. (39), to indicate that the matrix on the
LHS must be positive semidefinite. Similarly, finding the minimum allowed
value of $\hat{I}^{d}_{2,1}$ can be formulated as

maximize $\displaystyle\hat{I}^{d}_{2,1}\,,$ (50) subject to
$\displaystyle\mathbb{M}_{1}+\hat{I}^{d}_{2,1}\cdot\mathbb{M}_{2}\succcurlyeq
0\,.$

Finally, to implement our prescription of maximizing the smallest eigenvalue
to find the “central value” of the undetermined master integrals, we introduce
an additional undetermined parameter $\lambda$ and formulate the problem as

maximize $\displaystyle\lambda\,,$ (51) subject to
$\displaystyle\mathbb{M}_{1}+\hat{I}^{d}_{2,1}\cdot\mathbb{M}_{2}-\lambda\mathbb{I}\succcurlyeq
0\,.$

This is again in the form of a semidefinite program, where both
$\hat{I}^{d}_{2,1}$ and $\lambda$ are undetermined parameters whose values
will be fixed to satisfy the optimization objective, namely to maximize
$\lambda$. Note that in this case, finding a optimal solution to the
semidefinite program does not guarantee
$\mathbb{M}_{1}+\hat{I}^{d}_{2,1}\cdot\mathbb{M}_{2}\succcurlyeq 0$, _unless_
the value of $\lambda$ in the solution is non-negative. The value of
$\hat{I}^{d}_{2,1}$ in the solution is then taken as the central value for
this undetermined free parameter.

There exist many computer codes that specialize in solving semidefinite
programs. Wolfram Mathematica has supported semidefinite programming since
version 12 with the SemidefiniteOptimization function, and the default backend
(which can be changed by the Method option) is the open source library CSDP
borchers1999csdp  at least for the problems we deal with, working with double-
precision floating numbers, i.e. the standard machine precision on current
hardware. The SDPA family yamashita2012latest  of computer programs support
computation at double precision as well as a variety of extended precisions,
e.g. double-double precision with SDPA-DD, quadruple-double precision with
SDPA-QD, and arbitrary precision with SDPA-GMP. The SDPB solver by Simmons-
Duffin Simmons-Duffin:2015qma  specializes in polynomial programming problems
in the conformal bootstrap and works in arbitrary precision. Most of the work
in this paper will make use of the SDPA family, while some results from
Mathematica / CSDP will also be shown for the purpose of comparison.

We go on to discuss how to achieve higher numerical precision for the one-loop
bubble integral. We enlarge the ansatz for positive integrals in Eqs. (36) and
(47) to have more parameters, as

$0<\int\frac{d^{d}\bm{l}\,e^{\gamma_{E}\epsilon}}{i\pi^{d/2}}\frac{1}{\rho_{1}^{2}\rho_{2}}P(1/\rho_{1},1/\rho_{2})^{2}\,,$
(52)

where P is an arbitrary polynomial with a maximum degree $N$, i.e. an
arbitrary linear sum of all monomials in $1/\rho_{1}$ and $1/\rho_{2}$, with
each monomial multiplied by a free parameter $\alpha_{i}$. The $N=1$ and $N=3$
cases are shown previously in Eqs. (36) and (47), respectively. Generally, the
number of possible monomials in two variables with maximum degree $N$ is equal
to $(N+1)(N+2)/2$. For each value of $N$ from 1 to 10, we again set
$\hat{I}^{d}_{2,1}=I_{3,0}/I_{2,1}$ at $d=4,m=1,M^{2}=-p^{2}=2$ and solve the
semidefinite programs in Eqs. (49) to (51) to obtain the lower bound
$(\hat{I}^{d}_{2,1})_{\rm min}$, upper bound $(\hat{I}^{d}_{2,1})_{\rm max}$,
and central value $(\hat{I}^{d}_{2,1})_{\rm central}$ for the undetermined
parameter $\hat{I}^{d}_{2,1}$. To ensure numerical stability, we use the SDPA-
DD solver working at double-double precision. Then we compare with the exact
result $(\hat{I}^{d}_{2,1})_{\rm exact}$ to find the relative error in the
best estimate value, defined as

$\left|\frac{(\hat{I}^{d}_{2,1})_{\rm central}}{(\hat{I}^{d}_{2,1})_{\rm
exact}}-1\right|$ (53)

as well as the relative error of the rigorous bounds, defined as

$\left|\frac{(\hat{I}^{d}_{2,1})_{\rm max}-(\hat{I}^{d}_{2,1})_{\rm
min}}{2(\hat{I}^{d}_{2,1})_{\rm exact}}\right|$ (54)

In Fig. 7, we plot the relative errors against the cutoff degree $N$ of the
polynomial $P$ in Eq. (52).

Figure 7: Relative errors in numerical results for $\hat{I}^{d=4}_{2,1}$ at
the kinematic point Eq. (30) from solving positivity constraints Eq. (52)
using semidefinite programming. The horizontal axis is the cutoff degree of
the polynomial $P$. The relative errors are defined by Eq. (54) for the
rigorous bounds and Eq. (53) for the central values.

The plot is on a log scale, and we can see that the numerical results appear
to converge exponentially to the true value, reaching a precision as high as
$10^{-14}$ with a cutoff degree of 10. The central value from our somewhat
arbitrary prescription is seen to be consistently more precise than the
rigorous bounds.

We revisit the issue of numerical stability. In Fig. 8, we compare the
accuracy obtained for the central values obtained by SDPA-DD with double-
double precision, which was used above, and Mathematica / CSDP with double
precision. The two results visibly deviate from each other once the cutoff
degree is 5 or above, and only the double-double precision computation
continues to exhibit exponential reduction in errors as the cutoff degree is
increased. This signals that numerical instability has occurred if one
computes at double precision only. We have checked that the accuracies do not
improve further when computing with quadruple-double precision using SDPA-QD.

Figure 8: Relative errors in numerical results for the central values of
$\hat{I}^{d=4}_{2,1}$ at the kinematic point Eq. (30), obtained with
semidefinite programming solvers working at two different numerical
precisions, namely Mathematica / CSDP working at double precision and SDPA-DD
working at double-double precision. As the plot shows, double-double precision
is needed for the relative errors to improve exponentially beyond a cutoff
degree of 5.

###  2.3 Positivity constraints in Feynman parameter space

In Section 2.2, we used positivity constraints in loop momentum space to
evaluate bubble integrals defined in Eq. (1) in the case $p^{2}<0$ when it is
possible to Wick-rotate the integrals into Euclidean spacetime with real-
valued external momenta. We now use positivity constrains in Feynman parameter
space instead to evaluate bubble integrals for any value of $p^{2}$ less than
$4m^{2}$, which is what is commonly referred to as the “Euclidean region”
where the bubble integrals have no imaginary parts. We write down the Feynman
parameter representation of bubble integrals defined in Eq. (1),

$\displaystyle I_{a_{1},a_{2}}^{d}$
$\displaystyle\equiv\int\frac{d^{d}l\,e^{\gamma_{E}\epsilon}}{i\pi^{d/2}}\frac{1}{(-l^{2}+m^{2})^{a_{1}}[-(p+l)^{2}+m^{2}]^{a_{2}}}$
$\displaystyle=\frac{\Gamma(a_{1}+a_{2}-d/2)e^{\gamma_{E}\epsilon}}{\Gamma(a_{1})\Gamma(a_{2})}\int_{0}^{\infty}dx_{1}\int_{0}^{\infty}dx_{2}\,\delta(1-x_{1}-x_{2})$
$\displaystyle\quad\times
x_{1}^{a_{1}-1}x_{2}^{a_{2}-1}\frac{\mathcal{U}(x_{1},x_{2})^{a_{1}+a_{2}-d}}{\mathcal{F}(x_{1},x_{2})^{a_{1}+a_{2}-d/2}}$
(55)
$\displaystyle=\frac{\Gamma(a_{1}+a_{2}-d/2)e^{\gamma_{E}\epsilon}}{\Gamma(a_{1})\Gamma(a_{2})}\int_{0}^{1}dx\,x^{a_{1}-1}(1-x)^{a_{2}-1}\frac{1}{\mathcal{F}(x)^{a_{1}+a_{2}-d/2}}\,,$
(56)

where the graph polynomials were already given in Eq. (57) for the Schwinger
parametrization, printed again here:

$\mathcal{U}(x_{1},x_{2})=x_{1}+x_{2},\quad\mathcal{F}(x_{1},x_{2})=m^{2}(x_{1}+x_{2})^{2}-p^{2}x_{1}x_{2}-i0^{+}\,.$
(57)

In the last line of Eq. (56), we integrated $x_{2}$ over the delta function in
Eq. (55) to arrive at Eq. (56) with

$\mathcal{F}(x)\equiv\mathcal{F}(x,1-x)=m^{2}-p^{2}x(1-x)-i0^{+}\,,$ (58)

and

$\mathcal{U}(x)\equiv\mathcal{U}(x,1-x)\equiv 1\,,$ (59)

which appears as a unit numerator in Eq. (56).

_Aside:_ we note that Eq. (55) has the property that when ignoring the Dirac
delta function $\delta(1-\sum_{i}x_{i})$, the rest of the expression (with the
integration measure taken into account) is invariant under the rescaling

$x_{i}\to\lambda x_{i},$ (60)

where $\lambda$ is an arbitrary nonzero real number. This is called
_projective invariance_ and holds for the Feynman parameter form of arbitrary
Feynman integrals written down in Eq. (80). The deeper reason is that Feynman
parameter integrals can generally be written as integrals in real projective
space $\mathbb{R}\mathbb{P}^{N-1}$ (see e.g. Section 2.5.3 of Ref.
Weinzierl:2022eaz ). $\mathbb{R}\mathbb{P}^{N-1}$ is the space of $N$ real
coordinates $x_{i}$, excluding the origin, where any ray, i.e. a set of points
related to each other by a rescaling Eq. (60), is identified as the same
point. Abusing the language of gauge theory, Eq. (60) is a gauge symmetry and
$\delta(1-\sum_{i}x_{i})$ in Eq. (55) is a gauge-fixing term that restricts
the integration to one of the infinitely many gauge-equivalent slices. The
Fadeev-Popov Jacobian associated with this gauge-fixing term is unity and
therefore does not appear explicitly. Projective invariance is not a
prerequisite for following the rest of the paper, though it helps motivate
some of the developments.

Now we specialize to the following kinematic region for bubble integrals,

$0<p^{2}<4m^{2},$ (61)

i.e. with the value of $p^{2}$ below the Cutkosky cut threshold but cannot be
trivially Wick-rotated into Euclidean spacetime. We have

$m^{2}-p^{2}x(1-x)\geq m^{2}-p^{2}/4>0,$ (62)

so the $-i0^{+}$ prescription in Eq. (58) is negligible and can be dropped,
and the integral is real. For the rest of the paper, we will adopt the common
terminology of the Euclidean region to be the kinematic region in which all
graph polynomials are non-negative and the Feynman integral is real-valued due
to the lack of Cutkosky cuts. Generally such integrals cannot be embedded into
Euclidean spacetime. This is e.g. the working definition when the literature
refers to the Euclidean region of Feynman integrals with massless external
legs, since nonzero massless momenta cannot be literally embedded into
Euclidean spacetime.

If we set $a_{2}=1$, we can invert Eq. (56) to obtain

$\displaystyle\quad\int_{0}^{1}dx\,x^{a_{1}-1}\frac{1}{\mathcal{F}(x)^{1+a_{1}-d/2}}$
$\displaystyle=\int_{0}^{1}dx\,x^{a_{1}-1}\frac{1}{\left[m^{2}-p^{2}x(1-x)\right]^{1+a_{1}-d/2}}$
$\displaystyle=\frac{\Gamma(a_{1})}{e^{\gamma_{E}\epsilon}\,\Gamma(a_{1}+1-d/2)}I^{d}_{a_{1},1}\,,$
(63)

It will be more useful to have a version of the above equation with a fixed
exponent for $\mathcal{F}(x)$ on the LHS, even when the value of $a_{1}$
changes. Below is a version with a fixed exponent $d/2-3$ for
$\mathcal{F}(x)$, obtained by replacing $d\rightarrow d+2(a_{1}-2)$ in Eq.
(63) and multiplying by a constant prefactor,

$\displaystyle\quad
2(m^{2})^{3-d/2}\int_{0}^{1}dx\,x^{a_{1}-1}\frac{1}{\mathcal{F}(x)^{3-d/2}}$
$\displaystyle=2\int_{0}^{1}dx\,x^{a_{1}-1}\frac{1}{\left[1-p^{2}x(1-x)/m^{2}\right]^{3-d/2}}$
$\displaystyle=\frac{2\Gamma(a_{1})(m^{2})^{1+\epsilon}}{e^{\gamma_{E}\epsilon}\,\Gamma(3-d/2)}I^{d+2a_{1}-4}_{a_{1},1}$
$\displaystyle=I^{d+2a_{1}-4}_{a_{1},1}/I^{d}_{3,0}\,,$ (64)

where the last line used the explicit result for $I^{d}_{3,0}$ in Eq. (155).
We define

$\hat{F}(x)=\mathcal{F}(x)/m^{2}=1-p^{2}x(1-x)/m^{2}\,,$ (65)

and rewrite Eq. (64) as

$2\int_{0}^{1}dx\,x^{a_{1}-1}\frac{1}{\hat{F}(x)^{3-d/2}}=I^{d+2a_{1}-4}_{a_{1},1}/I^{d}_{3,0}\,.$
(66)

The RHS of Eq. (66) can be simplified further, as IBP identities and
dimension-shifting identities can be applied to reduce
$I^{d+2a_{1}-4}_{a_{1},1}$ to a linear combination of the two master
integrals, $I^{d}_{2,1}$ and $I^{d}_{3,0}$. We will only use Eq. (66) in the
case $a_{1}\geq 2$, when the RHS involves a bubble integral in spacetime
dimension greater than or equal to $d=4-2\epsilon$. As an example, consider
the case $a_{1}=3$, and Eq. (66) becomes

$2\int_{0}^{1}dx\,x^{2}\frac{1}{\hat{F}(x)^{3-d/2}}=I^{d+2}_{3,1}/I^{d}_{3,0}\,.$
(67)

Then we simplify the expression using the IBP reduction result Eq. (11) with
$d$ replaced by $d+2$, obtaining

$2\int_{0}^{1}dx\,x^{2}\frac{1}{\hat{F}(x)^{3-d/2}}=\frac{1}{2p^{2}(4m^{2}-p^{2})}\Big{(}[(2-d)p^{2}+4m^{2}]I^{d+2}_{2,1}+2(p^{2}-2m^{2})I^{d+2}_{3,0}\Big{)}/I^{d}_{3,0}\,.$
(68)

Finally, applying dimension-shifting identities Eqs. (22) and (23), the above
equation becomes

$\displaystyle 2\int_{0}^{1}dx\,x^{2}\frac{1}{\hat{F}(x)^{3-d/2}}$
$\displaystyle=\left(\frac{(d-2)p^{2}-4m^{2}}{4(d-3)p^{2}}I^{d}_{2,1}+\frac{m^{2}}{(d-3)p^{2}}I^{d}_{3,0}\right)/I^{d}_{3,0}$
$\displaystyle=\frac{(d-2)p^{2}-4m^{2}}{4(d-3)p^{2}}\hat{I}^{d}_{2,1}+\frac{m^{2}}{(d-3)p^{2}}\,.$
(69)

This concludes our example for simplifying the RHS of Eq. (66) in the case of
$a_{1}=3$. We now formulate a first version of positivity constraints for any
$d<6$, to be improved upon later, as,

$0\leq 2\int_{0}^{1}dx\,xP(x)^{2}\frac{1}{\hat{F}(x)^{3-d/2}}\,,$ (70)

where $P(x)$ is an arbitrary polynomial in $x$, which is analogous to the
arbitrary polynomial $P(1/\rho_{1},1/\rho_{2})$ in Eq. (52) used to construct
positive integrals in momentum space. In the special case $P(x)=1$, the RHS of
Eq. (70) is simply proportional to $I^{d}_{2,1}$ in unshifted spacetime
dimension $d$, according to Eq. (66). Since $x$ is non-negative in the range
of integration $0\leq x\leq 1$, $xP(x)^{2}$ is non-negative. We have chosen to
use $xP(x)^{2}$ instead of just $P(x)^{2}$ to ensure that each monomial in the
expanded expression contains at least one power of $x$ and is related to an
ultraviolet convergent integral in Eq. (66) as discussed above. Another valid
choice is $(1-x)P(x)^{2}$, but this will not give more constraints for the
bubble integral, because in Eq. (56), $\mathcal{F}(x)=p^{2}-m^{2}x(1-x)$ is
invariant under the exchange $x\leftrightarrow 1-x$, owing to a reflection
symmetry of the bubble diagram.

It is possible to slightly refine the inequality Eq. (70) to make the
constraint stronger. As we assume $p^{2}>0$, we have

$\mathcal{F}(x)=m^{2}-p^{2}x(1-x)<m^{2},\quad\hat{F}(x)=\mathcal{F}(x)/m^{2}<1\,.$
(71)

So for any $d<6$, we can modify Eq. (70) with an extra term,

$0\leq
2\int_{0}^{1}dx\,xP(x)^{2}\left(\frac{1}{\hat{F}(x)^{3-d/2}}-1\right)\,.$ (72)

Eq. (72) can also be written in a form that manifests the projective
invariance discussed around Eq. (60),

$\displaystyle 0$ $\displaystyle\leq
2\int_{0}^{\infty}dx_{1}\int_{0}^{\infty}dx_{2}\,\delta(1-x_{1}-x_{2})$
$\displaystyle\quad\times\frac{x_{1}}{U(x_{1},x_{2})}P\left(\frac{x_{1}}{\mathcal{U}(x_{1},x_{2})}\right)^{2}\left(\frac{\mathcal{U}(x_{1},x_{2})^{3-d}}{\hat{F}(x)^{3-d/2}}-\frac{1}{\mathcal{U}(x_{1},x_{2})^{3}}\right)\,,$
(73)

but we will use the form Eq. (72) below.

To be more concrete, in Eq. (72), we use

$P(x)=\alpha_{1}+\alpha_{2}x^{2}+\dots+\alpha_{N}x^{N}\,,$ (74)

where $N$ is the cutoff degree of the polynomial and $\alpha_{i}$ with $1\leq
i\leq N$ are free parameters, and Eq. (72) must hold for any values of the
$\alpha_{i}$ parameters. For each monomial from expanding $xP(x)^{2}$, the
first term in the curly bracket of Eq. (72) gives a bubble integral in a
shifted dimension, normalized against the tadpole integral $I^{d}_{3,0}$,
according to the formula Eq. (66), while the second term in the curly bracket
of Eq. (72) contributes to an integral of a monomial in $x$ over $0\leq x\leq
1$ which can be evaluated trivially. Using both dimension shifting identities
and IBP identities, the bubble integrals produced above are rewritten as
linear combinations of finite master integrals Eq. (9). Therefore Eq. (72) is
turned into the form

${\vec{\alpha}}^{T}\,\mathbb{M}\,\vec{\alpha}\geq 0\,,$ (75)

similar to the momentum-space version Eq. (36), with

$\mathbb{M}=\mathbb{M}_{1}+\hat{I}^{d}_{2,1}\mathbb{M}_{2}\,,$ (76)

where we use the definition $\hat{I}^{d}_{2,1}=I^{d}_{2,1}/I^{d}_{3,0}$ as
before, and the “inhomogeneous” term $\mathbb{M}_{1}$ receives contribution
from both constant terms in Eq. (72) and tadpole integrals coming from
dimension-shifting and IBP identities. Following analogous developments in
Sections 2.2.2 and 2.2.3, we solve the constraint

$\mathbb{M}_{1}+\hat{I}^{d}_{2,1}\mathbb{M}_{2}\succcurlyeq 0$ (77)

while minimizing or maximizing $\hat{I}^{d}_{2,1}$ to find rigorous bounds for
$\hat{I}^{d}_{2,1}$, or alternatively maximizing the smallest eigenvalue of
$\mathbb{M}_{1}+\hat{I}^{d}_{2,1}\mathbb{M}_{2}$ to find a central value for
$\hat{I}^{d}_{2,1}$ using the same prescription as described before.

As an example, we pick numerical values

$p^{2}=2,m=1\,,$ (78)

for bubble integrals in $d=4$, and compare our numerical results against the
exact result. As $p^{2}$ is positive, though below the Cutkosky cut threshold
$4m^{2}$, the Euclidean momentum space treatment in Section 2.2 is not
applicable. SDPA-QD working at quadruple-double precision is used to compute
central values and SDPA-GMP working at 8 times the double precision is used to
compute rigorous bounds. In Fig. 9, we plot the relative error of the central
value for $\hat{I}^{d}_{2,1}$ as well as the relative error of the rigorous
bounds for $\hat{I}^{d}_{2,1}$.

Figure 9: Relative errors in numerical results for $\hat{I}^{d=4}_{2,1}$ at
the kinematic point Eq. (78) from solving positivity constraints Eq. (72) in
Feynman-parameter space using semidefinite programming. The horizontal axis is
the cutoff degree of the polynomial $P$. The relative errors are defined by
Eq. (54) for the rigorous bounds and Eq. (53) for the central values.

We can see on the log-scale plot that the numerical result again converge
exponentially to the exact result as the cutoff degree is increased. In
particular, with cutoff degree $N=14$, the numerical result is

$\hat{I}_{2,1}^{d=4}\Big{|}_{p^{2}=2,\,m=1}\approx 1.57079632679413\,,$ (79)

which is slightly smaller than the exact result, with a relative error of
$4.9\times 10^{-13}$.

###  2.4 Constraints for expansions in dimensional regularization parameter

The methods described in Sections 2.2 and 2.3 are applicable to fixed
spacetime dimensions, i.e. $d=4-2\epsilon$ with fixed values of $\epsilon$,
which can be 0 if we target the 4-dimensional case, or any value larger than
$(-1)$ which will preserve the ultraviolet convergence properties of the
integrals involved in the text above. However, for practical applications,
Feynman integrals typically need to be evaluated as a Laurent expansion in
$\epsilon$. In the examples given in this paper, we can choose master
integrals which are finite as $\epsilon\to 0$, so the task is to calculation
their Taylor expansions in $\epsilon$. Any divergent integral can be reduced
to rational-linear combinations of the master integrals, with all divergences
absorbed into $\epsilon$ poles of the coefficients.111In fact, it is believed
that in general, one can choose “quasi-finite” master integrals
vonManteuffel:2014qoa  which are convergent as $\epsilon\to 0$ except for a
possible $1/\epsilon$ pole that appears as an overall prefactor in the Feynman
parameter representation.

We now two strategies for calculating the $\epsilon$ expansion, using the one-
loop bubble integral example. The first strategy presented below is directly
formulating positivity constraints for the $\epsilon$ expansion terms, and the
second strategy presented is numerical differentiation of the results with
respect to $\epsilon$ around $\epsilon=0$.

####  2.4.1 Generic constraints

We will take a break from the bubble integrals and write down the general form
of the Feynman parametrization for an $L$-loop integral with $n$ propagators,

$\displaystyle\quad
I^{d}_{a_{1},a_{2},\dots,a_{n}}\equiv\left(\prod_{i=1}^{L}\int\frac{d^{d}l_{i}\,e^{\gamma_{E}\epsilon}}{i\pi^{d/2}}\right)\frac{1}{\rho_{1}^{a_{1}}\rho_{2}^{a_{2}}\dots\rho_{n}^{a_{n}}}$
$\displaystyle=\frac{\Gamma(a-Ld/2)e^{\gamma_{E}\epsilon}}{\Gamma(a_{1})\Gamma(a_{2})\dots\Gamma(a_{n})}\int_{x_{i}\geq
0}d^{n}x_{i}\,\delta\left(1-\sum
x_{i}\right)\,\left(\prod_{i}x_{i}^{a_{i}-1}\right)\frac{\mathcal{U}(x_{i})^{a-(L+1)d/2}}{\mathcal{F}(x_{i})^{a-Ld/2}}\,,$
(80)

where $a\equiv\sum a_{i}$, and $\mathcal{U}$ and $\mathcal{F}$ are graph
polynomials. In the Euclidean region, i.e. when external kinematics do not
allow any Cutkosky cuts, $\mathcal{U}$ and $\mathcal{F}$ are non-negative in
the range of integration. To make it easier to formulate positivity
constraints, we adjust constant prefactors and slightly rewrite Eq. (80) as

$\displaystyle\quad\tilde{I}^{d}_{a_{1},a_{2},\dots,a_{n}}\equiv\frac{\Gamma(a_{1})\Gamma(a_{2})\dots\Gamma(a_{n})}{\Gamma(a-Ld/2)e^{\gamma_{E}\epsilon}}I^{d}_{a_{1},a_{2},\dots,a_{n}}$
$\displaystyle=\frac{\Gamma(a_{1})\Gamma(a_{2})\dots\Gamma(a_{n})}{\Gamma(a-Ld/2)e^{\gamma_{E}\epsilon}}\left(\prod_{i=1}^{L}\int\frac{d^{d}l_{i}\,e^{\gamma_{E}\epsilon}}{i\pi^{d/2}}\right)\frac{1}{\rho_{1}^{a_{1}}\rho_{2}^{a_{2}}\dots\rho_{n}^{a_{n}}}$
$\displaystyle=\int_{x_{i}\geq 0}d^{n}x_{i}\,\delta\left(1-\sum
x_{i}\right)\,\left(\prod_{i}x_{i}^{a_{i}-1}\right)\frac{\mathcal{U}(x_{i})^{a-(L+1)d/2}}{\mathcal{F}(x_{i})^{a-Ld/2}}\,.$
(81)

We will restrict our attentions to values of $d$ and $a_{i}$ under which the
RHS of Eq. (81) is convergent. Since there are otherwise no restrictions on
$a_{i}$, generally the integrals are not master integrals which are usually
chosen to have small values of $a_{i}$. We set

$d=d_{0}-2\epsilon\,,$ (82)

where $d_{0}$ is usually an integer spacetime dimension such as 4. The Taylor
expansion of the LHS of Eq. (81) is written as

$\tilde{I}_{a_{1},a_{2},\dots,a_{n}}=\tilde{I}_{a_{1},a_{2},\dots,a_{n}}\big{|}_{\epsilon^{0}}+\epsilon\cdot\tilde{I}_{a_{1},a_{2},\dots,a_{n}}\big{|}_{\epsilon^{1}}+\epsilon^{2}\cdot\tilde{I}_{a_{1},a_{2},\dots,a_{n}}\big{|}_{\epsilon^{2}}\,\dots$
(83)

The only $\epsilon$ dependence of the RHS of Eq. (81) is in the exponents on
the graph polynomials, so the $\mathcal{O}(\epsilon^{k})$ term in the Taylor
expansion is

$\tilde{I}_{a_{1},a_{2},\dots,a_{n}}\Big{|}_{\epsilon^{k}}=\int_{x_{i}\geq
0}d^{n}x_{i}\left(1-\sum
x_{i}\right)\,\left(\prod_{i}x_{i}^{a_{i}-1}\right)\frac{\mathcal{U}(x_{i})^{a-(L+1)d_{0}/2}}{\mathcal{F}(x_{i})^{a-Ld_{0}/2}}\frac{1}{k!}\log^{k}\frac{\mathcal{U}^{L+1}}{\mathcal{F}^{L}}\,.$
(84)

Note that $\mathcal{U}^{L+1}/\mathcal{F}^{l}$ in the equation above is a
quantity that is invariant under the rescaling symmetry Eq. (60), since
$\mathcal{U}$ is a homogeneous polynomial of degree $l$ and $\mathcal{F}$ is a
homogeneous polynomial of degree $l+1$. For the Euclidean region, as
$\mathcal{U}$ and $\mathcal{F}$ are positive in the range of integration, no
branch-cut singularities from the logarithm are encountered. Now we write down
a positivity constraint using our usual trick of constructing non-negative
integrands from squares of polynomials,

$0\leq\int_{x_{i}\geq 0}d^{n}x_{i}\left(1-\sum
x_{i}\right)\,\left(\prod_{i}x_{i}^{a_{i}-1}\right)\frac{\mathcal{U}(x_{i})^{a-(L+1)d_{0}/2}}{\mathcal{F}(x_{i})^{a-Ld_{0}/2}}P^{2}\left(\log\frac{\mathcal{U}^{L+1}}{\mathcal{F}^{L}}\right)\,,$
(85)

where $P$ is an arbitrary polynomial (of the argument in the bracket) under a
cutoff degree $N$, as a sum of monomials each multiplied by a free parameter.
After expanding the square of the polynomial, each monomial term is identified
with a term in the Taylor expansion over $\epsilon$ using Eq. (84). Using the
same manipulations as in Sections 2.2 and 2.3, Eq. (85) implies that a certain
symmetric matrix is positive semidefinite. We first define an auxiliary
notation

$H_{k}\equiv(k!)\tilde{I}^{d}_{a_{1},a_{2},\dots,a_{n}}\Big{|}_{\epsilon^{k}}\,.$
(86)

Then we have

$\displaystyle\begin{pmatrix}H_{0}&H_{1}&H_{2}&\dots&H_{N}\\\
H_{1}&H_{2}&H_{3}&\dots&H_{N+1}\\\ H_{2}&H_{3}&H_{4}&\dots&H_{N+2}\\\
\vdots&\vdots&\vdots&\ddots&\vdots\\\
H_{N}&H_{N+1}&H_{N+2}&\dots&H_{2N}\end{pmatrix}\succcurlyeq 0\,,$ (87)

where we again used the notation $\succcurlyeq 0$ to indicate that a matrix is
positive semidefinite. The matrix is in a special form called a Hankel matrix,
where all matrix entries are defined through a sequence $H_{0}$, $H_{1}$, …,
$H_{N}$. Hankel matrices have also appeared in the context of EFT positivity
bounds, e.g. in Ref. Arkani-Hamed:2020blm .

Eq. (87) is extremely general and applies to any convergent Feynman integral
(or quasi-finite Feynman integral vonManteuffel:2014qoa  after dropping an
overall divergent prefactor) in the Euclidean region with arbitrary powers of
propagators and no numerators. This tells us that the $\epsilon$ expansion of
such Feynman integrals are not arbitrary but are constrained by positivity
constraints which, to our best knowledge, have not been previously revealed in
the literature.

####  2.4.2 Taylored constraints for specific Feynman integrals

From Eq. (84), it is not hard to anticipate that more specialized positivity
constraints exist if we focus on a particular family of Feynman integrals if
$\log(\mathcal{U}^{L+1}/\mathcal{F}^{l})$ has either an upper bound or lower
bound, or both, in the range of integration. For example, for one-loop bubble
integrals, the $\mathcal{U}$ polynomial is equal to $x_{1}+x_{2}$ and is set
to 1 by the Dirac delta function in Eq. (80). So we recover the Feynman
parametrization for the one-loop bubble integral, Eq. (56), with
$x_{1}=x,\,x_{2}=1-x,\,\mathcal{U}(x)=1,\,\mathcal{F}(x)=m^{2}-p^{2}x(1-x)$.
We have, for $0<p^{2}<4m^{2}$ under consideration, in the integration range
$0\leq x\leq 1$,

$\log\frac{\mathcal{U}^{L+1}}{\mathcal{F}^{L}}=\log\frac{1}{\mathcal{F}}=\log\frac{1}{m^{2}-p^{2}x(1-x)}\leq\log\frac{1}{m^{2}-p^{2}/4}\equiv\log\max\frac{\mathcal{U}^{L+1}}{\mathcal{F}^{L}}\,.$
(88)

Therefore

$\log\max\frac{\mathcal{U}^{L+1}}{\mathcal{F}^{L}}-\log\frac{\mathcal{U}^{L+1}}{\mathcal{F}^{L}}$
(89)

is a positive quantity. Similarly, if a minimum of
$\log\mathcal{U}^{L+1}/\mathcal{F}^{L}$ exists over the range of integration,
then

$\log\frac{\mathcal{U}^{L+1}}{\mathcal{F}^{L}}-\log\min\frac{\mathcal{U}^{L+1}}{\mathcal{F}^{L}}$
(90)

is a positive quantity. In the bubble integral example, again under
$0<p^{2}<4m^{2}$ and $0\leq x\leq 1$,

$\log\frac{\mathcal{U}^{L+1}}{\mathcal{F}^{L}}=\log\frac{1}{\mathcal{F}}=\log\frac{1}{m^{2}-p^{2}x(1-x)}\geq\log\frac{1}{m^{2}}\equiv\log\min\frac{\mathcal{U}^{L+1}}{\mathcal{F}^{L}}\,.$
(91)

Now we show an example of using Eq. (88) to contrain the $O(\epsilon)$ term in
the expansion of the finite bubble integral $I_{2,1}$, in the style of an ad
hoc constraint as was done for the $\mathcal{O}(\epsilon^{0})$ part in Section
2.2.1.222It is also possible to use Eq. (91) instead, or in combination. Using
the definition $\hat{F}(x)=\mathcal{F}(x)/m^{2}$ in Eq. (65), Eq. (88) is
rewritten as

$\log\frac{1}{\hat{F}(x)}\leq\log\max\frac{1}{\mathcal{F}}=\log\frac{1}{1-p^{2}/(4m^{2})}\,,$
(92)

i.e.,

$\log\frac{1}{1-p^{2}/(4m^{2})}-\log\frac{1}{\hat{F}(x)}\geq 0\,,$ (93)

We expand both the LHS and RHS of Eq. (69), with $d=4-2\epsilon$, as a Taylor
series in $\epsilon$. Equating the $\epsilon^{0}$ terms gives

$2\int_{0}^{1}dx\,x^{2}\frac{1}{\hat{F}(x)}=\frac{p^{2}-2m^{2}}{2p^{2}}\left(\hat{I}_{2,1}\big{|}_{\epsilon^{0}}\right)+\frac{m^{2}}{p^{2}}\,,$
(94)

while equating the $\epsilon^{1}$ terms gives

$2\int_{0}^{1}dx\,x^{2}\frac{1}{\hat{F}(x)}\cdot\log\left(\frac{1}{\hat{F}(x)}\right)=\frac{p^{2}-4m^{2}}{2p^{2}}\left(\hat{I}_{2,1}\big{|}_{\epsilon^{0}}\right)+\frac{p^{2}-2m^{2}}{2p^{2}}\left(\hat{I}_{2,1}\big{|}_{\epsilon^{1}}\right)+\frac{2m^{2}}{p^{2}}\,.$
(95)

Now we use Eq. (93) to write down the positivity constraint

$2\int_{0}^{1}dx\,x^{2}\frac{1}{\hat{F}(x)}\left(\log\frac{1}{1-p^{2}/(4m^{2})}-\log\frac{1}{\hat{F}(x)}\right)\geq
0\,.$ (96)

Then applying Eqs. (94) and (95) leads to a constraint on
$\hat{I}_{2,1}\big{|}_{\epsilon^{1}}$,

$\displaystyle 0$
$\displaystyle\leq\left(\frac{p^{2}-2}{2p^{2}}\log\frac{1}{1-p^{2}/(4m^{2})}+\frac{4m^{2}-p^{2}}{2p^{2}}\right)\left(\hat{I}_{2,1}\big{|}_{\epsilon^{0}}\right)+\frac{2m^{2}-p^{2}}{2p^{2}}\left(\hat{I}_{2,1}\big{|}_{\epsilon^{1}}\right)$
$\displaystyle+\left(\log\frac{1}{1-p^{2}/(4m^{2})}-2\right)\frac{m^{2}}{p^{2}}\,.$
(97)

We plot the RHS of Eq. (97) versus $p^{2}/m^{2}$ in Fig. 10. We can see that
it is indeed non-negative in the range $0<p^{2}/m^{2}<4$.

Figure 10: The RHS of Eq. (97) versus $p^{2}/m^{2}$ in the range of validity
$0<p^{2}<4$.

####  2.4.3 Bubble integral up to second order in $\epsilon$

We proceed to combine positivity constraints for $\epsilon$ expansion
coefficients covered in Sections 2.4.1 and 2.4.2 with the semidefinite
programming technique already used in Sections 2.2 and 2.3, in order to obtain
high-precision results for the $\epsilon$ expansion of the bubble integral.

Recall that we evaluated the $\mathcal{O}(\epsilon^{0})$ part, i.e. the $d=4$
result, for the bubble integral starting from Eq. (70) and the refined version
Eq. (72). For simplicity, we will build upon the first version, Eq. (70), and
use the additional positive building block Eq. (93) to write down the
following positivity constraint at $d=d_{0}-2\epsilon=4-2\epsilon$,

$\displaystyle 0$
$\displaystyle\leq\int_{0}^{1}dx\,xP(x)^{2}\frac{1}{\hat{F}(x)^{3-d_{0}/2}}\left(\log\max\frac{1}{\hat{F}}-\log\frac{1}{\hat{F}(x)}\right)$
(98)
$\displaystyle=\int_{0}^{1}dx\,xP(x)^{2}\frac{1}{[1-x(1-x)p^{2}/m^{2}]}\left(\log\frac{1}{1-p^{2}/(4m^{2})}-\log\frac{1}{1-x(1-x)p^{2}/m^{2}}\right)\,,$
(99)

where $P(x)$ is again an arbitrary polynomial in $x$, and we are free to
choose the maximum degree of monomials that are included, depending on the
accuracy we would like to attain. After expanding $P(x)^{2}$ into a sum of
monomials, the contribution from each monomial can be evaluated following the
same procedure as used in the example in Section 2.4.2. In particular, the
result will be a sum of terms proportional to
$\hat{I}_{2,1}\big{|}_{\epsilon^{0}}$, terms proportional to
$\hat{I}_{2,1}\big{|}_{\epsilon^{1}}$, and constant terms. The remaining
calculation steps are very similar to those of Section 2.3. Re-using the
parametrization Eq. (74) for the polynomial $P$ with a cutoff degree $N=14$,
we again arrive at

${\vec{\alpha}}^{T}\,\mathbb{M}\,\vec{\alpha}\geq 0\,,$ (100)

stating that an appropriate matrix $\mathbb{M}$ is positive semidefinite, i.e.

$\mathbb{M}\succcurlyeq 0\,.$ (101)

In this case, $\mathbb{M}$ is a sum of three terms,

$\mathbb{M}=\mathbb{M}_{1}+I_{2,1}\big{|}_{\epsilon^{0}}\cdot\mathbb{M}_{2}+I_{2,1}\big{|}_{\epsilon^{1}}\cdot\mathbb{M}_{3}\,.$
(102)

where the three matrices $\mathbb{M}_{1,2,3}$, with rational dependence on
$p^{2}$ and $m$, are obtained from dimension shifting, IBP, and finally
$\epsilon$-expansion as in the example of Section 2.4.2. We approximate
$I_{2,1}\big{|}_{\epsilon^{0}}$ to be the central value Eq. (79) obtained in
the previous $d=4$ calculation with the same cutoff degree $N=14$. At this
point, $I_{2,1}\big{|}_{\epsilon^{1}}$ remains the only unknown parameter on
the RHS of Eq. (102), and its allowed range as well as the central value can
be determined by semidefinite programming solvers as covered in previous
sections. We again use SDPA-QD to produce the “central value” for the
numerical result, defined by the same prescription as before, obtaining

$I_{2,1}\big{|}_{\epsilon^{1}}\approx 0.74313814320586\,,\ \text{at
}{p^{2}=2,m=1}\,,$ (103)

which is slightly larger than the exact result with a relative error of
$4.3\times 10^{-12}$.

We continue to present how the $\mathcal{O}(\epsilon^{2})$ term of the bubble
integral is calculated. We use the positivity constraint,

$0\leq\int_{0}^{1}dx\,xP\left(x,\log\frac{1}{\hat{F}}\right)^{2}\left[m^{2}-p^{2}x(1-x)\right]^{d_{0}/2-3}\,,$
(104)

where $P$ is a polynomial in $x$ and $\log 1/\hat{F}(x)$ with maximum degrees
$N_{1}$ and $N_{2}$ in the two variables, parametrized as

$P\left(x,\log\frac{1}{\hat{F}(x)}\right)=\sum\limits_{0\leq i_{1}\leq
N_{1}}\sum\limits_{0\leq i_{2}\leq
N_{2}}\alpha_{i_{1},i_{2}}\,x^{i_{1}}\left(\log\frac{1}{\hat{F}(x)}\right)^{i_{2}}\,.$
(105)

We use $N_{1}=14$ and $N_{2}=1$, so that there are at most 14 power of $x$ in
$P(x)$ and at most one power of $\log[1/\hat{F}(x)]$. After expanding the
$P^{2}$, there are at most two powers of the aforementioned logarithm,
therefore there are $\epsilon$ expansions coefficients at orders
$\epsilon^{0}$, $\epsilon^{1}$, and $\epsilon^{2}$. We obtain an expression of
the form

${\vec{\alpha}}^{T}\,\mathbb{M}\,\vec{\alpha}\geq 0\,,$ (106)

where the column vector $\vec{\alpha}$ groups together all the
$\alpha_{i_{1},i_{2}}$ parameters. The matrix $\mathbb{M}$ is the sum of a
constant term, a term proportional to $I_{2,1}\big{|}_{\epsilon^{0}}$, a term
proportional to $I_{2,1}\big{|}_{\epsilon^{1}}$, and finally a term
proportional to $I_{2,1}\big{|}_{\epsilon^{2}}$. We use previous numerical
results for $I_{2,1}\big{|}_{\epsilon^{0}}$ and
$I_{2,1}\big{|}_{\epsilon^{1}}$, and solve a semidefinite programming problem
to find the central value of $I_{2,1}\big{|}_{\epsilon^{2}}$ to be

$I_{2,1}\big{|}_{\epsilon^{2}}\approx 0.208108744452\,,\ \text{at
}{p^{2}=2,m=1}\,,$ (107)

which is slightly larger than the exact result with a relative error of
$1.9\times 10^{-11}$.

There is no obstruction to obtaining results at even higher orders in the
$\epsilon$ expansion for the bubble integral example. Generally, we calculate
iteratively to higher and higher orders in $\epsilon$, at each step taking
previous numerical results as known input. The positivity constraint is Eq.
(104) with an appropriate cutoff degree for $\log[1/\hat{F}(x)]$ depending on
the desired order in the $\epsilon$ expansion. The RHS of Eq. (104) can be
optionally multiplied by either Eq. (89) or Eq. (90), with $\mathcal{U}=1$ and
$L=1$ in the case of bubble integrals, to give more constraints. This will
produce constraints for bubble integrals to any desired order in the
$\epsilon$ expansion.

###  2.5 $\epsilon$ expansion from numerical differentiation w.r.t. spacetime
dimension

Here we present an alternative method for numerically evaluating the
$\epsilon$-expansion of the normalized bubble integral $\hat{I}^{d}_{2,1}$
defined in Eq. (13). Instead of formulating constraints for the $\epsilon$
expansion coefficients, we calculate the $\hat{I}^{d}_{2,1}$ at numerical
values of the spacetime dimension near 4, and use finite-difference
approximations to obtain derivatives w.r.t. $\epsilon$. The derivatives are
related to the terms in the $\epsilon$ expansion via

$\hat{I}_{2,1}\big{|}_{\epsilon^{k}}=\frac{1}{k!}\frac{d^{k}}{d\epsilon^{k}}\hat{I}^{d=4-2\epsilon}_{2,1}\Big{|}_{\epsilon=0}\,.$
(108)

For an arbitrary function $f(\epsilon)$, we use 4th-order finite-difference
approximations,

$\displaystyle\frac{d}{d\epsilon}f(\epsilon)\bigg{|}_{\epsilon=\epsilon^{0}}$
$\displaystyle\approx\frac{1}{\Delta\epsilon}\bigg{[}\frac{1}{12}f(\epsilon_{0}-2\Delta\epsilon)-\frac{2}{3}f(\epsilon_{0}-\Delta\epsilon)$
$\displaystyle\qquad+\frac{2}{3}f(\epsilon_{0}+\Delta\epsilon)-\frac{1}{12}f(\epsilon_{0}+2\Delta\epsilon)\bigg{]}\,,$
(109)
$\displaystyle\frac{d^{2}}{d\epsilon^{2}}f(\epsilon)\bigg{|}_{\epsilon=\epsilon^{0}}$
$\displaystyle\approx\frac{1}{\Delta\epsilon^{2}}\bigg{[}-\frac{1}{12}f(\epsilon_{0}-2\Delta\epsilon)+\frac{4}{3}f(\epsilon_{0}-\Delta\epsilon)-\frac{5}{2}f(\epsilon_{0})$
$\displaystyle\qquad+\frac{4}{3}f(\epsilon_{0}+\Delta\epsilon)-\frac{1}{12}f(\epsilon_{0}+2\Delta\epsilon)\bigg{]}\,,$
(110)

where $\Delta\epsilon$ is the step size. As the name suggests, these formulas
are exact when $f$ is a polynomial with a degree up to 4. Note that the method
of Section 2.3 can be used to evaluate the normalized bubble integral
$\hat{I}^{d}_{2,1}$ in any spacetime dimension $d<6$, i.e. $\epsilon>-1$, so
Eqs. (109) and (110) can be readily used with $\epsilon_{0}=0$ and a small
$\Delta\epsilon$, chosen to be

$\Delta\epsilon=10^{-3}\,.$ (111)

We again choose kinematic parameter values $p^{2}=2$ and $m=1$. The final
results from numerical differentiation, up to the second order in $\epsilon$,
are

$\hat{I}_{2,1}\big{|}_{\epsilon^{1}}\approx 0.7431381432049\,,$ (112)

which is slightly larger than the exact result with a relative error of
$3.2\times 10^{-12}$, and

$\hat{I}_{2,1}\big{|}_{\epsilon^{2}}\approx 0.20810874450\,,$ (113)

which is slightly larger than the exact result with a relative error of
$2.8\times 10^{-10}$.

##  3 Three-loop banana integrals with unequal masses

###  3.1 Definitions and conventions

Here we present a three-loop example, the so called banana integrals. Banana
integrals are various loop orders have received intense interest from an
analytic perspective due to connections with Calabi-Yau manifold Groote:2005ay
; Klemm:2019dbm ; Bonisch:2020qmm ; Bonisch:2021yfw ; Kreimer:2022fxm ;
Forum:2022lpz ; Pogel:2022ken . We apply our numerical method developed in the
previous section 2 on one-loop bubble integrals, with some minor adaptations,
to evaluate $11$ nontrivial master integrals for the banana diagram. We assume
the readers to be familiar with the previous section as many shared techniques
will not be introduced again. The diagram for the integrals is shown in Fig.
11.

Figure 11: Three-loop banana family of integrals, with internal and external
squared masses labeled.

Due to dimension-shifting identities reviewed in Section 2.1, the $\epsilon$
expansions of integrals in $d=4-2\epsilon$ and $d=2-2\epsilon$ can be related
to each other. We will always use

$d=2-2\epsilon\,,$ (114)

for three-loop banana integrals, which is convenient as the scalar integral
has no ultraviolet divergence in this spacetime dimension. The banana family
of integrals is defined as

$\displaystyle I_{a_{1},a_{2},a_{3},a_{4}}$
$\displaystyle\equiv\left(\prod_{i=1}^{3}\int\frac{d^{d}l_{i}\,e^{\gamma_{E}\epsilon}}{i\pi^{d/2}}\right)\frac{1}{(-l_{1}^{2}+m_{1}^{2})^{a_{1}}}\frac{1}{(-l_{2}^{2}+m_{2}^{2})^{a_{2}}}\frac{1}{(-l_{3}^{2}+m_{3}^{2})^{a_{3}}}$
$\displaystyle\qquad\times\frac{1}{[-(p+l_{1}+l_{2}+l_{3})^{2}+m_{4}^{2}]^{a_{4}}}\,,\qquad\text{with
}d=2-2\epsilon\,.$ (115)

We have suppressed the $d$ dependence on the LHS of the above equation, unlike
the one-loop case Eq. (1), since we will not make use of dimension-shifting in
the treatment of three-loop banana integrals and will exclusively work with
$d=2-2\epsilon$. If any one of the four indices $a_{i}$ is non-positive in Eq.
(115), the remaining propagators have the structure of three one-loop massive
tadpole integrals. For example, if $a_{4}=0$, Eq. (115) clearly factorizes
into the product of three scalar tadpole integrals. For facilitating the
discussion of positivity constraints, it will be convenient to define a
variant of Eq. (115) with slightly adjusted constant factors,

$\displaystyle\hat{I}_{a_{1},a_{2},a_{3},a_{4}}$
$\displaystyle\equiv\left(\prod_{i=1}^{3}\int\frac{d^{d}l_{i}}{i\pi^{d/2}}\right)\frac{1}{\Gamma(4-3d/2)}\frac{1}{(-l_{1}^{2}+m_{1}^{2})^{a_{1}}}\frac{1}{(-l_{2}^{2}+m_{2}^{2})^{a_{2}}}\frac{1}{(-l_{3}^{2}+m_{3}^{2})^{a_{3}}}$
$\displaystyle\qquad\times\frac{1}{[-(p+l_{1}+l_{2}+l_{3})^{2}+m_{4}^{2}]^{a_{4}}}\,.$
(116)

By IBP reduction, all integrals of the banana family, with integer values of
$a_{i}$ in Eq. (116), can be expressed as linear sums of 15 master integrals.
Publicly available software, such as those presented Refs. Anastasiou:2004vj ;
vonManteuffel:2012np ; Smirnov:2008iw ; Smirnov:2019qkx ; Lee:2012cn ;
Lee:2013mka ; Maierhofer:2017gsa , can be used to give a list of master
integrals as well as performing the actual IBP reduction of integrals. There
are 11 nontrivial “top-level” master integrals that are not products of
tadpole integrals, shown in three groups below according to the total number
of the four indices,

$\displaystyle\hat{I}_{1,1,2,2},\,\hat{I}_{1,2,1,2},\,\hat{I}_{1,2,2,1},\,\hat{I}_{2,1,1,2},\,\hat{I}_{2,1,2,1},\,\hat{I}_{2,2,1,1},$
(117)
$\displaystyle\hat{I}_{1,1,1,2},\,\hat{I}_{1,1,2,1},\,\hat{I}_{1,2,1,1},\,\hat{I}_{2,1,1,1},$
$\displaystyle\hat{I}_{1,1,1,1}\,.$

In addition, there are 4 master integrals that are trivial products of tadpole
integrals. These master integrals are chosen as

$\hat{I}_{0,2,2,2},\,\hat{I}_{2,0,2,2},\,\hat{I}_{2,2,0,2},\,\hat{I}_{2,2,2,0}\,,$
(118)

where we raised every propagator to a 2nd power to make the integral UV finite
in $d=2-2\epsilon$. Each of these four master integrals is a product of three
one-loop tadpole integrals given in Eq. (157) with $n=2$, with some adjustment
of the overall factor according to Eq. (116),

$\displaystyle\hat{I}_{0,2,2,2}=\frac{\Gamma^{3}(1+\epsilon)}{\Gamma(1+3\epsilon)}\left(\frac{1}{m_{2}^{2}m_{3}^{2}m_{4}^{2}}\right)^{1+\epsilon},\quad\hat{I}_{2,0,2,2}=\frac{\Gamma^{3}(1+\epsilon)}{\Gamma(1+3\epsilon)}\left(\frac{1}{m_{1}^{2}m_{3}^{2}m_{4}^{2}}\right)^{1+\epsilon},$
$\displaystyle\hat{I}_{2,2,0,2}=\frac{\Gamma^{3}(1+\epsilon)}{\Gamma(1+3\epsilon)}\left(\frac{1}{m_{1}^{2}m_{2}^{2}m_{4}^{2}}\right)^{1+\epsilon},\quad\hat{I}_{2,2,2,0}=\frac{\Gamma^{3}(1+\epsilon)}{\Gamma(1+3\epsilon)}\left(\frac{1}{m_{1}^{2}m_{2}^{2}m_{3}^{2}}\right)^{1+\epsilon}\,.$
(119)

The values of the 11 remaining master integrals in Eq. (117) will be
calculated numerically from positivity constraints. We will work with
kinematic variables in the range

$p^{2}<(m_{1}+m_{2}+m_{3}+m_{4})^{2}\,,$ (120)

i.e. below the particle production threshold, which ensures that the integrals
are real-valued. Similar to the case of one-loop bubble integrals, when the
chosen value of $p^{2}$ is non-negative, we cannot embed the integrals into
Euclidean momentum space and need to use Feynman-parameter space to formulate
positivity constraints. The Feynman parametrization follows from the general
formula Eq. (80) with adjustment of constant factors according to Eq. (116),

$\displaystyle\hat{I}_{a_{1},a_{2},a_{3},a_{4}}$
$\displaystyle=\frac{\Gamma(a-3d/2)/\Gamma(4-3d/2)}{\Gamma(a_{1})\Gamma(a_{2})\Gamma(a_{3})\Gamma(a_{4})}\int_{x_{i}\geq
0}dx_{1}dx_{2}dx_{3}dx_{4}\,\delta(1-x_{1}-x_{2}-x_{3}-x_{4})$
$\displaystyle\quad\times\left(\prod_{i=1}^{4}x_{i}^{a_{i}-1}\right)\frac{\mathcal{U}(x_{i})^{a-2d}}{\mathcal{F}(x_{i})^{a-3d/2}}\,,$
(121)

where we used the definition $a\equiv a_{1}+a_{2}+a_{3}+a_{4}$. The two graph
polynomials $\mathcal{U}$ and $\mathcal{F}$ are, for the banana family of
integrals,

$\displaystyle\mathcal{U}(x_{1},x_{2},x_{3},x_{4})$
$\displaystyle=x_{2}x_{3}x_{4}+x_{1}x_{3}x_{4}+x_{1}x_{2}x_{4}+x_{1}x_{2}x_{3},$
$\displaystyle\mathcal{F}(x_{1},x_{2},x_{3},x_{4})$
$\displaystyle=p^{2}x_{1}x_{2}x_{3}x_{4}+(m_{1}^{2}x_{1}+m_{2}^{2}x_{2}+m_{3}^{2}x_{3}+m_{4}^{2}x_{4})\mathcal{U}(x_{1},x_{2},x_{3},x_{4})\,.$
(122)

Note that we have

$\mathcal{U}(x_{1},x_{2},x_{3},x_{4})\geq
0,\quad\mathcal{F}(x_{1},x_{2},x_{3},x_{4})\geq 0\,,$ (123)

in the integration region of Eq. (121), i.e. $x_{i}\geq 0,\,\sum_{i}x_{i}=1$.
This will help us formulate positivity constraints. As in the one-loop bubble
case, except for the “gauge fixing” Dirac delta function, the rest of Eq.
(121) has the projective invariance Eq. (60), as the $\mathcal{U}$ and
$\mathcal{F}$ polynomials are homogeneously of degree 3 and 4, respectively.

###  3.2 Positivity constraints

We define $\tilde{x}_{i}$ variables

$\tilde{x}_{i}=\frac{\mathcal{U}(x_{1},x_{2},x_{3},x_{4})}{\mathcal{F}(x_{1},x_{2},x_{3},x_{4})}x_{i},\quad
i=1,2,3,4\,,$ (124)

which are invariant under the scaling Eq. (60). We rewrite Eq. (121) as

$\displaystyle\quad\frac{\Gamma(4-3d/2)}{\Gamma(a-3d/2)}\Gamma(a_{1})\Gamma(a_{2})\Gamma(a_{3})\Gamma(a_{4})\,\hat{I}_{a_{1},a_{2},a_{3},a_{4}}$
$\displaystyle=\int_{x_{i}\geq
0}dx_{1}dx_{2}dx_{3}dx_{4}\,\delta(1-x_{1}-x_{2}-x_{3}-x_{4})\left(\prod_{i=1}^{4}\tilde{x}_{i}^{a_{i}-1}\right)\frac{\mathcal{U}(x_{i})^{4-2d}}{\mathcal{F}(x_{i})^{4-3d/2}}\,,$
(125)

again with $a\equiv a_{1}+a_{2}+a_{3}+a_{4}$. Note that if $a_{i}\geq 1$,
$a\geq 4$,

$\frac{\Gamma(4-3d/2)}{\Gamma(a-3d/2)}=\frac{1}{(4-3d/2)(5-3d/2)\dots(a-1-3d/2)}$
(126)

is a rational function in $d$. With any non-negative polynomial
$Q(\tilde{x}_{i})$, we formulate a positivity constraint,

$\displaystyle 0$ $\displaystyle\leq\int_{x_{i}\geq
0}dx_{1}dx_{2}dx_{3}dx_{4}\,\delta(1-x_{1}-x_{2}-x_{3}-x_{4})Q(\tilde{x}_{i})\frac{\mathcal{U}(x_{i})^{4-2d}}{\mathcal{F}(x_{i})^{4-3d/2}}\,,$
(127)

which is compatible with the projective invariance Eq. (60). After expanding
the polynomial $Q(\tilde{x}_{i})$ into a sum of monomials, the contribution of
each monomial $\prod_{i}\tilde{x}_{i}^{a_{i}-1}$ can be written as some
$\hat{I}_{a_{1},a_{2},a_{3},a_{4}}$ multiplied by a prefactor that is rational
in $d$, according to Eq. (125). All such integrals are UV convergent by power
counting and also IR convergent due to internal masses. No change of spacetime
dimensions is involved, unlike the treatment of one-loop bubble integrals in
Section 2.3.

We consider the following choices of $Q(\tilde{x}_{i})$, with the help of an
arbitrary polynomial $P(\tilde{x}_{i})$ under a chosen maximum degree,

$\displaystyle\text{choice 1:}\quad Q(\tilde{x}_{i})$
$\displaystyle=P(\tilde{x}_{i})^{2},$ (128) $\displaystyle\text{choice
2:}\quad Q(\tilde{x}_{i})$ $\displaystyle=\tilde{x}_{1}P(\tilde{x}_{i})^{2},$
(129) $\displaystyle\text{choice 3:}\quad Q(\tilde{x}_{i})$
$\displaystyle=\tilde{x}_{2}P(\tilde{x}_{i})^{2},$ (130)
$\displaystyle\text{choice 4:}\quad Q(\tilde{x}_{i})$
$\displaystyle=\tilde{x}_{3}P(\tilde{x}_{i})^{2},$ (131)
$\displaystyle\text{choice 5:}\quad Q(\tilde{x}_{i})$
$\displaystyle=\tilde{x}_{4}P(\tilde{x}_{i})^{2}\,.$ (132)

With any of the above five choices for $Q(\tilde{x}_{i})$ and with any choice
of $P(\tilde{x}_{i})$, the inequality Eq. (127) must hold. The general form of
$P(\tilde{x}_{i})$ is a sum of all monomials under a chosen cutoff degree,
each multiplied by an arbitrary coefficient. For example, if the cutoff degree
is 1, then $P(\tilde{x}_{i})$ is parametrized as

$\text{cutoff degree 1:}\quad
P(\tilde{x}_{i})=\alpha_{0,0,0,0}+\alpha_{1,0,0,0}\tilde{x}_{1}+\alpha_{0,1,0,0}\tilde{x}_{2}+\alpha_{0,0,1,0}\tilde{x}_{3}+\alpha_{0,0,0,1}\tilde{x}_{4}\,.$
(133)

With cutoff degree $N$, the parametrization is

$\text{cutoff degree $N$:}\quad
P(\tilde{x}_{i})=\sum_{i_{1},i_{2},i_{3},i_{4}\geq
0}^{i_{1}+i_{2}+i_{3}+i_{4}\leq
N}\alpha_{i_{1},i_{2},i_{3},i_{4}}\tilde{x}_{1}^{i_{1}}\tilde{x}_{2}^{i_{2}}\tilde{x}_{3}^{i_{3}}\tilde{x}_{4}^{i_{4}}\,,$
(134)

where the number of free parameters $\alpha_{i_{1},i_{2},i_{3},i_{4}}$ is
equal to ${N+4\choose 4}$ by combinatorics.

Since we have already discussed how to set up semidefinite optimization
programs in the context of one-loop bubble integrals, we will be brief in
covering the analogous steps here. Grouping the
$\alpha_{i_{1},i_{2},i_{3},i_{4}}$ parameters into a column vector
$\vec{\alpha}$ of length ${N+4\choose 4}$, the five choices of $Q$, Eqs. (128)
to (132), lead to

$\displaystyle(\vec{\alpha})^{T}\mathbb{M}^{(1)}\vec{\alpha}\geq
0,\quad(\vec{\alpha})^{T}\mathbb{M}^{(2)}\vec{\alpha}\geq
0,\quad(\vec{\alpha})^{T}\mathbb{M}^{(3)}\vec{\alpha}\geq 0,$
$\displaystyle(\vec{\alpha})^{T}\mathbb{M}^{(4)}\vec{\alpha}\geq
0,\quad(\vec{\alpha})^{T}\mathbb{M}^{(5)}\vec{\alpha}\geq 0,$ (135)

respectively, for any values of the vector $\alpha$. For reasons we do not
fully understand, the first constraint
$(\vec{\alpha})^{T}\mathbb{M}^{(1)}\vec{\alpha}\geq 0$ leads to poor numerical
convergence and is discarded. The remaining four constraints are rewritten as
requiring the matrices to be positive semidefinite using the notation Eq.(39),

$\mathbb{M}^{(2)}\succcurlyeq 0,\quad\mathbb{M}^{(3)}\succcurlyeq
0,\quad\mathbb{M}^{(4)}\succcurlyeq 0,\quad\mathbb{M}^{(5)}\succcurlyeq 0\,.$
(136)

For convenience, this can be rephrased as the positive semidefiniteness of a
single matrix which contains the above four matrices as diagonal blocks,

$\mathbb{M}=\begin{pmatrix}\mathbb{M}_{2}&0&0&0\\\ 0&\mathbb{M}_{3}&0&0\\\
0&0&\mathbb{M}_{4}&0\\\ 0&0&0&\mathbb{M}_{5}\\\ \end{pmatrix}\succcurlyeq
0\,.$ (137)

Analogous to the case of one-loop bubble integrals, IBP reduction expresses
$\mathbb{M}$ as a linear combination of the 15 master integrals in Eq. (117)
and (118), each multiplied by a matrix of rational functions in
$p^{2},m_{1}^{2},m_{2}^{2},m_{3}^{2},m_{4}^{2}$. It is necessary to perform
IBP reduction for banana integrals with up to 13 “dots”, i.e. additional
powers of propagators beyond the standard first power, since the positive
polynomial $Q$ in Eqs. (129)-(132) bring 13 powers of $x_{i}$ when $P$ has
degree 6. The four master integrals in Eq. (118) are known analytically in Eq.
(119), and the values of the remaining 11 master integrals are unknown
parameters to be constrained by Eq. (137).

Before presenting numerical results, we also formulate positivity constraints
for the $\epsilon$ expansion of banana integrals. Recall that banana integrals
in $d=2-2\epsilon$, normalized according to Eq. (116), has a Feynman-parameter
representation Eq. (125) using redefined Feynman parameters in Eq. (124). For
any integers $a_{i}\geq 1$, Taylor-expanding both sides of Eq. (125) and
equating the coefficients of the $\epsilon^{k}$ term for any integer $k$, we
have

$\displaystyle\quad\left[\frac{\Gamma(4-3d/2)}{\Gamma(a-3d/2)}\Gamma(a_{1})\Gamma(a_{2})\Gamma(a_{3})\Gamma(a_{4})\,\hat{I}_{a_{1},a_{2},a_{3},a_{4}}\right]\Bigg{|}_{\epsilon^{k}}$
(138) $\displaystyle=\int_{x_{i}\geq
0}dx_{1}dx_{2}dx_{3}dx_{4}\,\delta(1-x_{1}-x_{2}-x_{3}-x_{4})\left(\prod_{i=1}^{4}\tilde{x}_{i}^{a_{i}-1}\right)\frac{1}{\mathcal{F}(x_{i})}$
$\displaystyle\quad\times\frac{1}{k!}\log^{k}\frac{\mathcal{U}^{4}(x_{i})}{\mathcal{F}^{3}(x_{i})}\,.$

By integration-by-parts reduction, the LHS of the above equation can be
written as linear combinations of the $\epsilon$ expansions of the 15 master
integrals up to the $\epsilon^{k}$ order, assuming that the coefficients of
the master integrals (from IBP reduction) are finite as $\epsilon\to 0$, which
is the case here. Now we are ready to write down positivity constraints for
the $\epsilon$ expansion by extending Eq. (127),

$\displaystyle 0$ $\displaystyle\leq\int_{x_{i}\geq
0}dx_{1}dx_{2}dx_{3}dx_{4}\,\delta(1-x_{1}-x_{2}-x_{3}-x_{4})Q\left(\tilde{x}_{i},\log\frac{\mathcal{U}^{4}(x_{i})}{\mathcal{F}^{3}(x_{i})}\right)\frac{\mathcal{U}(x_{i})^{4-2d}}{\mathcal{F}(x_{i})^{4-3d/2}}\,,$
(139)

where $Q$ is now a positive polynomial in its two arguments above. To build
the most general form of $Q$, we follow Section 2.4.2 and use the building
block

$\log\max\frac{\mathcal{U}^{4}}{\mathcal{F}^{3}}-\log\frac{\mathcal{U}^{4}(x_{i})}{\mathcal{F}^{3}(x_{i})}\geq
0\,.$ (140)

The value of $\max(\mathcal{U}^{4}/\mathcal{F}^{3})$ will be found numerically
once $p^{2}$ and $m_{i}^{2}$ parameters are specified in Eq. (143), in the
next subsection on numerical results. However, we find no minimum of
$\log(\mathcal{U}^{4}/\mathcal{F}^{3})$ at the same parameter values, as
$\mathcal{U}^{4}/\mathcal{F}^{3}$ can become arbitrarily close to zero (from
above) in the range of integration. We will use the following choices of $Q$,

$\displaystyle
Q\left(\tilde{x}_{i},\log\frac{\mathcal{U}^{4}(x_{i})}{\mathcal{F}^{3}(x_{i})}\right)=\tilde{x}_{k}P^{2}\left(\tilde{x}_{i},\log\frac{\mathcal{U}^{4}(x_{i})}{\mathcal{F}^{3}(x_{i})}\right),$
(141) or $\displaystyle
Q\left(\tilde{x}_{i},\log\frac{\mathcal{U}^{4}(x_{i})}{\mathcal{F}^{3}(x_{i})}\right)=\tilde{x}_{k}\left(\log\max\frac{\mathcal{U}^{4}}{\mathcal{F}^{3}}-\log\frac{\mathcal{U}^{4}(x_{i})}{\mathcal{F}^{3}(x_{i})}\right)P^{2}\left(\tilde{x}_{i},\log\frac{\mathcal{U}^{4}(x_{i})}{\mathcal{F}^{3}(x_{i})}\right)\,,$
(142)

where $k$ can be 1, 2, 3, or 4, and $P$ is an arbitrary polynomial with a
maximum total degree $N_{1}$ for the four $\tilde{x}_{i}$ variables and
maximum degree $N_{2}$ for $\log(\mathcal{U}^{4}/\mathcal{F}^{3})$. Similar to
the bubble integral case in Section 2.4, to constrain the
$\mathcal{O}(\epsilon)$ part of the master integrals, we will use Eq. (142)
with $N_{2}=0$, and to constrain the $\mathcal{O}(\epsilon^{2})$ part, we will
use Eq. (141) $N_{2}=1$. For both $\mathcal{O}(\epsilon^{1})$ and
$\mathcal{O}(\epsilon^{2})$ parts, $N_{1}$ will be chosen to be the same as
the cutoff degree used for the $\mathcal{O}(\epsilon^{0})$ calculation.

###  3.3 Numerical results

We present numerical results for the 11 nontrivial master integrals of the
banana family in Eq. (117) at the following numerical values for kinematic
variables,

$p^{2}=2,\quad m_{1}^{2}=2,\quad m_{2}^{2}=3/2,\quad m_{3}^{2}=4/3,\quad
m_{4}^{2}=1\,.$ (143)

We remind readers that the spacetime dimension is set to

$d=2-2\epsilon\,.$ (144)

As this paper is aimed at illustrating a new method, we have chosen example
integrals that are known to high precision in the existing literature. For
three-loop banana integrals, high precision results from series solutions to
differential equations are available from the DiffExp package Hidding:2020ytt
. In fact, we have chosen the same values for the masses in Eq. (143) as the
example in the aforementioned paper, though we chose a different value of
$p^{2}$ as we restrict to the Euclidean region Eq. (120). DiffExp is used to
compute the master integrals to a precision of about $10^{-54}$, which can be
taken as exact values for the purpose of validating our numerical results.

For the values of the master integrals at $\mathcal{O}(\epsilon^{0})$, i.e. in
exactly $d=2$, we will use cutoff degrees of up to 6 in Eq. (134). Unlike the
one-loop case in Section 2, we will not fully characterize the allowed
parameter region which is a sub-region of an 11-dimensional parameter space
and cannot be described by just a lower bound and an upper bound. Instead, we
will only compute the central values for the 11 undetermined master integrals.
Following the prescription laid out in Section 2, the central values are
defined to maximize the lowest eigenvalue of the matrix $\mathbb{M}$ in Eq.
(137). With the largest cutoff degree 6, there are ${{6+4}\choose 4}=210$ free
parameters in $\vec{\alpha}$. Therefore, each of the four diagonal blocks in
Eq. (137) has size $210\times 210$, and the full matrix $\mathbb{M}$ has size
$840\times 840$. We use SDPA-QD as the semidefinite programming solver working
at quadruple-double precision. The solver is able to take advantage of the
block diagonal structure of the matrix $\mathbb{M}$ to improve efficiency.

In Fig. 12, we plot the relative errors of the central values of three
representative master integrals against the cutoff degree, for the
$\mathcal{O}(\epsilon^{0})$ order only. The actual results are given later in
Eq. (147) together with further terms in the $\epsilon$ expansion.

Figure 12: Relative errors of the central values of three representative
master integrals of the banana family, versus the cutoff degree in the
calculation.

The vertical axis of the plot is on a logarithmic scale, and we can see that
the results converge rapidly, in a apparently exponential fashion, as the
cutoff degree is raised. With the largest cutoff degree 6, each of the 11
master integrals is evaluated to an accuracy of at least $10^{-9}$.

For the values of master integrals at $\mathcal{O}(\epsilon^{1})$. We will
need the result for the parameter values Eq. (143),

$\max\left(\mathcal{U}^{4}/\mathcal{F}^{3}\right)\approx 5000/229059\,,$ (145)

which we found by numerical maximization in Mathematica. To be conservative,
we have slighted rounded up the numerical result to a larger nearby rational
number to guarantee that the inequality Eq. (140) is true when using Eq.
(145). This maximum value occurs at

$x_{1}\approx 0.12222,\quad x_{2}\approx 0.22592,\quad x_{3}\approx
0.26701,\quad x_{4}\approx 0.38485\,.$ (146)

While Eq. (145) will be used directly in calculations, Eq. (146) is only
included for completeness. The normalization in Eq. (146) does not matter
since $\mathcal{U}^{4}/\mathcal{F}^{3}$ is invariant under the scaling
transformation Eq. (60).

Then the calculation is similar to the calculation of one-loop bubble interals
to $\mathcal{O}(\epsilon^{1})$ in Section 2.4.3. We use the positivity
constraint Eq. (139) with Eq. (142) for the positive polynomial $Q$, taking
the values $k=1,2,3,4$ and combining the constraints from the four different
choices. For the $P$ polynomial in Eq. (142), we use a maximum total degree
$N_{1}=6$ for the $\tilde{x}_{i}$ variables and a maximum degree of $0$ for
the logarithm, i.e. dropping any terms involving the logarithm. The logarithm
still appears in the bracket preceding $P^{2}$ in Eq. (142) and contributes to
$O(\epsilon)$ parts of integrals by Eq. (138). Using the
$\mathcal{O}(\epsilon^{0})$ results as known inputs, we again solve a
semidefinite programming problem involving a $840\times 840$ matrix with four
diagonal blocks, each of size $210\ \times 210$, to obtain values for the
$\mathcal{O}(\epsilon^{1})$ terms of the master integrals.

For the values of master integrals at $\mathcal{O}(\epsilon^{2})$, we use the
positivity constraint Eq. (139) with Eq. (141) for the positive polynomial
$Q$. We again use a maximum total degree $N_{1}=6$ for the $\tilde{x}_{i}$
variables but now uses a maximum degree of $1$ for the logarithm. Since the
logarithm can appear in a monomial in $P$ with either power 0 or power 1, the
size of the matrix in the semidefinite programming problem is doubled to
$1680\times 1680$, with four diagonal blocks each of size $420\times 420$.
Taking both $\mathcal{O}(\epsilon^{0})$ and $\mathcal{O}(\epsilon^{1})$
results as known inputs, we run SDPA-QD to find the central values for the
$\mathcal{O}(\epsilon^{2})$ results. For brevity, we show results for 3
representative master integrals out of the 11 top-level master integrals, with
kinematic variables taking values of Eq. (143),

$\displaystyle\hat{I}_{1122}$ $\displaystyle\approx
0.31328353052153-0.12137516161264\epsilon-1.5577062442336\epsilon^{2},$ (147)
$\displaystyle\hat{I}_{1112}$ $\displaystyle\approx
1.3758733318476-3.5451169250640\epsilon+0.61363537070259\epsilon^{2},$
$\displaystyle\hat{I}_{1111}$ $\displaystyle\approx
5.9437542439912-33.914772364319\epsilon+106.87640125797\epsilon^{2}\,.$

For documenting the computational outputs, we have kept each number to 14
significant figures, even their actual accuracies are lower as shown in plots
in this section.

We have also calculated both $\mathcal{O}(\epsilon^{1})$ and
$\mathcal{O}(\epsilon^{2})$ results using numerical differentiation of
integrals evaluated at fixed values of dimensions, following the same strategy
of Section 2.5 for one-loop bubble integrals. The calculations are identical
to the $d=2$, i.e. $\mathcal{O}(\epsilon^{0})$ case and are based on Eq. (127)
without any Taylor expansion in $\epsilon$, with the only change being that
$\epsilon$ is set to small numerical values different from 0, i.e. $d$ is set
to numerical values that slightly deviate from $2$. The 4th-order numerical
differentiation formulas, Eqs. (109) and Eq. (110) are applied with
$\epsilon_{0}=0$ and $\Delta\epsilon=10^{-3}$, with the spacetime dimension
$d=2-2\epsilon$. Example results from this alternative method are, again
keeping each number of 14 significant figures,

$\displaystyle\hat{I}_{1122}$ $\displaystyle\approx
0.31328353052153-0.12137519105424\epsilon-1.5576503067221\epsilon^{2},$ (148)
$\displaystyle\hat{I}_{1112}$ $\displaystyle\approx
1.3758733318476-3.5451170400199\epsilon+0.61369255775305\epsilon^{2},$
$\displaystyle\hat{I}_{1111}$ $\displaystyle\approx
5.9437542439912-33.914771261794\epsilon+106.87318272740\epsilon^{2}\,.$

Note that the $\mathcal{O}(\epsilon^{0})$ results are copied from Eq. (147) as
they are not re-calculated.

Figure 13: Log-scale plot of relative errors of numerical results for 11
nontrivial master integrals of the three-loop banana family, Eq. (117), with
the normalization Eq. (116), up to second order in the dimension
regularization parameter $\epsilon$. The four numbers under each bar indicates
the subscript indices of $\hat{I}_{a_{1},a_{2},a_{3},a_{4}}$, with commas
omitted as all four indices are single-digit numbers (either 1 or 2). For each
master integral, the five vertical bars, from left to right, show the relative
errors for the $\epsilon^{0}$ term, the $\epsilon^{1}$ term calculated from
direct positivity constraints (abbreviated as cons in the legend), the
$\epsilon^{1}$ term calculated from numerical differentiation (abbreviated as
diff in the legend), the $\epsilon^{2}$ term calculated from direct positivity
constraints, and the $\epsilon^{2}$ term calculated from numerical
differentiation.

These numerical results for the $\epsilon$ expansion of master integrals Eq.
(117) are obtained with the normalization of Eq. (116). The reference results
from DiffExp are have the normalization of Eq. (115) and additional factors
for individual master integrals. The reference results have been converted to
use our normalizations for comparison. DiffExp results for the three sample
integrals, truncated to 14 significant digits, are

$\displaystyle\hat{I}_{1122}$ $\displaystyle\approx
0.31328353056677-0.121375191032390\epsilon-1.5577067713048\epsilon^{2},$ (149)
$\displaystyle\hat{I}_{1112}$ $\displaystyle\approx
1.37587333189510-3.5451170391547\epsilon+0.61363351857945\epsilon^{2},$
$\displaystyle\hat{I}_{1111}$ $\displaystyle\approx
5.9437542414259-33.914771263107\epsilon+106.876390717227\epsilon^{2}\,.$

The final numerical accuracy for the 11 master integrals, with values of
kinematic parameters chosen in Eq. (143), is shown in Fig. 13. The $\epsilon$
expansion results from “direct positivity constraints” are labeled cons and
results from numerical differentiation are labeled diff in the plot legend. We
can see that numerical differentiation gives very good accuracy for
$\mathcal{O}(\epsilon^{1})$ terms, comparable with the accuracy of
$\mathcal{O}(\epsilon^{0})$ terms, while for the $\mathcal{(}\epsilon^{2})$
terms, direct positivity constraints yield more accurate results. In any case,
both methods for the $\epsilon$ expansion have demonstrated their potentials
in this initial investigation, as all results for $\mathcal{O}(\epsilon^{1})$
terms have relative errors below $10^{-6}$ and all results for
$\mathcal{O}(\epsilon^{2})$ terms have relative errors below $10^{-3}$.

We briefly comment on computational resources used. IBP reduction takes a few
CPU-hours with FIRE6 Smirnov:2019qkx  with numerical kinematics Eq. (143). The
IBP reduction results are obtained with analytic dependence on $d$ and can be
subsequently expanded in $\epsilon$, so no extra IBP reduction is needed for
obtaining the $\epsilon$ expansions of master integrals beyond the zeroth
order. Running the semidefinite programming solver SDPA-QD takes a few CPU-
hours for every run, including one run for solving positivity constraints for
the $\mathcal{O}(\epsilon^{i})$ part for each $i=0,1,2$, and for the
alternative method based on numerical differentiation, several runs at
different numerical values of $\epsilon$ to generate the data needed to feed
into finite-difference approximations.

##  4 Discussions

We have demonstrated a new method for evaluating Feynman integrals which, to
our best knowledge, is the first method based on inequality constraints, while
previously exploited consistency conditions for Feynman integrals are based on
equality constraints (such as the vanishing of the coefficient of a certain
spurious singularity). Our calculation strategy is writing down an infinite
class of convergent integrals with non-negative integrands and reducing them
to linear sums of a set of master integrals. This constraints an infinite
number of linear sums of master integrals to be non-negative. A truncated set
of the constraints can be solved as a semidefinite programming problem in
mathematical optimization. Surprisingly, the constraints appear strong enough
to determine the integrals to any desired precision since the bounds appear to
converge exponentially as the truncation cutoff is increased. Like the method
of differential equations Kotikov:1990kg ; Bern:1993kr ; Remiddi:1997ny ;
Gehrmann:1999as ; Henn:2013pwa , our method relies on integration-by-parts
(IBP) identities, but instead of using differential equations to transport the
values of the integrals across kinematic space, we only use IBP information at
a single point in kinematic space.

Though our study is preliminary, the numerical results are promising. We have
demonstrated the applicability of our methods to a nontrivial example, namely
three-loop banana integrals with four unequal internal masses in
$d=2-2\epsilon$ dimensions. With modest computational resources, we evaluated
the $\mathcal{O}(\epsilon^{0})$ part of all the 11 nontrivial master integrals
to a relative accuracy of at least $10^{-9}$. The accuracies for
$\mathcal{O}(\epsilon^{1})$ and $\mathcal{O}(\epsilon^{2})$ terms are lower,
though only slightly so for $\mathcal{O}(\epsilon^{1})$ terms when the
numerical differentiation method is used. For all but the smallest problems,
extended-precision floating point arithmetic is needed to ensure numerical
stability in the semidefinite programming solver, similar to what was
encountered in the conformal bootstrap Simmons-Duffin:2015qma  and the quantum
mechanics bootstrap Berenstein:2022unr . We note that extended precision is
also generally need in the evaluation of Feynman integrals by series solutions
of differential equations as observed in e.g. Refs. Moriello:2019yhu ;
Hidding:2020ytt .

We have also revealed hidden consistency relations that link different terms
in the $\epsilon$ expansions of Feynman integrals. As explained in Section
2.4.1, for any (quasi-) finite Feynman integral without numerators, the
$\epsilon$ expansion terms (appropriately normalized) must give rise to a
positive-semidefinite Hankel matrix. This is an extremely general statement
which can be checked against a huge number of Feynman integral computations in
the literature, because many Feynman integrals have Euclidean regions and it
is believed that a quasi-finite basis exist for any family of integrals
vonManteuffel:2014qoa . This result is an elementary consequence of our
analysis but has not been previously exposed in the literature. Such
constraints have been solved numerically in our paper to predict the
$\epsilon$ expansion terms to high accuracy. We have also formulated an
alternative method to obtain the $\epsilon$ expansion terms by numerical
differentiation of semidefinite programming solutions with respect to the
spacetime dimension. The above two methods for calculating $\epsilon$
expansion terms are complementary and we have found cases in which either of
them outperforms the other in accuracy.

Our new method for calculating Feynman integrals is analogous to recent
developments in bootstrapping quantum mechanics systems and lattice models
Lin:2020mme ; Han:2020bkb ; Berenstein:2021dyf ; Berenstein:2022unr ;
Anderson:2016rcw ; Kazakov:2022xuh ; Cho:2022lcj . For example, the role of
IBP and dimensional-shifting identities in our work is analogous to the role
of moment recursion relations in the quantum mechanics bootstrap. Analogous
identities also appear in EFT bounds as “null constraints” from crossing
symmetry Caron-Huot:2020cmc .333We thank Francesco Riva for pointing out this
connection. While our work has imported techniques developed in non-
perturbative contexts to perturbative physics, in the reverse direction, the
differential equation method of perturbative calculations has been applied to
non-perturbative lattice correlation functions in Refs. Gasparotto:2022mmp ;
Weinzierl:2020nhw ; Cacciatori:2022mbi , also exploiting identities similar to
those from IBP. Therefore, we expect a fruitful exchange of techniques between
perturbative and non-perturbative calculations.

Finally, we speculate on possible future work. Except for the generic
constraints on the $\epsilon$ expansion, this paper has mainly treated massive
Feynman integrals, and for integral families involving massless internal
lines, it would be necessary to identify non-negative integrals free of not
only ultraviolet but also infrared divergences to generate the positivity
constraints. To extend our method to integrals outside the Euclidean region,
it remains to be seen how positivity constrains can be formulated, possibly
for real and imaginary parts separately after a suitable deformation of the
integration contour. Connections with other notions of positivity relevant for
Feynman integrals Abreu:2019wzk ; Yang:2022gko  remain to be explored. Another
interesting question is whether positivity constraints can be used to
understand complete scattering amplitudes (rather than individual Feynman
integrals) at a fixed order in perturbation theory, in light of numerical
observations in $\mathcal{N}=4$ super-Yang-Mills theory amplitudes in Ref.
Dixon:2016apl .

## Acknowledgements

M.Z.’s work is supported in part by the U.K. Royal Society through Grant
URF\R1\20109. For the purpose of open access, the author has applied a
Creative Commons Attribution (CC BY) license to any Author Accepted Manuscript
version arising from this submission.

##  Appendix A Analytic results for bubble integrals

The material here is well known in the literature but included for
completeness. We work in spacetime dimension $d=4-2\epsilon$. Using Feynman
parametrization, Eq. (1) is written as, for $a_{1}\geq 1,a_{2}\geq 1$.

$I^{d}_{a_{1},a_{2}}=\frac{\Gamma(a_{1}+a_{2}-d/2)e^{\gamma_{E}\epsilon}}{\Gamma(a_{1})\Gamma(a_{2})}\int_{0}^{1}dx\,x^{a_{1}-1}(1-x)^{a_{2}-1}\left[m^{2}-p^{2}x(1-x)-i0^{+}\right]^{d/2-a_{1}-a_{2}}\,.$
(150)

This allows us to evaluate the first of the master integrals in Eq. (9) as

$I^{d=4-2\epsilon}_{2,1}=\Gamma(1+\epsilon)(m^{2})^{-1-\epsilon}e^{\gamma_{E}\epsilon}\int_{0}^{1}dx\,\frac{x}{\left[1-x(1-x)p^{2}/m^{2}-i0^{+}\right]^{1+\epsilon}}$
(151)

Since the denominator in the above integrand is symmetric under
$x\leftrightarrow 1-x$, we can symmetrize the numerator as
$x\rightarrow[x+(1-x)]/2=1/2$, obtaining

$I^{d=4-2\epsilon}_{2,1}=\Gamma(1+\epsilon)(m^{2})^{-1-\epsilon}e^{\gamma_{E}\epsilon}\int_{0}^{1}dx\,\frac{1}{2\left[1-x(1-x)p^{2}/m^{2}-i0^{+}\right]^{1+\epsilon}}$
(152)

For $p^{2}<0$, this evaluates to

$I^{d=4-2\epsilon}_{2,1}=\Gamma(1+\epsilon)(-p^{2}+i0^{+})^{-1-\epsilon}e^{\gamma_{E}\epsilon}\left[\frac{1}{\beta}\log\frac{\beta+1}{\beta-1}+\mathcal{O}(\epsilon)\right]\,,$
(153)

where we defined

$\beta=\sqrt{1-\frac{4m^{2}}{p^{2}}-i0^{+}}\,.$ (154)

The omitted $\mathcal{O}(\epsilon)$ term in Eq. (152) is easy to evaluate and
we do not include the explicit result here. The second master integral in Eq.
(9) is a trivial tadpole integral,

$I^{d=4-2\epsilon}_{3,0}=\frac{1}{2}\Gamma(3-d/2)e^{\gamma_{E}(4-d)/2}(m^{2})^{-3+d/2}=\frac{1}{2}\Gamma(1+\epsilon)e^{\gamma_{E}\epsilon}(m^{2})^{-1-\epsilon}\,.$
(155)

So the ratio defined in Eq. (13) is equal to, when $d=4-2\epsilon$,

$\displaystyle\hat{I}^{d=4-2\epsilon}_{2,1}$
$\displaystyle=I^{d=4-2\epsilon}_{2,1}/I^{d=4-2\epsilon}_{3,0}=2\int_{0}^{1}dx\,\frac{x}{\left[1-x(1-x)p^{2}/m^{2}-i0^{+}\right]^{1+\epsilon}}$
$\displaystyle=-\frac{2m^{2}}{p^{2}}\frac{1}{\beta}\log\frac{\beta+1}{\beta-1}+\mathcal{O}(\epsilon)\,.$
(156)

This result is real and positive for any $p^{2}<4m^{2}$, after taking into
cancellation of imaginary parts with the definition Eq. (154). The corrections
at higher order in the dimensional regularization parameter $\epsilon$ are
also readily obtained but we do not include the results here.

The tadpole integral with an arbitrary power for the propagator is well known,

$\int\frac{d^{d}l\,e^{\gamma_{E}\epsilon}}{i\pi^{d/2}}\frac{1}{(-l^{2}+m^{2})^{n}}=\frac{1}{2}e^{\gamma_{E}(4-d)/2}\frac{\Gamma(n-d/2)}{\Gamma(n)}\frac{1}{(m^{2})^{n-d/2}}\,.$
(157)

## References

  * (1) T. Binoth and G. Heinrich, _An automatized algorithm to compute infrared divergent multiloop integrals_ , _Nucl. Phys. B_ 585 (2000) 741 [hep-ph/0004013].
  * (2) C. Bogner and S. Weinzierl, _Resolution of singularities for multi-loop integrals_ , _Comput. Phys. Commun._ 178 (2008) 596 [0709.4092].
  * (3) T. Kaneko and T. Ueda, _A Geometric method of sector decomposition_ , _Comput. Phys. Commun._ 181 (2010) 1352 [0908.2897].
  * (4) S. Borowka, G. Heinrich, S. P. Jones, M. Kerner, J. Schlenk and T. Zirke, _SecDec-3.0: numerical evaluation of multi-scale integrals beyond one loop_ , _Comput. Phys. Commun._ 196 (2015) 470 [1502.06595].
  * (5) S. Borowka, G. Heinrich, S. Jahn, S. P. Jones, M. Kerner, J. Schlenk et al., _pySecDec: a toolbox for the numerical evaluation of multi-scale integrals_ , _Comput. Phys. Commun._ 222 (2018) 313 [1703.09692].
  * (6) A. V. Smirnov, _FIESTA 3: cluster-parallelizable multiloop numerical calculations in physical regions_ , _Comput. Phys. Commun._ 185 (2014) 2090 [1312.3186].
  * (7) A. V. Smirnov, _FIESTA4: Optimized Feynman integral calculations with GPU support_ , _Comput. Phys. Commun._ 204 (2016) 189 [1511.03614].
  * (8) F. Moriello, _Generalised power series expansions for the elliptic planar families of Higgs + jet production at two loops_ , _JHEP_ 01 (2020) 150 [1907.13234].
  * (9) M. Hidding, _DiffExp, a Mathematica package for computing Feynman integrals in terms of one-dimensional series expansions_ , _Comput. Phys. Commun._ 269 (2021) 108125 [2006.05510].
  * (10) 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 [1711.09572].
  * (11) X. Liu and Y.-Q. Ma, _Multiloop corrections for collider processes using auxiliary mass flow_ , _Phys. Rev. D_ 105 (2022) L051503 [2107.01864].
  * (12) X. Liu and Y.-Q. Ma, _AMFlow: A Mathematica package for Feynman integrals computation via auxiliary mass flow_ , _Comput. Phys. Commun._ 283 (2023) 108565 [2201.11669].
  * (13) T. Armadillo, R. Bonciani, S. Devoto, N. Rana and A. Vicini, _Evaluation of Feynman integrals with arbitrary complex masses via series expansions_ , _Comput. Phys. Commun._ 282 (2023) 108545 [2205.03345].
  * (14) M. Hidding and J. Usovitsch, _Feynman parameter integration through differential equations_ , 2206.14790.
  * (15) N. I. Usyukina and A. I. Davydychev, _An Approach to the evaluation of three and four point ladder diagrams_ , _Phys. Lett. B_ 298 (1993) 363.
  * (16) N. I. Usyukina and A. I. Davydychev, _Exact results for three and four point ladder diagrams with an arbitrary number of rungs_ , _Phys. Lett. B_ 305 (1993) 136.
  * (17) V. A. Smirnov, _Analytical result for dimensionally regularized massless on shell double box_ , _Phys. Lett. B_ 460 (1999) 397 [hep-ph/9905323].
  * (18) J. B. Tausk, _Nonplanar massless two loop Feynman diagrams with four on-shell legs_ , _Phys. Lett. B_ 469 (1999) 225 [hep-ph/9909506].
  * (19) M. Czakon, _Automatized analytic continuation of Mellin-Barnes integrals_ , _Comput. Phys. Commun._ 175 (2006) 559 [hep-ph/0511200].
  * (20) A. V. Smirnov and V. A. Smirnov, _On the Resolution of Singularities of Multiple Mellin-Barnes Integrals_ , _Eur. Phys. J. C_ 62 (2009) 445 [0901.0386].
  * (21) J. Gluza, K. Kajda and T. Riemann, _AMBRE: A Mathematica package for the construction of Mellin-Barnes representations for Feynman integrals_ , _Comput. Phys. Commun._ 177 (2007) 879 [0704.2423].
  * (22) A. V. Belitsky, A. V. Smirnov and V. A. Smirnov, _MB tools reloaded_ , _Nucl. Phys. B_ 986 (2023) 116067 [2211.00009].
  * (23) S. Laporta, _High precision calculation of multiloop Feynman integrals by difference equations_ , _Int. J. Mod. Phys. A_ 15 (2000) 5087 [hep-ph/0102033].
  * (24) R. N. Lee and V. A. Smirnov, _The Dimensional Recurrence and Analyticity Method for Multicomponent Master Integrals: Using Unitarity Cuts to Construct Homogeneous Solutions_ , _JHEP_ 12 (2012) 104 [1209.0339].
  * (25) R. N. Lee and A. F. Pikelner, _Four-loop HQET propagators from the DRA method_ , _JHEP_ 02 (2023) 097 [2211.03668].
  * (26) S. Catani, T. Gleisberg, F. Krauss, G. Rodrigo and J.-C. Winter, _From loops to trees by-passing Feynman’s theorem_ , _JHEP_ 09 (2008) 065 [0804.3170].
  * (27) R. Runkel, Z. Szőr, J. P. Vesga and S. Weinzierl, _Causality and loop-tree duality at higher loops_ , _Phys. Rev. Lett._ 122 (2019) 111603 [1902.02135].
  * (28) Z. Capatti, V. Hirschi, D. Kermanschah and B. Ruijl, _Loop-Tree Duality for Multiloop Numerical Integration_ , _Phys. Rev. Lett._ 123 (2019) 151602 [1906.06138].
  * (29) M. Borinsky, _Tropical Monte Carlo quadrature for Feynman integrals_ , 2008.12310.
  * (30) M. Borinsky, H. J. Munch and F. Tellander, _Tropical Feynman integration in the Minkowski regime_ , 2302.08955.
  * (31) S. Weinzierl, _Feynman Integrals_. 1, 2022, 10.1007/978-3-030-99558-4, [2201.03593].
  * (32) V. A. Smirnov, _Analytic tools for Feynman integrals_ , vol. 250. 2012, 10.1007/978-3-642-34886-0.
  * (33) R. Rattazzi, V. S. Rychkov, E. Tonni and A. Vichi, _Bounding scalar operator dimensions in 4D CFT_ , _JHEP_ 12 (2008) 031 [0807.0004].
  * (34) M. Kruczenski, J. Penedones and B. C. van Rees, _Snowmass White Paper: S-matrix Bootstrap_ , 2203.02421.
  * (35) A. Adams, N. Arkani-Hamed, S. Dubovsky, A. Nicolis and R. Rattazzi, _Causality, analyticity and an IR obstruction to UV completion_ , _JHEP_ 10 (2006) 014 [hep-th/0602178].
  * (36) B. Bellazzini, J. Elias Miró, R. Rattazzi, M. Riembau and F. Riva, _Positive moments for scattering amplitudes_ , _Phys. Rev. D_ 104 (2021) 036006 [2011.00037].
  * (37) S. Caron-Huot and V. Van Duong, _Extremal Effective Field Theories_ , _JHEP_ 05 (2021) 280 [2011.02957].
  * (38) A. J. Tolley, Z.-Y. Wang and S.-Y. Zhou, _New positivity bounds from full crossing symmetry_ , _JHEP_ 05 (2021) 255 [2011.02400].
  * (39) N. Arkani-Hamed, T.-C. Huang and Y.-t. Huang, _The EFT-Hedron_ , _JHEP_ 05 (2021) 259 [2012.15849].
  * (40) S. Caron-Huot, D. Mazac, L. Rastelli and D. Simmons-Duffin, _Sharp boundaries for the swampland_ , _JHEP_ 07 (2021) 110 [2102.08951].
  * (41) S. Caron-Huot, Y.-Z. Li, J. Parra-Martinez and D. Simmons-Duffin, _Causality constraints on corrections to Einstein gravity_ , 2201.06602.
  * (42) Z. Bern, D. Kosmopoulos and A. Zhiboedov, _Gravitational effective field theory islands, low-spin dominance, and the four-graviton amplitude_ , _J. Phys. A_ 54 (2021) 344002 [2103.12728].
  * (43) H. Chen, A. L. Fitzpatrick and D. Karateev, _Nonperturbative bounds on scattering of massive scalar particles in d $\geq$ 2_, _JHEP_ 12 (2022) 092 [2207.12448].
  * (44) B. Bellazzini, G. Isabella and M. M. Riva, _Classical vs Quantum Eikonal Scattering and its Causal Structure_ , 2211.00085.
  * (45) H. W. Lin, _Bootstraps to strings: solving random matrix models with positivity_ , _JHEP_ 06 (2020) 090 [2002.08387].
  * (46) X. Han, S. A. Hartnoll and J. Kruthoff, _Bootstrapping Matrix Quantum Mechanics_ , _Phys. Rev. Lett._ 125 (2020) 041601 [2004.10212].
  * (47) D. Berenstein and G. Hulsey, _Bootstrapping Simple QM Systems_ , 2108.08757.
  * (48) D. Berenstein and G. Hulsey, _A Semidefinite Programming algorithm for the Quantum Mechanical Bootstrap_ , 2209.14332.
  * (49) P. D. Anderson and M. Kruczenski, _Loop Equations and bootstrap methods in the lattice_ , _Nucl. Phys. B_ 921 (2017) 702 [1612.08140].
  * (50) V. Kazakov and Z. Zheng, _Bootstrap for lattice Yang-Mills theory_ , _Phys. Rev. D_ 107 (2023) L051501 [2203.11360].
  * (51) M. Cho, B. Gabai, Y.-H. Lin, V. A. Rodriguez, J. Sandor and X. Yin, _Bootstrapping the Ising Model on the Lattice_ , 2206.12538.
  * (52) 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.
  * (53) L. Vandenberghe and S. Boyd, _Semidefinite programming_ , _SIAM review_ 38 (1996) 49.
  * (54) D. Poland, D. Simmons-Duffin and A. Vichi, _Carving Out the Space of 4D CFTs_ , _JHEP_ 05 (2012) 110 [1109.5176].
  * (55) A. von Manteuffel, E. Panzer and R. M. Schabinger, _A quasi-finite basis for multi-loop Feynman integrals_ , _JHEP_ 02 (2015) 120 [1411.7392].
  * (56) O. V. Tarasov, _Connection between Feynman integrals having different values of the space-time dimension_ , _Phys. Rev. D_ 54 (1996) 6479 [hep-th/9606018].
  * (57) S. Laporta and E. Remiddi, _The Analytical value of the electron (g-2) at order alpha**3 in QED_ , _Phys. Lett. B_ 379 (1996) 283 [hep-ph/9602417].
  * (58) C. Anastasiou and A. Lazopoulos, _Automatic integral reduction for higher order perturbative calculations_ , _JHEP_ 07 (2004) 046 [hep-ph/0404258].
  * (59) A. von Manteuffel and C. Studerus, _Reduze 2 - Distributed Feynman Integral Reduction_ , 1201.4330.
  * (60) A. V. Smirnov, _Algorithm FIRE – Feynman Integral REduction_ , _JHEP_ 10 (2008) 107 [0807.3243].
  * (61) A. V. Smirnov and F. S. Chuharev, _FIRE6: Feynman Integral REduction with Modular Arithmetic_ , _Comput. Phys. Commun._ 247  (2020) 106877 [1901.07808].
  * (62) R. N. Lee, _Presenting LiteRed: a tool for the Loop InTEgrals REDuction_ , 1212.2685.
  * (63) R. N. Lee, _LiteRed 1.4: a powerful tool for reduction of multiloop integrals_ , _J. Phys. Conf. Ser._ 523 (2014) 012059 [1310.1145].
  * (64) P. Maierhöfer, J. Usovitsch and P. Uwer, _Kira—A Feynman integral reduction program_ , _Comput. Phys. Commun._ 230 (2018) 99 [1705.05610].
  * (65) B. Borchers, _Csdp, ac library for semidefinite programming_ , _Optimization methods and Software_ 11 (1999) 613.
  * (66) M. Yamashita, K. Fujisawa, M. Fukuda, K. Kobayashi, K. Nakata and M. Nakata, _Latest developments in the sdpa family for solving large-scale sdps_ , _Handbook on semidefinite, conic and polynomial optimization_ (2012) 687\.
  * (67) D. Simmons-Duffin, _A Semidefinite Program Solver for the Conformal Bootstrap_ , _JHEP_ 06 (2015) 174 [1502.02033].
  * (68) S. Groote, J. G. Korner and A. A. Pivovarov, _On the evaluation of a certain class of Feynman diagrams in x-space: Sunrise-type topologies at any loop order_ , _Annals Phys._ 322 (2007) 2374 [hep-ph/0506286].
  * (69) A. Klemm, C. Nega and R. Safari, _The $l$-loop Banana Amplitude from GKZ Systems and relative Calabi-Yau Periods_, _JHEP_ 04 (2020) 088 [1912.06201].
  * (70) K. Bönisch, F. Fischbach, A. Klemm, C. Nega and R. Safari, _Analytic structure of all loop banana integrals_ , _JHEP_ 05 (2021) 066 [2008.10574].
  * (71) K. Bönisch, C. Duhr, F. Fischbach, A. Klemm and C. Nega, _Feynman integrals in dimensional regularization and extensions of Calabi-Yau motives_ , _JHEP_ 09 (2022) 156 [2108.05310].
  * (72) D. Kreimer, _Bananas: multi-edge graphs and their Feynman integrals_ , 2202.05490.
  * (73) A. Forum and M. von Hippel, _A Symbol and Coaction for Higher-Loop Sunrise Integrals_ , 2209.03922.
  * (74) S. Pögel, X. Wang and S. Weinzierl, _Taming Calabi-Yau Feynman Integrals: The Four-Loop Equal-Mass Banana Integral_ , _Phys. Rev. Lett._ 130 (2023) 101601 [2211.04292].
  * (75) A. V. Kotikov, _Differential equations method: New technique for massive Feynman diagrams calculation_ , _Phys. Lett. B_ 254 (1991) 158.
  * (76) Z. Bern, L. J. Dixon and D. A. Kosower, _Dimensionally regulated pentagon integrals_ , _Nucl. Phys. B_ 412 (1994) 751 [hep-ph/9306240].
  * (77) E. Remiddi, _Differential equations for Feynman graph amplitudes_ , _Nuovo Cim. A_ 110 (1997) 1435 [hep-th/9711188].
  * (78) T. Gehrmann and E. Remiddi, _Differential equations for two loop four point functions_ , _Nucl. Phys. B_ 580 (2000) 485 [hep-ph/9912329].
  * (79) J. M. Henn, _Multiloop integrals in dimensional regularization made simple_ , _Phys. Rev. Lett._ 110 (2013) 251601 [1304.1806].
  * (80) F. Gasparotto, A. Rapakoulias and S. Weinzierl, _Nonperturbative computation of lattice correlation functions by differential equations_ , _Phys. Rev. D_ 107 (2023) 014502 [2210.16052].
  * (81) S. Weinzierl, _Correlation functions on the lattice and twisted cocycles_ , _Phys. Lett. B_ 805 (2020) 135449 [2003.05839].
  * (82) S. L. Cacciatori and P. Mastrolia, _Intersection Numbers in Quantum Mechanics and Field Theory_ , 2211.03729.
  * (83) S. Abreu, R. Britto, C. Duhr, E. Gardi and J. Matthew, _From positive geometries to a coaction on hypergeometric functions_ , _JHEP_ 02 (2020) 122 [1910.08358].
  * (84) Q. Yang, _Schubert problems, positivity and symbol letters_ , _JHEP_ 08 (2022) 168 [2203.16112].
  * (85) L. J. Dixon, M. von Hippel, A. J. McLeod and J. Trnka, _Multi-loop positivity of the planar $\mathcal{N}$ = 4 SYM six-point amplitude_, _JHEP_ 02 (2017) 112 [1611.08325].