Landau Singularities Revisited: Computational Algebraic Geometry for Feynman Integrals

We reformulate the analysis of singularities of Feynman integrals in a way that can be practically applied to perturbative computations in the Standard Model in dimensional regularization. After highlighting issues in the textbook treatment of Landau singularities, we develop an algorithm for classifying and computing them using techniques from computational algebraic geometry. We introduce an algebraic variety called the principal Landau determinant, which captures the singularities even in the presence of massless particles or UV/IR divergences. We illustrate this for 114 example diagrams, including a cutting-edge 2-loop 5-point non-planar QCD process with multiple mass scales. The algorithms introduced in this work are implemented in the open-source Julia package PLD.jl available at https://mathrepo.mis.mpg.de/PLD/.


INTRODUCTION
Modern high-precision computations in particle physics rely on the knowledge of the singularity structure of perturbative scattering amplitudes.This emphasizes the importance of being able to predict it without the explicit evaluation of Feynman integrals.A pioneering step in this direction was taken by Bjorken, Landau, and Nakanishi [1][2][3], who demonstrated that poles and branch points of finite Feynman integrals can be determined by pinches of on-shell hypersurfaces in the loop-momentum space.They are known as Landau singularities; see [4,Sec. 1] for a historical overview and more complete references.Nevertheless, it is not known how to systematically apply this analysis to Standard Model processes beyond one-loop level, due to many problems outlined below.The goal of this letter is to fix these issues and give a practical algorithm for determining singularities (physical or not) in scenarios with massless particles and UV/IR divergences.
Landau analysis is particularly important in the context of finding analytic expressions, differential equations, and symbol calculus for Feynman integrals, see, e.g., [5].As a simple example, if the amplitude has a term of the form 1  P log[(Q + √ R)/(Q − √ R)] for some polynomials P, Q, R in the external kinematics, then its Landau singularities are located at P = 0 (simple pole), R = 0 (square root), and Q 2 − R = 0 (logarithmic) on some Riemann sheets.By definition, Landau singularities coincide with zeros and singularities of symbol letters in polylogarithmic cases [6,7].By basic elimination theory, their positions can be always expressed as ∆ = 0, where ∆ is a polynomial in the Mandelstam invariants and masses-squared with coefficients in Q.The challenge is to determine all such ∆'s without having to compute Feynman integrals.
According to textbooks, positions of singularities of Feynman integrals can be computed by solving Landau equations [8][9][10][11][12].The leading singularity is computed by setting all propagators on-shell and imposing that all interaction vertices are located at space-time points.Subleading singularities are leading singularities of reduced diagrams obtained by contracting any subset of edges.In addition, second-type and mixed singularities are those corresponding to one or more loop momenta diverging.All of these are often summarized as solving the equations where q µ i , m i , and α i are respectively the momentum, mass, and the Schwinger parameter of the i-th propagator, and the sum goes over all edges belonging to the a-th loop with signs denoting orientations.
Our work starts with an observation that the above classification is not complete.Indeed, there are known examples where a naive application of (1) does not detect all singularities and more careful blow-ups are needed [13][14][15].There are two fundamental problems with the above textbook lore: (i) it is not clear why only α i ̸ = 0 or α i = 0 solutions are considered, in contrast with more complicated scaling patterns; and (ii) the prescription is ambiguous in the presence of UV/IR divergences, i.e., when (1) has solutions for any external kinematics.In this work, building on [16], we introduce tools from nonlinear algebra [17], which tackle these issues as follows.
(i) Beyond the standard classification.We find that, in general, a singularity can come from scalings Each region is given by the set of rational exponents (w 1 , w 2 , . ..) called weights.The internal momenta q µ i will also have a specific scaling as a consequence of (2), which will be given later in (10).Their combinatorics is captured by a lattice polytope P described below.
Hence, it is not enough to study reduced diagrams, but instead a whole suite of configurations obtained by shrinking and expanding edges at specific relative rates is needed.Below, we will explain that scalings beyond (2) are often possible as well.This fact may not be surprising to the readers familiar with the method of regions [18], which also exploits similar scalings, though only for α i ⩾ 0, while we consider any α i ∈ C.
(ii) Interplay with UV/IR divergences.The second issue is best illustrated on an example of a reduced diagram obtained by contracting a number of edges, say in a given QCD process with all external legs on-shell: In this case, the diagram on the right-hand side is always singular due to a collinear divergence between gluons (curly) and possibly quarks (solid) if they are taken to be massless.Hence, the Landau equations in the form (1) have solutions regardless of the remaining external kinematics, i.e., they evaluate to 0 = 0.However, it does not mean we can simply discard them, because the same region of phase-space will in general lead to other genuine kinematic singularities.Moreover, it is clear one could not find them solely by considering the reduced diagram in (3) which does not depend on any Mandelstam invariants.Instead, the limit (3) needs to be taken more carefully.
More generally, we explain how to disentangle the divergences that are always present (UV/IR) from singularities that exist only for special values of the kinematics (Landau) using a combination of elimination theory and numerical irreducible decomposition.
The goal of this letter is to summarize the main construction in a way accessible to physicists and introduce the software.The follow-up paper [4] will provide mathematical foundations and algorithmic details.

