Cutkosky representation and direct integration

We present a new method of direct integration of Feynman integrals based on the Cutkosky representation of the integrals. In this representation we are able to explicitly compute the integrals which yield square root singularities and leave only the integrals which yield logarithmic singularities, thus making the transcendentality weight manifest. The method is elementary, algorithmic, does not introduce spurious non-physical singularities and does not require a reduction to a basis of pure integrals.


Introduction
This paper grew out of trying to understand two basic facts about Feynman integrals.The first fact is that a large class of Feynman integrals at L loops and in D dimensions can be written as iterated integrals of length ⌊ DL 2 ⌋.This is less than or equal to half of the number of integrals in momentum space.It seems that, outside of a few experts, this fact was not widely appreciated and, in fact, before the work of ref. [1] we lacked even the language to properly discuss transcendental weight in the Feynman integral literature. 1 It is hard to see who should get the credit of this simple but important observation, but Nima Arkani-Hamed has forcefully made this point to me in a private discussion.This fact is not only surprising, but actually very useful practically since it reduces in half the number of integrals one needs to perform.
The second basic fact is that the iterated integrals are a sequence of one-dimensional integrals.This also looks very surprising.How can one rewrite the original integral as a sequence of one-dimensional integrals?And what should these one-dimensional integration variables be?
A first clue to answering these questions arose from an important dichotomy of singularities, noticed already by Landau in his original paper [4] on singularities of integrals.As Landau showed, the singularities broadly divide into two categories: square root2 and logarithmic.These singularities have very different nature.Square root singularities are algebraic and do not contribute to transcendental weight, while logarithmic singularities are transcendental and do contribute to transcendental weight.Furthermore, ref. [5] showed in examples that about half of the singularities are of square root type and half are of logarithmic type.It is then reasonable to guess that it is only the logarithmic singularities which contribute to the transcendental weight and in turn this gives a way to explain Arkani-Hamed's observation.
But the second question remains.What should be the one-dimensional integration variables?A second clue towards answering this come from Cutkosky's work [6] where he described discontinuities across branch cuts in terms of cut integrals.In that reference Cutkosky introduced a new way to write the Feynman integrals, designed with the explicit purpose of making both the singularities and the discontinuities manifest.This representation, as we will review below in sec. 2 and sec.3, is a sequence of one-dimensional integrals potentially followed by a higher-dimensional integral.As a general rule, when this last integral is non-trivial, an example being that of a non-singular elliptic integral, the Feynman integral can not be written as a polylogarithm (we discuss an example in sec.6).
In order to make contact with the iterated integral form of the answer, it should then be possible to explicitly do the integrals producing square root singularities.This would then constitute an alternative way of performing direct integration to the methods of refs.[7,8].
We explicitly show how to do this in a few simple examples.The method relies on nothing more than basic complex analysis and Cauchy theorem.As we will show, in order for this to work, there should be only one pair of square root branch points at each step involving square root singularities.This does indeed happen, sometimes via some non-trivial algebraic identities.The examples presented in this paper are not meant to challenge the state of the art computations, but rather to showcase how the method works in simple examples.We hope to tackle more complicated cases in future work.
In sec.4.3 we discuss the case of a reducible integral and show that applying the same ideas to this case poses no difficulty.Therefore, this method, unlike the differential equation method (see ref. [9]), can avoid a potentially expensive integral reduction (see ref. [10]) step.Clearly, applying the integration algorithm to each Feynman integral separately is not economical.Instead, as in the unitarity method (see refs.[11,12]), one can group together all the diagrams which contribute to a given Landau singularity (at a given order in perturbation theory), compute the on-shell state sums, take the internal momenta offshell and compute the resulting integral.The full answer can be obtained by merging (not adding!) different contributions which account for all potential singularities.
The integration method we introduce has the following advantages over approaches in Feynman parameter space.First, it can apply to a large variety of mixed types of iϵ conditions (advanced, retarded, Feynman, anti-Feynman) which are sometimes required.This is so because, unlike in Feynman parameter space, we have several denominators and we can choose contours for each independently.Second, the momentum space method applies effortlessly to cut integrals, which are more difficult in Feynman parameter language (see ref. [13]).A third advantage is that the on-shell varieties arising in momentum space are typically less singular and do not require as many blow-ups (see refs.[14,15] for examples of blow-ups required in Feynman parameter space).A more technical difference is that (the properly compactified) contours of integration in momentum space are not relative homology classes so are easier to deal with.But most importantly, a huge advantage of our method is that we only need to think about one variable at a time and in principle all the singularities in that variable are visualizable in the complex plane of that variable.
The polylogarithmic integrals are fairly well understood and the next frontier is that of integrals which are not polylogarithmic.If the on-shell space of the leading singularity has a non-trivial topology, such as that of a (non-singular) elliptic curve or a Calabi-Yau variety, then the integrals can not be computed in terms of polylogarithms.Sometimes, as in the case of the bubble integral in three dimensions, the on-shell space has the topology of a circle.When complexified, the circle becomes non-compact and adding two points to achieve a compactification amounts to adding two singularities of pole type.
In refs.[16,17] a formalism for defining a coaction has been developed by using cuts instead of differentials.In principle this can be applied to integrals of elliptic or Calabi-Yau type.However, the entries of the ensuing symbol will not be as simple as in the polylogarithmic case.It is therefore not clear yet how to use this method for writing the answer in a canonical form or how effective this method can be in that case.Indeed, unlike for the polylogarithmic case, even the notion of a prefactor for the integral does not seem to be well-defined (see ref. [18]).Other approaches have been proposed in refs.[19][20][21].

