Unveiling Regions in multi-scale Feynman Integrals using Singularities and Power Geometry

We introduce a novel approach for solving the problem of identifying regions in the framework of Method of Regions by considering singularities and the associated Landau equations given a multi-scale Feynman diagram. These equations are then analyzed by an expansion in a small threshold parameter via the Power Geometry technique. This effectively leads to the analysis of Newton Polytopes which are evaluated using a Mathematica based convex hull program. Furthermore, the elements of the Gr\"{o}bner Basis of the Landau Equations give a family of transformations, which when applied, reveal regions like potential and Glauber. Several one-loop and two-loop examples are studied and benchmarked using our algorithm which we call ASPIRE.


Introduction
Over the last couple of decades, Effective Field Theories (EFTs) have become mainstream in problems that involve separation of scales. While there are several applications in low-energy QCD, its foray into high-energy processes is still relatively new and recent years have seen the application of EFT ideas to problems in particle physics that involve several scales. Collider physics provides a perfect example of a multi-scale problem involving particles with high energies as well as comparatively low-mass particles such as protons. Multi-scale Feynman diagrams arise naturally in many branches of elementary particle physics. The analysis of such diagrams has led to new EFTs such as the soft-collinear effective theory, heavy quark effective theory and so on. For an accessible introduction to the subject of effective field theories, see, e.g., Refs. [1][2][3]. In a multi-loop problem, one Electronic supplementary material The online version of this article (https://doi.org/10.1140/epjc/s10052-019-6533-x) contains supplementary material, which is available to authorized users. a e-mail: ratansarkar@iisc.ac.in typically encounters several mass and kinematic scales. A strategy that has been very effective in these problems is the Method of Regions (MoR). This allows one to carry out asymptotic expansions of Feynman integrals within dimensional regularization. The Feynman integral for a process is expanded as a sum of simpler integrals which can then be done term by term. Further, the diagrams obtained in MoR can also be obtained from an appropriate EFT. The different regions identified correspond to different EFTs characterized by a threshold parameter.
An accessible example of separation of scales can be found in three flavor chiral perturbation theory, where one finds instances of diagrams with a hierarchy of masses, namely the masses of the pion, the kaon and the eta. Recent analytical progress can be found in Refs. [4,5]. The MoR was applied by Kaiser and Schweizer [6] for studying π − k scattering processes in the context of ChPT. In this work, we also study this process within our framework to identify the associated regions.
The application of the MoR to multi-scale problems has been studied now for nearly two decades, starting with the fundamental work of Beneke and Smirnov [7]. Subsequently, this approach has been used for Drell-Yan Processes [8], studying massless fermionic processes at next to leading order [9], investigating processes in the Sudakov limit [10] and so on. To find the leading order contribution of Feynman integrals at threshold, one needs to sum contributions from several regions which span the entire loop momentum space. The contributing regions result from the presence of a hierarchy of masses, or from components of some momenta becoming small or large compared to others.
In a recent work, Pak and Smirnov [11] proposed a geometric algorithm which could be automated to find the regions of a given Feynman diagram. By going to the Alpharepresentation, they show that the problem can be turned into 0123456789().: V,-vol a geometric one and can be solved by finding the convex hull of a set of carefully chosen points. The Mathematica package developed that automates the finding of regions, called asy.m, which we will refer to as ASY, uses a C++-based QuickHull algorithm [12].
In its original implementation, ASY fails to identify potential and Glauber regions. These regions usually manifest as differences between the Alpha-parameters. In the updated asy2.1 [13,14] a new feature called PreResolve was introduced that eliminates the differences of Alpha-parameters from the Symanzik polynomials using linear transformations. We will refer to the upgraded version of the code as ASY2.
In the present work, we approach the problem from another perspective by looking at the singular structure of the Feynman integral in the Alpha-representation. The analysis of the Pinched Singular Surfaces, in the momentum space, in connection to the regions [15][16][17][18] is well understood in terms of the Landau Equations. We set up and study the set of Landau Equations, in the alpha parameter space, for a given process using the Gröbner Basis and derive a criterion for the determination of the transformations required to reveal the regions within the framework of Power Geometry [19][20][21], thus, demonstrating a new way of solving the problem of finding regions.
The paper is organized as follows. In Sect. 2, we introduce elements essential to our algorithm. We review the strategy of the MoR in Sect. 2.1 and then give an overview of the geometric approach discussed by Pak and Smirnov in Sect. 2.2, by considering an example from [11]. In Sect. 2.3, we revisit the problem of finding the regions through an alternative approach that links the analytic structure of the Feynman integral to the regions. In fact, the contents of Sect. 2.3 constitute the important theoretical progress being reported in this work. We conclude Sect. 2 by summarizing the steps of our algorithm, "Algebro-geometric analysis of Singular Polynomials for Identification of REgions (ASPIRE)". We apply the algorithm to one-loop and two-loop examples in Sect. 3. Finally, in Sect. 4, we discuss our results and present our future goals. The details of our Mathematica notebooks and external packages used in this work are given in the appendix.

Formalism
In this section, we set up the formalism for identifying different regions using the singular structure of the Feynman integral. This process can be automated using ideas from power geometry. However, for the sake of completeness, we also summarize the technique of Pak and Smirnov in the subsequent sub-section. The technique of MoR was proposed in an attempt to analytically approximate various processes within perturbation theory [7,[22][23][24]. The idea of MoR is to provide an expansion of the integrand in ratio of the scales involved, usually in the form of low-energy scale to high-energy scale. This results in expressing the original Feynman integral as a sum over simpler integrals, all of which need to be integrated over their corresponding domains, which are called regions.
The original integral is the sum of the contributions of the above three regions each of which is evaluated within dimensional regularization. While it is easy to identify the regions in this particular example by looking at the different scales in the Feynman integral, the procedure is a non-trivial task in general, especially when one wants to evaluate multi-loop processes. Even at the one-loop level, one encounters nontrivial regions, that involve a multitude of scales. Another obvious difficulty in identifying regions is when the components of the momenta scale differently or when scalings of the difference of momenta are involved. In the following, we review a specific implementation by Pak and Smirnov that allows to isolate the regions.

ASY implementation
Pak and Smirnov proposed an algorithm to automate the process of finding regions, which was documented in [11,13] together with their codes "asy.m" and "asy2.m" (referred to here as ASY and ASY2). The second version adds crucial improvement to the first, as we summarize in the ensuing discussions. The basic idea is to parameterize the Feynman integral using the Alpha-parameters and then carry out the integration over the loop momenta using dimensional regularization to obtain its Alpha-representation. Expanding a process in the momentum space in scalings of momenta (or its components) is equivalent to expanding the Alpha-representation in scalings of combinations of Alpha-parameters.
ASY starts with expressing the integral in the Alpharepresentation. For the process in Fig. 1 this yields, where U and F are the Symanzik polynomials given by: and which are homogeneous in the Alpha-parameters, x 1 and x 2 . Furthermore, Pak and Smirnov, build a new polynomial U ·F that allows for a combined analysis of both the polynomials. All the terms in the leading order Symanzik polynomials have the same scaling in terms of the expansion (or threshold) parameter, ρ ≡ |m 2 / p 2 |. More precisely, if one constructs a set of vector exponents {r i } for each monomial and scalings {v i } such that x i scales as ρ v i then x r i i will scale as ρ r i v i . Hence, the monomials in the Symanzik polynomials will scale as where r 0 is the exponent of the threshold parameter appearing in the prefactor of the monomials and v 0 = 1. We can now construct an n + 1 dimensional vector with components r = (r 0 , r 1 , . . . r n ) such that the vector exponents obtained from each monomial of the leading order expansion, v leading , lie on a plane described by the equation where c is a constant. All the other terms which do not appear at leading order, i.e v subleading , will lie above the surface described by Eq. (6). These points satisfy the condition In general, it can be seen that if we construct the collection of vector exponents with n + 1 components and plot them in an n-dimensional sub-space then the leading order terms corresponding to a region will be points lying on the same surface and all the other points will lie above it. This immediately leads to the interpretation of the surface as a bottom facet of the convex hull of the set of vector exponents. Thus, finding the regions amounts to finding the convex hull of the set of vector exponents and then finding the normals of the lower facets of the convex hull.
For the Symanzik polynomials in Eqs. (3) and (4), ASY first calculates the product from which one can extract the set of vector exponents(using threshold parameter ρ ≡ |m 2 /q 2 |) Each vector exponent, corresponds to a monomial in Eq. (8) and the components give the exponents of the expansion parameter ρ, followed by the set { x i }. ASY projects this set of vector exponents onto a lower dimensional subspace. The new set of exponents, after performing the projection, is now r = {(1, 3), (1, 2), (0, 2), (1, 1), (0, 1), (1, 0)}. Following the discussion earlier in this section, ASY finds the convex hull of this projected set of exponents as shown in Fig. 2 and then finds the normals of the bottom facets of the convex hull.
For example, the leading order term for the hard region (defined by x 1 ∼ ρ 0 , x 2 ∼ ρ 0 ), is q 2 x 2 1 x 2 + q 2 x 1 x 2 2 corresponding to the projected points (0, 2), (0, 1) lie on the plane with the normal vector (1,0). This is also seen from the fact that the leading order terms are independent of ρ and thus scale as ρ 0 . All the other terms have scalings larger than ρ 0 and thus will lie above the plane.
The first version of the code ASY could identify regions, except the Glauber and the potential, in several instances. The second version of the code ASY2 fixes this shortcoming by linearly transforming the Alpha-parameters, eliminating any term in U · F that appears with a difference between the Alpha-parameters. This is done by an option PreResolve in the code. The main feature that distinguishes the two versions of the ASY codes are the implementation of linear transformations which allows the code to now identify the potential and Glauber regions. This approach has been very successful in determining regions and has been applied to several examples. To this extent, ASY also provides a very useful crosscheck on prior studies besides a useful benchmark for comparison. In the following sub-sections we will demonstrate a new way of solving the problem of finding regions based on the singularities of Feynman integral in Alpha representation.

Determination of the regions using the analytic structure of the propagator
The different regions, where a particular mass or kinematical scale becomes important can be linked to the underlying singularities of the Feynman integral. In the following, we will introduce the main concepts and motivate the ideas that will lead to the development of the final algorithm. We first give an overview of the singularities that are of interest for our problem, followed by a review of the basic understanding of particle thresholds as pinched singularities in momentum space. This interpretation is well understood and can be mathematically expressed using a set of equations called the Landau equations. Since expansions in the neighborhood of the singular surfaces give us the leading order behavior of Feynman amplitudes, we perform similar expansions in the Alpha-parameter space in carefully chosen neighborhoods of the singular points. This requires us to use techniques from the field of power geometry. We then motivate the use of Gröbner basis for the identification of all neighborhoods of the singular points.

Singularities and threshold processes
Understanding the analytic structure of the amplitude is crucial to identifying the different regions. The poles in the integrand of the amplitude for a given process are functions of kinematical invariants, loop momenta etc. Therefore, when these parameters vary, the poles in the integration domain move. In the case of isolated singularities, it is always possible to deform the contour of integration to avoid these singular points. However, sometimes, the poles migrate so as to pinch the contour of integration (pinch singularities) or move to one of the end point of the integration (end point singular- Fig. 3 Types of singularities: P1 is a simple pole, P2 and P3 are Pinched Singularities and P4 is an End Point Singularity. While the contour between the points A and B can be deformed so as to avoid the simple pole P1, the same is not true for the pinch and end-point singularities ities) as illustrated in Fig. 3. In such cases, these singularities cannot be avoided by contour deformations. The condition for a point to be one of these unavoidable singular points is the usual condition for establishing a singularity for a polynomial. For an arbitrary polynomial g({ x i }) that appears in the denominator of the Feynman integral, the point x i is singular point iff Therefore, at these unavoidable singular points, hereby referred to as just singular points, the integrand diverges. We will now adopt the approach of Coleman and Norton [25], also discussed in [26], to explain the connection between the singularities and physical events.
Consider a general Feynman amplitude in the Alpha representation and k i are the loop momenta, q j are internal momenta which are linear functions in loop and external momenta, m j are the masses and α j are the Alpha-parameters. For a singularity, corresponding to Eq. (9) we have the conditions: and Now, since each q i is a linear combination of the loop momenta, we have the condition with the constraint Equations (12) and (13) are the Landau equations [27]. If some of the α i = 0, then the corresponding internal lines get contracted to a vertex. Such a contraction leads to formation of effective vertices which can then be described via an EFT. Further, in the Alpha-parameter space, the end points are 0 and ∞. α i = 0 ∀ i corresponds to an end-point singularity. The relevant polynomial to analyze in this case is the D polynomial in Eq. (11) and its singular points. Given the Feynman graph of the process, one can define a separation between the vertices in terms of the momentum carried by the connecting lines as If α i = 0, Eq. (12) requires that q 2 i = m 2 i and we see that an on shell particle propagates from one vertex to the other, that is the reason for divergence of the integrand is if some of the internal lines become onshell. This process of setting some internal lines on-shell, referred to as performing unitary cuts, produces a set of sub-integrals called cut integrals. The resulting cut diagrams must now describe processes where on-shell particles propagate from one vertex to the other [28]. By finding an appropriate subset of these cut diagrams, it is possible to evaluate the original integral [29]. The parameter α i is identified with the proper time divided by the mass of the particle. The fact that α i > 0 then implies that the particle is propagating forward in time. An immediate corollary of the above for a closed loop is Therefore, it is evident that the Landau equations are statements for finding the singular points of a polynomial derived from the Feynman Integral. The analysis of Norton and Coleman followed by the work of Libby and Sterman [16,17] shows the connection between singularities and threshold processes which are described by EFTs.

Regions from Feynman graph in the alpha parametric form
A Feynman graph having l-loops, m-denominators, and rexternal momenta ( p 1 , . . . , p r ) in d-dimension has the form [30], where A, B are respectively l × l, l × r matrices and C j are constants.
In Alpha-parametric form, I (n) can be written as, where U, F are the Symanzik polynomials (of degree l and (l + 1) respectively) and |n| = n 1 + n 2 + · · · + n m .
In terms of the Symanzik polynomials, the Landau equations can be written as [15], Therefore, we seek approximate solutions, of the type α i ∼ cρ v i , of the Landau equations, near the singular surfaces. In our notation ρ is the threshold expansion parameter and c is a constant. The sets of {v i } corresponding to each of the solution (at leading order) branches, represent all the regions associated with the integral. We will extract these leading order solution in the neighborhood of the singular points using the techniques of power geometry, which we discuss next.

Newton polytope and power geometry
We are interested in obtaining the leading order scaling of the Alpha-parameters with respect to the expansion parameter ρ. This can be achieved using ideas from power geometry, developed by Bruno and Bathkin, which allows for obtaining solutions of a polynomial [19][20][21] in certain limits. Consider a generic polynomial in n-variables, where Q i are the exponents of the variables x i for each monomial, i.e., of the form x . . . and the Q i s are a set of natural numbers. Let χ = {X 0 } be a set of points such that g(X 0 ) = 0. If it turns out that ∂g(X 0 ) = 0, whereχ = {X 0 } andχ ⊆ χ , then the setχ contain the singular points of the polynomial. Power geometry allows one to obtain solutions to polynomials and is particularly useful around singular points. Before we outline the procedure for obtaining solutions to polynomials using the techniques developed by Bruno [19][20][21], we briefly summarize the basic definitions and concepts which we have used in our subsequent analysis.
Let us now define the following: (i) Support: The support S(g) is defined as the set of all vector exponents. For example, given a polynomial in two variables (x, y) (iv) Normal Cone: Let S ∈ R n be a compact convex set which is the support S(g), with faces {S } = {S 1 , S 2 . . . S r }. Let |ξ be a vector in R n . Then, lim sup{ ξ |η |η ∈ S} is attained by a vector |v belonging to a face S i for some i. It is easy to see that if |v is another vector in {S }, then ξ |v = ξ |v . The set of all |ξ such that ξ |v ≥ ξ |η for |η ∈ S is defined as the normal cone U d i , where i denotes the ith face and d its dimension.
(v) Cone of the problem: The Cone of the problem is a convex cone of vectors K = (s 1 , . . . , s n ) such that curves of the form where t parametrizes the polynomial, fill those regions of the X -variables space that we are interested in. (vi) Truncated polynomial: The truncation of the sum on the boundary subset is defined as Such a truncated polynomial should be quasi-homogenous, that is for the polynomialĝ d j (X ), there exists n integers {w 1 , . . . , w n }, called weights of the variables, such that the sum w = w 1 Q 1 + · · · + w n Q n is the same for all nonzero monomials ofĝ d j (X ), where Q = {Q 1 , . . . , Q n } is the vector exponent of terms in the polynomialĝ d j (X ). w is then the degree of the polynomial.
The truncated polynomials corresponding to the faces of the convex hull in Fig. 4, as well as its weights are listed in Table 1.
Finally, the algorithm to obtain the leading order solution in terms of a parameter t, given a polynomial g(X ) can be summarized as follows: 1. The set of singular points (or regular zeros) once identified, can lie on a generalized surface. It turns out that it is very convenient to change variables such that the singular points now lie at the origin, coordinate axis or on a coordinate plane. Such a choice allows one to define the cone of the problem, K , easily. As a result of these transformations, the polynomial g(X ) → g (X ), where X = T X under the map T . The problem now reduces to analyzing the g (X ) polynomial. 2. The support of g (X ), denoted as S(g (X )) is determined and convex hull of the support defines the Newton polyhedron. 3. For each of the two dimensional faces, 2 j of the Newton polytope, the normal vectors are determined, leading to the construction of the normal cone.
4. Only those normal vectors that intersect with the cone of the problem are retained. 5. For each normal cone lying in the cone of the problem, the truncated polynomial is determined, where the truncated polynomial is defined on the faces of the convex hull. 6. Finally, the vector P which gives us the scaling of the variables with respect to the parameter t is obtained, using the following theorem [20]: where a, b, c and p i are constants, belongs to the set g = {X : g(X ) = 0} and the vector P = ( p 1 , p 2 , p 3 ) belongs to U d j , then the first approximation x = at p 1 , y = bt p 2 , z = ct p 3 of the curve satisfies the truncated equationĝ d j (X ) = 0. As a result of the theorem, we get the leading order behavior of the solution. While the method is particularly useful for solutions around singular points where the implicit function theorem fails, this method can also be applied to obtain solutions about any zero of the polynomial, including regular zeros and has the advantage of being easily automated on a computer. Therefore, we will not make the distinction between the singular and the regular zeros of the polynomial. An important step in this algorithm is to determine the set of transformations that map the zeros of the polynomial to either the origin or the coordinate axes or the coordinate plane. The Gröbner basis [31] of the Landau equations conveniently gives us the required set of transformations as well as the appropriate neighborhoods of singular points where one needs to perform a leading order expansion of the polynomial.
With all these key definitions and ideas in place, we go on to enumerate the steps in our algorithm, ASPIRE, that tailors the techniques of power geometry to determining regions in Sect. 2.3.4 and we discuss the corresponding Mathematica package that automates the determination of region.

Algorithm: ASPIRE
The algorithm proposed by Bruno and Batkhin [19][20][21] is a very powerful method for obtaining the asymptotic behavior of algebraic curves near singular points. For our purposes however, we need to extract only the leading order scaling behavior of the alpha parameters with respect to ρ which is the expansion parameter. Therefore, using Bruno's theorem [20], we can conclude that the truncated polynomials must have a solution of the form where: where ρ = 1/t so that as t → ∞, ρ → 0. In cases when multiple Alpha parameters scale differently one has to use an approach that systematically finds the leading order expansions of solutions of the Landau equations and thus revealing all the regions. Using the results of Lee and Pomeransky [30], we obtain the following parametric form of a Feynman integral, 1 where Since G is not quasi-homogeneous in general, one can obtain faces of the convex hull of the support of G that correspond to different planes, which in the end leads to different regions. However, we note that the singular points of F now become regular zeros of G.
We now enumerate the steps of our algorithm.
2. Find the Gröbner Basis for the set of Landau equations for F. 4. Using the definition of the small threshold expansion parameter (x), in terms of the kinematic invariants, reexpress all the constant coefficients like mass and external momenta appearing in the above equations. 5. For every transformation applied to the Alpha-parameters, find the support of G which has the structure, S(G) = (Q 0 , Q), where Q 0 is the vector exponent of the small expansion parameter x and Q = (q 1 , q 2 , . . . , q n ) are vector exponents of the Alpha-parameters. 6. Find the convex hull of the support. 7. Find the boundary subset for every facet of every Newton polytope. 8. Find the normal cone for each of the facets. This amounts to finding the normal vector to the surfaces. 9. Using the theorem from Bruno [20] we conclude that the above truncated polynomials are satisfied by the following expressions for the alpha-parameters Here a i ∈ C and the set of P = {p i } defines the region. The normal to the surface corresponding to the truncated polynomial gives us P.
The scaling of the threshold parameter with respect to itself, which we call the zeroth component of the normal, is by definition 1. This is ensured by simply rescaling the normal which is possible as long as the zeroth component is not zero.
Jantzen et al. [13] attribute the fact that the potential and the Glauber regions were missing in ASY, to the cancellations amongst the Alpha-parameters themselves and thus try to resolve it by performing transformations in the Alphaparameter space which eliminate all differences between the Alpha-parameters. In our algorithm we identify these transformations by studying the Gröbner basis elements. In Sect. 3, we demonstrate the working of our algorithm ASPIRE via examples.
It is worth noting here that the set of scalings we obtain correspond to asymptotic expansions near the singular points. If the expansion corresponds to regions where the Alpha parameters are far from zero then Norton-Coleman analysis tells us that the Landau equations can be satisfied only if the internal lines are put on-shell which corresponds to a physical particle traveling from one vertex to another. However, if the expansion is in a region where the Alpha-parameters are zero (or close to zero) then the Landau equations can be easily satisfied and such scenarios correspond to shrinking the internal lines and creating effective vertices, which can be interpreted as the emergence of effective field theories. According to our analysis, the scalings obtained from the bottom facets of the Newton Polytope correspond to the latter case. The physi-

One-loop examples
We first consider some one-loop examples already discussed in Ref. [13].

Example 1: Two-point one loop diagram
Consider the integral (Fig. 5), Here q is the external momentum, k, the loop momentum and the threshold expansion parameter is defined as y = m 2 − q 2 4 .
The energy scales involved in this integral are set by: q, y q and √ y. From [7,13], the contributing regions are the hard and the potential regions whose scaling in the momentum space with respect to the threshold parameter are given as, Potential region : We will reproduce all those contributing regions using the algorithm developed in Sect.
yielding an output The first and second elements are the U and F polynomials respectively, while the third element of the output is the number of loops, which is 1 in this example. There are two Alpha parameters, denoted as x 1 and x 2 corresponding to the two denominators in Eq. (31). For ease of reference, we list U and F polynomials below: We wish to find the locations of the singularities in the Alphaparameter space with the help of Landau equations, Next we find the Gröbner basis of the set of Landau Equations for which, we use the Mathematica function GroebnerBasis via the command which gives the elements, The F polynomial given in Eq. (37) can be written in terms of the elements in Eq. (41). The simultaneous zeros of F and its first order derivatives define the singular points which in general coincide with the zeros of the Gröbner basis elements. As seen in Fig. 6, the solution curves of the Gröbner basis elements partition the Alpha-parameter space and so one can now choose different neighborhoods for studying the leading behavior of F or equivalently G. To perform the expansion in the neighborhood of the solution curve of third element of Gröbner Basis i.e. q 2 (x 1 − x 2 ), we define a set of linear transformations x 1 → ax 1 , x 2 → x 2 + ax 1 and x 1 → x 1 + ax 2 , x 2 → ax 2 respectively. Under these transformations x 1 − x 2 → x 2 = 0 and x 1 − x 2 → x 1 = 0 respectively. We can now expand F or equivalently G in the variable x 2 or x 1 . In all of the above calculations we need to keep in mind the constraint x i ≥ 0. These transformations are analogous to the approach in ASY2, where linear transformations were performed when the Alpha-parameters appeared with a negative sign between them. Such transformations reveal the potential and the Glauber regions. We now list all distinct transformations: • Identity transformation: • Non-trivial transformations: In the above list of transformations, we have fixed the constant a by ensuring that the transformations leave the delta function in the integral, unchanged. Now we go on to compute G with all the above transformations T = {T 1 , T 2 , T 3 }, where the G = {G 1 , G 2 , G 3 } corresponding to the three transformations. Therefore: and Here, we have substituted y → x and q 2 → x 0 (i.e. q 2 → 1). We next find the support of the G polynomials. Here we consider the threshold parameter, x, as an independent coordinate, and therefore, while extracting the vector exponents of the Alpha-parameters x i , we extract the exponents of x as well.
The support S i of G i , where i enumerates the three polynomials coming from the three transformations are, Eur. Phys. J. C (2019) 79:57

Fig. 7
From left to right -Newton polytopes for G 1 , G 2 and G 3 Each row is a point in R 3 and we denote each row as P i . The next step in the algorithm is to determine the convex hull of the support in R 3 for which, we use the function CHNQuickHull [34] as follows Finally, we can obtain the scaling of the Alpha parameters and hence the regions. For G 1 , the vertices are, 0, 1, 0), P 2 (0, 2, 0), P 3 (1, 2, 0), P 4 (0, 0, 1), as seen in Fig. 7 in the left-most panel labelled (a). We then go on to find the normal vectors corresponding to each of the surfaces using the function genNormalCoordinates and obtain (2) are the components of the normal vector of the facets of the Newton Polytope and c is a constant. The presence of the element Null implies that the code was not able to find any normal vector based on the conditions 1 and 2. In our code, we determine the normal vector corresponding to the facets of the Newton polytope based on following conditions: 1. r . v = c and r . v < c, where r belongs to a boundary subset of Newton polytope and r does not. We name the surface giving the normal vector depending on this condition "the top facet" and assign a label "sur f → 1" for that surface. 2. r . v = c and r . v > c. For this condition, we call the surface giving the normal vector "the bottom facet" and label the surface by "sur f → −1".
We use the function UniqueRegions (explained in the Appendix) to select only the unique normal vectors and hence the unique regions. In the above list, we only have one region, which is the hard region {0, 0}.
For the polynomial G 2 , we have six vertices, out of which we get two unique normal vectors from the bottom facets. Therefore, we have two regions are {0, 0} and {− 1 2 , −1}. Similarly, for the support of G 3 Finally, taking the union of all of the above regions we have the following set of regions: ⎛ where the first entry with (0, 0) scaling corresponds to the hard region, while the other two − 1 2 , −1 and −1, − 1 2 are the potential regions, in agreement with ASY/ASY2. As mentioned earlier, we defer the discussion of the components of the normal vectors obtained from the top region to a future publication.

Example 2: A five point one loop diagram
In this example, we consider the following one loop five point integral seen in Fig. 8, which has also been discussed in [13].
Here we have p 1 = p 2 = p and q 1 = q 2 = q. At threshold, p 2 → 0, q 2 → 0, 2 p q → Q 2 where Q 2 is the hard scale and m 2 Q 2 . The threshold expansion parameter is, x = m 2 /Q 2 .
The Symanzik polynomials are, and the Gröbner basis elements for the Landau equations are given by, From the elements of the Gröbner basis we can immediately conclude that we will need transformations listed in Table 2.
Next we apply these transformations to the G polynomial that we construct from the Symanzik polynomials to get four versions for each of the four transformations. We find the Newton Polytope of the support for each of the G polynomials and in each case, we determine the normal vectors listed as follows: Eur. Phys. J. C (2019) 79:57  [13,35].
It needs to be noted here that the regions we obtained above correspond to choosing a specific orientation of the Newton Polytope (for a discussion on rotation of the Newton Polytope in the alpha parametric space the reader is referred to Appendix A.2). It is possible, by arbitrary rotations, to generate an infinite set of scaleless/scaleful regions. However, for a given orientation there is only a finite set of regions. In our implementation, we determine the unique set of regions for a fixed orientation.

A two-loop fish diagram
We apply our technique for the two-loop fish diagram, seen in Fig. 9, that has been discussed in [6] in the context of expansion by regions for the π −K scattering at the threshold. This process can be studied as the scattering of a pion having mass m and a kaon having mass M. The momenta of the pion and kaon are p and P respectively.
For this type of diagram, there are two mass scales (m and M) based on which the scalings of two loop momenta, k and l, one can have the following five possible regions as discussed: In [6], it has been discussed that the h-h region starts contributing at order one while the s-s region at order m M . The h-h' and h-s regions contribute at order m 2 M 2 . The integral having an internal pion loop is given by At threshold, p 2 → m 2 , P 2 → M 2 and ( p + P) 2 = (m + M) 2 . The expansion parameter is x = m M . For the integral in Eq. (59), we get the following Symanzik polynomials, As mentioned in our algorithm, we find the Gröbner basis elements of F and its derivatives with respect to We use the elements of the Gröbner basis, together with the constraint x i ≥ 0, to obtain five transformations out of which one is trivial and others are non-trivial and are listed below: • Identity transformation: • Non-trivial transformations: As before, we compute G polynomials by applying all of these transformations, determine the support in each case and the corresponding normal vectors. This leads us to the following list of unique normal vectors: ⎛ The regions are (once again from the bottom facet) are: The region {0, −2, 0, 0} is also isolated by ASY. This region has not been identified in the existing literature to the best of our knowledge, hence, we give them a generic name of "scaleful region".

One loop scalar triangle diagrams in Sudakov limits
As a final illustration of our algorithm, we discuss the case of one-loop Sudakov integrals in the following limits, illustrated in Fig. 10. Sudakov limits appear in processes where the square of the momentum transfer is large compared to the squares of the masses. A detailed discussion may be found in the chapter on Sudakov limits in Ref. [24]. See also Ref. [1].

Limit (a)
In this limit the integral becomes, We define the threshold expansion parameter for this limit as x = m 2 /Q 2 , where Q 2 is the large scale.
The U and F polynomials in terms of the alpha parameters are respectively given by (73) and the Gröbner basis elements for the Landau equations that F satisfy are, Clearly, we have only the Identity transformation : } and the normal vectors we get in this limit are ⎛ which in turn correspond to the following regions: ASY also reports the same regions. In [24], it is confirmed that only the hard and collinear contributions suffice for the evaluation of this diagram at leading order.