LANDAU EQUATIONS REVISITED
A family of Feynman integrals with n external legs in D dimensions can be parametrized by where the P i 's correspond to m denominators and irreducible scalar products, each raised to a (possibly noninteger) power ν i , and ℓ a 's are the L loop momenta.
We collectively call all the kinematic parameters (Mandelstam invariants and masses-squared) z ∈ E, where E is the specific subspace of the kinematics we consider.The overall normalization does not affect the singularity structure.Introducing m Schwinger parameters α = (α 1 , α 2 , . . ., α m ), we write m i=1 where L and c can depend on the internal masses and external kinematics, but are independent of the loop momenta.The matrix Q and vector L give rise to the Symanzik polynomials: Hence, F depends on external kinematic invariants and masses, while U does not.Following standard steps, the integrals in ( 4) are then proportional to [19] where is often called the graph polynomial.We fix the auxiliary mass scale µ = 1 from now on.Singularities can also be analyzed in other representations, which leads to equivalent equations [4,Sec. 3.4].One can show that propagators with ν a ⩽ 0 do not give rise to new singularities and can be excluded from the analysis [4, App.B].
The simplest Landau singularity (known as leading second-type) is obtained by solving the critical point equations away from the boundaries: where It is the condition for the denominator in (7) to be singular in a way that cannot be cured by contour deformations.The Feynman integral I(z) possibly develops a singularity whenever a solution to (9) exists for a given value of z.Back in the loop-momentum space, completing the square in ( 5) is equivalent to solving the loop Landau equations (1b).Hence, one can read off the behavior of the loop momenta in terms of Schwinger parameters for every solution.Similarly, solving (9) imposes the on-shell conditions (1a).Most of the singularities are hiding on the boundaries.They can be classified using the Newton polytope: The same polytope appears in other recent tropical approaches to Feynman integrals [20][21][22][23][24].Its faces can be parametrized by weights w = (w 1 , w 2 , . . ., w m ) and are in one-to-one correspondence with scalings α i → ε wi α i with ε → 0. In particular, the behavior of the loop momenta ℓ a in terms of ε can be read off by plugging in these scalings into (10).On each face f of P, we keep only the terms of G leading in ε.The result is called the initial form and denoted by G f .The Landau equations on this face are The system ( 9) is recovered in the special case when f is the full-dimensional face, i.e., the entire polytope.Mathematically, the faces of P correspond to torus orbits in a toric compactification of (C * ) m .The incidence variety Y f (also called "pinch surface") is the variety defined by (12).Some of the faces can be matched onto the textbook classification, though others are new.For example, the facet with weight w = (−1, −1, . . ., −1) computes the leading first-type singularity.Replacing −1's with 0's and 1's corresponds to various subleading singularities, see [4,Sec. 3.5].All other faces go beyond the standard classification (note that they also appear in the geometric approach to the expansion by regions [25]).By definition, solutions with all α i ⩾ 0 are on the physical sheet of the kinematic space [26]; we are interested in singularities on any sheet where α i 's are generically complex.We can further categorize the solutions depending on which value of the external kinematics they occur as follows.
(a) Any value of the external kinematics, z ∈ E. These are the UV/IR divergences, which below are called dominant components (also known as "permanent pinches").After identifying from which part of the integrand they arise, we want to strip them away since they are instead taken care of by dimensional regularization.
(b) Special values of the external kinematics, z ∈ {∆ = 0} ⊂ E. These are the kinematic singularities we want to identify.Since analytic functions can only have codimension-1 singularities in E, i.e., the ones that can be written as ∆ = 0, we focus on them here.
A given face can have none, either, or both types of solutions which need to be distinguished.We now move on to studying the geometry of this problem.