Cutkosky's argument
In ref. [6], Cutkosky described a change of variable from the usual loop momentum integration variables to q 2 e , where q e are internal (not necessarily independent) loop momenta.One can change variables from k i , independent loop momenta to q 2 e and other "angular" (in Cutkosky's terminology) variables ξ.
After this change of variables the integral reads where J is a Jacobian factor.It can potentially contain numerator factors of the original integral as well.
In favorable cases, the last form

J
has further residues and its integral can be computed in simple terms.In more complicated cases, this last form is a holomorphic form on elliptic curves (or hyper-elliptic curves) or Calabi-Yau manifolds and γ in eq.(2.1) is a real homology cycle.Despite much study, a general theory of integrals of elliptic or Calabi-Yau type is not available yet.
Then, Cutkosky describes the integration limits as the solutions to a modified form of Landau equations, j≤i β j q j = 0 where q j have norms fixed by the values of the outer integrals.This relies on some (arbitrary) ordering of the propagators.Obviously, a judicious choice of ordering can simplify the calculations.
It is worth pointing out that the q 2 e integrals in Cutkosky's representation (2.1) have a superficial resemblance to G-functions However, the G-function representation is more restricted since the boundaries of integration depend in a much simpler way on the previous integration variables.In Cutkosky's representation the integration boundaries have a potentially complicated functional dependence on previous integration variables a s (q 2 1 , . . ., q 2 s−1 ) and b s (q 2 1 , . . ., q 2 s−1 ) instead.Nevertheless, Cutkosky's representation has one important qualitative similarity to the G-functions: in both cases we are dealing with one-dimensional integrals on Riemann spheres P 1 .
In this representation of the integral, the singularities arise as follows (see ref. [6, eq. 9]).If we denote the result of doing all integrals except the outer one by F (1) (q 2  1 , p) where p are external kinematics, then the full answer is If there is a singularity when p → p 0 , then this means that the integration contour in the q 2 1 complex plane is pinched between q 2 1 = m 2 1 and a singularity of F (1) .More precisely, we must have that F (1) (m 2  1 , p) is singular when p → p 0 .By a contour deformation we can pick up a residue at q 2 1 = m 2 1 and this is the only part of the integral which is singular.Therefore, the singularity is given by ±2πiF (1) (m 2 1 , p).As remarked above, this indeed becomes singular when p → p 0 .

Cuts and spectral densities
Given a function ρ(x), we can build a function F (x) with a branch cut between a and b whose discontinuity is ρ(x).This construction is well-known and Indeed, we have Kinematics of the bubble integral.
where we have used 1 x∓iϵ = pv 1 x ± iπδ(x).Notice the obvious similarity between eq.(2.3) and eq.(3.1).The function F (1) in eq.(2.3) is the cut of the function defined by the integral.The function F (1) itself can be represented by a similar integral and the function F (2) it contains corresponds to another cut, where more propagators are set on-shell.In turn, this writing looks similar to the writing in terms of G-functions, with one major difference: the number of integrals in the G-function representation is equal to ⌊ LD 2 ⌋, while the number of integrals in the Cutkosky representation is at least the number of propagators.We will show below that by explicitly integrating square root singularities and the "angular integrals" one can make the number of integrals match.
Clearly in the representation of eq.(3.1) z = a and z = b are branch points.If ρ(a) or ρ(b) are non-vanishing finite constants then we have logarithmic branch points.It is also possible to have singularities such as ρ(z) ∼ (z − a) γ when z → a, for ℜγ > −1 (to ensure convergence) and similarly for z → b.In this case the value of ρ at the branch points is either zero or infinity, depending on the sign of ℜγ.
If γ is a half-integer then we have a square root branch point at z = a.Obviously, the type of branch cut (i.e.square root versus logarithmic) must match at z = a and z = b.A common form for ρ which we will encounter in the following is ρ(z) = This is closely connected to the Mandelstam representation (see ref. [22]).The Mandelstam representation has been the subject of many studies (see also ref. [23] for an application in a similar context to the present one).Our proposal amounts to building some kind of spectral densities in perturbation theory.However, unlike in Mandelstam's approach, we explicitly integrate the square root cuts and keep only the logarithmic singularities.For example, the double-spectral function for a box diagram is of square root type (see ref. [24, eq.B-42, p. 215]) which in our approach would not survive the integration.

The bubble integral
Consider the bubble integral where p = q 1 + q 2 and p is the external momentum (see fig. 1).This can be rewritten as The integration region for the bubble integral in the Euclidean region is inside the curve.It is outside the curve, including for negative values of q 2 1 and q 2 2 in the Lorentzian region (for p 2 > 0).
where (q 2 2 ) min and (q 2 2 ) max are the minimum and maximum values of q 2 2 , subject to the constraints that q 2 1 is fixed and p = q 1 + q 2 .The integration domain is in fig. 2.

Euclidean signature
In Euclidean signature we have where ϵ(v, w) = v 0 w 1 − v 1 w 0 .We have Then, the integral to compute becomes where ∥p∥ = p 2 and we have used the fact that the minimal value of q 2 1 in Euclidean signature is zero and its maximal value is infinity.Once the value of q 2 1 is fixed, the minimal value of q 2  2 is (∥p∥ − ∥q 1 ∥) 2 , obtained when q 1 and p are aligned and the maximal value is (∥p∥ + ∥q 1 ∥) 2 when q 1 and p are anti-aligned.
The inner integral reads where with a < b and c ̸ ∈ (a, b) can be computed as follows.We introduce a curve y 2 = (b − x)(x − a) which can be rationally parametrized by t as follows Then, the integrand can be written as The value x = a corresponds to t = 1 while x = b corresponds to t = −1.Then we obtain b a dx (4.11)The logarithm contributes π, so the transcendental weight is purely numerical and does not have any dependence on the kinematics.
The same integral can be done by contour integration.Let us briefly describe the method since it will generalize to more complicated cases.We want to define a function Then we have where γ is a contour as in fig. 3. We have We therefore have Finally, we continue with the evaluation of the outer integral with The integral is along the real axis where the quantity under the square root is positive so, with the usual convention for the square root branch cut, the integration path does not intersect the cut.
To compute the integral we introduce a curve y 2 = (x − a) 2 + b 2 .This curve can be parametrized rationally by Then we have where To finish the computation of the integral we need to compute the values of t at the boundary of the integration region.At x = 0 we have y = √ a 2 + b 2 , where we need to pick the positive sign for y.This implies that whence . For x → ∞, we also have y → ∞.The condition x → ∞ follows from t → ∞ or from t → 0 − , but in the second case we obtain y → −∞.
Finally, the integral we wanted to compute becomes Plugging in the values for a, b and c we find that the bubble integral in Euclidean signature reads where

Lorentzian signature
In two-dimensional Lorentz signature we have where we have used dq 2 = −dq 1 (since p is considered constant) and fact that (v Here we have defined and It can be checked by a simple calculation that where Let us first compute the stationary points of q 2 2 subject to the constraints mentioned above.Using the Lagrange multiplier λ 1 , we find which reads Using momentum conservation this implies that q 1 = 1 1+λ 1 p, q 2 = λ 1 1+λ 1 p.We can then determine . Using this value of λ 1 we find for the stationary points.At this stage we don't know yet if these are true minima or maxima.To decide the nature of the stationary points and find the minima and maxima one can follow the general procedure described in sec.6 which involves computing a bordered Hessian (see eq. (6.14)).Their nature depends on the signs of p 2 and q 2 1 .In this case we do not need the full power of a general theorem since a direct analysis of the inequalities suffices.Going back to the equation eq. ( 4.27), we see that for real Lorentz kinematics we have (q Let us assume p 2 > 0.Then, if q 2 1 < 0 the inequality is satisfied for all values of q 2 2 .The same holds if p 2 < 0 and q 2 1 > 0. But if p 2 > 0 and q 2 1 > 0, then we have It follows that the possible values for q 2 2 are either If instead p 2 < 0 and q 2 1 < 0, then we have Then, either The end-points of the integration domain (except the ones at infinity) are the same as obtained by the stationary point study.
The q 2 2 integral has two forms where the roots a 2 and b 2 are not real, and where the roots a 2 and b 2 are real with a 2 < b 2 .These integrals can be seen as integrals along contours in the curve , which is a double cover of the complex x 2 plane, branched at two points x 2 = a 2 and This curve can be rewritten as This provides a rationalization of the curve and is a useful change of variables in the integral.
In terms of the coordinate t 2 we have Then, we have Then we have The square root prefactor is now independent on t 2 and can be combined with the outer differential form We will also set x 1 = q 2 1 for brevity.Replacing these in the square root we find with This can be treated in the same way as before.Indeed, we can introduce a variable y 1 defined by ) and a uniformizing variable t 1 .In the end we get 1 The square root in front can be written as which is the familiar Källén function.
The integration domain is outside the curve in fig. 2.More precisely, for q ), ∞).For q 2  1 ∈ (−∞, 0] we can do the inner integral and find that, as function in q 2 1 it has no singularities in the upper half plane (it has a logarithmic branch cut along the positive real axis and there is also a pole at q 2 1 = m 2 1 − iϵ).In particular there is no pole at infinity and the contour in q 2 1 along the negative real axis can be rotated clockwise to sit on the positive real axis.Changing the direction of integration introduces a minus sign and combining with the previous integration along q 2  1 ∈ [0, ∞) produces This has the same form as in Euclidean signature, but this form is a result of a cancellation between different regions.
Here we have two integrals but we expect only a single logarithm.Therefore, it should be possible to do one of the integrals and get a rational multiple of 2πi.Since the integration contour in x 2 variable goes between a 2 and b 2 , then in the t 2 variable it goes between The points x 2 = a 2 and x 2 = b 2 have a unique point in the double cover and therefore each corresponds to a unique value of t 2 .Then, we have . This cross-ratio is very special, since ( t+ ) where we have used the expressions in eq.(4.43).
Then, the integral becomes where the sign depends on the determination of the logarithm.The final answer for the integral resembles the one of the Euclidean signature.The biggest difference is that we need to replace to obtain the square root in Lorentzian signature.This replacement is the same as the one arising from a Wick rotation.