Limit (b)
In this limit, we have the leading order integral The threshold expansion parameter is The Symanzik polynomials and the Gröbner basis for the Landau equation satisfied by F are respectively, The Gröbner basis of Landau equations are, Once again, we have only the Identity transformation and the normal vectors identified are ⎛ We see in this case also the complete agreement of our result with ASY and also with the contributions, reported in [24].

Limit (c)
For this limit, the integral is, with the threshold expansion parameter x = m 2 /Q 2 .
The Symanzik polynomials U and F are, and and the Gröbner basis elements for the Landau equation satisfied by F is given by, using which, we once again see that only Identity transformations are required. Finally, the distinct normal vectors corresponding to the facets are ⎛ The regions are:

Limit (d)
In this limit, the integral is The threshold expansion parameter here is x = m 2 /M 2 . The Symanzik polynomials are The Gröbner basis elements that generate the same ideal as the Landau equation for F is given by, Once more with only the Identity transformation, The scalings above correspond to the following regions: (i) The hard region : {0, 0, 0} (ii) The 1-collinear region: {0, −1, −1}.

Discussion and conclusion
The MoR is a powerful technique for obtaining expression of Feynman amplitudes at any order. The field however needs a robust algorithm for systematically finding all the regions.
The study of such algorithms may provide insight into ideas crucial for validity of the MoR. In this work, we present an algorithm to identify the regions using ideas from power geometry, which is a powerful technique for analyzing properties of algebraic polynomials. Our algorithm, ASPIRE, has allowed us to develop an implementation in Mathematica where we have also used external programs including UF.m and NDConvexHull.m. We have benchmarked the code by reproducing one and two loop examples from the literature. The salient steps in developing this algorithm can be summarised as follows: 1. We translate the traditional problem of finding the regions in the momentum and alpha parameter space to purely in the Alpha-parameter space by integrating out the loop momenta and then expanding the resulting integral in the new regions. 2. Using the form of the Landau Equations in the Alphaparameter space, we reduce the problem to finding the leading order behavior of solution of the Landau equations in the different regions. Instead of working with F or U, which are homogeneous, we use the polynomial G = U + F, as defined by [30], which has the advantage of not being homogeneous, while the truncated polynomials are quasi-homogeneous. 3. These regions lie in the neighborhood of the origin and we systematically expand the integral in these neighborhoods via simple transformations which are obtained from the analysis of the Gröbner basis. These transformations are strongly constrained via the delta function in the parametric integral which forces α i = 1, where α i are the Alpha-parameters. 4. The expansion of the integrals in these neighborhoods is obtained via the use of an old but unvisited topic of Power Geometry. We use only a small portion of this powerful technique to speed up our algorithm for finding the regions in the Alpha-parameter space. This technique allows one to obtain the leading order expansion of the integral in a particular neighborhood by analyzing the support of a G polynomial, and constructing its Newton Polytope.
5. The regions are then defined as scalings of the Alphaparameters which lead to the leading order expansion of the integral in that particular neighborhood of the origin. 6. The scalings are obtained by finding the normals to the surfaces of the Newton Polytope with a constraint on the first component of the normal vector. 7. We demonstrate that linear transformations in the Alphaparameters lead to non-trivial transformations of the Newton Polytope, which leads to the uncovering of previously hidden regions. 8. We present a stand-alone documented Mathematica implementation of the above developed algorithm and provide several examples covering one loop and two loop amplitudes.
Our results are in agreement with previous work in this field by Jantzen et al. [13] and Pak and Smirnov [11]. We also settle the issue of PreResolve as raised in the work of Jantzen, Smirnov and Smirnov [13] and provide a Mathematically justified way of unveiling the full set of transformations that one needs, to successfully identify all the regions in the Alpha parameter space.
Future extensions of this work includes looking into the connection between the sub-leading contributions of regions and the resulting geometry in the Alpha-parameter space as well as predicting the number of distinct regions a priori by studying the topology of amplitudes. For the moment, this work, as indeed is the case with the work of Pak and Smirnov, Jantzen, Smirnov and Smirnov, is based on the U and F polynomials. Having more than one approach based on these may be profitable in the sense of eliminating the possibility of missing some regions when one or the other is used at one given time. The work here is based on the general properties of Landau singularities and closer to the spirit of the classical analyses of Feynman amplitudes to analyze their analyticity in the past, is now being employed to identify the regions associated with Feynman diagrams.
The investigation of multi-loop non-planar vertex diagrams encountered in references, e.g., [36] and also in the book of Smirnov [24] are an interesting class of amplitudes where our initial analysis suggests a rich family of solutions of the Gröbner Basis elements. Treatment of such non-trivial integrals would require automation of the analysis of Gröbner Basis elements for determining the complete set of transformations which will lead to the identification of all contributing regions. As a preliminary example we solve Landau equations and the Gröbner basis of the non-planar vertex diagram considered by [24] in the Appendix A.5. It may be readily seen that this is highly complicated and we do not yet have the full solution for this system, and work is in progress.
It may also be noted here that the present work, as in work of Pak and Smirnov, Jantzen, Smirnov and Smirnov is limited to dimensionally regularized Feynman diagrams. We only provide the regions that maximally partition the Alpha parametric space at leading order. With the regions in hand one still needs, for a complete evaluation of the original integral at leading order, to use appropriate phase space regulators, as discussed in [35], and follow the procedure of Jantzen [22] to perform zero-bin subtractions in the framework of dimensional regularization. In the past, analytic regulators above and beyond dimensional regularization have been employed in the MoR approach to separate soft and collinear regions, see e.g. [35]. The presence of the analytic regulators has been automated in ASY code as well as in FIESTA [11,37]. Such an extension to the ASPIRE algorithm would be very interesting. although beyond the scope of the current work.
An important question to ask is whether the ASPIRE algorithm can find regions that the ASY algorithm cannot? In order to answer this question, we have visited some examples that have been studied using ASY. These include the diagrams found in Ref. [6] all of which have been analyzed by us using ASY, correspond to (a) the J type one-loop integral, (b) the fish diagram with the kaon loop, and (c) the fish diagram with the π − K loop. When implemented on ASPIRE the results have been confirmed. In addition we have considered the two loop Master Integral vertex diagram in eq. (7.30) of Ref. [24] and find agreement between the results from ASY and ASPIRE. To this extent, at the present level of investigation we find complete agreement. It may yet be that ASPIRE has the potential to probe new regions in the setting of nonplanar diagrams not necessarily in the Sudakov limit. For the moment, this analysis has not been performed partly due to the highly non-trivial Gröbner basis elements. Such investigations are deferred to the future.
All the work reported here has used the Alpha parametrization. While it has been convenient to analyze the regions here, the connection to the actual scaling behavior in the momentum space is less transparent. In order to actually assign the nature of the scalefulness to an isolated region, one has to go back to the momentum scaling behavior. This identification has been carried out by hand. This part of the algorithm needs to be automated as well. Rotation of the Newton polytope corresponds to making a linear transformation of the support which translates to redefinition of the alpha parameters themselves. For the above let us make the transformation In the above we performed the transformation The new polynomial we get is In the above it needs to be noted that the matrix is not orthogonal. The only condition that needs to be satisfied by a matrix, M representing a rotation of the Newton polytope is M T M = cI, c ∈ R and I is the identity.
The new support of polynomial G (x, y) is The new Newton polytope for the new support is (Fig. 12): It is evident that rotation of the Newton polytope does change the normal vector of the surfaces. However, the boundary subset of the Newton polytope remains the same and thus the rotated and unrotated normal vectors would give the same boundary subset and hence the same truncated polynomial which implies that they both correspond to the same region. This can be seen Mathematically. Let us take the boundary subset of the Newton polytope S 1 j . Now for a normal vector n we have, Rotating the polytope corresponds to Thus, we now have the condition, Thus, we see that at the end of the transformation the boundary subset does not change. Hence the truncated polynomial still does not change and the truncated polynomial integrates to the same expression as the unrotated one. After making such a transformation, of course one needs to also include the jacobian in the integral and also change the limits of integration if necessary. One interesting thing to note here is that one can generate an infinite set of regions by rotating the polytope. These look different from each other at the outset but in fact represent the same process. The purpose of this discussion is to show that the invariances of the Newton polytope under the rotations above, does not change the regions identified in a substantive manner, in the case of the analysis of the U and F polynomials. Rather the effect of the rotations is to reexpress the rotated regions in terms of the original ones. This does not give any further information on scaleless regions, are required to be eliminated in any event.

