title: Tropical Monte Carlo quadrature for Feynman integrals
authors: Michael Borinsky
text: 
  
  
  
  
  Tropical Monte Carlo quadrature for Feynman integrals
      Michael BorinskyNikhef, Science Park 105, Amsterdam 1098 XG, The Netherlands — Preprint nr. 2020-027
  
      May 2020
  =========================================================================================================
  
  
  
   
  We introduce a new method to evaluate algebraic integrals over the simplex numerically. This new approach employs techniques from tropical geometry and exceeds the capabilities of existing numerical methods by an order of magnitude.  The method can be improved further by exploiting the geometric structure of the underlying integrand. As an illustration of this, we give a specialized integration algorithm for a class of integrands that exhibit the form of a generalized permutahedron. This class includes integrands for scattering amplitudes and parametric Feynman integrals with tame kinematics. A proof-of-concept implementation is provided with which Feynman integrals up to loop order 17 can be evaluated.
  
  
  
  
  
  
  
  § INTRODUCTION
  
  
  
   §.§ Motivation
  
  
  Feynman integrals are ubiquitous in various branches of theoretical physics. They are hard to evaluate and predictions for particle physics experiments rely heavily on them. Their evaluation even poses a bottleneck for the analysis of the data from some high accuracy experiments <cit.>. This situation has fostered the development of extremely sophisticated and specialized technologies aimed to obtain a manageable analytic expression for a given Feynman integral. The state of the art technique is the differential equation method <cit.>. A slightly less powerful method, which is amendable to more general algebraic integrals, is systematic algebraic integration <cit.>. See also <cit.> for an overview on other methods.
  
  The rapid development of these technologies in the last decades has been driven largely by new deep insights into the underlying mathematical structures. 
  Recent advances in the differential equation method were inspired by the simple analytic expressions which can be obtained in supersymmetric quantum field theories via generalized unitarity and recursion relations <cit.>.
  A program to study the arithmetic properties of parametric Feynman integrals <cit.> led to the development of systematic algebraic integration algorithms. This arithmetic understanding of the relevant function classes was also one of the driving forces of the differential equation method <cit.> and is still driving new developments in the especially challenging elliptic regime <cit.>.
  
  All these technologies aim to obtain a closed form analytic expression for the Feynman integral and they all fail once the underlying graph and the associated physical parameters exceed a certain complexity. In these cases a numerical approach is the only way to proceed <cit.>. 
  
  
  
  The most established numerical approach to tackle such integrals is sector decomposition. Sector decomposition as a tool for numerical evaluation of Feynman integrals has been developed by Binoth and Heinrich <cit.>. It was subsequently improved by Bogner and Weinzierl <cit.>. Another conceptual innovation of the overall method was achieved by Kaneko and Ueda <cit.> who brought sector decomposition on a geometric footing. Today, geometric sector decomposition is still the most powerful method for the numerical evaluation of Feynman integrals. It lies at the heart of two popular software tools <cit.>. 
  Another promising numerical technique for Feynman integration is loop-tree duality <cit.>, which is in an active development phase (see for instance <cit.> and the references therein).
  
  In contrast to analytic evaluation methods the mathematical structures exhibited by Feynman integrals are an essentially untapped resource in the context of numerical evaluation. Most numerical techniques are completely oblivious to the rich specific structure of the integrals as they are designed to be applicable to arbitrary algebraic integrands.
  
  
  For this reason, the major objective of this paper is to use some of these mathematical structures to improve the numerical evaluation of Feynman graphs and to show that these dormant resources can be harnessed.
  The overall endeavour behind this objective consists of making progress towards the following two goals:
  
  The first goal is to make numerical evaluation techniques more applicable to real world phenomenology. There are integrals which contribute to interesting measurable processes, but cannot be calculated analytically with available methods. For these integrals numerical evaluation is currently the only way to make predictions for experiments. 
  Numerical methods naturally come with a caveat: they suffer from long evaluation times or they are limited in accuracy. 
  Feynman integrals usually need to be evaluated a large number of times in a big parameter space. This is not difficult if an analytic expression for the integral is known, which can be evaluated sufficiently fast, but poses a tough challenge for numerical methods which sacrifice evaluation speed and accuracy for generality. The task for this goal is therefore to increase the performance of numerical methods. 
  
  
  The second goal is to obtain reliable data in the large-order regime where analytic methods hopelessly fail. There are many indications that the large-order behaviour of perturbation theory is deeply intertwined with non-perturbative phenomena <cit.>. The analysis of the large-order behaviour of perturbation theory in quantum mechanics by Bender and Wu <cit.> has sparked an extremely fruitful branch of research in theoretical and mathematical physics <cit.>. Non-perturbative analytic calculations in quantum field theory are plagued with various gaps in our understanding of the underlying mathematics <cit.>. A repetition of an explicit Bender-Wu like numerical analysis for perturbative quantum field theory is very desirable as it would shed some light into a highly unexplored territory. Unfortunately, this is extremely challenging as the evaluation of large numbers of Feynman integrals of order ∼ 100 would be necessary. 
  It is hopeless to approach this task using the naive method of evaluation Feynman integrals one by one. New methods with which whole classes of diagrams can be evaluated at once need to be developed. The growing understanding of the geometry of amplitude integrals could lead the way to such methods. It is also necessary that these methods are computationally efficient: The demand for computing resources shall at most depend polynomially on the size of the problem (i.e. the respective order in the perturbative expansion).
  
  
  
  
  This paper achieves some progress towards both these goals.  The strategy is to employ tropical geometry <cit.> for numerical quadrature. Panzer <cit.> recently showed that a tropical version of a Feynman graph's period behaves similar to the period itself and anticipated that this tropical version may be used for explicit numerical evaluation. Tropical geometry has also recently been used in the context of string theory and scattering amplitudes <cit.>. 
  
  We will introduce a new Monte Carlo algorithm with which the numerical evaluation of Feynman integrals can be significantly accelerated. It is based on the established geometric sector decomposition principle and can be applied to general algebraic integrals such as the one below in eq. (<ref>). The improvement comes from a stratified sampling approach to Monte Carlo quadrature, which is driven by the (tropical) geometric structure of the algebraic integrand. We will call this method tropical sampling. Tropical sampling effectively decouples the complexity of the underlying integral from the achievable accuracy within the Monte Carlo approach. 
  
  Even though this new tropical-geometric technique offers a significant improvement over the traditional procedure already for general algebraic integrals, a lot more can be achieved if further information on the geometric structures of the integral is used. 
  
  As an example of this, we will give a specialized algorithm for cases where the integrand at hand exhibits the form of a generalized permutahedron <cit.> in a certain sense which will be defined later. Many integrals in quantum field theory and string theory fall under this category. For instance, integrands for complete amplitudes in various theories <cit.> and Feynman integrals with generic Euclidean kinematics are of this kind <cit.>. 
  
  A proof-of-concept implementation of the resulting algorithm is provided. With this implementation high dimensional Feynman integrals can be numerically integrated using widely available hardware. High dimensional explicitly means that integrals corresponding to Feynman graphs with around 20 edges can be estimated up to 10^-3 relative accuracy in a couple of CPU-seconds and integrals for graphs with up to 30 edges in about half an hour. Ultimately, this approach is not CPU but memory constrained when the complexity increases. For instance, for a graph with 30 edges already 16  GB of computer memory are required to run the algorithm. To integrate an 18 loop φ^4-theory four-point graph with the implementation 1  TB of memory would be necessary. 
  
  Although both algorithms are not efficient in the strong sense, as exponential runtime and memory requirements start to dominate at some point, there is hope for the existence of an algorithm that evaluates a Feynman graph of loop order n up to a given accuracy with runtime and memory demands bounded by a polynomial in n (see Section <ref> (1)).
  
  
  
   §.§ Algebraic integrals over the simplex
  
  
  The central object of study in this article is the integral
  
      I = ∫__>0^n-1∏_i a_i(x)^ν_i/∏_j b_j(x)^ρ_jΩ   with
  
  
  
    * the positive orthant of real projective space as integration domain,
  
      _>0^n-1 = { x= [x_1:…:x_n] ∈^n-1(): x_k > 0  for all  k =1,…,n},
  
  
    * the differential form 
  
      Ω = ∑_k=1^n (-1)^n-k x_1/x_1∧…∧ x_k/x_k∧…∧ x_n/x_n,
  
  
    * the sets of homogeneous polynomials {a_1,a_2,…},{b_1,b_2,…}⊂[x_1, …, x_n], s.t. 
  
      ∑_i ν_i  a_i = ∑_j ρ_j  b_j,
  
  
    * where the coefficients ν_i,ρ_j ∈ have non-negative real part: ν_i, ρ_j ≥ 0 and
  
    * a fixed branch choice for the each of the non-integer powers, e.g. a_i( x)^ν_i↦ e^ν_i log a_i( x).
  
  It follows that I is a projective integral over the projective simplex. The differential form Ω, which is homogeneous degree 0, is also called the canonical form on this simplex <cit.>.  
  
  Parametric Feynman integrals in quantum field theory are of the same type as the integral in eq. (<ref>) <cit.>. In string theory this type of integrals plays an equally  important role <cit.>. Integrals over positive geometries, which appear in the theory of scattering amplitudes can also be brought into this form <cit.>.
  
  The integral in eq. (<ref>) can be written as an integral over the positive orthant of ^n: ^n-1_>0 = {(x_1,…,x_n-1)∈^n-1: x_k > 0} by picking an affine chart for projective space, for instance 
  
      I = ∫__>0^n-1. ∏_i a_i(x)^ν_i/∏_j b_j(x)^ρ_j|_x_n=1∏_k=1^n-1 x_k/x_k,
  
  by pulling back the diffeomorphism ^n-1_>0→^n-1_>0, (x_1,…,x_n-1) ↦ [x_1, …, x_n-1, 1].
  
  Such an integral is called a generalized Euler-Mellin integral. Continuing a program started by Gelfand, Kapranov and Zelevinsky (GKZ) <cit.>, these integrals have been studied extensively by Nilsson and Passare and others <cit.>. 
  This analysis is compatible with the geometric sector decomposition approach as was shown by Schultka <cit.>, who studied these integrals using toric geometry.  Recently, generalized Euler-Mellin integrals gained new attention in the context of positive geometries, scattering amplitudes and string theory <cit.>. Along these lines also the analysis of Feynman integrals as GKZ-type hypergeometric functions has recently gained a lot of attention <cit.>. 
  
  As every non-homogeneous polynomial in n-1 variables can be homogenized by introducing a new variable, every generalized Euler-Mellin integral such as the one in eq. (<ref>) is equivalent to an integral of the projective form in eq. (<ref>). 
  For our considerations it will be more convenient to work with the projective form.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
   §.§ Outline of the paper
  
  After introducing the necessary preliminaries from polyhedral geometry and numerical integration in Section <ref>, we will establish the most important tool in this paper in Section <ref>: 
  an approximation of a multivariate polynomial, which is obtained by setting all its coefficients to 1 and replacing + by max. Starting for instance with the polynomial p(x_1, x_2, x_3) = a x_1^2 x_2 + b x_1 x_2 x_3 + c x_3^3, we get the `approximation' p^(x_1,x_2,x_3) = max(x_1^2 x_2, x_1 x_2 x_3, x_3^3). 
  This procedure is inspired from and closely related to tropical geometry. Hence, p^ will be called the tropical approximation of p. 
  This tropical approximation will be the subject of the main Theorem <ref> of this article, where it will be proven that p^ can indeed by used to approximate the polynomial p in a certain sense, as long as p is completely non-vanishing. In the rest of the article we will apply this property in various contexts.
  
  In Section <ref>, we will reformulate Kaneko-Ueda geometric sector decomposition in a tropical geometric framework. This reformulation will enable us to introduce the new tropical sampling algorithm in Section <ref>. This algorithm is significantly more efficient than traditional Monte Carlo methods as the achievable accuracy is effectively decoupled from the complexity of the integral. Only the runtime of a preprocessing step still depends heavily on the complexity of the integral. 
  
  This new algorithm can be improved further if more is known about the structure of the integrand polynomials {a_i} and {b_j}. As an example for this, we will specialize to the case where the Newton polytopes of these polynomials are generalized permutahedra in Section <ref>. Generalized permutahedra are a family of polytopes with a rich combinatorial structure. Many polytopes are from this family including associahedra, the Newton polytopes of Symanzik polynomials and other polytopes at play in the theory of scattering amplitudes <cit.>. We will use results by Postnikov <cit.>, Aguiar, Ardila <cit.> and Fujishige, Tomizawa <cit.> to formulate a specialized algorithm. This second new algorithm has more favorable runtime and memory requirements and is easier to implement. 
  
  Even though the improved algorithm can be applied to all generalized permutahedra integrands, as for instance the ones for complete scattering amplitudes introduced by Arkani-Hamed, Bai, He and Yan <cit.>, we will specify to Feynman integrals in Section <ref> for illustrative purposes. 
  
  The first, general tropical sampling algorithm can always be applied to Feynman integrals regardless of their explicit form. The second algorithm can only be applied if the Newton polytopes of the Symanzik polynomials are generalized permutahedra. We will use results by Brown <cit.> and Schultka <cit.> which ensure that this is the case as long as we are in a non-exceptional Euclidean kinematic region. Subsequently, we will discuss some experimental results which have been obtained using a proof-of-concept implementation of the second algorithm. 
  
  We conclude with a selection of future research directions resulting from this project in the last Section <ref>.
  
  
  
   §.§ Acknowledgements
  
  I am indebted to David Broadhurst, Iain Crump, Alejandro Morales and Erik Panzer for helpful discussions during the workshops at SFU in 2016 and at the University of Waterloo in 2018 both organized by Karen Yeats, the Programme on `Algorithmic and Enumerative Combinatorics' at the Erwin-Schrödinger Institute in 2017 and the `Summer school on structures in local quantum field theory' at the École de physique des Houches in 2018 organized by Dirk Kreimer.
  This work also greatly benefited from discussions with Francis Brown, Gudrun Heinrich, Franz Herzog, Dirk Kreimer, Oliver Schnetz, Konrad Schultka, Jos Vermaseren and Karen Yeats.
  I also wish to thank Erik Panzer for comments and suggestions on an early version of this manuscript and Jos Vermaseren for granting me generous access to his computing resources.
  
  This work has been supported by the NWO Vidi grant 680-47-551 `Decoding Singularities of Feynman graphs'.
  
  
  § PRELIMINARIES
  
  
  
  
   §.§ Notation for polytopes and multivariate polynomials
  
  The integral in eq. (<ref>) is convergent if the sets of polynomials {a_i} and {b_j} fulfill certain properties, which essentially have been determined by Nilsson and Passare <cit.>. In this section, we will briefly review these properties and introduce the necessary vocabulary from polyhedral and tropical geometry.
  
  To keep the notation simple, we will identify the space of linear forms on ^n with ^n via the usual scalar product v w = ∑_k=1^n v_k w_k. A polytope is the intersection of a finite number of half-spaces in ^n. A subset F ⊂𝒫 of a polytope 𝒫⊂^n is a face of 𝒫 if there is a vector y ∈^n and a scalar ξ∈ such that 𝒫 is contained in the half-space { v∈^n:  y v≤ξ} and F is the intersection of 𝒫 with the hyperplane { v∈^n:  y v = ξ}. Equivalently, a face of a polytope is a subset of 𝒫 which maximizes a given linear functional y ∈^n, F = { v ∈𝒫:  y v = max_ w ∈𝒫 y w}. We will assume polytopes to be bounded, i.e. max_ w ∈𝒫 y w < ∞ for all y∈^n.
  
  For a pair of non-negative real numbers λ, μ≥ 0 the weighted Minkowski sum of two polytopes 𝒫, 𝒬⊂^n is λ𝒫 + μ𝒬 = {λv + μw : v∈𝒫, w∈𝒬}⊂^n.  The relative interior, (𝒫), of a polytope 𝒫 is the interior of 𝒫 determined in the subtopology of the affine hull of 𝒫, which is the affine subspace of minimal dimension that contains 𝒫.
  
  We can write a generic multivariate polynomial p in the variables x_1, …, x_n as
  
      p(x_1,…, x_n) = p( x) = ∑_ℓ∈(p) c_ℓ∏_k=1^n x_k^ℓ_k = ∑_ℓ∈(p) c_ℓ x^ℓ,
  
  where (p), the support of p, is the set of all multi-indices (ℓ_1, …, ℓ_n) = ℓ∈^n such that c_ℓ≠ 0. We will make regular use of the multiplicative multi-index notation x^ℓ = ∏_k=1^n x_k^ℓ_k as above.
  The Newton polytope of p is the convex hull of the elements in (p) interpreted as vectors in ^n:
  
      _p= {∑_ℓ∈(p)λ_ℓℓ : ∑_ℓ∈(p)λ_ℓ = 1  and λ_ℓ≥ 0 }⊂^n.
  
  The Newton polytope _p of a homogeneous polynomial p in n variables is at most (n-1)-dimensional as it is contained in the hyperplane {v∈^n : v = p}⊃_p, where  is the only-ones-vector =(1,…, 1)∈^n.
  
  
  
  For each face F of _p associated to a polynomial p, we define the truncated polynomial p_F by
  
      p_F( x) = ∑_ℓ∈ F ∩(p) c_ℓ x^ℓ.
  
  
  
  To ensure convergence of integrals such as the one in eq. (<ref>), the following property of polynomials is useful:
  
  A polynomial p ∈[x_1, …, x_n] is completely non-vanishing on a domain X if for each face F ⊂_p, the truncated polynomial p_F does not vanish on X.
  
  
  
  With this terminology at hand we can give a convergence criterion for the integral in eq. (<ref>).
  
  
  We define the polytopes 𝒜,ℬ⊂^n as  the weighted Minkowski sums
  
      𝒜   = ∑_i (ν_i) _a_i   ℬ   = ∑_j (ρ_j) _b_j
  
  of the Newton polytopes of the numerator and denominator polynomials {a_i} and {b_j}.
  
  The integral in eq. (<ref>) is convergent if
  
  
    R1
  the denominator polytope ℬ is (n-1)-dimensional,
  
    R2
  the numerator polytope 𝒜 is contained in the relative interior of ℬ: 𝒜⊂ℬ,
  
    R3
  all the denominator polynomials {b_j} are completely non-vanishing on ^n-1_>0.
  
  
  
  
  
  
  Because each of the {a_i} and {b_j} polynomials is homogeneous, neither 𝒜 nor ℬ is full dimensional in ^n. 
  The condition in eq. (<ref>), which implies that ∑_i ν_i  a_i = ∑_j ρ_j  b_j, guarantees that 𝒜 and ℬ both lie in the same hyperplane 𝒜, ℬ⊂{v∈^n : v = ξ}, where ξ = ∑_i ν_i  a_i = ∑_j ρ_j  b_j in which ℬ is required to be full-dimensional by requirement R1.
  
  
  A similar theorem in the equivalent context of Euler-Mellin integrals was proven in <cit.> (see also <cit.>). The tropical approximation that we will introduce later will lead to an alternative proof of Theorem <ref> which we postpone to Section <ref>. 
  
  
  
  
  
  
  
  
   §.§ Monte Carlo quadrature
  
  We will be interested in situations where the dimension n of the integral in eq. (<ref>) is `not small'. The dimension of an integral is small from the perspective of numerical quadrature if fast-converging deterministic quadrature methods are feasible. The computational demands of deterministic quadrature methods such as Gauss-quadrature grow exponentially with the dimension. For this reason, it is necessary to use non-deterministic methods which do not suffer from an exponential slow-down if n is not small. Monte Carlo quadrature is the most elementary of these. 
  The working principle behind it is the following fact:
  
  
  If x^(1),…, x^(N) are independent random variables with probability density measure μ, i.e. 1 = ∫_Γμ and μ>0 on the domain Γ and
  
      G^(N) = 1/N∑_ℓ=1^N f( x^(ℓ)), then[ G^(N) ] =  [ f( x) ] = ∫_Γ f( x) μ
      and [ G^(N) ] = 1/N[ f( x) ]  where [f( x)] = ∫_Γ | f( x) - [ f( x) ] |^2 μ,
  
  provided that the integrals in the last two lines exist.
  
  This theorem may be applied to approximate the integral ∫_Γ f( x) μ as long as we have a way of generating samples from the distribution μ.  The condition that the integral for the variance shall exist effectively restricts the set of functions f, which can be integrated numerically, to the set of square integrable functions under the measure μ over the domain Γ. Note that in our convention of the statement of Theorem <ref> the expectation value [ · ] may be complex, but the variance [ ·] is always real and non-negative. 
  
  
  
  
  The integral in eq. (<ref>) is not directly amendable to Monte Carlo quadrature as the differential form Ω over the domain _>0^n-1 as defined for eq. (<ref>) is not a probability distribution: it is not normalizable. Even if we use the affine representation in eq. (<ref>) and map _>0^n onto a bounded domain via a variable transformation (for instance by x ↦ x/(1+x) which maps _>0→ (0,1) smoothly and injectively), the resulting integral will, in the general case, not be square integrable. A pragmatic solution to this problem is sector decomposition <cit.>, where the integral I is expressed as a sum of integrals, which are each individually directly amendable to Monte Carlo integration. 
  
  
  
  
  
  
  
  
   §.§ Sector decomposition
  
  
  
  In the context of quantum field theory, sector decomposition goes back to Hepp and Speer who used the technique to prove the finiteness of renormalized Feynman integrals <cit.>. Even though Hepp/Speer sector decomposition can be employed to deal with the singularities of Euclidean Feynman integrals, it turned out to be insufficient to handle more general singularities which appear in Minkowski space Feynman integrals (see <cit.> and <cit.> for reviews on the topic). 
  A more general approach was pioneered by Binoth and Heinrich <cit.>, who introduced a recursive algorithm that decomposes general integrals of the form in eq. (<ref>) (or equivalently as in eq. (<ref>)) into a set of sector integrals: 
  
      I    = ∑_s ∈ S I_s     I_s    = C_s ∫_[0,1]^n-1 x^ m^(s)∏_i p^ν_i_s,i( x)/∏_j q^ρ_j_s,j( x)∏_k=1^n-1 x_k/x_k,
  
  such that the auxiliary polynomials p_s,i( x) and q_s,j( x) do not vanish inside the integration domain [0,1]^n-1 (at least as long as all the coefficients of the initial denominator polynomials b_j are positive, which implies that these polynomials are completely non-vanishing) and C_s is a prefactor for each sector s∈ S. If all components of the vector m^(s) are positive, i.e. m^(s)_k > 0 for all k∈{1,…, n-1}, a simple reparametrization ξ =  x^ m^(s) produces an integral with a bounded integrand 
  
      I_s = C_s ∫_[0,1]^n-1∏_i p^ν_i_s,i( x(ξ))/∏_j q^ρ_j_s,j( x(ξ))∏_k=1^n-1ξ_k,
  
  to which basic Monte Carlo as described in Theorem <ref> can immediately be applied using the uniform measure μ = ∏_k=1^n-1ξ_k on the unit (n-1)-cube Γ = [0,1]^n-1. 
  
  
  Although a problem of the method which impeded the recursion from terminating was solved by Bogner and Weinzierl <cit.>, this class of algorithms suffers from a proliferation in the numbers of sectors with rising complexity of the underlying polynomials. Moreover, a new set of polynomials p_s,i, q_s,i is associated to each sector. This means that we are not necessarily dealing with a partition of the integration domain, but a non-trivial distribution of the volume of the integral into each of the sectors s∈ S. 
  
  A both conceptual and practical innovation was achieved by Kaneko and Ueda, who reinterpreted this decomposition as a geometric problem <cit.>. This geometric viewpoint results in more economical decompositions, in terms of the total number of sectors (see <cit.> for a comparison of different methods), while also arguably being conceptually more elegant. 
  
  
  
  
  
   §.§ Analytic continuation
  
  By Theorem <ref>, the convergence of the integral in eq. (<ref>) depends on the values of the parameters {ν_i} and {ρ_j}.  Provided that there is an extended domain of such parameters where the integral is convergent, we can interpret it as a function of these parameters and perform an analytic continuation. It turns out that this analytic continuation is a meromorphic function in these parameters <cit.>. 
  
  A sector decomposition as in eq. (<ref>) provides a pragmatic way to perform this analytic continuation. 
  A violation of the condition 𝒜⊂ℬ in Theorem <ref> corresponds to a component of m^(s) in the integral in eq. (<ref>) being non-positive, i.e. m^(s)_k ≤ 0 for some sector s∈ S. Hence, if we assume that the polynomials p_s,i, q_s,i are non-vanishing on the integration domain, then the associated integral I_s is divergent. Performing an analytic continuation of this integral, interpreted as a function of the coefficients of m^(s) is a simple task. A standard approach is to integrate over a Pochhammer contour instead of the unit interval in eq. (<ref>). This avoids the singularity at the integration boundary and agrees with the integral over the unit interval if convergence is ensured (see for instance <cit.>).
  
  For explicit computations it is sufficient to compute a Taylor expansion of the rational function in eq. (<ref>) up to an appropriate order, integrate the analytically continued expansion terms analytically and the remainder term numerically. See for instance <cit.>, where this process is described in detail.
  
  A more sophisticated strategy to perform this analytic continuation is based on iteratively performing `directed integration by parts' on the integral in eq. (<ref>) and thereby extending the domain of {ν_i}, {ρ_j} parameters in which the integral converges. The inner workings of this procedure are of geometric nature and make use of the structure of the polytopes 𝒜, ℬ. See <cit.> and thereafter and also <cit.> for a description of this method. This approach gives, after being applied to a given integral as the one in eq. (<ref>), a sum of integrals of the same type where each integral has a larger region of convergence than the original one. A similar procedure has been developed independently in <cit.> for the special case of parametric Feynman integrals. 
  
  In this work, we will therefore assume that the integral has been subjected to such a procedure and we can assume that we are within the region of convergence in terms of the {ν_i}, {ρ_j} parameters.
  
  
  § THE TROPICAL APPROXIMATION
  
  
  
  For the considerations in this article the following tropical approximation of a polynomial will be central:
   
  
  For a polynomial p ∈[x_1, …, x_n] define 
  p^( x) = max_ℓ∈(p) x^ℓ.
  
  Such an object has been defined by Panzer <cit.> for the Kirchhoff polynomial to study the Hepp-bound, a graph invariant relevant for Feynman period integral calculations. We adopt Panzer's notation and denote tropically approximated polynomials with a superscript ^.
  
  To give some additional motivation to consider this `tropical approximation' suppose that a polynomial p has only real and positive coefficients and interpret it as a function p:_> 0^n →_> 0. The tropical limit is
  
      lim_ξ→∞ p(x_1^ξ,…, x_n^ξ)^1/ξ = lim_ξ→∞( ∑_ℓ∈(p) c_ℓ x^ξℓ)^1/ξ = max_ℓ∈(p) x^ℓ = p^( x).
  
  This way, p^ can be seen as a deformed version of p: the function p( x^ξ)^1/ξ interpolates between p and p^ with ξ between 1 and ∞. A limit as ξ→∞ with the associated phenomenon of transforming a very smooth object—in this case a polynomial—into a function with non-differentiable singularities, is something commonly encountered in physics. For instance, the thermodynamical limit is of similar nature. These kind of limits give rise to numerous critical phenomena. Also the weak string coupling limit α' → 0 shows this behavior <cit.>.
  
  In our case p^ is of `simpler' nature as the original polynomial p. Information is lost while going from a polynomial to its tropical approximation, as p^ only depends on the support of p. In fact, p^ is nothing but a realization of a geometric object: the Newton polytope of the polynomial p.
  
  
  
  To make this explicit, change to logarithmic coordinates y_k = log x_k in Definition <ref> and use the fact that the Newton polytope is the convex hull of the support of the underlying polynomial. We find that
  
      log p^(x) = max_ℓ∈(p)log x^ℓ = max_ℓ∈(p)∑_k=1^n y_k ℓ_k = max_ℓ∈(p)yℓ = max_v∈_pyv,
  
  which is a piece-wise linear function ^n →,  y ↦max_v∈_pyv. 
  This function is the support function of the Newton polytope _p <cit.>.
  
  
  
  
  
  
  Often it is useful to write p^ in exponential form using the support function:
  
  
  p^(e^y) = e^max_ v ∈_p y v,
  
  where we used the notation e^y=(e^y_1, …, e^y_n)  to denote the component-wise exponential.
  
  
  
  
  
  
  
  
  This support function is also a tropicalization of the polynomial p, which uses the trivial valuation on  to tropicalize. Because much of the algebro-geometrical information of the polynomial p carries over to its tropicalization, tropical geometry developed into a fruitful branch of algebraic geometry in the recent years <cit.>.  
  It is tempting to call p^ the tropicalization of p. Unfortunately, this name is reserved for y ↦log p^(e^ y) and we will use the name tropical approximation instead. The motivation for this is that besides the fact that p^ is a simplification of p, it can also be used to approximate p.
  
  
  
   §.§ The approximation property
  
  The main theorem of this article is the following approximation property of p^ with respect to the polynomial p: 
  theorem
  
  
  
  For every polynomial p ∈[x_1, …, x_n] there is a constant C>0 such that
  
      | p(x) |    ≤ C p^ (x)  for all x∈^n_> 0.
  
  
  
  C = ∑_ℓ∈(p) |c_ℓ|.
  
  
  
  
  If p ∈[x_1, …, x_n] is completely non-vanishing on ^n_>0, then there is a constant C>0 such that
  
      C p^(x) ≤ | p(x) |  for all x∈^n_> 0.
  
  
  
  
  This property trivially extends to homogeneous polynomials which are naturally considered to be functions on projective space:
  
  
  For every homogeneous polynomial p ∈[x_1, …, x_n] there is a constant C>0 such that
  
      | p(x) |    ≤ C p^ (x)  for all x∈^n-1_> 0.
  
  and if p is additionally completely non-vanishing on ^n-1_>0, then there is a constant C>0 such that
  
      C p^(x) ≤ | p(x) |  for all x∈^n-1_> 0.
  
  
  Strictly speaking p(x) and p^(x) are not well-defined objects for x∈^n-1_> 0. The quotient p(x)/p^( x), on the other hand, makes sense for all projective x∈^n-1_> 0. Inequalities as the one above, which can be written as quotients of homogeneous objects, shall be interpreted accordingly, in this obvious sense as statements on these quotients. 
  
  If p is homogeneous and completely non-vanishing on ^n-1_>0, it is also completely non-vanishing on _>0^n. The inequality in Theorems <ref> and <ref> is homogeneous, therefore it trivially extends to ^n-1_>0.
  
  
  By Theorem <ref>, p^ can indeed be used to `approximate' p, i.e. it provides a lower and an upper bound of p with appropriate prefactors, as long as p is completely non-vanishing. 
  
  Theorem <ref> is substantially harder to prove than Theorem <ref>. This proof of Theorem <ref> will be given in the next Section <ref>. Only a special case is also trivial: If the polynomial p has only positive coefficients (which implies that p is completely non-vanishing on _>0^n), there is a simple lower bound for |p( x)|: take for instance C = min_ℓ∈(p) c_ℓ. 
  For such a lower bound to exist it is not necessary for the polynomial to have only positive coefficients; it is sufficient for the polynomial to be completely non-vanishing. In fact, the existence of such a lower bound is also necessary for a polynomial to be completely non-vanishing, which can be proven using a similar argument as in the proof of Theorem <ref> below. 
  
  
  
  
  
  
   §.§ Cones and normal fans
  
  
  
  A (polyhedral) cone is a subset of ^n that is closed under linear combinations with only non-negative scalars, e.g. 𝒞 = {∑_kλ_k  u^(k) : λ_k ≥ 0 } for some set of given vectors u^(1),  u^(2), …∈^n.
  A fan in ^n is a family ℱ = {𝒞_1, 𝒞_2, …} of cones with the property that every face of a cone in ℱ is also in ℱ and that the intersection of two cones 𝒞_1, 𝒞_2 ∈ℱ is a face of both 𝒞_1 and 𝒞_2. The normal cone associated to a face F of the polytope 𝒫 is the set of all linear functionals that are maximal on the respective face, 
  
      𝒞_F = { y ∈^n :  y v = max_ w ∈𝒫 y w for all  v ∈ F }         .
  
  The set of normal cones of a polytope is its normal fan: ℱ_N = {𝒞_F : F ∈(𝒫), F ≠∅}. The normal fan is always complete, that means ^n = _𝒞∈ℱ_N𝒞, where ⊎ denotes the disjoint union. If a face F has dimension d, then the associated normal cone 𝒞_F has dimension n-d, where n is the dimension of the ambient space. The maximal cones in the fan are the cones of maximal dimension. Figure <ref> and <ref> depict a polytope and its normal fan. Note that the maximal cones can be associated to the vertices of the polytope. 
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  For the proof of Theorem <ref>, it is convenient to have three further properties of the tropical approximation p^ and the associated polytopes at hand:
  
  
  
  For each face F of the Newton polytope _p of a polynomial p, the truncated polynomial p_F fulfills
  
      p_F(e^ s+ t) = p^(e^ s) p_F(e^ t)  for all  s ∈𝒞_F  and  t ∈^n.
  
  
  
  Use Definition <ref>, eq. (<ref>) and Proposition <ref>.
  
  
  
  
  p^(e^ s+ t) ≤ p^(e^ s) p^(e^ t)  for all  s,  t ∈^n.
  
  
  This follows from Proposition <ref> and 
  
      max_ v ∈𝒫( s +  t) v≤max_ v ∈𝒫 s v +max_ v ∈𝒫 t v.
  
  
  
  
  
  
  
  
  If 𝒫⊂^n is a polytope, then there exists a constant C > 0 such that
  
      max_ v ∈𝒫 y v -  y k≥ C R 
       for all  R ≥ 0  and all faces  F ⊂𝒫 with  k ∈ (𝒫) ∖ F 
       and  y ∈𝒞_F∖⋃_ F' >  F( 𝒞_F' + B_R ),
  
  where (𝒫) is the set of vertices of 𝒫, the union is over all faces of 𝒫 of higher dimension than F and B_R ⊂^n is a ball of size R.
  
  A set of modified cones 𝒞_F = 𝒞_F∖⋃_ F' >  F( 𝒞_F' + B_R ) is depicted in Figure <ref>.
  Note that the size of the `gaps' between the modified cones is the diameter 2R of the ball B_R.
  
  For a given face F ⊂𝒫 and k ∈(𝒫) ∖ F, choose some vertex v ∈(F) and consider the hyperplane H_ v = { y ∈^n :  y( k -  v) = 0 }.
  By definition of the normal cone, this hyperplane will not intersect with the interior of 𝒞_F as y v≥ y w for all w ∈𝒫 and y ∈𝒞_F.
  We can project any point y ∈𝒞_F onto H_ v using the orthogonal projection 
  y^⊥ =  y - ( k - v)  y( k - v)/ k - v^2∈ H_ v. The line segment from y to y^⊥ will intersect a face of 𝒞_F.
  
  Let 𝒞_F' be this face. Clearly, F' >  F.
  If y ∉𝒞_F' + B_R, then the vector y must have a larger distance than R from all points in 𝒞_F'.
  
  By construction, the orthogonal projection y^⊥ is at least as far away from y as the closest point in 𝒞_F'.
  
  
  Hence, y^⊥ - y = |  y( v -  k) | / v -  k > R and it follows that 
  
      y( v -  k) > R  v -  k for all  y ∈𝒞_F ∖ (𝒞_F' + B_R),  v ∈(F)  and  k ∈ (𝒫) ∖ F,
  
  where we used that y( v -  k)≥ 0 for all v∈ F, k ∈𝒫 and y ∈𝒞_F.
  To prove the statement, choose C = min_ v≠ w ∈(𝒫)  v -  w.
  
  
  Theorem <ref> follows now as a Corollary from the following 
  
  
  If p ∈[x_1,…,x_N] is completely non-vanishing on _>0^n and R ≥ 0, then
  there is a constant C > 0 such that 
  
      C p^(e^s+t) ≤ | p(e^ s +  t  ) |  for all  s ∈𝒞_F  and  t ∈ B_R  for each face F ⊂_p,
  
  where B_R⊂^n is a ball of radius R.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  We are going to prove this by induction in the codimension of F. Starting with codimension 0, i.e. F = _p, we have by Definition <ref>, Lemma <ref> and Lemma <ref>,
  
      | p(e^s + t)| = | p__p(e^s + t)| = p^( e^s) | p(e^t) | ≥ p^(e^ s +  t) |p(e^t)| /p^(e^ t) for all  s ∈𝒞__p and  t ∈^n.
  
  Because p(e^ t) is non-vanishing on the compact domain t ∈ B_R, we can choose the constant C = min_ t ∈ B_R |p(e^t)|/p^(e^ t) > 0 to get the desired bound. 
  
  Suppose F is of codimension d and eq. (<ref>) holds for all faces up to codimension d-1. 
  By the induction hypothesis, for each R' > 0 there exists a constant C >0 such that eq. (<ref>) is fulfilled in a ball of size R' around all cones 𝒞_F' of lower dimension, 𝒞_F' < 𝒞_F. We therefore only need to prove the existence of such a constant for the smaller domain,
  𝒞_F = 𝒞_F ∖⋃_ F' >  F ( 𝒞_F' + B_R' ), where the union is over all faces of _p of higher dimension than F. See Figure <ref> for an illustration of this smaller domain.
  
  
  
  By definition of the truncated polynomial, we can write p(e^s + t) as,
  
      p(e^s + t)    = p_F(e^ s + t) + ∑_k∈(p) ∖ F c_ k e^(s + t) k for all  s ∈𝒞_F and  t ∈^n,
  
  and estimate using Proposition <ref>, Lemma <ref>, Lemma <ref> and Lemma <ref>,
  
      |p(e^s + t)|    ≥ p^(e^ s) ( |p_F(e^t)| - ∑_k∈(p) ∖ F |c_ k| e^t k e^s k - max_ v ∈_p s v) 
         ≥ p^(e^ s+  t) ( |p_F(e^t)|/p^(e^ t) - e^-R' C' ∑_k∈(p) ∖ F |c_ k| e^t k/p^(e^ t)) 
         ≥ p^(e^ s+  t) ( |p_F(e^t)|/p^(e^ t) - e^-R' C' ∑_k∈(p) ∖ F |c_ k| )  for all  s ∈𝒞_F and  t ∈^n,
  
  where C' is the constant we obtained from Lemma <ref>.
  We can choose a constant C such that 2 C = min_ t ∈ B_R |p_F(e^t)|/p^(e^ t) > 0 and 
  R' > 1/C'log (∑_k∈(p) ∖ F |c_ k|/C). This gives the desired bound.
  
  
   
  Use Proposition <ref> and the completeness of the normal fan. 
  
  
  
  In the proofs of Lemma <ref> and Proposition <ref>, we actually constructed explicit bounds for the constants in Theorem <ref> which depend on the geometry of the relevant polytopes and polynomials. These explicit bounds might be useful for further considerations, but we will not make use of them in this article.
  
  
  
  
  § GEOMETRIC SECTOR DECOMPOSITION
  
  
  
  
  The rough overall plan of our take on the integral in eq. (<ref>) is as follows: we can trivially `factorize' its integrand and write it as
  
      I = ∫__>0^n-1∏_i a_i^(x)^ν_i/∏_j b_j^(x)^ρ_j( ∏_i a_i( x)^ν_i/a_i^( x)^ν_i/∏_j b_j( x)^ρ_j/b_j^( x)^ρ_j) Ω .
  
  The second factor is bounded by Corollary <ref> as long as the {b_j} polynomials are completely non-vanishing. The first term has a geometric interpretation in terms of the polytopes 
  𝒜 = ∑_i (ν_i) _a_i and ℬ = ∑_j (ρ_j) _b_j as defined in Theorem <ref>.
  
  As before it will be handy to change to logarithmic coordinates to expose this geometric interpretation. The component-wise exponential : ^n→^n_>0,  y ↦ e^ y = (e^y_1, …, e^y_n), extends to a smooth bijective map : ^n/→^n-1_>0, as  respects the respective equivalence relation. That means if x = e^ y and x' = e^ y' with y,  y' ∈^n,
  then y' =  y + μ for some μ∈ if and only if
  x' = λ x  for some λ∈_>0.
  For this reason the quotient ^n/ is also called tropical projective space.
  
  By Proposition <ref> and the definition of the weighted Minkowski sum with x = e^ y
  
      ∏_i a_i^( x)^ν_i/∏_j b_j^( x)^ρ_j   = exp( ∑_iν_i max_ v∈_a_i y v - ∑_jρ_j max_ v∈_b_j y v) 
         = exp(max_ v∈𝒜 y v - max_ v ∈ℬ y v).
  
  If 𝒜 and ℬ fulfill the
  requirements R1 and R2 of Theorem <ref>, this exponent is falling sufficiently fast for large y for the integral in eq. (<ref>) to be convergent.
  
  
  
  
  Let 𝒜 and ℬ be the polytopes defined in Theorem <ref>. If 𝒜 and ℬ fulfill the
  requirements R1 and R2 of Theorem <ref>,
  
  
  then there is a constant ε >0 such that 
  
      max_ v ∈ℬ y v - max_ v ∈𝒜 y v≥ε y_^n/ for all  y ∈^n/,
  
  where y _^n/=inf_μ∈ y+μ is the norm on the quotient space ^n/ induced from the standard norm · on ^n.
  
  
  
  First note that the inequality in the statement is well-defined for y ∈^n/, as 𝒜 and ℬ both lie in the same hyperplane 𝒜, ℬ⊂ H_ξ = {v∈^n : v = ξ} for some ξ∈. Therefore, y v =  y w for all v ∈𝒜, w ∈ℬ and y ∈.
  
  As ℬ is full-dimensional in H_ξ (see Remark <ref>) 
  and 𝒜⊂ℬ,
  we can Minkowski add a ball B_ε to 𝒜 such that 𝒜 + B_ε⊂ℬ, provided that this ball only extends in the subspace orthogonal to  which is parallel to H_ξ and ε is sufficiently small. Let B_ε = { v ∈^n :  v = 0  and  v≤ε} be such a ball.  The resulting convex set 𝒜 + B_ε is the outer parallel body of 𝒜 restricted to its affine hull.
  
  
  Observe that max_ v ∈ B_ε y v = εinf_μ∈ y+μ= ε y _^n/. By the definition of the Minkowski sum and because 𝒜 + B_ε⊂ℬ, 
  
      max_ v ∈𝒜 y v + ε y _^n/ = max_ v ∈𝒜 y v + max_ v ∈ B_ε y v
       = max_ v ∈𝒜 + B_ε y v≤max_ v ∈ℬ y v for all  y ∈^n.
  
  
  
  To be able to handle Kaneko-Ueda geometric sector decomposition with our tropical approach, we will need additional tools from convex geometry. 
  
  A cone is called pointed if it contains no 1-dimensional subspace. A fan is pointed if all its cones are pointed. If a polytope 𝒫⊂^n is full-dimensional, its normal fan ℱ is pointed. For lower dimensional polytopes 𝒫⊂^n, the normal cone 𝒞_𝒫 associated to the polytope itself is non-trivial: 𝒞_𝒫 = 𝒫^⊥ = { y ∈^n :  y v =  y w for all  v,  w ∈𝒫}. It consists of all linear functionals that are constant on 𝒫. This subspace is contained in each cone of the normal fan. Taking the quotient with respect to this subspace within each cone in the normal fan 𝒞∈ℱ results in a pointed fan ℱ/𝒫^⊥ on the quotient vector space ^n/𝒫^⊥. This fan is the reduced normal fan.
  
  Given a fan ℱ, another fan ℱ' refines ℱ if every cone in ℱ is a union of cones in ℱ'.  If ℱ and 𝒢 are both fans, then their common refinement is defined as ℱ∧𝒢 = {𝒞∩𝒞' : 𝒞∈ℱ, 𝒞' ∈𝒢}.  
  Let ℱ_𝒜ℬ be the common refinement of the normal fans of the polytopes 𝒜 and ℬ which were defined in Theorem <ref>.
  Recall that the polytopes 𝒜, ℬ are not full-dimensional, because they are weighted Minkowski sums of the Newton polytopes of homogeneous polynomials (see Remark <ref>). As ℬ is required to be full-dimensional in a (n-1)-dimensional hyperplane which is orthogonal to the -vector, we have ℬ^⊥ =. We will therefore consider the reduced refined normal fan ℱ_𝒜ℬ /, which is pointed.
  The following lemma identifies the exponentiated cones of 𝒞∈ℱ_𝒜ℬ / as the domains where the function ∏_i a_i^( x)^ν_i/∏_j b_j^( x)^ρ_j behaves like a monomial.
  
  
  
  
  
  Let 𝒜 and ℬ be the polytopes defined in Theorem <ref>. If 𝒜 and ℬ fulfill the
  requirements R1 and R2 of Theorem <ref>
  
  and 𝒞 is a cone in the reduced common refinement ℱ_𝒜ℬ/, then
  
      ∏_i a_i^( x)^ν_i/∏_j b_j^( x)^ρ_j =  x^- w for all  x ∈ (𝒞),
  
  where 
  w =  w_ℬ -  w_𝒜 
  and w_𝒜∈𝒜,  w_ℬ∈ℬ 
  
  such that 
  y w_𝒜 = max_ v ∈𝒜 y v
  and y w_ℬ = max_ v ∈ℬ y v
  for all y ∈𝒞. Moreover, w = 0 and y w > 0 for all y ∈𝒞∖{0}.
  
  
  
  
  
  As 𝒞 is a refinement of the normal fans of 𝒜 and ℬ, there must be normal cones 𝒞_F_𝒜^𝒜 and 𝒞_F_ℬ^ℬ associated to respective faces F_𝒜⊂𝒜 and F_ℬ⊂ℬ such that 𝒞⊂𝒞_F_𝒜^𝒜∩𝒞_F_ℬ^ℬ.
  Hence,
  max_ v∈ℬ y v - max_ v ∈𝒜 y v =  y( w_ℬ -  w_𝒜 ) for all  y ∈𝒞,
  where we can choose arbitrary w_𝒜∈ F_𝒜 and w_ℬ∈ F_ℬ by definition of the normal cone in eq. (<ref>).
  
  Since 𝒜 and ℬ are required to lie in the same hyperplane orthogonal to the -vector, we also have ( w_ℬ -  w_𝒜 )  =  w = 0.
  Due to Lemma <ref>, 
  max_ v∈ℬ y v - max_ v ∈𝒜 y v =  y( w_ℬ -  w_𝒜 )≥ε y _^n/> 0  for all  y ∈𝒞∖{0} and the statement follows from eq. (<ref>).
  
  
  
  
  A cone 𝒞 is simplicial if it is generated as, 𝒞 = {∑_k=1^d λ_k  u^(k) : λ_k ≥ 0 } where u^(1), … u^(d) are linear independent.  
  
  For a given cone 𝒞, we can always find a set of simplicial cones 𝒞^Δ_1, 𝒞^Δ_2, … such that 𝒞 = _i𝒞^Δ_i and 𝒞^Δ_i ∩𝒞^Δ_j is a simplicial cone which is a face of both 𝒞^Δ_i and 𝒞^Δ_j. Such a set of simplicial cones is called a triangulation of 𝒞. 
  Let ℱ_𝒜ℬ^Δ/ be a simplicial refinement of ℱ_𝒜ℬ/, i.e. a refinement of ℱ_𝒜ℬ/ such that each cone in ℱ_𝒜ℬ^Δ/ is simplicial. 
  
  A feature of simplicial cones is that there are convenient coordinates describing points in their interior. This fact is important while proving the following lemma:
  
  
  If a pointed simplicial cone 𝒞⊂^n/ is generated by
  linear independent vectors u^(1), …,  u^(n-1)∈^n/, i.e. 𝒞 = {∑_k=1^n-1λ_k  u^(k): λ_k ≥ 0 }, w ∈^n with w = 0 and
  y w > 0 for all y ∈𝒞∖{0}, then 
  
      ∫_(𝒞) x^- w f( x) Ω = |( u^(1), …,  u^(n-1),) |/∏_k=1^n-1 u^(k) w∫_[0,1]^n-1 f(  x(ξ) ) ∏_k=1^n-1ξ_k,
  
  where 
  f:_>0^n-1→ is a measurable homogeneous function of degree 0 and
  x(ξ) ∈(𝒞) is given component-wise by x_k = ∏_i = 1^n-1ξ_i^-u_k^(i)/( u^(i) w). 
  
  
  
  
  By slightly abusing the notation, we identified the vectors u^(1), …,  u^(n-1) with appropriate representatives in ^n in the statement of this lemma.  The value of the integral does not depend on the specific choice of representatives, because w = 0 and f is homogeneous of degree zero. Hence, the expression on the right hand side is invariant under shifts u^(k)→ u^(k) + μ_k for all μ_k ∈.
  It is also invariant under rescalings of the vectors u^(k)→λ_k  u^(k) for all λ_k > 0 as it should be due to the equivalence of the cone representation. 
  Even though changing the representatives of the u-vectors modifies the vector x(ξ), it only does so by an overall scaling, which does not modify the point in ^n-1_>0, which x(ξ) represents.
  
  
  
  Start by changing to logarithmic coordinates x = e^ y,
  
      ∫_(𝒞) x^- w f( x) Ω = ∫_𝒞 e^- y w f(e^ y) Ω
  
  where
  Ω = ^* Ω= ∑_k=1^n (-1)^n-k y_1∧…∧ y_k ∧…∧ y_n is the pullback of Ω under .
  Using barycentric coordinates y = ∑_k=1^n-1 u^(k)λ_k shows that this is equal to
  
      = |( u^(1), …,  u^(n-1),)| ∫__> 0^n-1 e^-∑_k=1^n-1λ_k  u^(k) w f(e^∑_k=1^n-1λ_k  u^(k)) ∏_k=1^n-1λ_k.
  
  The form of the determinant follows from the form Ω, Laplace's expansion and 
  
      y_1 ∧…∧ y_k∧…∧ y_n =  u_1^(𝒞,1)       u_k^(𝒞,1)       u_n^(𝒞,1)
      ⋮   ⋱   ⋮   ⋱   
       u_1^(𝒞,n-1)       u_k^(𝒞,n-1)       u_n^(𝒞,n-1)λ_1 ∧…∧λ_n-1.
  
  As y w > 0 for all y ∈𝒞∖{0}, it follows that u^(k) w > 0 for all k ∈{ 1,…, n-1}.
  We can therefore change variables via λ_k = -1/ u^(k) wlogξ_k which proves the statement.
  
  
  With these tools at hand, we are ready to give our tropical formulation of geometric sector decomposition:
  
  
  
  
  Let 𝒜 and ℬ be the polytopes defined in Theorem <ref>. If 𝒜 and ℬ fulfill the
  requirements R1 and R2 of Theorem <ref>
  
  and ℳ_𝒜ℬ^Δ⊂ℱ_𝒜ℬ^Δ/ is the set of maximal cones, i.e. the cones of maximal dimension, in a simplicial refinement of the reduced common normal fan of 𝒜 and ℬ, then we can write the integral 
  
      I[f] = ∫__>0^n-1∏_i a_i^( x)^ν_i/∏_j b_j^(x)^ρ_j f( x) Ω.
  
  as a sum 
  I[f] = ∑_𝒞∈ℳ_𝒜ℬ^Δ I_𝒞[f] with
  
      I_𝒞[f]    = |( u^(𝒞,1), …,  u^(𝒞,n-1),) | /∏_k=1^n-1 u^(𝒞,k) w^(𝒞)∫_[0,1]^n-1 f (  x^(𝒞)(ξ) ) ∏_k=1^n-1ξ_k,
  
  where 
  
  
    * f:_>0^n-1→ is a measurable homogeneous function of degree 0,
  
    * w^(𝒞) =  w^(𝒞)_ℬ -  w^(𝒞)_𝒜 with
  some w^(𝒞)_𝒜∈𝒜 and w^(𝒞)_ℬ∈ℬ such that
  y w^(𝒞)_𝒜 = max_ v ∈𝒜 y v and 
  y w^(𝒞)_ℬ = max_ v ∈ℬ y v for all y ∈𝒞,
  
    * the vectors u^(𝒞,1), …,  u^(𝒞,n-1)∈^n/ span the simplicial cone 𝒞 such that 𝒞 = {∑_k=1^n-1λ_k  u^(𝒞,k): λ_k ≥ 0 }⊂^n/,
  
    * x^(𝒞)(ξ) ∈(𝒞) is given component-wise by x_k^(𝒞) = ∏_i = 1^n-1ξ_i^-u_k^(𝒞,i)/( u^(𝒞,i) w^(𝒞)) and
  
    * the prefactor |( u^(𝒞,1), …,  u^(𝒞,n-1),) | /∏_k=1^n-1 u^(𝒞,k) w^(𝒞) is finite and positive for each 𝒞∈ℳ_𝒜ℬ^Δ.
  
  
  
  
  The fan ℱ_𝒜ℬ^Δ/ is complete, i.e. it corresponds to a partition of ^n/ = _𝒞∈ F_𝒜ℬ^Δ/𝒞. 
  Because : ^n/→^n-1_>0 is smooth and bijective this partition gives also a partition of ^n-1_>0= _𝒞∈ F_𝒜ℬ^Δ/ (𝒞). 
  
  
  Since we would like to integrate over ^n-1_>0 or (equivalently over ^n/) it is enough to only consider the cones of maximal dimension ℳ_𝒜ℬ^Δ⊂ℱ_𝒜ℬ^Δ/ as other cones in ℱ_𝒜ℬ^Δ/ only describe measure zero subsets of ^n-1_>0. Hence,
  
      I    = ∑_𝒞∈ℳ_𝒜ℬ^Δ I_𝒞    I_𝒞   = ∫_(𝒞)∏_i a_i^( x)^ν_i/∏_j b_j^(x)^ρ_j f( x) Ω.
  
  Because each cone 𝒞∈ℳ_𝒜ℬ^Δ refines a cone in ℱ_𝒜ℬ/ and because of Lemma <ref>, 
  
      I_𝒞   = ∫_(𝒞) x^- w^(𝒞) f( x) Ω for all 𝒞∈ℳ_𝒜ℬ^Δ,
  
  where 
  w^(𝒞) =  w^(𝒞)_ℬ -  w^(𝒞)_𝒜 with
  some w^(𝒞)_𝒜∈𝒜 and w^(𝒞)_ℬ∈ℬ such that
  y w^(𝒞)_𝒜 = max_ v ∈𝒜 y v and 
  y w^(𝒞)_ℬ = max_ v ∈ℬ y v for all y ∈𝒞. By Lemma <ref>, we also have w^(𝒞) = 0 and y w^(𝒞) > 0 for all y ∈𝒞∖{0}.
  
  Eq. (<ref>) follows from Lemma <ref>, because  we can always pick a set of generators u^(𝒞,1), …,  u^(𝒞,n-1)∈^n/ for every simplicial cone 𝒞∈ℳ_𝒜ℬ^Δ.
  As y w^(𝒞) > 0 for all y ∈𝒞∖{0}, we also have u^(𝒞,k) w^(𝒞) > 0 for all k ∈{ 1,…,n-1}. The positivity of the determinant is obvious because of the linear independence of the vectors u^(𝒞,k).
  
  
  If we specify f( x) = R_a/b( x) given by 
  
      R_a/b( x) = ∏_i a_i( x)^ν_i/a_i^( x)^ν_i/∏_j b_j( x)^ρ_j/b_j^( x)^ρ_j,
  
  in Theorem <ref>, we recover the integral in eq. (<ref>).
  
  
  We only need to prove that each sector integral in the geometric sector decomposition of Theorem <ref> with f( x) = R_a/b( x) from eq. (<ref>) is finite. As all the denominator polynomials {b_j} are completely non-vanishing, Corollary <ref> implies that |R_a/b( x)| is bounded on _>0^n-1. Hence, each integral I_𝒞[R_a/b] is finite.
  
  
  
  
  
  Theorem <ref> provides a sector decomposition as it was formulated in eq. (<ref>), because the sector integrands in Theorem <ref> are bounded as long as the function f is bounded on _>0^n-1. This way, Theorem <ref> not only ensures finiteness of the integral in eq. (<ref>) under appropriate conditions, but also allows to evaluate the integral via Monte Carlo quadrature. 
  
  
  If we have triangulated the reduced normal fan ℱ_𝒜ℬ/, i.e. we have computed a simplicial refinement ℱ^Δ_𝒜ℬ/ and stored the vectors u^(𝒞,1), …,  u^(𝒞,n-1) and w^(𝒞) for each maximal cone 𝒞∈ℳ^Δ_𝒜ℬ⊂ℱ^Δ_𝒜ℬ/ in a table, then we can estimate the integral using Algorithm <ref>.
  
  
  
  
  If the conditions of Theorem <ref> are fulfilled, then the random value I^(N) returned by Algorithm <ref> has expectation value equal to the integral in eq. (<ref>), I = [I^(N)] and [I^(N)] = C/N with some constant C ≥ 0.
  
  
  Algorithm <ref> is an application of Theorem <ref> on the integral I_𝒞[f] for each cone 𝒞∈ℳ_𝒜ℬ^Δ from Theorem <ref> with f( x) = R_a/b( x): 
  
      ∫_[0,1]^n-1 R_a/b(  x^(𝒞)(ξ) ) ∏_k=1^n-1ξ_k.
  
  As |R_a/b( x)| is bounded on _>0^n-1 and x^(𝒞)(ξ) ∈_>0^n-1 by construction, the integrand is bounded and therefore also square integrable. Hence, there is a constant C_𝒞≥ 0 for each cone integral such that [I_𝒞^(N)] = C_𝒞/N and [I^(N)] = ∑_𝒞∈ℳ_𝒜ℬ^Δ [ I_𝒞^(N)] = C/N.
  
  
  Effectively, Proposition <ref> ensures that we can consider the random variable I^(N) as an approximation for I with relative accuracy δ = 1/I√(C/N). Estimating the constant C is usually easy in practice: as long as sufficiently high powers of the integrand f( x) are integrable, we can also use Theorem <ref> to estimate [ f( x)].
  
  Variants of Algorithm <ref> are implemented e.g. as  <cit.> and as  <cit.>. 
  Both these implementations provide a variety of different ways to perform the preprocessing triangulation step which computes ℳ_𝒜ℬ^Δ. A dedicated tool to perform such a triangulation is  <cit.> which is also used internally in . 
  Subsequently, both programs use a version of the VEGAS algorithm <cit.> to numerically integrate I_𝒞[f] in eq. (<ref>) over the unit hypercube [0,1]^n-1 or, equivalently, to execute the inner loop of Algorithm <ref> for each individual cone 𝒞∈ℳ_𝒜ℬ^Δ. 
  
  We can estimate the computational complexity of the algorithm by counting the number of necessary evaluations of the function R_a/b( x). This is justified because we can assume that the runtime to evaluate R_a/b( x) overshadows the time it takes to compute a random vector ξ∈ [0,1]^n-1 and the value of x(ξ) ∈^n-1_>0 from it. 
  
  Therefore, the estimation step summarized in Algorithm <ref> needs N |ℳ_𝒜ℬ^Δ| evaluations to produce the estimate I^(N) for the integral in eq. (<ref>). Equivalently, as the relative accuracy δ≈ I^(N)/I of the resulting estimate is inverse proportional to √(N), the number of evaluations needed is proportional to δ^-2 |ℳ_𝒜ℬ^Δ| to achieve an estimate of δ accuracy.
  
  A severe bottleneck is the number of maximal cones |ℳ_𝒜ℬ^Δ| which tends to grow exponentially with growing dimension n of the problem. A particularly unsatisfying aspect of this bottleneck is that the value of the individual sector contributions I_𝒞 typically varies quite much in magnitude. Consequently, only a fraction of the geometric sector contributions in eq. (<ref>) are relevant for the overall integral I and much of the computational effort spent to estimate each of the integrals I_𝒞 is wasted. In the next section, we will explain how to overcome this bottleneck.
  
  
  
  § TROPICAL SAMPLING
  
  
  
  In summary, the strategy to overcome this problem is the following: instead of numerically integrating each of the sector integrals individually and eventually summing all the resulting numbers to obtain an estimate for the integral in eq. (<ref>), we can use a more `inclusive' Monte Carlo approach, where we evaluate both the individual integrals I_𝒞[R_a/b] in eq. (<ref>) and the sum over these integrals ∑_𝒞∈ℳ_𝒜ℬ^Δ I_𝒞[R_a/b] via Monte Carlo methods. This approach is much more efficient than the traditional one because there is a canonical way to perform importance sampling on the sum. That means that we can expose the individual sectors to our sampler `undemocratically' such that more important sectors are sampled more often than less important contributions.
  
  To do this it is convenient to define a tropically approximated version of the integral in eq. (<ref>):
  
      I^ = ∫__>0^n-1∏_i a^_i(x)^ν_i/∏_j b^_j(x)^ρ_jΩ.
  
  Such tropically approximated integrals have been considered as a simple avatar of period Feynman integrals <cit.> and identified to appear in the weak string coupling limit <cit.>. Moreover, this tropically approximated integral also gives rise to the canonical function of a polytope under certain conditions on the polynomials {a_i} and {b_j}, which has applications in the theory of scattering amplitudes <cit.>.
  
  It follows from Theorem <ref> with f( x) =1 that I^ is finite and that I^ > 0, provided that the conditions R1 and R2 on the polytopes 𝒜 and ℬ in Theorem <ref> are fulfilled. Moreover, the integrand in eq. (<ref>) is obviously positive for all x ∈^n-1_>0. Hence, we can define a probability distribution given by the differential form,
  
      μ^= 1/I^∏_i a_i^( x)^ν_i/∏_j b_j^(x)^ρ_jΩ,
  
  such that 1 =∫__>0^n-1μ^. The integral in eq. (<ref>) can now be written as,
  
      I = I^∫__>0^n-1 R_a/b( x) μ^,
  
  with R_a/b as defined in eq. (<ref>). 
  As μ^ is a properly normalized probability distribution on ^n-1_>0, we 
  can use Theorem <ref> to get a direct estimation algorithm for I from this, provided that we have a reasonably efficient way to sample from the distribution μ^. 
  
  
  
  
  
  
  
  
  
  
  Algorithm <ref> is not obviously simpler or more efficient than Algorithm <ref>, as the complicated part—generating a sample from the random distribution given by μ^—has been conveniently out-sourced. 
  
  A simple method to sample from μ^ is to again use a geometric sector decomposition. By Theorem <ref> the tropically approximated integral in eq. (<ref>) can be written as a sum,
  
      I^ = ∑_𝒞∈ℳ_𝒜ℬ^Δ I_𝒞^ with  I_𝒞^   = |( u^(𝒞,1), …,  u^(𝒞,n-1),) | /∏_k=1^n-1 u^(𝒞,k) w^(𝒞),
  
  where I_𝒞^ > 0 for all maximal cones 𝒞∈ℳ^Δ_𝒜ℬ. Hence, we can interpret I^_𝒞/I^ as a probability assigned to each cone 𝒞∈ℳ^Δ_𝒜ℬ and draw a random cone accordingly. Drawing a random sample from a finite discrete probability distribution is a classic problem. It can be solved in constant time independent of the number of possible outcomes if a table of the probabilities of the respective outcomes is appropriately preprocessed, for instance by using the alias method <cit.>. 
  Provided that we have generated such a table together with a table of appropriate values of w^(𝒞) and u^(𝒞,1), …,  u^(𝒞,n-1) we can execute the following algorithm:
  
  
  
  
  
  
  Algorithm <ref> generates a sample x ∈^n-1_>0, distributed as μ^ in eq. (<ref>).
  
  
  For any test function f:_>0^n-1→  and a random sample x ∈_>0^n-1 generated by Algorithm <ref>, we have
  
      [ f( x) ] = ∑_𝒞∈ℳ^Δ_𝒜ℬ I^_𝒞/I^∫_[0,1]^n-1 f(  x^(𝒞)(ξ)) ∏_k=1^n-1ξ_k.
  
  Using eq. (<ref>) and Theorem <ref> gives 
  
      [ f( x) ] = 1/I^∫_^n-1_>0∏_i a_i^( x)^ν_i/∏_j b_j^(x)^ρ_j f( x) Ω = ∫_^n-1_>0 f( x) μ^.
  
  
  
  
  
  To run both Algorithms <ref> and <ref> together we need N evaluations of the function R_a/b( x). Equivalently, we need proportional to δ^-2 evaluations to obtain an estimate I of δ accuracy. This is a significant improvement over Algorithm <ref> as the runtime is now independent of the number of sectors |ℳ^Δ_𝒜ℬ|. 
  
  It has to be stressed that this suggested direct comparison between Algorithm <ref> and the combination of the Algorithms <ref> and <ref> is flawed by the inherent difference in the respective proportionality factors for δ^-2 or equivalently, in the number of samples N that results in a given accuracy. In a situation, in which the sector integrals all contribute roughly the same value to the overall integral, Algorithms <ref> and <ref> offer no advantage over Algorithm <ref>. For practical applications the values of the sector integrals tend to differ heavily in magnitude, which makes Algorithms <ref> and <ref> favorable.
  
  Just as for Algorithm <ref> a preprocessing step needs to be performed for Algorithms <ref> and <ref>: the triangulation ℳ^Δ_𝒜ℬ and the associated table needs to be calculated. 
  
  This computation is also necessary to compute the normalization factor I^=∑_𝒞∈ℳ^Δ_𝒜ℬ I_𝒞^. In the best case, the time it takes to create such a table will be proportional to the number of sectors |ℳ^Δ_𝒜ℬ|, but we only need to compute this table once and can evaluate an arbitrary large number of samples afterwards. 
  
  Therefore, even though we are still effectively constrained by the dimension of the problem, which has to be small enough for the preprocessing step to be finished in a reasonable time, this constraint on the dimension is decoupled from the achievable accuracy. 
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  Recall that so far, we considered completely general integrals in eq. (<ref>). Although we already managed to accelerate the integration for the general case in comparison to the traditional approach, further improvements are possible if more specific properties of the integrand are used. Especially, integrals that come from physical applications are well-known to carry a very rich geometric structure, whose exploitation offers a whole new set of tools to improve numerical approximation methods. In the following, we will achieve a further improvement in runtime, memory requirement and overall complexity by using a specific structure which is exhibited by a large family of integrals. Integrals of this family appear in many contexts in high energy physics. This family consists of all integrals as in eq. (<ref>) where the Newton polytopes of the polynomials {a_i} and {b_j} are generalized permutahedra.
  
  
  
  § GENERALIZED PERMUTAHEDRA
  
  
  
  
  
  The permutahedron Π_n is an (n-1)-dimensional polytope in ^n. It can be defined as the convex hull of n! vertices determined by permutations in S_n:
  
      Π_n = {∑_σ∈ S_nλ_σ v^(σ) : ∑_σ∈ S_nλ_σ = 1  and λ_σ≥ 0 }⊂^n,
  
  where the vector v^(σ) = (σ(1), …, σ(n)) ∈^n encodes the permutation σ. The permutahedron is contained in the hyperplane Π_n ⊂{ v ∈^n:  v = n(n+1)/2 } and is full-dimensional within this hyperplane. The permutahedron Π_3 is depicted in Figure <ref>.
  The cones of maximal dimension in the reduced normal fan ℱ_Π_n/ of Π_n are labelled by permutations as well. They are of the form,
  
      𝒞_σ = { y ∈^n/: y_σ(1)≤…≤ y_σ(n)}
  
  such a domain is called a Weyl chamber. 
  It is not hard to see that these are simplicial cones as
  
      𝒞_σ = {∑_k = 1^n-1λ_k  u^(σ,k): λ_k ≥ 0 } with  u^(σ,k)_σ(i) =  -1    if k ≤ i
       0    else
  
  where we chose the set of representatives in ^n of the vectors in u^(σ,k)∈^n/ by fixing u^(σ,k)_σ(n) = 0 for all k∈{1,…,n-1}. The remaining cones of the reduced normal fan ℱ_Π_n can be constructed by taking arbitrary intersections of these cones. This fan is also called the braid arrangement fan. The reduced normal fan of Π_3 is illustrated in Figure <ref>.
  
  
  
  A polytope whose normal fan is a coarsening of ℱ_Π_n is a generalized permutahedron.
  
  Generalized permutahedra have a large number of remarkable properties <cit.>. E.g. 
  
  
  
  A generalized permutahedron 𝒢_z has the facet presentation 
  
      𝒢_z = { v ∈^n : ∑_i ∈ [n] v_i = z([n])  and ∑_i ∈ I v_i ≥ z(I)  for all  I ⊂ [n] },
  
  where [n] = {1,…, n} and z is a supermodular boolean function z:  2^[n]→ with z(∅) = 0. In fact every supermodular boolean function, that means z:  2^[n]→ with 
  
      z(A) + z(B) ≤ z(A ∩ B) + z(A ∪ B )  for all  A,B ⊂ [n],
  
  gives rise to a generalized permutahedron by the inequality description in eq. (<ref>)[We are using a different sign notation than <cit.>, but agree with <cit.>.].
  
  
  
  
  If both 𝒢_z_1 and 𝒢_z_2 are generalized permutahedra and 
  z_1(A) > z_2(A) for all non-empty A ⊊ [n] and z_1([n]) = z_2([n]), then
  𝒢_z_1⊂𝒢_z_2.
  
  
  By Theorem <ref> it follows immediately that 
  𝒢_z_1⊂𝒢_z_2.
  The inequalities in eq. (<ref>) are strict <cit.>. Therefore the statement follows.
  
  
  The Minkowski sum of two generalized permutahedra is again a generalized permutahedron (see for instance <cit.>):
  
  
  If both 𝒢_z_1 and 𝒢_z_2 are generalized permutahedra, then also their Minkowski sum 
  𝒢_z_12 = 𝒢_z_1 + 𝒢_z_2 is a generalized permutahedron with the boolean functions z_1,z_2,z_12: 2^[n]→ related by z_12(A) = z_1(A) + z_2(A) for all A ⊂ [n] = {1,…, n}.
  
  
  A vector v∈𝒢_z which maximizes all linear functionals in a Weyl chamber 𝒞_σ is a vertex of 𝒢_z. 
  This gives a canonical map from permutations σ∈ S_n to the vertices of a generalized permutahedron. We can use a result of Fujishige and Tomizawa to explicitly construct this map:
  
  
  
  
  If z is a supermodular function z:  2^[n]→, σ∈ S_n a permutation and w^(σ,z)∈^n is the vector given component-wise by
  
      w^(σ,z)_σ(k) = z(A_k^σ) - z(A_k-1^σ)  for all  k ∈ [n],
  
  where A_k^σ = {σ(1), …, σ(k)}⊂ [n]={1,…,n}, then
  w^(σ,z) is a vertex of the generalized permutahedron 𝒢_z and
  
      y w^(σ,z) = max_ v ∈𝒢_z y v for all  y ∈𝒞_σ,
  
  where 𝒞_σ is a Weyl-chamber in the braid arrangement fan as defined in eq. (<ref>). 
  
  
  
  
  
   §.§ Tropical sampling for generalized permutahedra
  
  In general, it is necessary to compute a triangulation of the reduced refined normal fans of the 𝒜 and ℬ polytopes to perform the procedure described in Section <ref>. This cumbersome computation can be circumvented if the 𝒜 and ℬ polytopes are generalized permutahedra. In this case, there is an especially simple way to sample from the associated tropical measure μ^ defined in eq. (<ref>) without the need for an explicit triangulation as required for Algorithm <ref>.
  
  From now on, we will therefore assume that the polytopes 𝒜 and ℬ are both generalized permutahedra. This implies by Theorem <ref> that there are unique boolean functions z_𝒜 :  2^[n]→ and z_ℬ analogously which describe these polytopes. Using these functions and the properties of generalized permutahedra introduced above, we can state
  
  
  
  Let 𝒜 and ℬ be the polytopes defined in Theorem <ref>. 
  If 𝒜 and ℬ are generalized permutahedra with associated boolean functions z_𝒜,z_ℬ:  2^[n]→ which 
  fulfill the requirements R1 and R2 of Theorem <ref>,
  
  then we can write the integral 
  
      I[f] = ∫__>0^n-1∏_i a_i^( x)^ν_i/∏_j b_j^(x)^ρ_j f( x) Ω
  
  as a sum 
  I[f] = ∑_σ∈ S_n I_σ[f] with
  
      I_σ[f]    = 1/∏_k=1^n-1 r(A^σ_k) ∫_[0,1]^n-1 f (  x^(σ)(ξ) ) ∏_k=1^n-1ξ_k,
  
  where 
  
  
    * f:_>0^n-1→ is a measurable homogeneous function of degree 0,
  
    * A_k^σ = {σ(1), …, σ(k)}⊂ [n]={1,…,n},
  
    * r(A) = z_𝒜(A) - z_ℬ(A),
  which fulfills r(A) > 0 for all non-empty proper subsets A ⊊ [n] and
  
    * x^(σ)(ξ) ∈(𝒞_σ) is given component-wise by x_σ(k) = ∏_i=k^n-1ξ_i^1/r(A^σ_i) and x_σ(n) = 1.
  
  
  
  This theorem is a specialization of Theorem <ref> to the generalized permutahedron case. The braid arrangement fan defined in eq. (<ref>) provides an appropriate reduced simplicial fan.
  By Lemma <ref> we have vertices w^(σ,z_𝒜)∈𝒜 and w^(σ,z_ℬ)∈ℬ such that
  y w^(σ,z_𝒜) = max_ v ∈𝒜 y v
  and 
  y w^(σ,z_ℬ) = max_ v ∈ℬ y v
  for all σ∈ S_n and y ∈𝒞_σ.
  Using the explicit representatives of the generators u^(σ,1), …, u^(σ,n-1) of the cone 𝒞_σ from eq. (<ref>) together with Lemma <ref> gives
  u^(σ,k) w^(σ,z_𝒜) = -z_𝒜(A_k^σ) and
  u^(σ,k) w^(σ,z_ℬ) = -z_ℬ(A_k^σ).  It follows from this and Lemma <ref> that u^(σ,k) w^(σ,z_ℬ) -  u^(σ,k) w^(σ,z_𝒜) = z_𝒜(A_k^σ) - z_ℬ(A_k^σ) > 0 for all σ∈ S_n and k∈{1,…,n-1} which implies r(A) > 0 for all non-empty A⊊ [n].
  From the form of the u^(σ,k) vectors in eq. (<ref>) it is obvious that 
  |( u^(σ,1), …, u^(σ,n-1),)|=1.
  
  
  Theorem <ref> ensures that we can proceed as above and perform the Monte Carlo Algorithms <ref> and <ref> just as in the general case. It is clear that the preprocessing step will be straightforward as generalized permutahedra come with an appropriately simplicial fan `built in'. Algorithm <ref> requires us to generate a table of size n! as we need one entry for each cone in the braid arrangement fan.  For this algorithm to be applicable in a computationally feasible way that table needs to be stored in the memory of the computer. Hence the naive algorithm is only practically applicable for relatively small values of n. 
  
  However, a further significant improvement can be achieved: it is not necessary to store an entry for each permutation in a table. If a small additional computation for each sampled point is performed, a table of size proportional to 2^n suffices. We will describe this specialized version of Algorithm <ref> in the rest of this section.
  
  First observe that the overall normalization factor needed to apply Algorithm <ref> is given by 
  
      I^   = ∑_σ∈ S_n I^_𝒞_σ with  I^_𝒞_σ = 1/∏_k=1^n-1 r(A^σ_k) ,
  
  where r(A) = z_𝒜(A) - z_ℬ(A) for all non-empty A ⊊ [n]. This equation is just eq. (<ref>) specified using Theorem <ref> to the generalized permutahedron case. 
  For the following considerations it will be convenient to declare r(∅) = 1, which opens the way towards the following generalization that promotes I^ to a boolean function on 2^[n]:
  
  
  For a boolean function r: 2^[n]→ with r(∅) = 1 and r(A) > 0 for all non-empty A ⊊ [n], we define the boolean function J_r:  2^[n]→_>0 recursively as
  
      J_r(A) = ∑_e∈ AJ_r(A∖ e)/ r(A∖ e)  for all non-empty  A ⊂ [n]  where  J_r(∅) = 1.
  
  
  
  
  If r(A) = z_𝒜(A) - z_ℬ(A) for all non-empty A ⊊ [n] and r(∅) =1, then
  I^ = J_r([n]).
  
  
  
  We will prove that 
  J_r(A) = ∑_σ : [m] → A1/∏_k=1^m-1 r(A^σ_k), where the sum is over all bijections σ: [m] → A.
  Fixing such a bijection is equivalent to fixing a pair (e,μ) of an element e ∈ A and a bijection μ:[m-1] → A∖ e. Decomposing the sum in this way and using eq. (<ref>) gives the statement.
  
  
  This recursive method to calculate the normalization factor I^ might also be useful in other contexts. For instance, this can be used to calculate the volume of the polar dual of a generalized permutahedron fairly efficiently. In the context of scattering amplitudes this recursion can also be used to calculate the canonical form of a generalized permutahedron.
  
  
  If we prepare a table of the values J_r(A) and r(A) for all A ⊂ [n], we can run the following algorithm:
  
  Note that the probability distribution p_e = 1/J_r(A)J_r(A∖ e)/r(A∖ e) over the elements e∈ A is properly normalized due to Definition <ref>.
  
  
  If r(A) = z_𝒜(A) - z_ℬ(A) for all non-empty A ⊊ [n], r(∅) =1 and J_r is the boolean function given in Definition <ref>, then
  Algorithm <ref> generates a sample x ∈^n-1_>0, distributed as μ^ in eq. (<ref>) in the generalized permutahedron case.
  
  
  For any test function f:_>0^n-1→ and a random sample x ∈_>0^n-1 generated by Algorithm <ref>, 
  
      [ f( x) ] = ∑_e_n ∈ A_n1/J_r(A_n)J_r(A_n∖ e_n)/r(A_n∖ e_n)…∑_e_1∈ A_11/J_r(A_1)J_r(A_1∖ e_1)/r(A_1∖ e_1)∫_[0,1]^n-1f( x(ξ)) ∏_k=1^n-1ξ_k,
  
  where we gave distinguished subscripts to the numbers e and sets A in the reverse order in which they appear in Algorithm <ref> and x(ξ) is component-wise x_e_k = ∏_i=k^n-1ξ_i^1/r(A_i). We identify A_k ∖ e_k = A_k-1. The terms J_r(A_k∖ e_k) telescope, J_r(∅) = r(∅) = 1 and we get
  
      [ f( x) ] = 1/J_r(A_n)∑_e_n ∈ A_n…∑_e_1∈ A_11/r(A_n) ⋯ r(A_1)∫_[0,1]^n-1f( x(ξ)) ∏_k=1^n-1ξ_k.
  
  The sum can be written as a sum over all permutations in σ∈ S_n and A_k = A_k^σ.
  The statement follows from Proposition <ref> and Theorem <ref>.
  
  
  Algorithm <ref> allows us to integrate any integral of the form in eq. (<ref>) via Monte Carlo quadrature without actually performing any complicated triangulation or non-trivial sector decomposition step if the polytopes 𝒜 and ℬ are generalized permutahedra. Compared to the naiver approach where a table of size n! is needed, only a table of size 2^n is required. The complexity of the preprocessing step is similarly reduced as the recursion in Definition <ref> gives an efficient way to calculate all the necessary constants: The table for J_r(A) can be calculated in (n 2^n) steps. 
  
  All this achieves not only a huge improvement in the required runtime and memory of the algorithm, but also significantly reduces the complexity of the overall algorithm. Triangulating an n-dimensional polytope is an involved algorithmic task. Circumventing this triangulation with the approach above makes it straightforward to implement an efficient integration algorithm. A detailed example is given in the following section. 
  
  
  
  § FEYNMAN INTEGRALS
  
  
  
  
  
  
  
  
  
  
  
  A scalar Feynman integral associated to a Feynman graph G with E edges and V vertices in parametric representation in D-dimensional Euclidean space can be written as,
  
      I_G    = ∫_^E-1_> 0∏_e x_e^ν_e/Ψ_G( x)^D/2( Ψ_G( x)/Φ_G( x))^ω(G)Ω,
  
  which depends on the edge weights ν_1, …, ν_E, which we will assume to be positive and real. The superficial degree of divergence ω(G) is given by ω(G) = ∑_eν_e - ℓ(G) D/2, where ℓ(G) is the number of loops of G (i.e. the first Betti number of G). The Kirchhoff-Symanzik polynomials Ψ_G and Φ_G are homogeneous of degree ℓ(G) and ℓ(G)+1 in the x_e variables. Obviously, the integral I_G is a specific instance of an integral of the form in eq. (<ref>). To simplify the notation, we omitted a prefactor of Γ(ω_G)/∏_e Γ(ν_e), which is usually included in the definition of scalar Feynman integrals. See for instance <cit.> for details on this representation of Feynman integrals.
  
  A complete finiteness proof of the Euclidean space Feynman integral I_G together with an analysis of its analytic continuation properties in the ν_1,…,ν_E parameters has been achieved by Speer <cit.>. More recently, Brown <cit.> showed that there is a canonical way to associate the integral I_G to a motivic avatar, which can be thought of as a specific representation of a conjectured cosmic Galois group. This group suggests the existence of a coaction principle which relates different Feynman integrals in a highly non-trivial way and it allows to analyse Feynman integrals with a whole new toolkit of homological methods and representation theory. 
  
  For our endeavour to merely evaluate the integrals I_G, we can make use of parts of Brown's analysis <cit.> to ensure that the relevant polytopes associated to the Ψ_G and Φ_G polynomials are generalized permutahedra. We will start by giving some additional details on these polynomials including an efficient way to evaluate them.
  
  
  
  
  
   §.§ Symanzik polynomials
  
  Explicitly, the polynomials can be expressed as sums over spanning trees T_1 and spanning 2-forests (spanning forests with two connected components) T_2 of the graph G:
  
      Ψ_G( x)    = ∑_T_1∏_e∉ T_1 x_e    Φ_G( x)    = ∑_T_2 p(T_2)^2 ∏_e∉ T_2 x_e + Ψ_G ∑_e x_e m_e^2,
  
  where  p(T_2) is the total momentum flowing between the two components of the 2-forest T_2. Only the Φ_G( x) polynomial depends on the external physical parameters: a set of momenta p^(1), …,  p^(V)∈^D incoming into each of the vertices and a set of masses m_1, …, m_E ∈ associated to the edges of the graph. These polynomials can also be written in terms of the weighted V× V Laplace matrix of the graph (see for instance <cit.>), which is component-wise
  
      L_v,w   =  -x_e^-1   if there is an edge e between v and w
      ∑_e  incident to  v x_e^-1   if v=w
       0    else
  
  This matrix is only positive semi-definite whereas the reduced Laplacian L_v,w( x), which is given by an arbitrary leading principle minor of the matrix L_v,w( x), is positive definite. The Symanzik polynomials can be written as
  
      Ψ_G( x)    = (∏_e x_e ) ( L )    Φ_G( x)    = Ψ_G( ( P^T L^-1 P) + ∑_e x_e m_e^2 ),
  
  where P is the (V-1) × D matrix, given row-wise by the incoming momenta, p^(v)∈^D: P_v,μ = p^(v)_μ where v=1,…,V-1 and μ is a D-dimensional spacetime index. Note that due to momentum conservation no information is lost when only V-1 of the V incoming momenta are used.
  
  The second representation of the Symanzik polynomials in eq. (<ref>) is more suitable for numerical evaluation than eq. (<ref>). The number of spanning trees of a graph grows exponentially with the number vertices V <cit.> and the evaluation of the expressions in eq. (<ref>) quickly becomes intractable when the graph gets large. The evaluation of the determinant with the other matrix operations in eq. (<ref>) is computationally much more favorable: With a Cholesky decomposition of the matrix L( x) both the value of Ψ_G( x) and Φ_G( x) can be immediately calculated. Computing the Cholesky decomposition of a (V-1)× (V-1)-matrix takes (V^3) time. Due to the special structure of the problem—the matrix L( x) being the reduced Laplace matrix of a graph—there even exists a nearly linear time approximation algorithm <cit.> to compute this decomposition.
  
  To give a precise account on the Newton polytopes of the Symanzik polynomials we need some additional notation from <cit.> for subgraphs of Feynman graphs. A subgraph γ⊂Γ is equivalent to a subset of edges of the graph Γ. The set of subgraphs is therefore isomorphic to the set 2^[E] and we will identify boolean functions 2^[E]→ with functions defined on the set of subgraphs of the graph G. Just as for G, we will denote the first Betti number of a subgraph (i.e. the number of loops) as ℓ(γ). A subgraph γ⊂ G is called mass-momentum-spanning (m.m.) in G if the second Symanzik polynomial of the contracted graph G/γ vanishes Φ_G/γ = 0. Mass-momentum-spanning graphs can also be defined combinatorially as subgraphs that contain all massive edges and one connected component which connects all vertices with non-zero incoming momentum. See <cit.> for details on these types of subgraphs.
  
  
  
  If we restrict to Euclidean and non-exceptional kinematics, then the Newton polytope of Ψ_G and Φ_G are generalized permutahedra. A facet presentation of these polytopes is given by the supermodular functions 
  
      z_Ψ_G(γ)    = ℓ(γ) 
       z_Φ_G(γ)    = ℓ(γ) +1     if γ is m.m. in  G 
      ℓ(γ)     else
  
  for all subgraphs γ⊂Γ.
  
  See <cit.> for a definition of non-exceptional or generic kinematics. Briefly, this condition ensures that there is 
  no non-trivial combination of the external momenta that adds up to 0. 
  It is worth remarking that this condition is not necessary if the combinatorial concept of mass-momentum-spanning is slightly generalized while keeping the equivalence Φ_G/γ = 0 ⇔mass-momentum-spanning. With this generalization it is sufficient to 
  require that any external momenta are non-zero. 
  
  Theorem <ref> has been proven by Schultka <cit.> using results from Brown <cit.>. 
  
  The first statement that the Newton polytopes of Ψ_G and Φ_G are generalized permutahedra can be traced back to Hepp <cit.> and Speer <cit.>, who realized that a complete ordering of the integration parameters in eq. (<ref>) is sufficient to capture the relevant singularities of parametric integrals in the Euclidean non-exceptional case. See also <cit.> for a comparison of this viewpoint with modern sector decomposition techniques.
  
  
  
  
  
  
  
  The form of the boolean functions z_Ψ_G and z_Φ_G follows directly from the factorization laws <cit.>, <cit.> and <cit.> of the Ψ_G and Φ_G polynomials. Their supermodularity follows from the argument in <cit.> after Corollary 4.12.
  
  
  
  It was implicitly proved by Panzer <cit.> that the Newton polytope of Ψ_G is a generalized permutahedron using an elegant argument based on Kruskal's algorithm <cit.>. This argument can also be generalized to the Φ_G polynomials by an extension of Kruskal's algorithm to minimal 2-forests. 
  
  
  
  Theorem <ref> is also of interest in a different context: Generalized permutahedra have a universal property with respect to their Hopf monoid structure. Feynman graphs carry a Hopf algebra structure which is deeply intertwined with renormalization <cit.> and encodes the singularity structure of the integrand <cit.>. The relationship between these two structures remains to be explored.
  
  
  
  For general non-Euclidean kinematics, the Newton polytope of Φ_G is not a generalized permutahedron. An explicit counterexample is given in <cit.>. We emphasize that the general tropical sampling algorithm introduced in Section <ref> still applies. The caveat is that an explicit triangulation has to be computed in contrast to the generalized permutahedron case where no explicit triangulation is necessary.  
  
  
  
  
   §.§ Tropical Monte Carlo quadrature of Euclidean Feynman integrals
  
  To perform the generalized permutahedron tropical Monte Carlo routine from Section <ref> on the parametric Feynman integral in eq. (<ref>), we still have to ensure that the numerator monomials ∏_e x_e^ν_e are generalized permutahedra. This is of course trivial, as the Newton polytope of a monomial is zero-dimensional and its normal fan is trivial. The braid arrangement fan is automatically a refinement of this fan and the conditions for Definition <ref> are fulfilled. The facet presentation of these 0-dimensional polytopes associated to the Newton polytope of the polynomial p_e( x) = x_e in the form of Theorem <ref> is given by the boolean function z_p_e(γ) = 1 if e ∈γ and z_p_e(γ) = 0 if e∉γ for all subgraphs γ.
  
  Because we assume that the edge weights ν_e, the dimension D and the superficial degree of divergence ω(G) are real, 
  we have
  
      𝒜   = ∑_eν_e _p_e + ω(G) _Ψ_G   ℬ   = 1/2 D _Ψ_G + ω(G) _Φ_G      if ω(G) ≥ 0 and 
      𝒜   = ∑_eν_e _p_e +  (-ω(G)) _Φ_G   ℬ   = 1/2 D _Ψ_G +  (-ω(G)) _Ψ_G      if ω(G) < 0.
  
  we can define the boolean function r_G: 2^[E]→ as in Theorem <ref>,
  
      r_G(γ)    = z_𝒜( γ ) - z_ℬ( γ ) 
         = ∑_eν_e z_p_e(γ) - D/2 z_Ψ_G(γ) + ω(G) (z_Ψ_G(γ) - z_Φ_G(γ) ) 
         = ∑_e ∈γν_e - D/2ℓ(γ) - ω(G) δ_m.m.(γ)  for all non-empty γ⊂ G,
  
  where we used Lemma <ref> and Theorem <ref> and where δ_m.m.(γ) = 1 if γ is mass-momentum-spanning and 0 otherwise. Note that up to the δ_m.m.-term the function r_G(γ) is equal to the superficial degree of divergence ω(γ) of a subgraph.
  
  For Euclidean kinematics, the polynomials Ψ_G and Φ_G have only positive coefficients. Therefore, they are completely non-vanishing on _>0^E-1. We can apply Theorem <ref> independently of the sign of ω(G) and find that the parametric integral in eq. (<ref>) is convergent if r_G(γ) > 0 for all non-empty proper γ⊊ G by using Corollary <ref> which implies 𝒜⊂ℬ in this case. In fact, it is sufficient that r_G(γ) > 0 holds for all proper motic subgraphs γ as defined in <cit.> for I_G to be convergent. 
  
  As defined in eq. (<ref>) the tropical differential form associated to I_G is
  
      μ^_G = 1/I_G^∏_e x_e^ν_e/Ψ_G^( x)^D/2( Ψ_G^( x)/Φ_G^( x))^ω(G)Ω,
  
  with an appropriate normalization factor I^_G such that 1 = ∫_^n-1_>0μ^_G. For ω(G)=0 this normalization factor is a certain invariant of the graph G which has been studied by Panzer <cit.>. This invariant is the Hepp-bound. The Hepp-bound is independent of the physical parameters encoded in the masses and external momenta. It mirrors many properties of the period, which is given by the integral in eq. (<ref>) in the same special case ω(G) = 0. The period is another graph invariant which has interesting number theoretical properties <cit.>.
  
  By Definition <ref>, the normalization factor can be generalized to a subgraph function J_G :  2^[E]→ which is determined by the recursion
  
      J_G(γ) = ∑_e∈γJ_G(γ∖ e)/ r_G(γ∖ e)  for all non-empty γ⊂ G  with  J_G(∅) = 1  and  r_G(∅) = 1.
  
  The actual normalization factor is recovered for γ = G, i.e. I_G^ = J_G(G)
   by Proposition <ref>. With a precalculated table of the values r_G(γ) and J_G(γ) for all γ⊂ G, Algorithm <ref> provides an efficient way to sample from the distribution given by the differential form μ^_G on ^n-1_>0. Using this sampling algorithm we can obtain estimates for the parametric Feynman integral eq. (<ref>) by the standard Monte Carlo procedure from Theorem <ref> or equivalently Algorithm <ref>.
  
  
  
   §.§ Expansions in regularization parameters
  
  Often not only the integral in eq. (<ref>) is of interest, but also the Taylor expansions of the parameters D and ν_e around specific points. Very important is the ε-expansion of the parametric Feynman integral in eq. (<ref>) in the context of dimensional regularization. Effectively, such an expansion results in integrals of the form
  
      I_G    = ∫_^E-1_> 0∏_e x_e^ν_e/Ψ_G( x)^D/2( Ψ_G( x)/Φ_G( x))^ω(G)(∏_elog^k_e (x_e) ) log^s (Ψ_G) log^t (Ψ_G/Φ_G) Ω,
  
  for some set of integers s,t∈ and k_1,…,k_E ∈.
  The estimation of this generalization is also possible using Algorithm <ref> or Algorithm <ref>. Using
  
      I_G = I^_G ∫_^E-1_> 01/(Ψ_G( x)/Ψ^_G( x))^D/2( Ψ_G( x)/Ψ_G^( x)/Φ_G( x)/Φ_G^( x))^ω(G)×
      ×(∏_elog^k_e (x_e) ) log^s (Ψ_G) log^t (Ψ_G/Φ_G) μ^,
  
  gives the desired estimate. A caveat is that the integrand is not bounded anymore, as the logarithms will exhibit singularities at the boundary of the integration domain. This is not a severe problem, as these singularities are square integrable and Theorem <ref> may still be applied. 
  
  
   §.§ Some experimental results
  
  
  
  A proof-of-concept  implementation of this algorithm, which evaluates general Euclidean Feynman integrals, is available on the author's personal web page[https://michaelborinsky.commichaelborinsky.com] and in the ancillary files to the arXiv version of this article. The algorithm has been tested on various graphs from φ^4-theory in four dimensions, which have been generated using tools from <cit.>. To illustrate the performance of the algorithm a benchmark is given in Table <ref>. The benchmark has been performed on a single core of an  processor. The columns E and ℓ(G) show the number of edges (equivalently the dimension of the integral +1) and the corresponding number of loops of the underlying φ^4-graph. The column σ_I/I gives the relative standard deviation of the samples, i.e. if δ^-2·σ_I/I samples are drawn, then a relative accuracy δ can be expected from the resulting estimate. Up to this expected accuracy, all obtained estimates are consistent with the available analytic results from <cit.>. The implementation has also been checked using numerical calculations of non-φ^4 graphs with non-trivial masses and kinematics performed with  <cit.>. 
  
  Recall that the algorithm can be applied to arbitrary D-dimensional scalar Feynman integrals with arbitrary kinematics in the Euclidean regime and the benchmark results can expected to be representative for the evaluation of all such graphs with the same number of edges. The choice for φ^4-theory and D=4 is practical because much analytic data is available even at high loop orders, which allows for convenient checks of the numerical estimates. 
  
  As can be seen from the table, the number of samples per second decreases slowly with the loop order or equivalently the dimension of the problem. The necessary time for the preprocessing step on the other hand depends exponentially on the dimension.  For example: it takes 2.5 CPU-seconds to evaluate a graph with 10 edges and general kinematics up to δ= 10^-3 relative accuracy. The necessary time for the preprocessing step of 6.6 · 10^-4 s is negligible and the memory requirements of 16  KB insignificant. It takes 20 CPU-seconds to evaluate a graph with 20 edges and general kinematics up to the same relative accuracy. The time for the preprocessing step is 1 second and the memory requirements of 16  MB are still very manageable.  Similarly, it takes about 2 CPU-minutes to evaluate a Feynman graph with 30 edges up to this accuracy, after the preprocessing step has been performed. At this point this preprocessing step unfortunately already takes about 30 minutes and 16  GB of RAM are necessary.
  
  
  
  
  The evaluation step of the algorithm is fully parallelizable and the preprocessing step partially. The memory requirements can be reduced in the special case ω(G)=0 or by using a more efficient storage of the relevant constants. The overall picture of exponentially growing memory demands and an exponential time for the preprocessing step will not change without modifying the algorithm significantly.
  
  
  An interesting example of a φ^4-graph in D=4, whose evaluation was not approachable by any previously existing techniques, is the graph in Figure <ref>. It is one of the smallest graphs in φ^4-theory whose period is not a linear combination of multiple zeta values or multiple polylogarithms at roots of unity. This has been proven in <cit.> for a graph which is equivalent with respect to its period by the completion identity <cit.>.
  Sampling 10^12 points in about 24 hours on 54-CPU-cores results in the following estimate for the period of this graph,
  
      I_G    = ∫_^E-1_> 0∏_e x_e/Ψ_G( x)^2Ω≈ 422.9610 ± 0.0009 .
  
  
  
  
  § FURTHER RESEARCH DIRECTIONS
  
  
  
  
    * 
  (Markov chain Monte Carlo based sampling)
  
  The tropical Monte Carlo algorithms are still very limited in terms of the complexity of the integrals to which they apply, because of the cumbersome preprocessing step that has to be performed for each integral. 
  To overcome this bottleneck without relying on special structures of the integrals, it would be necessary to find a more efficient way to sample from μ^ than Algorithm <ref> or Algorithm <ref>, while also having access to the normalization factor I^. Eventually, one has to settle with a still relatively slow algorithm for this task, as we have a `no-go Theorem' in a special case: if all the numerator polynomials are monomials, i.e. a_i( x) = x_i, then I^ corresponds to the volume of a certain polytope. Computing or approximating the volume of a general n-dimensional polytope is a task that cannot be performed deterministically in polynomial time <cit.>. A workaround is to use a non-deterministic algorithm  for both the computation of the normalization factor I^ and to obtain samples from μ^. There are many highly advanced Markov chain Monte Carlo algorithms that have been developed to perform exactly this task (see for instance <cit.> and the references therein). It is very plausible that adapting these polytope integration and sampling algorithms to our algebraic integral quadrature application should result in the sought after polynomial time algorithm for algebraic and Feynman integral evaluation.
  
  
  
    * (Physical integration regions and components of coamoeba)
  The last condition in Theorem <ref> is closely related to the coamoeba of the set of polynomials {b_j}. 
  If a polynomial p∈[x_1,…,x_n] has zero locus Z_p = { z ∈ (∖{ 0 })^n: p( z) = 0 }, then the coamoeba of p is the image of Z_p under the coordinate-wise complex -function: 𝒜_p' = Arg( Z_p ) ⊂ [0,2π]^n. The coamoeba is related to the amoeba which goes back to Gelfand, Kapranov and Zelevinsky <cit.> and has numerous applications in tropical geometry. By a result proven independently by Johansson <cit.> and Nisse, Sottile <cit.>, a polynomial is completely non-vanishing if the origin is not in the closure of its coamoeba 0 ∉𝒜̅_p'. In <cit.> it was shown via Cauchy's theorem that the integration cycle ^n-1_>0 of the integral in eq. (<ref>) can be replaced with the Arg^-1(θ) as long as θ and 0 lie in the same connected component of the intersection of the coamoeba of the denominator polynomials. A similar argument works for the projective version of generalized Euler-Mellin integrals which was considered here.
  
  A strikingly reminiscent procedure is necessary while evaluating Feynman integrals with kinematics in Minkowski space. The necessary analytic continuation in this case is governed by the i ε-prescription, which ultimately results from causality and unitarity constraints on the amplitude <cit.>. Formulating this procedure in terms of a canonical choice of a component in the respective coamoeba would result in a canonical analytic continuation procedure in the Minkowski case.
  See also 
  <cit.> where related observation regarding parametric Feynman integrals and coamoeba have been made.
  
  
    * (Further acceleration of the algorithms by using more structures)
  In the ω(G)=0 case the normalization factor of the μ^_G distribution for the parametric Euclidean Feynman integral in eq. (<ref>) reduces to the Hepp-bound studied by Panzer <cit.>. He gave more efficient ways to compute this normalization factor I^_G which likely can also be used to sample from the μ^_G distribution more efficiently. Moreover, these more elaborate ways to compute the Hepp-bound can probably be generalized to deal with the ω(G) ≠ 0 case using results from Brown <cit.>.
  
  In a broader sense, an extension of Algorithm <ref> beyond the generalized permutahedron case could be possible. 
  Especially attractive would be an extension which includes the interesting Minkowski space Feynman integral case. Further analysis of the relevant structures for the tropical geometric framework, starting for instance with the explicit counterexample in <cit.>, could lead to an appropriate refinement of the braid arrangement fan. Such a refinement could lead to a direct generalization of the generalized permutahedron sampling Algorithm <ref>, which would make the integration of high dimensional Feynman integrals (i.e. with ∼ 30 edges) also possible in the Minkowski regime.
  
    * (BPHZ renormalization)
  
  As mentioned above, Feynman integrals as the one in eq. (<ref>) with non-integrable singularities are often interesting. A common approach to deal with these singularities is to subject the integral to an analytic continuation procedure before any numerical integration is performed. Ultimately, these singularities have a well-studied physical origin and are handled via renormalization. The momentum BPHZ renormalization scheme takes care of these singularities before any integration is performed. This renormalization scheme can be implemented on the level of the parametric integrand <cit.>. Such an implementation would make the analytical continuation step in the ν_e edge weights and the dimension D unnecessary.
  
    * (Estimates for large loop order β-functions)
  The estimates that can be obtained using the proof-of-concept implementation of Feynman graph integrals up to loop order 17 in φ^4-theory can immediately be used for a numerical estimation of the β-function up to this loop order. This has phenomenological applications for the calculation of critical exponents for various complex systems and even works without the need for an analytic continuation as it has been observed that the non-primitive contributions to the β-function become negligible with sufficiently large loop order. The approach can be further amplified by making use of the observed Hepp-bound - period correlation which has been applied by Panzer and Kompaniets <cit.> to obtain estimates of the φ^4-theory β-function up to order 13. 
  
  Even if it is not possible to evaluate all necessary Feynman diagrams individually for the respective loop order, a numerical approach could be sufficient to gain enough insights on the distribution of the value of these integrals. If such statistical knowledge is available, the numbers of (renormalized) Feynman diagrams are sufficient to extrapolate values for the entire β-function contribution <cit.>.
  
  Such an approach naturally extends to a question for the inherently large-order regime: suppose that G is a random 1PI φ^4-graph without subdivergences and ℓ(G) loops. Is there a limiting distribution 
  
      lim_ℓ(G) →∞ C_ℓ(G)^-1∫_^E-1_> 0∏_e x_e^ν_e/Ψ_G( x)^D/2( Ψ_G( x)/Φ_G( x))^ω(G)Ω,
  
  with an appropriate normalization constant C_L for each loop order L=ℓ(G) and if yes, what does it look like?
  The analysis <cit.> gives some positive indication for the existence of such a distribution.
  
  An overall normalization constant for such a distribution, which normalizes its expectation value to one, can be calculated using instanton methods <cit.> and renormalized graph counting <cit.>, as was pointed out by Panzer <cit.>:
  
      C_L = 4 e^-3γ_E/√(2 π) A^6 L^5/2(3/2)^L+3,
  
  where γ_E is the Euler–Mascheroni constant and A is the Glaisher–Kinkelin constant.
  
  An exhaustive statistical analysis using the algorithms from this article should give further indication for or against the existence of such a limit. A combination with analytic combinatorial methods for Dyson-Schwinger equations might lead to an explicit form of a limit distribution of Feynman integrals <cit.>.
  
    * (Phase-space integration)
  Simple phase space integrals, which are another type of integrals necessary for particle physics phenomenology, also fall under the category of integrals in eq. (<ref>). For more elaborate phase space integrals more complicated non-simplicial integration domains are necessary. It is possible that the algorithms discussed in this article may be extended to these more complicated domains. Writing the phase space integrals in terms of kinematic variables as in <cit.> and using a geometric subtraction scheme for the infrared singularities <cit.> could be instrumental for this extension. 
  
    * (Tropical sampling applied to sums of Feynman diagrams)
  The general tropical sampling algorithm described in Section <ref> is made possible by a well-calculated emancipation from the rigid concept of sectors as parts of the integral, which each have to be attacked individually. 
  
  Following a line of thought from <cit.>, we can say that a similar but much stronger bias exists on the level of the amplitude. The time-honoured approach to amplitude calculation is to write it as a sum over Feynman graphs with the same number of loops L (pictorially in disregard of renormalization and the explicit form of the integrals),
  
      A_L = ∑_ℓ(G) = L1/| G|∫( …) Ω,
  
  and evaluate the indicated integrals one by one.
  There are promising indications that there is a superior structure which can be `triangulated' into Feynman integrals in an appropriate sense yielding a sum as the one above, similar to the sector decomposition approach, where an individual integral is decomposed in terms of the triangulation of the respective normal fan. 
  
  An especially suggestive candidate for such a superior object in the case of scalar quantum field theories is Outer space <cit.> and its quotient formed under the action of Out(F_n) which is the moduli space of graphs. For instance, unitarity and branch cut properties of Feynman integrals can be understood using Outer space <cit.>. This space can be seen as a tropical analogue of Teichmüller space <cit.> and the moduli space of curves, which holds a similar superior role in string theory. Recently, quantum field theory inspired techniques have been successfully applied in the theory of Outer space <cit.>.
  
  A problem to overcome for such an approach are the UV-divergences that naturally appear in renormalizable QFT calculations. It is well-known how such divergences can be handled both on the amplitude or on a per integral level <cit.> and also mathematically these divergences are quite well understood, even in the large-order regime <cit.>. These divergences would also appear in a geometric setting for the amplitude and dealing with them would mean to work on a certain compactification of Outer space and the moduli space of graphs. One such compactification has been constructed by Berghoff <cit.> (see also <cit.>), which might be usable for the numerical evaluation of amplitudes.
  
    * (Quasi Monte Carlo)
  
  State of the art implementations for numerical Feynman integral integration employ quasi Monte Carlo methods <cit.> for the actual integration of the sector integrals instead of traditional Monte Carlo methods. This has the simple and obvious advantage of a significantly increased rate of convergence. The disadvantage of the quasi Monte Carlo approach is that it is mathematically much more challenging to handle. It is plausible that the algorithms introduced in this article can be accelerated using quasi Monte Carlo methods. A challenge will be the handling of the mixture of discrete and continuous probability distributions in the sampling algorithms from Sections <ref> and <ref>.
  
    * (Systematic tropical expansions)
  The tropical approximation p^ can be interpreted as a certain limit as shown in Section <ref>. It is natural to ask if one can interpret p^ as the `zeroth' order in a systematic expansion and it is plausible that a systematic improvement of the technique can be obtained this way. The ultimate aim of such investigations would be an efficient approximation scheme that gets by without a final Monte Carlo step and immediately yields a deterministic result. 
  
  
  
  
  
  
  
  
  
  100
  
  aguiar2017hopf
  M. Aguiar and F. Ardila, Hopf monoids and generalized permutahedra,
    https://arxiv.org/abs/1709.075041709.07504.
  
  Arkani-Hamed:2017mur
  N. Arkani-Hamed, Y. Bai, S. He and G. Yan, Scattering forms and the
    positive geometry of kinematics, color and the worldsheet,
    https://doi.org/10.1007/JHEP05(2018)096JHEP 05
    (2018) 096 [https://arxiv.org/abs/1711.091021711.09102].
  
  Arkani-Hamed:2017tmz
  N. Arkani-Hamed, Y. Bai and T. Lam, Positive Geometries and Canonical
    Forms, https://doi.org/10.1007/JHEP11(2017)039JHEP
    11 (2017) 039
    [https://arxiv.org/abs/1703.045411703.04541].
  
  ArkaniHamed:2010kv
  N. Arkani-Hamed, J.L. Bourjaily, F. Cachazo, S. Caron-Huot and J. Trnka,
    The all-loop integrand for scattering amplitudes in planar N=4
    SYM, https://doi.org/10.1007/JHEP01(2011)041JHEP 01 (2011) 041 [https://arxiv.org/abs/1008.29581008.2958].
  
  Arkani-Hamed:2019mrd
  N. Arkani-Hamed, S. He and T. Lam, Stringy canonical forms,
    https://doi.org/10.1007/JHEP02(2021)069JHEP 02
    (2021) 069 [https://arxiv.org/abs/1912.087071912.08707].
  
  barany1987computing
  I. Bárány and Z. Füredi, Computing the volume is difficult,
    https://doi.org/10.1007/BF02187886Discrete Comput. Geom.
    2 (1987) 319.
  
  bender1969anharmonic
  C.M. Bender and T.T. Wu, Anharmonic oscillator,
    https://doi.org/10.1103/PhysRev.184.1231Physical Review
    184 (1969) 1231.
  
  Berghoff:2017dyq
  M. Berghoff, Feynman amplitudes on moduli spaces of graphs,
    https://doi.org/10.4171/AIHPD/84Annales de l'I.H.P. D
    7 (2020) 203 [https://arxiv.org/abs/1709.005451709.00545].
  
  Berghoff:2020bug
  M. Berghoff and D. Kreimer, Graph complexes and Feynman rules,
    https://arxiv.org/abs/2008.095402008.09540.
  
  berkesch2014euler
  C. Berkesch, J. Forsgård and M. Passare, Euler-Mellin integrals and
    A-hypergeometric functions,
    https://doi.org/10.1307/mmj/1395234361Michigan Math. J.
    63 (2014) 101 [https://arxiv.org/abs/1103.62731103.6273].
  
  Bern:1994zx
  Z. Bern, L.J. Dixon, D.C. Dunbar and D.A. Kosower, One loop n point
    gauge theory amplitudes, unitarity and collinear limits,
    https://doi.org/10.1016/0550-3213(94)90179-1Nucl. Phys. B
    425 (1994) 217
    [https://arxiv.org/abs/hep-ph/9403226hep-ph/9403226].
  
  Binoth:2000ps
  T. Binoth and G. Heinrich, An automatized algorithm to compute infrared
    divergent multiloop integrals,
    https://doi.org/10.1016/S0550-3213(00)00429-6Nucl. Phys. B
    585 (2000) 741
    [https://arxiv.org/abs/hep-ph/0004013hep-ph/0004013].
  
  Bloch:2005bh
  S. Bloch, H. Esnault and D. Kreimer, On motives associated to graph
    polynomials, https://doi.org/10.1007/s00220-006-0040-2Comm.
    Math. Phys. 267 (2006) 181
    [https://arxiv.org/abs/math/0510011math/0510011].
  
  bloch2015cutkosky
  S. Bloch and D. Kreimer, Cutkosky rules and Outer space,
    https://arxiv.org/abs/1512.017051512.01705.
  
  Bogner:2007cr
  C. Bogner and S. Weinzierl, Resolution of singularities for multi-loop
    integrals, https://doi.org/10.1016/j.cpc.2007.11.012Comput.
    Phys. Commun. 178 (2008) 596
    [https://arxiv.org/abs/0709.40920709.4092].
  
  Bogner:2010kv
  C. Bogner and S. Weinzierl, Feynman graph polynomials,
    https://doi.org/10.1142/S0217751X10049438Int. J. Mod. Phys. A
    25 (2010) 2585
    [https://arxiv.org/abs/1002.34581002.3458].
  
  Borinsky:2014xwa
  M. Borinsky, Feynman graph generation and calculations in the Hopf
    algebra of Feynman graphs,
    https://doi.org/10.1016/j.cpc.2014.07.023Comput. Phys. Commun.
    185 (2014) 3317
    [https://arxiv.org/abs/1402.26131402.2613].
  
  Borinsky:2017hkb
  M. Borinsky, Renormalized asymptotic enumeration of Feynman diagrams,
    https://doi.org/10.1016/j.aop.2017.07.009Annals Phys.
    385 (2017) 95
    [https://arxiv.org/abs/1703.008401703.00840].
  
  borinsky2018graphs
  M. Borinsky, Graphs in perturbation theory: Algebraic structure and
    asymptotics, Springer (2018),
    https://doi.org/10.1007/978-3-030-03541-910.1007/978-3-030-03541-9,
    [https://arxiv.org/abs/1807.020461807.02046].
  
  Borinsky:2019rtu
  M. Borinsky and K. Vogtmann, The Euler characteristic of
    Out(F_n),
    https://doi.org/10.4171/CMH/501Comment. Math. Helv. 95 (2020) 703 [https://arxiv.org/abs/1907.035431907.03543].
  
  Borowka:2018goh
  S. Borowka, G. Heinrich, S. Jahn, S. Jones, M. Kerner and J. Schlenk, A
    GPU compatible quasi-Monte Carlo integrator interfaced to pySecDec,
    https://doi.org/10.1016/j.cpc.2019.02.015Comput. Phys. Commun.
    240 (2019) 120
    [https://arxiv.org/abs/1811.117201811.11720].
  
  Borowka:2017idc
  S. Borowka, G. Heinrich, S. Jahn, S. Jones, M. Kerner, J. Schlenk et al.,
    pySecDec: a toolbox for the numerical evaluation of multi-scale
    integrals, https://doi.org/10.1016/j.cpc.2017.09.015Comput.
    Phys. Commun. 222 (2018) 313
    [https://arxiv.org/abs/1703.096921703.09692].
  
  Borowka:2015mxa
  S. Borowka, G. Heinrich, S. Jones, M. Kerner, J. Schlenk and T. Zirke,
    SecDec-3.0: numerical evaluation of multi-scale integrals beyond one
    loop, https://doi.org/10.1016/j.cpc.2015.05.022Comput. Phys.
    Commun. 196 (2015) 470
    [https://arxiv.org/abs/1502.065951502.06595].
  
  Borowka:2018dsa
  S. Borowka, T. Gehrmann and D. Hulme, Systematic approximation of
    multi-scale Feynman integrals,
    https://doi.org/10.1007/JHEP08(2018)111JHEP 08
    (2018) 111 [https://arxiv.org/abs/1804.068241804.06824].
  
  Broadhurst:1995km
  D.J. Broadhurst and D. Kreimer, Knots and numbers in φ^4 theory
    to 7 loops and beyond,
    https://doi.org/10.1142/S012918319500037XInt. J. Mod. Phys. C
    6 (1995) 519
    [https://arxiv.org/abs/hep-ph/9504352hep-ph/9504352].
  
  Broedel:2018iwv
  J. Broedel, C. Duhr, F. Dulat, B. Penante and L. Tancredi, Elliptic
    symbol calculus: from elliptic polylogarithms to iterated integrals of
    Eisenstein series,
    https://doi.org/10.1007/JHEP08(2018)014JHEP 08
    (2018) 014 [https://arxiv.org/abs/1803.102561803.10256].
  
  Brown:2008um
  F. Brown, The massless higher-loop two-point function,
    https://doi.org/10.1007/s00220-009-0740-5Comm. Math. Phys.
    287 (2009) 925
    [https://arxiv.org/abs/0804.16600804.1660].
  
  Brown:2009ta
  F. Brown, On the periods of some Feynman integrals,
    https://arxiv.org/abs/0910.01140910.0114.
  
  Brown:2015fyf
  F. Brown, Feynman amplitudes, coaction principle, and cosmic Galois
    group, https://doi.org/10.4310/CNTP.2017.v11.n3.a1Commun.
    Num. Theor. Phys. 11 (2017) 453
    [https://arxiv.org/abs/1512.064091512.06409].
  
  Brown:2011pj
  F. Brown and D. Kreimer, Angles, scales and parametric renormalization,
    https://doi.org/10.1007/s11005-013-0625-6Lett. Math. Phys.
    103 (2013) 933
    [https://arxiv.org/abs/1112.11801112.1180].
  
  Brown:2010bw
  F. Brown and O. Schnetz, A K3 in ϕ^4,
    https://doi.org/10.1215/00127094-1644201Duke Math. J.
    161 (2012) 1817
    [https://arxiv.org/abs/1006.40641006.4064].
  
  Brown:2015ztw
  F. Brown and O. Schnetz, Single-valued multiple polylogarithms and a
    proof of the zig-zag conjecture,
    https://doi.org/10.1016/j.jnt.2014.09.007J. Number Theor.
    148 (2015) 478.
  
  Brown:2009rc
  F. Brown and K. Yeats, Spanning forest polynomials and the
    transcendental weight of Feynman graphs,
    https://doi.org/10.1007/s00220-010-1145-1Comm. Math. Phys.
    301 (2011) 357
    [https://arxiv.org/abs/0910.54290910.5429].
  
  bruns2010normaliz
  W. Bruns and B. Ichim, Normaliz: algorithms for affine monoids and
    rational cones,
    https://doi.org/10.1016/j.jalgebra.2010.01.031J. Algebra
    324 (2010) 1098
    [https://arxiv.org/abs/0910.28450910.2845].
  
  Cachazo:2019ngv
  F. Cachazo, N. Early, A. Guevara and S. Mizera, Scattering equations:
    From projective spaces to tropical Grassmannians,
    https://doi.org/10.1007/JHEP06(2019)039JHEP 06
    (2019) 039 [https://arxiv.org/abs/1903.089041903.08904].
  
  Capatti:2019ypt
  Z. Capatti, V. Hirschi, D. Kermanschah and B. Ruijl, Loop-tree duality
    for multiloop numerical integration,
    https://doi.org/10.1103/PhysRevLett.123.151602Phys. Rev. Lett.
    123 (2019) 151602
    [https://arxiv.org/abs/1906.061381906.06138].
  
  Catani:2008xa
  S. Catani, T. Gleisberg, F. Krauss, G. Rodrigo and J.-C. Winter, From
    loops to trees by-passing Feynman's theorem,
    https://doi.org/10.1088/1126-6708/2008/09/065JHEP 09 (2008) 065 [https://arxiv.org/abs/0804.31700804.3170].
  
  chan2013tropical
  M. Chan, M. Melo and F. Viviani, Tropical Teichmüller and Siegel
    spaces, https://doi.org/10.1090/conm/589Contemp. Math.
    589 (2013) 45 [https://arxiv.org/abs/1207.24431207.2443].
  
  Collins:1984xc
  J.C. Collins, Renormalization, Cambridge University Press (1986).
  
  Connes:1999yr
  A. Connes and D. Kreimer, Renormalization in quantum field theory and
    the Riemann-Hilbert problem. 1. The Hopf algebra structure of graphs and the
    main theorem, https://doi.org/10.1007/s002200050779Comm.
    Math. Phys. 210 (2000) 249
    [https://arxiv.org/abs/hep-th/9912092hep-th/9912092].
  
  Courtiel:2019dnq
  J. Courtiel and K. Yeats, Next-to^k leading log expansions by chord
    diagrams, https://doi.org/10.1007/s00220-020-03691-7Comm.
    Math. Phys. 377 (2020) 469
    [https://arxiv.org/abs/1906.051391906.05139].
  
  Culler1986ModuliOG
  M. Culler and K. Vogtmann, Moduli of graphs and automorphisms of free
    groups, https://doi.org/10.1007/BF01388734Invent. Math.
    84 (1986) 91.
  
  deCalan:1981szv
  C. de Calan and V. Rivasseau, Local existence of the Borel transform in
    Euclidean ϕ^4 in four-dimensions,
    https://doi.org/10.1007/BF01206946Comm. Math. Phys. 82 (1981) 69.
  
  delaCruz:2019skx
  L. de la Cruz, Feynman integrals as A-hypergeometric functions,
    https://doi.org/10.1007/JHEP12(2019)123JHEP 12
    (2019) 123 [https://arxiv.org/abs/1907.005071907.00507].
  
  AIHPA_1999__71_1_1_0
  E. Delabaere and F. Pham, Resurgent methods in semi-classical
    asymptotics, Annales de l'I.H.P. Physique théorique 71
    (1999) 1.
  
  doker2011geometry
  J.S. Doker, Geometry of generalized permutohedra, Ph.D. thesis, UC
    Berkeley, 2011.
  
  Drummond:2010cz
  J.M. Drummond, J.M. Henn and J. Trnka, New differential equations for
    on-shell loop integrals,
    https://doi.org/10.1007/JHEP04(2011)083JHEP 04
    (2011) 083 [https://arxiv.org/abs/1010.36791010.3679].
  
  dyer1991random
  M. Dyer, A. Frieze and R. Kannan, A random polynomial-time algorithm for
    approximating the volume of convex bodies,
    https://doi.org/10.1145/102782.102783JACM 38
    (1991) 1.
  
  Eden:1966dnq
  R.J. Eden, P.V. Landshoff, D.I. Olive and J.C. Polkinghorne, The
    analytic S-matrix, Cambridge University Press (1966).
  
  Feng:2019bdx
  T.-F. Feng, C.-H. Chang, J.-B. Chen and H.-B. Zhang, GKZ-hypergeometric
    systems for Feynman integrals,
    https://doi.org/10.1016/j.nuclphysb.2020.114952Nucl. Phys. B
    953 (2020) 114952
    [https://arxiv.org/abs/1912.017261912.01726].
  
  fujishige1983note
  S. Fujishige and N. Tomizawa, A note on submodular functions on
    distributive lattices, J. Oper. Res. Soc. Japan 26
    (1983) 309.
  
  Gehrmann-DeRidder:2003pne
  A. Gehrmann-De Ridder, T. Gehrmann and G. Heinrich, Four particle phase
    space integrals in massless QCD,
    https://doi.org/10.1016/j.nuclphysb.2004.01.023Nucl. Phys. B
    682 (2004) 265
    [https://arxiv.org/abs/hep-ph/0311276hep-ph/0311276].
  
  gkz1994
  I. Gelfand, M. Kapranov and A. Zelevinsky, Discriminants, Resultants and
    Multidimensional Determinants, Birkhäuser Boston Inc. (1994).
  
  gelfand1990generalized
  I. Gelfand, M. Kapranov and A. Zelevinsky, Generalized Euler integrals
    and A-hypergeometric functions,
    https://doi.org/10.1016/0001-8708(90)90048-RAdv. Math.
    84 (1990) 255.
  
  Green:1987sp
  M.B. Green, J. Schwarz and E. Witten, Superstring Theory. Vol. 1:
    Introduction, Cambridge University Press (1988).
  
  Hahn:2004fe
  T. Hahn, CUBA: A library for multidimensional numerical integration,
    https://doi.org/10.1016/j.cpc.2005.01.010Comput. Phys. Commun.
    168 (2005) 78
    [https://arxiv.org/abs/hep-ph/0404043hep-ph/0404043].
  
  hammersley1964monte
  J.M. Hammersley and D.C. Handscomb, Monte Carlo methods, Springer
    (1964).
  
  He:2020onr
  S. He, Z. Li, P. Raman and C. Zhang, Stringy canonical forms and binary
    geometries from associahedra, cyclohedra and generalized permutohedra,
    https://doi.org/10.1007/JHEP10(2020)054JHEP 10
    (2020) 054 [https://arxiv.org/abs/2005.073952005.07395].
  
  Heinrich:2008si
  G. Heinrich, Sector decomposition,
    https://doi.org/10.1142/S0217751X08040263Int. J. Mod. Phys. A
    23 (2008) 1457
    [https://arxiv.org/abs/0803.41770803.4177].
  
  Heinrich:2020ybq
  G. Heinrich, Collider Physics at the Precision Frontier,
    https://arxiv.org/abs/2009.005162009.00516.
  
  henk200416
  M. Henk, J. Richter-Gebert and G.M. Ziegler, Basic properties of convex
    polytopes, Handbook of discrete and computational geometry (2018)
    383.
  
  Henn:2013pwa
  J.M. Henn, Multiloop integrals in dimensional regularization made
    simple, https://doi.org/10.1103/PhysRevLett.110.251601Phys.
    Rev. Lett. 110 (2013) 251601
    [https://arxiv.org/abs/1304.18061304.1806].
  
  henn2015lectures
  J.M. Henn, Lectures on differential equations for Feynman integrals,
    https://doi.org/10.1088/1751-8113/48/15/153001J. Phys. A
    48 (2015) 153001
    [https://arxiv.org/abs/1412.22961412.2296].
  
  Hepp:1966eg
  K. Hepp, Proof of the Bogolyubov-Parasiuk theorem on renormalization,
    https://doi.org/10.1007/BF01773358Comm. Math. Phys. 2 (1966) 301.
  
  Herzog:2018ily
  F. Herzog, Geometric IR subtraction for final state real radiation,
    https://doi.org/10.1007/JHEP08(2018)006JHEP 08
    (2018) 006 [https://arxiv.org/abs/1804.079491804.07949].
  
  Hu:2018liw
  S. Hu, O. Schnetz, J. Shaw and K. Yeats, Further investigations into the
    graph theory of ϕ^4-periods and the c_2 invariant,
    https://arxiv.org/abs/1812.087511812.08751.
  
  johansson2013argument
  P. Johansson, The argument cycle and the coamoeba,
    https://doi.org/10.1080/17476933.2011.592581Complex Var.
    Elliptic Equ. 58 (2013) 373.
  
  Kaneko:2009qx
  T. Kaneko and T. Ueda, A geometric method of sector decomposition,
    https://doi.org/10.1016/j.cpc.2010.04.001Comput. Phys. Commun.
    181 (2010) 1352
    [https://arxiv.org/abs/0908.28970908.2897].
  
  Klausen:2019hrg
  R.P. Klausen, Hypergeometric series representations of Feynman
    integrals by GKZ hypergeometric systems,
    https://doi.org/10.1007/JHEP04(2020)121JHEP 04
    (2020) 121 [https://arxiv.org/abs/1910.086511910.08651].
  
  knuth2014art
  D.E. Knuth, The Art of Computer Programming, Volume 2: Seminumerical
    Algorithms, Addison-Wesley (2014).
  
  Kompaniets:2017yct
  M.V. Kompaniets and E. Panzer, Minimally subtracted six loop
    renormalization of O(n)-symmetric ϕ^4 theory and critical exponents,
    https://doi.org/10.1103/PhysRevD.96.036016Phys. Rev. D
    96 (2017) 036016
    [https://arxiv.org/abs/1705.064831705.06483].
  
  Kotikov:1990kg
  A. Kotikov, Differential equations method: New technique for massive
    Feynman diagrams calculation,
    https://doi.org/10.1016/0370-2693(91)90413-KPhys. Lett. B
    254 (1991) 158.
  
  Kreimer:2018kah
  D. Kreimer, Multi-valued Feynman graphs and scattering theory,  in
    KMPB Conference: Elliptic Integrals, Elliptic Functions and Modular
    Forms in Quantum Field Theory, pp. 295–325, 2019,
    https://doi.org/10.1007/978-3-030-04480-0_13DOI
    [https://arxiv.org/abs/1807.002881807.00288].
  
  Kreimer:2006ua
  D. Kreimer and K. Yeats, An etude in non-linear Dyson-Schwinger
    equations,
    https://doi.org/10.1016/j.nuclphysbps.2006.09.036Nucl. Phys. B
    Proc. Suppl. 160 (2006) 116
    [https://arxiv.org/abs/hep-th/0605096hep-th/0605096].
  
  Kreimer:2006gm
  D. Kreimer and K. Yeats, Recursion and growth estimates in
    renormalizable quantum field theory,
    https://doi.org/10.1007/s00220-008-0431-7Comm. Math. Phys.
    279 (2008) 401
    [https://arxiv.org/abs/hep-th/0612179hep-th/0612179].
  
  kruskal1956shortest
  J.B. Kruskal, On the shortest spanning subtree of a graph and the
    traveling salesman problem,
    https://doi.org/10.2307/2033241Proc. Am. Math. Soc. 7 (1956) 48.
  
  le2012large
  J.-C. Le Guillou and J. Zinn-Justin, Large-order behaviour of
    perturbation theory, Elsevier (2012).
  
  Lepage:1977sw
  G. Lepage, A new algorithm for adaptive multidimensional integration,
    https://doi.org/10.1016/0021-9991(78)90004-9J. Comput. Phys.
    27 (1978) 192.
  
  lovasz1993random
  L. Lovász and M. Simonovits, Random walks in a convex body and an
    improved volume algorithm,
    https://doi.org/10.1002/rsa.3240040402Random structures &
    algorithms 4 (1993) 359.
  
  maclagan2015introduction
  D. Maclagan and B. Sturmfels, Introduction to tropical geometry, AMS
    (2015).
  
  von2015quasi
  A. Manteuffel, Andreasvon Manteuffel, E. Panzer and
    R.M. Schabinger, A quasi-finite basis for multi-loop Feynman
    integrals, https://doi.org/10.1007/JHEP02(2015)120JHEP
    02 (2015) 120 [https://arxiv.org/abs/1411.73921411.7392].
  
  McKane:2018ocs
  A.J. McKane, Perturbation expansions at large order: results for scalar
    field theories revisited,
    https://doi.org/10.1088/1751-8121/aaf768J. Phys. A 52 (2019) 055401 [https://arxiv.org/abs/1807.006561807.00656].
  
  mckay1983spanning
  B.D. McKay, Spanning trees in regular graphs,
    https://doi.org/10.1016/S0195-6698(83)80045-6Eur. J. Comb.
    4 (1983) 149.
  
  Nakanishi:110324
  N. Nakanishi, Graph theory and Feynman integrals, Gordon and Breach
    (1971).
  
  niederreiter1992random
  H. Niederreiter, Random number generation and quasi-Monte Carlo methods,
    SIAM (1992),
    https://doi.org/10.1137/1.978161197008110.1137/1.9781611970081.
  
  nilsson2013mellin
  L. Nilsson and M. Passare, Mellin transforms of multivariate rational
    functions, https://doi.org/10.1007/s12220-011-9235-7J. Geom.
    Anal. 23 (2013) 24
    [https://arxiv.org/abs/1010.50601010.5060].
  
  nisse2013
  M. Nisse and F. Sottile, The phase limit set of a variety,
    https://doi.org/10.2140/ant.2013.7.339Algebra Number Theory
    7 (2013) 339.
  
  panzer2020personal
  E. Panzer. personal communication, 2020-09-15.
  
  Panzer:2014caa
  E. Panzer, Algorithms for the symbolic integration of hyperlogarithms
    with applications to Feynman integrals,
    https://doi.org/10.1016/j.cpc.2014.10.019Comput. Phys. Commun.
    188 (2015) 148
    [https://arxiv.org/abs/1403.33851403.3385].
  
  Panzer:2019yxl
  E. Panzer, Hepp's bound for Feynman graphs and matroids,
    https://arxiv.org/abs/1908.098201908.09820.
  
  Panzer:2016snt
  E. Panzer and O. Schnetz, The Galois coaction on ϕ^4 periods,
    https://doi.org/10.4310/CNTP.2017.v11.n3.a3Commun. Num. Theor.
    Phys. 11 (2017) 657
    [https://arxiv.org/abs/1603.042891603.04289].
  
  postnikov2009permutohedra
  A. Postnikov, Permutohedra, associahedra, and beyond,
    https://doi.org/10.1093/imrn/rnn153Int. Math. Res. Not.
    2009 (2009) 1026
    [https://arxiv.org/abs/math/0507163math/0507163].
  
  Remiddi:1997ny
  E. Remiddi, Differential equations for Feynman graph amplitudes,
    Nuovo Cim. A 110 (1997) 1435
    [https://arxiv.org/abs/hep-th/9711188hep-th/9711188].
  
  Runkel:2019yrs
  R. Runkel, Z. Szőr, J.P. Vesga and S. Weinzierl, Causality and
    loop-tree duality at higher loops,
    https://doi.org/10.1103/PhysRevLett.122.111603Phys. Rev. Lett.
    122 (2019) 111603
    [https://arxiv.org/abs/1902.021351902.02135].
  
  Schnetz:2008mp
  O. Schnetz, Quantum periods: A census of ϕ^4-transcendentals,
    https://doi.org/10.4310/CNTP.2010.v4.n1.a1Commun. Num. Theor.
    Phys. 4 (2010) 1
    [https://arxiv.org/abs/0801.28560801.2856].
  
  schnetz2018numbers
  O. Schnetz, Numbers and functions in quantum field theory,
    https://doi.org/10.1103/PhysRevD.97.085018Phys. Rev. D
    97 (2018) 085018
    [https://arxiv.org/abs/1606.085981606.08598].
  
  schultka2018toric
  K. Schultka, Toric geometry and regularization of Feynman integrals,
    https://arxiv.org/abs/1806.010861806.01086.
  
  Smirnov:2013eza
  A.V. Smirnov, FIESTA 3: cluster-parallelizable multiloop numerical
    calculations in physical regions,
    https://doi.org/10.1016/j.cpc.2014.03.015Comput. Phys. Commun.
    185 (2014) 2090
    [https://arxiv.org/abs/1312.31861312.3186].
  
  Smirnov:2008aw
  A. Smirnov and V. Smirnov, Hepp and Speer sectors within modern
    strategies of sector decomposition,
    https://doi.org/10.1088/1126-6708/2009/05/004JHEP 05 (2009) 004 [https://arxiv.org/abs/0812.47000812.4700].
  
  Smirnov:2012gma
  V.A. Smirnov, Analytic tools for Feynman integrals, Springer (2012).
  
  Speer:1975dc
  E. Speer, Ultraviolet and Infrared Singularity Structure of Generic
    Feynman Amplitudes, Annales de l'I.H.P. Physique théorique
    23 (1975) 1.
  
  spielman2014nearly
  D.A. Spielman and S.-H. Teng, Nearly linear time algorithms for
    preconditioning and solving symmetric, diagonally dominant linear systems,
    https://doi.org/10.1137/090771430SIAM J. Matrix Anal. Appl.
    35 (2014) 835
    [https://arxiv.org/abs/cs/0607105cs/0607105].
  
  whittaker1963course
  E. Whittaker and G. Watson, A Course of Modern Analysis, Cambridge
    University Press (1963).
  
  