Bubble integral in three dimensions
As an example where the inner "angular" integral is not zero-dimensional, consider the case of a bubble integral in three dimensions.As we will see, this integral yields a simpler answer than the bubble integral in two dimensions.
For simplicity we compute this integral in Euclidean signature.We have which can be written as The inner integral is over the space of triangles with side lengths ∥p∥, ∥q 1 ∥ and ∥q 2 ∥ in a three-dimensional space, where the vector p is fixed.We have
We are left with the integral Notice that this integral is odd under ∥p∥ → − ∥p∥.The prefactor also contains a ∥p∥ which is also odd under this transformation so overall the integral is invariant under ∥p∥ → − ∥p∥.
The inner integral can be computed straightforwardly, which leaves us with We rewrite the integrand so that it has a cut along the integration region Here the square root has the principal determination √ z = |z| exp( i 2 arg z) with arg z ∈ (−π, π).This means that along the real axis −q 2 1 = i q 2 1 and ℜ √ z ≥ 0 for all complex z.
We define (4.60) In principle we could have chosen the function ρ in several other ways, but this form also has the property that ρ(0) = 0 and lim q 2 1 →∞ ρ(q 2 1 ) = 0. Then the original integral can be written as where the discontinuity is across the branch cut along the positive real axis and the factor of one half is due to the fact that computing the discontinuity across the branch cut doubles the value of the function ρ just above the cut.We have where γ 1 is a contour from R to ϵ slightly displaced from the real axis in the lower half plane, γ 2 is an arc of circle of radius ϵ and center 0 which continues the path γ 1 and γ 3 is a contour Integration contour for the last integral of the bubble integral in three dimensions.The horizontal cut is of square root type while the vertical cuts are of logarithmic type.There is also a pole at from ϵ to R slightly displaced in the upper half plane.This contour can be completed to a closed contour by adding a circle of radius R and center 0, but paying attention to the logarithmic branch cuts of the function ρ.The resulting contour is sketched in fig. 4.
Let us find the logarithmic branch cuts of ρ.The branch cut condition is The branch points in q 2 1 arise when either the numerator or the denominator of the argument of the logarithm vanish.This condition can be rewritten as which can be rewritten as Solving this we find the logarithmic branch points The equations for the logarithmic branch cuts are more complicated.We have schematically represented them by the vertical lines going to infinity in fig. 4. The precise location of the branch cuts is not essential, but we of course require that the integration contours do not cross them.The integral along the small circle γ 2 vanishes in the limit ϵ → 0. The integral along the pieces of the large circle of radius R also vanishes in the limit R → ∞; the term contributes a pole at infinity, but its residue vanishes (here we are using the crucial property that lim q 2 1 →∞ ρ(q 2 1 ) = 0).When q 2 1 = p 2 − m 2 2 + 2im 2 ∥p∥ we have −q 2 1 = ±(m 2 − i ∥p∥) but for the principal determination of the square root we should choose −q 2 1 = m 2 − i ∥p∥.In this case, it is the numerator of the argument of the logarithm in ρ that vanishes.Since the circles around the logarithmic branch points turn in the clockwise direction, we have that the discontinuity across the branch cut ending at q 2 1 = p 2 − m 2 2 + 2im 2 ∥p∥ is −2πi.Similarly, the discontinuity across the branch cut ending at q 2 1 = p 2 −m 2 2 −2im 2 ∥p∥ is 2πi.Therefore, the contribution of the logarithmic branch cuts to the contour integral is where R ± are some complex numbers whose norm is of order R.They are determined by the intersection of the large circle of radius R and the logarithmic branch cuts.In the limit R → ∞ this integral becomes 2πi Equating the value of the integral along the contour in fig. 4 with the contribution of the residue at A non-trivial check on the computation is that the answer should be symmetric under the exchange m 1 ↔ m 2 .This can be simplified to where we assumed m 1 > 0, m 2 > 0 and ∥p∥ > 0.
Including the prefactor of − π 2∥p∥ we obtain the final answer for the bubble integral in three dimensions in Euclidean signature The result is real and positive, but this is not completely manifest.It can be rewritten as This integral is not difficult to compute by other methods (see for example ref. [25, eq. 12]).See also ref. [26, eq. 3.19] for a recent occurrence of the same integral.
Curiously, this form of the answer is not invariant under m i → −m i , which is a symmetry of the original integral.