GEOMETRY OF SINGULARITIES
Euler discriminants.Intuitively, Feynman integrals (7) may develop singularities whenever the surface {G(α; z) = 0} becomes more singular than for a generic z ∈ E. As a way of probing the topology of this surface, we compute the (signed) Euler characteristic: While non-trivial, one can prove that there exists a generic value χ * , such that χ z = χ * for almost all z ∈ E and χ z < χ * otherwise [4,Thm. 3.1].This leads us to define the Euler discriminant variety ∇ χ (E) as the locus of such special kinematic points: If the coefficients of all monomials in G were independent of each other, ∇ χ (E) coincides with the singularity locus of generic Gelfand-Kapranov-Zelevinsky (GKZ) hypergeometric functions, which is known as the principal Adeterminant [27].We propose that this relationship extends to Feynman integrals: that their singularities are contained in ∇ χ (E).
This claim is physically intuitive: χ z computes the number of master integrals in the appropriate sense (see [28,29] for reviews).The above proposal equivalently says that the system of first-order differential equations satisfied by ( 4) drops rank when evaluated on a singularity.
There are two obstacles with applying the above tools to Feynman integrals.Firstly, beyond the simplest of cases, they are not generic in the GKZ sense: for example, the G polynomial of the diagram in (3) has 64 monomials, but depends only on 7 kinematic invariants.Hence, one cannot directly apply GKZ technology such as principal A-determinants to Feynman integrals (see [30,31] for previous attempts).Explicit criteria are given in [4, Sec.2.2].Secondly, while the Euler discriminant can be used to determine whether a given candidate point z is singular, it does not give a constructive algorithm for finding ∇ χ (E).Both of these problems lead us to the definition of the principal Landau determinant (PLD) [44].
Principal Landau determinants.The PLD aims to systematize Landau analysis based on the combinatorics of the polytope P introduced above.A complication is that each face f can produce multiple components and one has to distinguish if they belong to (a) or (b).
The incidence variety Y f defined by ( 12) is a union of subvarieties of (C * ) m ×E.In general, it can have multiple irreducible components with different dimensions: Note that the dimension of is in general different from the dimension of f .We are interested in eliminating α's, which is the same as projecting all components of Y f "downstairs" to the kinematic space E. In other words, for every point (α; ẑ) contained in Y f , its projection π(α; ẑ) = ẑ "forgets" about α and remembers only the coordinates ẑ ∈ E. We illustrate this with a simple example: Here, α and z = (z 1 , z 2 ) are coordinates on the onedimensional C * and two-dimensional E respectively.The incidence variety Y f has five irreducible components, Y f , Y f , . . ., Y f .Let us analyze them one by one.Below, when we write fiber dimension of a component, we mean the dimension of the pre-image π −1 (z) for a generic point z on π(Y (i) f ), i.e., it is the number of unconstrained Schwinger parameters remaining after solving (12).
The orange component Y (3) f is dominant because it projects down to the whole kinematic space E. It has fiber dimension 0, i.e., every point downstairs comes from a projection of a unique point upstairs.In physics terms, it corresponds to a UV and/or IR divergence.In other words, the Landau equations ( 12) have a solution for any z ∈ E. This, however, does not mean that we can simply disregard the (fictional) face f : there are other irreducible components of type (b).
The red components Y has fiber dimension 1 and projects to codimension 1.We hence keep all three of them in PLD.Finally, the green component Y (2) f has fiber dimension 1, but projects to codimension 2. We do not include it in PLD.See [4,Ex. 5.1] for more details.
Let us call the result of the projection π(Y For any Feynman diagram G, PLD is defined as the union of all those projections that have codimension 1: This definition can be used in practice and is implemented in PLD.jl (we refer to [4,Sec. 3.2] for a formulation in the language of algebraic geometry).

IMPLEMENTATION
Code and database.Classic symbolic elimination tools, such as Gröbner bases, can be used, but are not ef-ficient enough to handle multi-loop examples.We hence introduce a numerical algorithm for computing PLD, based on homotopy continuation techniques.It performs the above irreducible decomposition and projection to the kinematic space.We refer to [4,Sec. 5.2] for the algorithmic details.We implemented both symbolic and numerical elimination algorithms in an open-source Julia package PLD.jl available at [32].This website also contains documentation and a guided tutorial through the functionality of PLD.jl.Together with the package, we provide a database of 114 example diagrams with various graph topologies and mass assignments.
Comparison with HyperInt.The only competitive tool for finding Landau singularities is the compatibility graph reduction algorithm cgReduction implemented in HyperInt [33] and based on the formalism of Pham and Brown in integer dimensions D [34,35].It computes an upper bound on the set of singularities, i.e., finds many components that are not genuine singularities of the corresponding Feynman integral.However, they can be efficiently filtered out by the Euler characteristic criterion (14), see [4, App.A] for a practical example.We found diagrams for which this pipeline gave more singularities than those computed numerically by PLD, which means that PLD G (E) ⊊ ∇ χ (E) ⊊ cgReduction.Indeed, it is known that the blow-ups implicitly performed by looking at the Newton polytope might not suffice to detect all singularities of Feynman integrals, see [15,Sec. 6.4] and [4,Ex. 3.9].In other words, there are scalings of Schwinger parameters that go beyond (2) with genuinely new singularities.These are captured by the Euler discriminant but not by PLD.
On the other hand, in practice, cgReduction terminated only on the simplest 70 out of 114 diagrams we considered, so PLD remains the only practical tool for obtaining (a subset of) singularities of more complicated diagrams.At any rate, most of the entries in the database give new predictions for singularities of previously unstudied diagrams.

EXAMPLE
As a concrete example, let us consider the following Higgs + jet production process: We take the external Higgs (dashed) to have masssquared M 2 = p 2 5 , the quarks (solid) to have masssquared m 2 , and all the remaining particles massless.
The m = 0 case has been computed in [36], but m ̸ = 0 remains unknown.Here, we use PLD.jl to predict its singularity structure.
In this case, the kinematic space E is parametrized by z = (s 12 , s 23 , s 34 , s 45 , s 51 , m 2 , M 2 ) , where s ij = (p i + p j ) 2 are the Mandelstam invariants.The diagram is specified by providing the internal edges as a list of pairs of vertices they connect, as in (18), and similarly nodes listing the vertices to which the external momenta p 1 , p 2 , . . ., p 5 are attached.Likewise, the internal and external masses are also assigned in order of appearance.One computes the PLD by running: [2,3], [3,6], [4,6], 1 [5,6], [5,7], [1,7], [4,7]] getPLD(edges, nodes, internal_masses = [0,0,0,m2,m2,m2,0,m2], 6 external_masses = [0,0,0,0,M2]) The code scans through all faces in increasing order of dimensions (roughly speaking, easier to harder).The generic Euler characteristic is χ * = 330 and the number of faces of P in dimensions 0, 1, . . ., 8 is 56, 294, 681, 884, 699, 343, 101, for a total of 3075 separate systems of equations to consider.Note that a naive analysis based on reduced diagrams would include only 2 8 = 256 systems, many of which ambiguous as in (3).We solved the simplest 1270 ones symbolically and the remaining 1805 numerically [45], while cgReduction did not terminate.PLD.jl found the total of 71 distinct kinematic singularities {∆ i = 0} with degrees between 1 and 12 in z.In general, multiple different faces f contribute to the same singularity.The full list is provided in the G = Hj-npl-pentb entry in our online database.It is a new prediction for the singularities of ( 18), though we do not claim that the list is exhaustive.Here, we give a few examples of ∆ i for illustration.
We can make a comparison with the m = 0 case.As emphasized throughout [4], specializing kinematic parameters and computing PLD in general do not commute, but one can obtain a subset of the m = 0 singularities by substituting m = 0 in each ∆ i .E.g., (23)  (24) Indeed, each factor in this polynomial was already identified as a singularity of the m = 0 specialization [36] and similarly for ∆ 47 | m=0 .This can also be verified with the EulerDiscriminantQ function in PLD.jl.The fact that ( 23) is so much larger than (24) suggests that the m ̸ = 0 Feynman integral may have a vastly more complicated analytic form.

DISCUSSION
In this letter, we reformulated Landau analysis in a way that can be applied to multi-loop computations in the Standard Model and made it practical by introducing PLD.jl.It was important to realize that singularities can arise from more complicated patterns of shrinking/expanding edges than just reduced diagrams.Previous results relying on this assumption should be revisited.Likewise, our work emphasizes the need to study asymptotics going beyond (2) more carefully.
This progress should be viewed as a necessary step in a more ambitious program of using Landau singularities to systematically constrain Feynman integrals relevant to collider and gravitational-wave physics.For example, expanding around incidence varieties would give information about the local nature of the singularity and their discontinuities in the kinematic space, see, e.g., [26,42] for progress on 0-dimensional fibers without UV/IR subdivergences.We leave these important questions for future work.
Data Statement.The supporting data for this article are openly available from the research data repository MathRepo [32].