UF.m
UF is a Mathematica based package designed for extracting the Symanzik polynomials, U and F, from the alpha represen-tation of any multiloop Feynman integral. The code is openly available for redistribution [33]. In our codes, UF function of the package UF.m takes as input, for a particular Feynman integral, the set of loop momenta, set of all propagators and the set of kinematical substitutions as input. Example: gives the output [1]x [2], 1}, where x [1] and x [2] are the Alpha-parameters. The output has three parts: first and second elements represent the U and F polynomials respectively and the third entry is the number of loops.

NDConvexHull.m
This package includes implementations of several algorithms like Chan's algorithm, gift wrapping, quick hull, incremental convex hull, for finding the convex hull of a set of points in multi-dimensional space. They take a set of points as input and return a sorted list of vertex points and a sorted list of simplexes as output. This package was developed by Loren Petrich [34].
In this work, we use the function CHNQuickHull for finding the convex hull of a set of points using the Quick Hull algorithm. Example: Consider the following set of points Output: where the first list is the set of point ids of the given points and the second matrix is the surfaces of the Newton polytope in terms of those ids.

A.4 Description of mathematica functions implemented for this work getMul
We use this function for finding the support of a given polynomial. This function returns the vector exponents of the corresponding variables in a given monomial term.

UniqueRegions
It takes as input the set of normal vectors and the length of the normal vectors, removes any instance of Null from the set and then eliminates all normals related by a constant shift. This gives us a unique set of regions Output: Scaleless.
From the above list, let us analyze the Gröbner basis element: In the above element the factor (x − 1)x (x 4 − x 5 ) (x 4 + x 6 ), becomes zero if x 4 − x 5 = 0 or x 4 + x 6 = 0. The positivity of the Alpha-parameters implies that x 4 + x 6 = 0 gives x 4 = x 6 = 0. Thus, this solution branch corresponds to a trivial transformation.
More complicated Gröbner basis elements have multiple solution branches. All of these branches look distinct at the beginning but as one does a careful study of these, it might be revealed that these seemingly distinct branches lead to similar transformations. To find the complete set of transformations needed to detect all the regions an automated approach to analyze and identify all the distinct solution branches is required and is deferred to future versions of our algorithm. While analyzing this system in our scheme by considering only the Identity transformation (i.e. the transformation which leaves the Alpha parameters unaltered), we get the same regions as obtained from ASY. However, the complicated nature of the Gröbner bases may allow for other nontrivial transformations which can be studied in the future. A notebook for the preliminary study has been provided.

A.6 List of mathematica notebooks used in this work
We provide brief description of the Mathematica notebooks used in this discussion. It may be noted that care has been taken to produce the U polynomial with the correct positive sign denoted by the extension "P" at the end of the name of the notebook.