Triangle integral in two dimensions
Consider the triangle integral in two dimensions.This integral is usually computed by reducing it to three bubble integrals.Each of the bubble integrals has a different prefactor so the integral is not "pure" in the language used in the literature.Such integrals can not be computed by the usual means of taking differentials and they are usually reduced to "pure" integrals using a procedure of integral reduction.In this section, we use this example to illustrate how this new method of direct integration copes with this difficulty.
At the same time, we have (5.6) Hence, when q 2 2 = a 2 or q 2 2 = b 2 both the numerator β 2 3 − 4α 3 γ 3 and the denominator factor det G(q 1 , q 2 ) vanish.An explicit computation reveals that Therefore, we can pull out a factor dependent only on external kinematics and we are left with the integrals (5.9) The integral in q 2 2 can be done by straightforward partial fractioning in 2 ) 2 + β 2 q 2 2 + γ 2 and this polynomial in q 2 2 has roots x ± , then An explicit calculation yields We also have that (5.12) To do the final integral, we proceed as before.First, notice that when ∥q 1 ∥ = 0 or ∥q 1 ∥ → ∞ we have that a 2 and b 2 coincide.This means that the logarithms vanish.In fact, the integrand of the q 2 1 integral has a square root branch point at q 2 1 = 0 and q 2 1 → ∞.It looks like the result of the integration in eq.(5.10) (and therefore the integrand for the q 2 1 integration) also has square root branch points at β 2 2 − 4α 2 γ 2 = 0. Interestingly, these square root branch points are actually canceled in the combination in eq.(5.10).Indeed, when taking q 2 1 along a path such that β 2 2 − 4α 2 γ 2 goes once around the origin, we have -17 - Then the expression in eq.(5.10) is sent to itself.
The integrand has some pole singularities.Indeed, α 2 m 4 2 − β 2 m 2 2 + γ 2 = α 1 (q 2 1 ) 2 + β 1 q 2 1 + γ 1 = 0 when q 2 1 = (q 2 1 ) ± .We have (5.13)There is also a pole when 1 ) − ) when x + + m 2 2 = 0 we have that either q 2 1 = (q 2 1 ) + or q 2 1 = (q 2 1 ) − .Let us assume for definiteness that we have q 2 1 = (q 2 1 ) + .Then, where we have used m 2 2 = −x + .The final integral to do is over q 2  1 and runs from 0 to ∞.This along a square root branch cut so the integral can be written as one half the integral around the branch cut.This contour can be deformed to a sum of contours around the poles described above and around the logarithmic branch cuts.The integrals around the logarithmic branch cuts can be written as integrals along a contour connecting the two logarithmic branch points where the new integrand is obtained by replacing the logarithm by 2πi.After this replacement the integral becomes easy to perform, by using identities such as eq.(A.20).
We will not go through all the trivial (but tedious) steps in detail (see sec.4.3 for more details on the method), except to describe the determination of the locations of logarithmic branch points.From eq. (5.10) we have that the first logarithmic branch points appear for values of q 2 1 where a 2 = −m 2 2 and b 2 = −m 2 2 .Using a 2 = (∥p 2 ∥ − ∥q 1 ∥) 2 and b 2 = (∥p 2 ∥ + ∥q 1 ∥) 2 we find logarithmic branch points at q 2 1 = (∥p 3 ∥ ± im 2 ) 2 .From the second and third terms we have x ± = b 2 which implies that α 2 b 2 2 + β 2 b 2 + γ 2 = 0.When written as an equation in ∥q 1 ∥ this equation is of degree four, so one might worry that we need to deal with roots of order four.Fortunately, it turns out that this degree four equation is very special and its roots can be written using a single square root where (5.16) Each one of these roots appears with multiplicity two.A similar analysis applies for a 2 .
In this section we have demonstrated an algorithm for computing a reducible Feynman integral with multiple prefactors, without performing an integral reduction.Integral reductions (see ref. [27]) are often resource-intensive and require a non-canonical and often symmetry-breaking choice of basis.
One interesting fact we can notice is the occurrence of higher order equations, but so far such that their roots only require quadratic field extensions.It is plausible that by this method of integration no unnecessary field extensions will be required.When computing a particular integral in ref. [28] using HyperInt (see ref. [8]) we encountered field extensions of degree 16 while the final answer was completely rational.Higher order equations are expected to occur (see refs.[29,30]) in general for more complicated integrals.

Triangle integral in three dimensions
It is instructive to consider the triangle integral in three dimensions. 3In this case we expect the answer to contain one logarithm, so we should be able to compute two integrals as rational multiples of 2πi.
We do the integral in Euclidean signature.The integral reads . (5.17) Putting the integral in the Cutkosky form we find (5.18) We will determine a 1 , a 2 , a 3 and b 1 , b 2 , b 3 in the following, but first we compute (5.19) 3 Often the integrals are studied in Feynman parameter space, not in momentum space as we do here.
In Feynman parameter space the integrals in odd dimensions may contain square roots, which complicates their analysis.One of the advantages of the momentum space approach is that it can be done in even or odd dimensions with no differences.
Using momentum conservation we find q 2 = q 1 +p 3 and q 3 = q 1 +p 1 +p 3 which implies that dq 2 2 = 2q 2 • dq 1 and dq 2 3 = 2q 3 • dq 1 since the external momenta are taken to be constant.Next, we find dq 2 1 ∧ dq 2 2 ∧ dq 2 3 = 8ϵ(q 1 , q 2 , q 3 )d 3 q 1 .We have (5.20) Therefore, we have . (5.21) To compute the boundary values, a i and b i we proceed as follows.First, we have a 1 = 0 and b 1 = ∞.Next, have the same problem as for the bubble case in sec.4. That is, given the fixed value of q 2 1 and p 3 = q 2 − q 1 , find the extremal values of q 2 2 .Finally, for the extremal values of q 2 3 there are several possibilities.First, there is a constraint arising from q 2 3 = (q 1 − p 2 ) 2 , or from the bubble with momenta q 1 , q 3 .Second, there is a constraint arising from the bubble with momenta q 2 , q 3 .Taken together, these imply the triangle Landau equations.
Thus, we see the hierarchical principle of ref. [31] arise in a quite concrete way.What is not so clear are singularities in the bubble with momenta q 2 , q 3 , for example.
Let us denote α 3 (q 2 3 ) 2 + β 3 q 2 3 + γ 3 = det(q i • q j ) 1≤i,j≤3 , where α 3 , β 3 and γ 3 are complicated polynomials which we will not need to spell out.We define c 3 = β 3 α 3 , d 3 = γ 3 α 3 and x 3 = q 2 3 , y 2 3 = x 2 3 + c 3 x 3 + d 3 .This curve can be parametrized by where ∆ 3 = c 2 3 − 4d 3 .This implies that (5.24) The boundaries of the integral in x 3 are for values of t 3 where y 3 vanishes.This means that which follows from and similarly for the denominator.This integral therefore will produce a constant transcendental factor of π.The final answer is This integral is positive if the quantity under the square root in the integrand is positive for q 2 3 ∈ (a 3 , b 3 ) and −m 2 3 ̸ ∈ (a 3 , b 3 ) (this is certainly the case if m 2 3 > 0 since a 3 , b 3 ≥ 0).Then the quantity under the square root in the answer is positive.Now we want to do the second integral, (5.28) The quantity under the square root is minus the Gram determinant of q 1 , q 2 , q 3 , evaluated at the "Euclidean on-shell" condition q 2 3 = −m 2 3 .We now have −α where The first term can be rewritten as (5.31) The second term, up to a prefactor, is the Gram determinant of p 1 and p 2 .Indeed, This Gram determinant, representing a volume in Euclidean space, is positive hence the second term in the factorization of ∆ 2 is negative.Therefore, ∆ 2 < 0 and the roots of α 2 (q 2 2 ) 2 + β 2 q 2 2 + γ 2 = 0 are complex.In particular, α 2 (q 2 2 ) 2 + β 2 q 2 2 + γ 2 > 0 in the region of integration.
If P is a quadratic polynomial, the integral can be rewritten as (see eq. (A.20) with p(x) = P (x)/lc(P ) and lc(P ) is the leading coefficient of the polynomial P .p q 1 q 2 q 3 p Figure 6.The kinematics of the two-loop sunrise integral.

The two-loop sunrise integral
This integral is famously elliptic, see refs.[14,27,32].We will see the elliptic curve and the holomorphic differential appear explicitly.We will do the integrals in Euclidean signature.We start with the integral (see fig. 6) Using Cutkosky's change of variables we have where p is the external momentum and p = q 1 + q 2 + q 3 , where q 1 , q 2 and q 3 are the momenta through the internal lines of the sunrise diagram.The values (q 2 3 ) min and (q 2 3 ) max are obtained by computing the stationary points of q 2 3 subject to the constraints that p = q 1 + q 2 + q 3 and q 2 1 and q 2 2 take fixed values.The contour γ is a real cycle on an elliptic curve which we will describe in more detail in the following.
Indeed, it is clear that q 2 1 and q 2 2 can take any values, but the values of q 2 3 are constrained.For clarity, let us denote the fixed values for q 2 1 and q 2 2 by q 2 1 = z 1 and q 2 2 = z 2 .Then, we look for the stationary points of q 2 3 = (p − q 1 − q 2 ) 2 subject to the constraints q 2 1 = z 1 and q 2 2 = z 2 where z 1 , z 2 are some fixed values.We define a function where λ 1 and λ 2 are Lagrange multipliers.The stationary point conditions read while the derivatives with respect to the Lagrange multipliers reproduce the constraints.
Taking the derivatives we find p − q 1 − q 2 + λ 1 q 1 = 0, ( p − q 1 − q 2 + λ 2 q 2 = 0. ( In particular, this implies that λ 1 q 1 = λ 2 q 2 .By squaring we have Let us first assume that λ 1 , λ 2 are non-vanishing.The system of equations above has a solution By squaring the first equation we find Using with (± 1 )(± 2 )(± 3 ) = 1.Then, we have In conclusion, we find four different critical points for q 2 3 .Let us now study their nature in more detail.In order to decide the nature of the critical points, we compute the bordered Hessian matrix In our example we have We have four variables (the components of q 1 and q 2 ) and two constraints.Therefore, we need to look at the signs of the 5 × 5 and 6 × 6 principal minors.Computing these minors we find The conditions for a minimum are where 2 is the number of constraints.The conditions for a maximum are Let us consider the conditions for the minimum: There are four cases 1. + 1 + 2 + 3 : in this case the minimum conditions become ∥p∥ + √ z 2 < 0 and ∥p∥ + √ z 1 + √ z 2 > 0. The first condition never holds since ∥p∥ ≥ 0 and Let us now assume that λ 1 = 0 or λ 2 = 0. Since λ 1 q 1 = λ 2 q 2 , then if λ 1 = 0, either λ 2 = 0 or q 2 = 0. We will not study the region q 2 = 0 anymore since in this case it does not yield a contribution to the integral.This conclusion should be re-evaluated when studying the case m 2 = 0 or in general when the integrals are divergent.We conclude that when λ 1 = λ 2 = 0 the two analogs of the Landau loop equations become a single equation p = q 1 + q 2 .This is consistent with the constraints q 2 1 = z 1 and q 2 2 = z 2 if the triangle inequalities are satisfied z 1 + z 2 > ∥p∥, ∥p∥ + z 1 > z 2 and ∥p∥ + z 2 > z 1 .
In conclusion, we have For the maximum, we have the following four cases: This is actually pretty obvious geometrically.The innermost integral contains a one-form We can compute this ratio as follows.We have q 1 +q 2 +q 3 = p and we take p to be constant (or dp = 0).Then we have where we have used dq 3 = −dq 1 − dq 2 .Then, taking v = p and using where x 2 = (p − q 1 ) 2 = (q 2 + q 3 ) 2 and y 2 = (p − q 3 ) 2 = (q 1 + q 2 ) 2 .The variables x 2 and y 2 are the lengths squared of the diagonals of the quadrilateral whose sides are the vectors q 1 , q 2 , q 3 and −p.
Since the quadrilateral is in a plane, we have that the volume of the simplex it generates vanishes.In other words, the following Gram determinant should vanish (see also ref. [33, lemma 4.1]) where with 3 ), (6.29) 2 ), (6.30) d 00 = (p 2 − q 2 1 + q 2 2 − q 2 3 )(p 2 q 2 2 − q 2 1 q 2 3 ).(6.32) The polynomial P can be made homogeneous of degree three so it describes an elliptic curve embedded in P 2 with homogeneous coordinates (u : v : w).However, an alternative compactification, described in sec.7 is more natural.
Using these results we have that .33)This one-form can be obtained from the two-form udvdw−vdudw+wdudv P by taking a residue at P = 0.
Similarly, we have where As before, taking into account the triangle inequality, the domain for v = y 2 is max (∥p∥ − ∥q 3 ∥) 2 , (∥q We know that (∥p∥ − ∥q 3 ∥) 2 < (∥p∥ + ∥q 3 ∥) 2 , (6.41) so the possible orderings of these four roots are obtained by shuffling the two sets of ordered roots.In total there are six possibilities.The interplay between these conditions and the boundary conditions for integrating over q 2 3 yield a large number of regions.
To make progress, we follow a bit of a different route that Cutkosky's.Instead of integrating over q 2 1 and q 2 2 first, we integrate over q 2 1 and the diagonal v = y 2 (see fig. 8).The triangle inequalities now imply that q 2  2 ∈ [(y − ∥q 1 ∥) 2 , (y + ∥q 1 ∥) 2 ], (6.43) Then the integral becomes Next, a short calculation reveals that Then, the integrals over q 2 2 and q 2 3 can be done as in eq.(4.14) and have the effect of introducing a factor of π each and replacing q 2 2 → −m 2 2 and q 2 3 → −m 2 3 .After doing these integrals, we find ) .
(6.47)Note that the quantity under the square root is always positive in the integration domain.
Note also that when the "angular" integral was the innermost integral it was a complete elliptic integral, once we pulled it through the q 2 2 and q 2 3 integrals it became an incomplete elliptic integral (see sec.B for a discussion of such integrals).
This result can be rewritten in several ways, but the number of integrals can only be reduced at the cost of introducing transcendental functions, such as logarithms.For example, we can also do the integral over q 2 1 which will produce a logarithm and will replace ∥q 1 ∥ → im 1 in the quartic under the square root.Similar results can be more quickly obtained by using Feynman parametrization.The quartic in v can be symmetrically reduced as in eq.(B.9).
At the Euclidean pseudo-threshold ∥p∥ = i(m 1 + m 2 − m 3 ) the integration can be done explicitly in terms of dilogarithms as shown in ref. [14].
Hyperelliptic integrals can also occur, see ref. [34].In that case a similar analysis applies and the angular integrals should yield a distinguished holomorphic form on the hyperelliptic curve along with a distinguished cycle.
The parametrization of the integral in terms of momenta squared may open up new possibilities for regularization.The "angular" integral, being along a compact real cycle and not meeting any singularities does not itself produce divergences, however this inner integral is the only one affected by dimensional regularization (see the discussion in sec.8).Divergences arise from integrals along the variables q 2 e .It is therefore more rational to regularize them instead, since they are producing divergences.One idea that immediately comes to mind are hard cut-offs (IR and/or UV) in q 2 e in Euclidean signature.This will undoubtedly make the integrals more complicated, but possibly not much more than dimensional regularization.We should point out that this type of regularization is better than the usual textbook cutoff regularization, which depends on the choice of loop momentum.
However, the mathematical literature has other types of regularizations which have been applied to polylogarithms and multiple zeta values (see ref. [35] for a longer discussion).Such regularizations, such as tangential basepoint regularization have already been used in refs.[7,8].The regularizations used in the mathematical literature have been designed to preserve various identities satisfied by quantities which did not require regularization.Similarly, in physics, we want to preserve various properties satisfied by physical quantities, which are often broken by traditional regularization choices.

Configurations of quadrilaterals as elliptic curves
It will prove convenient to compactify and complexify the integration domain.This is actually necessary for applying mathematical theorems such as those of refs.[36][37][38][39] and also refs.[40,41] for reviews.The compactification is essential if we want to study second type or mixed second type singularities.For the variables q 2 e an obvious choice of compactification is RP 1 with complexification CP 1 (see ref. [33] where the same compactification is used).
It may happen that, after complexification, the "angular" variables in Cutkosky's terminology (see sec. 2) do not parametrize a compact space.Sometimes, as in the case of the bubble integral in three dimensions (see sec.4.3), a compactification can be performed at the cost of introducing a pole for the innermost differential form.We should note that this rather natural (partial) compactification does not seem to have been discussed in the physics literature before.The type of compactifications that have been considered, see ref. [40, p. 107] involve representing the complexified compactified Minkowski space as a quadric in P 5 , a compactification familiar to twistor theorists.Here, instead, we are proposing to use a compactification to a product of P 1 , times an ad hoc compactification for the "angular" variety.Curiously, an embedding in a product of P 1 was discussed in a different context in ref. [42] but there the interpretation of the coordinate on P 1 was different from here.
In the compactification the integration path q 2 e ∈ (−∞, ∞) is closed since we have a single point at infinity.Indeed, on the complex projective line CP 1 or the Riemann sphere we can choose coordinates so that the origin is at the North pole and the infinity is at the South pole.Then the integration contour (−∞, ∞) is along a meridian.
The interpretation of quadrilateral configurations as points of elliptic curve was described in ref. [33].As explained in this reference, there are two moduli spaces of quadrilaterals; oriented and unoriented.For our purposes the unoriented moduli space is relevant.The oriented moduli space of quadrilaterals makes sense over the real numbers only.
In order to obtain an elliptic curve, we need to compactify the space, which involves taking the lengths of the sides of the quadrilateral to be valued in P 1 .p q 1 q 2 q 3 Figure 7.A quadrilateral configuration, which corresponds to a kinematic point in the sunrise integral.
p q 1 q 2 q 3 y x Figure 8.A quadrilateral configuration, together with diagonals, which has the same lengths of the sides, and same length for one of the diagonals as in fig.7, while the other diagonal has a different length.
The dual space parametrization used here goes back the initial studies of Landau singularities in refs.[4,[43][44][45].This representation is really useful since some non-obvious algebraic properties have simple geometric causes.An equation for the elliptic curve can be obtained by setting to zero the volume of a (degenerate) tetrahedron in fig.8.This yields an equation in u = x 2 and v = y 2 of bi-degree (2, 2) and is, coincidentally, also naturally embedded in P 1 × P 1 .

A peek at integrals in dimensional regularization
In dimensional regularization the Cutkosky representation is usually called Baikov representation (see ref. [46] for the original paper and ref. [47] for an introduction).This has a loop-by-loop version worked out in ref. [48].
Let us now do a sample computation in dimensional regularization.We will first attempt a simple integral, the massless bubble in dimension d = 4 − 2ϵ.Our computation will not use the Wick rotation, which is less natural for massless particles since it clashes with the on-shell conditions.The computation is more complicated than the textbook computation using Feynman parameters and Wick rotation.
We have As before, we have a 2 = (∥p∥ − ∥q 1 ∥) 2 and b 2 = (∥p∥ + ∥q 1 ∥) 2 .Recall that ∥p∥ = p 2 and if p 2 < 0, then ∥p∥ ∈ iR.In this case, when the boundaries of integration are not real, value z 1 = −1 provides a natural boundary between these two regions so one can canonically separate the integral into an IR and a UV region, to be studied separately.
Next, we will seek to establish that this hypergeometric integral has the following behavior  This was an woeful derivation, which we hope to improve later.What is needed is a way to systematically expand around z 1 = 0 in the first region and z 1 → −∞ in the second region.
Finally let us briefly discuss the massless box integral.We have (8.24) Let us put this integral in Cutkosky form.The range of k 2 is R.The range of (k + p 2 ) 2 can be determined by extremizing (k + p 2 ) 2 at fixed p 2 and subject to the constraint that k 2 is fixed at k 2 = z 1 .Using Lagrange multipliers we find Since p 2 2 = 0 and z 1 = k 2 ̸ = 0 we have that an extremum is never realized.Hence, the range of (k + p 2 ) 2 is also R. Note that this is not what happens in Cutkosky's approach of proving his theorem.It remains to be seen if and how his proof would have to be modified to cover this case.
The range of (k +p 23 ) 2 is determined in a similar way.Using the same idea of Lagrange multipliers we find the equation The compatibility condition is the existence of a tetrahedron with sides k, k+p 2 , k+p 23 , p 3 , p 14 and p 2 and the extrema arise when the tetrahedron becomes degenerate which is when its volume vanishes.Using the Cayley-Manger formula for the volume of the tetrahedron one can solve for the stationary point of z 3 = (k + p 23 ) 2 and we find a unique solution . It is remarkable that there is a unique solution which is another departure from the case analyzed by Cutkosky.We defer a more detailed study of the nature of this stationary point.
The extrema of (k − p 1 ) 2 can be analyzed similarly and in that case one finds two solutions, as usual.
The remaining integrals can be written as
(b − x)(x − a) to have a branch cut along the segment [a, b].We pick b > a and c ̸ ∈ [a, b].With the usual definition of the square root for complex numbers we have that (b − x)(x − a) has a branch cut (−∞, a] and another branch cut [b, ∞).Since we want to have a branch cut along the segment [a, b], we split the square root as (x − a)(x − b) → √ x − a √ x − b and we use the definition √ z = √ ρe i θ 2 where z = ρe iθ with θ ∈ [0, 2π) for the first square root and the definition √ z = √ ρe i θ 2 where z = ρe iθ with θ ∈ (−π, π] for the second square root.If we define z − a = ρ 1 e iθ 1 with θ 1 ∈ [0, 2π) and z − b = ρ 2 e iθ 2 with θ 2 ∈ (−π, π], then we have that the of the square root above the cut is i √ ρ 1 ρ 2 .Since for x ∈ [a, b] we have b − x = ρ 2 and x − a = ρ 1 have that if defined such that the branch cut is along [a, b], the function √ x − b √ x − a (where the two square roots are defined as above) is equal by continuity from above the cut to i (b − x)(x − a).

Figure 3 .
Figure 3.A contour around a square root branch cut.