An algorithm for the symbolic reduction of multi-loop Feynman integrals via generating functions

Multi-loop Feynman integrals with many external legs and several kinematic scales remain a central bottleneck in state-of-the-art perturbative calculations. Regarding master-integral evaluation, the last decade has witnessed striking progress driven by differential-equation methods, particularly in canonical form, alongside increasingly automated basis transformations and complementary numerical strategies. These developments now support genuinely high-multiplicity systems: automated transformations to canonical bases are available, auxiliary-mass-flow methods provide an efficient route to numerical evaluation, the complete reduction of two-loop five-light-parton amplitudes has been achieved, all two-loop five-point one-mass and six-point planar massless master integrals are known analytically, and analytic results have recently reached the three-loop five-point level [1−10]. However, for frontier calculations, these advances can only be fully exploited if the accompanying reduction problem can be handled with comparable efficiency and transparency.

To obtain physical amplitudes, one must first reduce the large space of integrals generated by Feynman rules to a finite basis of master integrals. The modern framework for this step originated with the integration-by-parts (IBP) method introduced by Tkachov and Chetyrkin, while Laporta subsequently transformed IBP identities into a systematic elimination algorithm suitable for computer implementation [11−13]. Around this core framework, several complementary approaches have been developed, including dimensional shifts, generalized recurrence relations, and difference-equation methods [13−15]. In this sense, the history of reduction is not merely the history of solving larger linear systems; it is also the history of learning how to encode and organize relations among integrals in a form adapted to high-loop, multi-leg, and multi-scale problems. Even today, for sufficiently complex families, reduction is often the dominant algebraic bottleneck in the entire calculation [16, 17].

The practical success of IBP reduction is inseparable from the development of public software packages. Early automated implementations such as $\mathrm{AIR}$, $\mathrm{FIRE}$, and $\mathrm{Reduze}$ established the viability of large-scale Laporta-style reductions, while $\mathrm{LiteRed}$ emphasized the heuristic search for symbolic sector reduction rules alongside symmetry detection, differential equations, and dimensional recurrences [18−22]. A modified algorithm from $\mathrm{LiteRed}$ appears in the recent program $\mathrm{SPIDER}$ [23]. Later generations of software, including $\mathrm{Reduze}$ 2, $\mathrm{FIRE}$, and $\mathrm{Kira}$, advanced the field through distributed reduction, modular arithmetic, finite-field methods, and improved equation management [24−28]. More recently, $\mathrm{Blade}$ and $\mathrm{NeatIBP}$ have pushed toward block-triangular decompositions and trimmed IBP systems [29−31]. A concise recent snapshot of this software-centered landscape from the viewpoint of symbolic reduction is provided in Ref. [32], while broader historical overviews can be found in Refs. [16, 17]. At the same time, these developments highlight a fundamental limitation of the traditional paradigm: in demanding multi-loop, multi-leg problems, one still often confronts extremely large sparse systems whose solution is expensive in both time and memory.

For this reason, recent work has increasingly focused on reformulating the reduction problem itself rather than merely accelerating Laporta elimination. One line of development constructs shorter or better-adapted IBP systems using generalized generators, unitarity cuts, algebraic geometry, module intersections, small-size IBP systems, and transverse integration identities [30, 33−36]. A second line seeks symbolic rules that remain valid on the lattice of indices, leading to syzygy-constrained symbolic reduction rules, diagonalized or triangularized IBP equations, seedless lowering operators, and covariant-differentiation-based reductions [37−40]. A third line explores data-driven strategies, where machine-learning, reinforcement-learning, and related heuristic approaches are used to optimize equation generation, equation selection, or even online reduction itself [41−43]. Recently, an optimized IBP reduction method has been proposed based on the fundamental filtration structure of Feynman integrals [44−46]. Alternative non-IBP reduction frameworks based on intersection theory [47−56] are also becoming increasingly practical for multi-leg problems [57]. Taken together, these developments make clear that the modern question is not only how to solve the reduction equations faster, but also how to represent them in a form where their hidden structure becomes manifest.

Among these new developments, the generating function method possesses unique characteristics at both conceptual and algorithmic levels. Instead of treating each lattice point of indices separately, one packages an entire sector into a single object, such that shifts of propagator powers or numerator degrees become differential operators, and IBP identities become differential equations for the generating function. Symbolic recurrence relations then arise naturally at the operator level. This approach makes the dependence on indices explicit, suggests completeness tests that do not rely on manually chosen seed sets, and opens the possibility of globally valid lowering relations. From this viewpoint, the generating-function program sits naturally at the intersection of the older recurrence-relation tradition and the modern search for seed-independent symbolic reductions [13−15].

The recent history of this approach is already rich enough to be instructive. The early use of generating functions appeared in [58, 59], where specific reduction tasks were carried out. Generating functions for one-loop tensor reduction first appeared in [60], after which explicit closed-form constructions were developed and the method was extended to arbitrary one-loop reduction [61, 62] and, more recently, to the multi-loop case [63]. In parallel, generating functions were proposed as a means to explore the linear space of Feynman integrals and to derive reduction rules from differential equations satisfied by those generating functions [64]. The next step was to connect this language directly to IBP relations in a systematic one-loop framework [65], then to Baikov-space generating functions for multi-loop reduction coefficients [66], and finally to a general proposal for the complete symbolic reduction of multi-loop integrals via generating functions [32]. Closely related developments, such as syzygy-constrained symbolic rules and seedless constructions of lowering operators, reinforce the same lesson: reduction may be organized most efficiently at the level of operator identities rather than by solving large seeded systems case by case [37, 39].

Our goal is to develop precisely such a formulation for sector-wise generating functions. The starting point is to package all integrals in a given sector into a single generating function. In this language, IBP identities become differential equations for the generating function, while symbolic reduction rules arise from relations among differential operators. The emphasis therefore shifts from isolated recurrence relations at specific lattice points to operator identities that can be generated, simplified, compared, and organized systematically.

This perspective leads naturally to an iterative algorithm. Starting from the differential equations induced by IBP identities, we derive reduction rules for selected operator structures by solving small linear systems, simplify newly generated equations with the rules already obtained, and test whether the remaining irreducible lattice points match the expected master-integral content. The resulting workflow separates the problem into four tightly linked ingredients: the generating-function formalism, the operator algebra underlying the reduction rules, the iterative procedure that constructs and updates those rules, and a lattice-based completeness criterion.

The paper is organized as follows. In Section II, we introduce the generating-function description of an integral family and the operator-algebra language used throughout the paper. In Section III, we present the iterative algorithm for generating and solving differential equations and for checking completeness. The remaining sections develop a sequence of examples. We begin with the sunset topology, which serves as the main pedagogical case, and then turn to planar and non-planar double-box topologies, representative subsectors, and finally a degenerate topology in which the top sector has no master integrals. Taken together, these examples demonstrate how the same framework handles top sectors, subsectors, and sector elimination within a unified formulation.

In this section, we establish the foundational notation and terminology used throughout the paper. Rather than viewing integral reduction as a set of linear relations among individual Feynman integrals, we encode this information in generating functions that satisfy differential equations. Under this reformulation, the reduction problem can be studied directly in terms of differential operators and their non-commutative algebra.

The central idea is to encapsulate all integrals within a fixed sector into a single generating function. Integration-by-parts identities then become differential equations for the generating function, and symbolic reduction rules emerge from operator identities obtained from these equations. This perspective is advantageous for two reasons. First, it clarifies the structural content of the reduction problem. Second, it provides the algebraic framework required by the iterative algorithm developed in the subsequent section.

An L-loop Feynman integral family with E independent external momenta can generally be represented by

where $ D = 4-2\epsilon $ is the spacetime dimension, $ \ell_i $ are the loop momenta, and $ \vec{\nu} = (\nu_1, \dots, \nu_N) $ is a vector of integer indices. We denote the number of inverse propagators by $ N_{\mathrm{prop}} $; these inverse propagators, denoted $ \mathcal{D}_1,\dots,\mathcal{D}_{N_{\mathrm{prop}}} $, may possess positive or negative powers. The remaining denominators, $ \mathcal{D}_{N_{\mathrm{prop}}+1},\dots,\mathcal{D}_N $, constitute the irreducible scalar products (ISPs) required for a complete basis of scalar denominators, and their corresponding powers $ \nu_{i} $ are non-positive integers. The total number of denominators is $ N = L(L+1)/2 + LE $.

An integral family comprises multiple sectors, each defined by the presence of a specific subset of denominators $ \{\mathcal{D}_i\} $. We encode this structure using a binary sector vector $ \vec{\mu} = (\mu_1, \dots, \mu_N) $, where $ \mu_i=1 $ if $ \mathcal{D}_i $ is a propagator in this sector and $ \mu_i=0 $ otherwise ¹. For each sector $ \vec{\mu} $, we define a generating function $ G_{\vec{\mu}} $ in terms of a set of formal variables $ \vec{\eta} = (\eta_1, \dots, \eta_N) $:

Here, $ s_0 $ is an arbitrary momentum scale introduced to render the variables $ \eta_i $ dimensionless.

Expanding $ G_{\vec{\mu}}(\vec{\eta}) $ in a Taylor series around $ \vec{\eta}=\vec{0} $ reveals the connection to the individual Feynman integrals $ I(\vec{\nu}) $. To illustrate, consider a sector where the first $ N_{\mathrm{prop}} $ denominators are propagators ($ \mu_i=1 $ for $ i=1,\dots,N_{\mathrm{prop}} $), while the remaining terms appear in the exponential ($ \mu_j=0 $ for $ j=N_{\mathrm{prop}}+1,\dots,N $). The expansion yields:

Within a fixed sector, we suppress the sector label and write $ F_{\vec n}\equiv F_{\vec\mu,\vec n} $; the full notation $ F_{\vec\mu,\vec n} $ is retained whenever the sector must be specified.

For later convenience, we summarize our multi-index conventions. For N-dimensional vectors $ \vec a $ and $ \vec b $, we define

Thus, these expressions are products rather than sums. For example, when $ N=2 $, $ \vec\eta^{(a_1,a_2)}\vec\partial^{(b_1,b_2)} =\eta_1^{a_1}\eta_2^{a_2}\partial_1^{b_1}\partial_2^{b_2} $. Let $ \vec e_i $ denote the ith standard basis vector, whose only nonzero component is $ (e_i)_i=1 $. Vector inequalities are understood componentwise:

By comparing the terms, we establish a direct mapping between the expansion coefficients $ F_{\vec{\mu}, \vec{n}} $ and the integrals $ I(\vec{\nu}) $:

This mapping is the foundation of our formalism. It identifies Feynman integrals with the expansion coefficients of the generating function, which are indexed by lattice points $ \vec{n} $. We define the degree of the coefficients $ F_{\vec{n}} $ (or, equivalently, of the corresponding Feynman integrals I) as $ {\displaystyle \sum}_{i=1}^N n_i $. The reduction problem for Feynman integrals is thus reformulated as finding relations among these coefficients, induced by the differential equations (DEs) satisfied by the generating function.

Although there are potentially $ 2^N $ different generating functions corresponding to different choices of $ \vec{\mu} $, some sectors vanish identically, allowing the corresponding generating functions to be set to zero from the outset. Furthermore, if global symmetries relate different sectors, it is sufficient to compute a single representative generating function and derive the others by symmetry. Consequently, the number of generating functions that must be explicitly computed is significantly smaller than $ 2^N $ in practice.

The well-known IBP identities arise because the integral of a total derivative with respect to the loop momenta vanishes; schematically, $\displaystyle \int \mathrm{d}^D\ell \, \dots \, \dfrac{\partial}{\partial \ell^\mu} (\dots) = 0 $. When applied to the integrand in Eq. (2), these identities yield a system of linear partial differential equations for $ G_{\vec{\mu}}(\vec{\eta}) $. The building blocks of these equations are differential operators—specifically, polynomials in the variables $ \eta_i $ and the partial derivatives $ \partial_i \equiv \partial/\partial\eta_i $. Using the multi-index products defined in Eq. (4), any such operator can be expressed as a linear combination of operators in standard form, $ \vec{\eta}^{\vec{a}} \vec{\partial}^{\vec{b}} $. For example, the operator $ \partial_1 \eta_1 \partial_1 $ can be rewritten in standard form as $ \eta_1 \partial_1^2 + \partial_1 $. In this form, an operator is uniquely specified by a pair of vectors $ (\vec{a},\vec{b}) $, which we call the finer index of the operator.

The general differential equations (DEs) for the generating function of a given sector $ \vec{\mu} $ encountered in this work take the form

where $ \widehat{O}_t $ are differential operators in standard form, and the coefficients $ c_t $ are polynomials in the spacetime dimension D and the kinematic invariants. The boundary term B receives contributions from the generating functions of subsectors (i.e., sectors with fewer propagators), possibly with additional operator insertions. When performing the reduction for a given sector, we assume that the reductions for all its subsectors are already known, thus treating B as a known source term. We expand this source as

where $ \widetilde B_{\vec n} $ denotes the coefficient of $ \vec\eta^{\vec n} $ following the required subsector reductions. When the lattice point is clear from the context, we abbreviate this coefficient as $ \widetilde B $.

Applying the operator $ \widehat{O}_t = \vec{\eta}^{\vec{a}_t} \vec{\partial}^{\vec{b}_t} $ to the series expansion $ G_{\vec{\mu}}(\vec{\eta}) = {\displaystyle \sum}_{\vec{n}} F_{\vec{n}} \vec{\eta}^{\,\vec{n}} $ yields:

where, for simplicity, we adopt the convention ²

Thus, Eq. (7) yields

This translates the differential Eq. (7) into a recurrence relation. From (11), it is evident that the role of an operator in the recurrence relation is largely determined by the combination $ \vec{o} \equiv \vec{b} - \vec{a} $, which we term the index of the operator. Similarly, we define the degree of an operator as $ |\vec{o}|\equiv \displaystyle \sum_{i=1}^N o_i $, which corresponds to the degree of the lattice points via (11). A differential equation is said to have degree M when the highest degree among all operators appearing in that equation is M.

Together, (7) and (11) recast integral reduction as an operator problem. Expressing higher-degree coefficients in (3) in terms of lower-degree ones is equivalent to expressing a higher-degree operator as a linear combination of lower-degree operators using the DEs. We refer to such expressions as symbolic reduction rules, as they generate recurrence relations that reduce the Feynman integrals. The primary objective of this work is to develop a systematic algorithm to determine the complete set of symbolic reduction rules.

Before proceeding, we emphasize that we do not seek to solve the differential Eq. (7) for the generating function—that is, to express the generating function as an analytic function of the external data and η's. Doing so would be exceedingly difficult. As previously noted, our aim is to solve for the highest-degree operators from a set of differential equations, which merely requires solving a system of linear equations. The number of linear equations involved is significantly smaller than that in the IBP system.

Before developing the operator algebra in full generality, we illustrate the main idea using the two-loop sunset family. This example is revisited in detail in Section 4; here, we present only the short chain from a sector generating function to a symbolic reduction rule.

The inverse propagators and auxiliary scalar products are

where K denotes the external momentum. The non-vanishing generating functions required for the top sector and its subsectors are

Using Eq. (2), they are

For later use, we define

One of the six fundamental top-sector IBP equations is

After the six equations are assembled, Gaussian elimination isolates a single highest-degree operator in each row. Schematically, the preceding row yields a rule featuring a single highest-degree operator (classified below as T1A) of the form

where $\widehat{O}_{ <2}$ contains lower-degree operators and B denotes known subsector terms. At the coefficient level, this yields

Consequently, all lattice points satisfying $n_3,n_5\geq1$ are reducible by this rule. A complete sunset analysis derives all six first-round rules, identifies the remaining lattice subsets, and generates descendant equations until the reduction terminates.

The preceding discussion demonstrates that the reduction problem is governed by a non-commutative algebra of differential operators. To exploit this structure algorithmically, we must identify which properties of the algebra remain stable under simplification and which dictate the shifts on the coefficient lattice. We therefore compile the operator properties that will be used repeatedly in the construction of reduction rules.

We first consider a special class of operators, namely, those with finer index $ (\vec{a},\vec{a}) $. Since their index is $ \vec{o}=\vec{0} $, we refer to them as zero-index operators. Hereafter, the term "zero-index operator" may denote either a single such operator or a sum thereof. These operators exhibit several key properties that substantially simplify the algebra.

$ {(A)} $ Commutativity: Any two zero-index operators commute, i.e., $ [\widehat{O}_{0,a}, \widehat{O}_{0,b}] = 0 $. This can be proved by showing that every operator of the form $ \eta^k \partial^k $ can be expressed as a polynomial in the fundamental operator $ \widehat{\eta} \equiv \eta \partial $. Since these operators are all functions of the single operator $ \widehat{\eta} $, they must mutually commute.

$ {(B)} $ Braiding Property: For any operator $ \widehat{O} $ and any zero-index operator $ \widehat{O}_1 $, there exists a second zero-index operator $ \widehat{O}_2 $ such that:

Essentially, a zero-index operator can be "braided," or moved through any other operator, altering its form while preserving its zero-index status. Consequently, all zero-index operators can be collected on one side of an expression. An analogous property holds when moving an operator from right to left. Relation (23) will facilitate subsequent simplifications; its explicit construction is provided in the following subsubsection.

Equation (11) shows that, within a recurrence relation, all contributions sharing the same operator index can be collected into a single term. Expressed in operator notation, the sum of operators with the same index can be compactly written as

where $ \widehat{O}^F_{\vec{o}} $ denotes the representative operator associated with index $ \vec{o} $, and $ \widehat{O}_0 $ denotes a zero-index operator. The representative operator is defined as follows:

where $ \vec{o}_{+} $ contains the positive components of $ \vec{o} $ (with the remaining components set to zero) and $ \vec{o}_{-} $ contains the absolute values of the negative components, such that $ \vec{o} = \vec{o}_{+} - \vec{o}_{-} $. Zero-index operators are linear combinations of operators of the form $ \vec{\eta}^{\,\vec{a}}\vec{\partial}^{\,\vec{a}} $. The identity operator is a special zero-index operator with $ \vec{a}=\vec{0} $. For instance, the combination of operators sharing the same index, $ \widehat{O} = \partial_1 + \eta_2\partial_1\partial_2 $, can be decomposed into a zero-index operator $ \widehat{O}_0 = (1+\eta_2\partial_2) $ and a representative operator $ \widehat{O}_{\vec{o}}^F = \partial_1 $. This decomposition provides an effective method for handling combinations of operators with the same index.

Using the decomposition in (24) and (25), we can explicitly construct the relation between the zero-index operators $ \widehat{O}_1 $ and $ \widehat{O}_2 $ in (23). First, we compute the commutator $ \left[\widehat {O}^F_{\vec{o}}, \widehat{O}_1\right] $. By matching the operator indices, this commutator must take the form $ \left[\widehat {O}^F_{\vec{o}}, \widehat{O}_1\right]=\widehat{O}_{a}\widehat {O}^F_{\vec{o}} $, where $ \widehat{O}_a $ is a zero-index operator. Multiplying both sides by $ \widehat{O}_0 $, we obtain $ \widehat{O} \widehat{O}_1-\widehat{O}_1\widehat{O}=\widehat{O}_{a} \widehat{O} $, yielding $ \widehat{O}_2= \widehat{O}_1+\widehat{O}_{a} $.

The system of differential equations exhibits a hierarchy among operators of different degrees. Many operators are not fundamentally distinct; rather, they can be generated from simpler ones. This motivates the following definition: an operator $ \widehat Q_2 $ is termed a descendant of an operator $ \widehat Q_1 $ if there exists another operator (or a combination of operators) $ \widehat{O}_1 $ such that

where $ \widehat{O}_0 $ is a zero-index operator (or a sum of such operators). According to (26), $ \widehat Q_1 $ is also referred to as the mother of $ \widehat Q_2 $. Because $ \widehat{O}_\delta $ may be the identity operator, an operator can be its own descendant (or mother).

A necessary and sufficient condition for $ \widehat{Q}_2 = \vec{\eta}^{\,\vec{a}_q} \vec{\partial}^{\,\vec{b}_q} $ to be a descendant of $ \widehat{Q}_1 = \vec{\eta}^{\,\vec{a}_o} \vec{\partial}^{\,\vec{b}_o} $ is that

Thus, we define $ \widehat{O}_\delta=\eta^{\vec{a}_q - \vec{a}_o}\vec{ \partial}^{\,\vec{b}_q - \vec{b}_o} $, and $ \widehat{O}_0 $ can be readily solved. Notably, $ \widehat{O}_0 $ contains the identity operator. This structure is fundamental to the algorithm, as it enables the simplification of newly generated differential equations by referencing a controlled set of previously solved "mother" rules.

To facilitate subsequent discussions, we present the construction of $ \widehat{O}_\delta $ and $ \widehat{O}_0 $ for two specific cases:

● Case A: In this case, $ \widehat{Q}_1 = \vec{\partial}^{\,\vec{b}_o} $. It follows that

● Case B: Here, $ \widehat{Q}_1 = \widehat{O}_{0;m} \vec{\partial}^{\,\vec{b}_o} $ and $ \widehat{Q}_2 = \widehat{O}_{0;d} \vec{\eta}^{\,\vec{a}_q} \vec{\partial}^{\,\vec{b}_q} $, where both $ \widehat{O}_{0;m} $ and $ \widehat{O}_{0;d} $ are zero-index operators. Both $ \widehat{Q}_1 $ and $ \widehat{Q}_2 $ are, in fact, combinations of operators with the same index and are written in the form (24). For this more general case, we consider

Now,

where we have used the result (23). From this computation, we obtain

With the operator algebra established, we can now precisely define a symbolic reduction rule. Such a rule isolates a highest-degree operator or operator class and expresses it by solving linear equations in terms of lower-degree operators and possible boundary contributions. Consequently, it is the operator-level analogue of a recurrence relation for the coefficients of the generating function.

For the recurrence relation (11) to be applicable for all values of $ \vec{n}\geqslant 0 $ (i.e., the left-hand side is non-zero), the operator to be solved for, termed a good operator, must satisfy the following two conditions:

1. The index of the target operator, $ \vec{o} = \vec{b}-\vec{a} $, must be non-negative (i.e., $ \vec{o} \geqslant 0 $).

2. As discussed in the previous subsection, if a sum of operators of the same highest degree is evaluated, it can be expressed as $ \widehat{O}_0 \widehat{O}^F_{\vec{o}} $. We then require $ \widehat{O}_0 $ to contain the identity operator with a non-zero coefficient.

If an equation does not contain a suitable operator, it cannot yield a reduction rule. With this established, we classify the differential equations and their corresponding rules according to their operator content:

● T1-type (Type 1): The equation contains a single highest-degree operator index. As the simplest case, it enables a direct solution for this operator index. For example, the toy equation

possesses the unique maximal operator index $ \vec e_1+\vec e_2 $ and is therefore of type T1.

– T1A-type: The equation involves a single operator of the highest degree, yielding a rule of the form:

– T1B-type: The equation involves multiple operators, all sharing the same highest-degree index $ \vec{o} $. These operators can be combined to derive a rule for the representative operator $ \widehat{O}^F_{\vec{o}} $:

where $ \widehat{O}_{0} $ is a zero-index operator.

The T1-type reduction rules are perfect in that their application to Feynman integrals always reduces the degree by at least one.

● T2-type (Type 2): The equation involves multiple operator indices of the same maximal degree. This presents a more complex scenario, requiring the selection of which operator index to solve for. For instance, when $ c\neq0 $, the toy equation

has two distinct maximal indices, $ 2\vec e_1 $ and $ \vec e_1+\vec e_2 $, each of degree two. A global ordering is therefore required to select the target index. Once a target index (denoted $ \vec{o}_1 $) is selected, the resulting rule is further classified:

– T2A-type: The selected target index $ \vec{o}_1 $ corresponds to a single operator in the equation:

– T2B-type: The selected target index $ \vec{o}_1 $ corresponds to multiple operators within the equation:

In both types, the letter A denotes that the selected maximal index is carried by a single standard-form monomial, whereas B indicates that several standard-form monomials share this index and combine to form a zero-index coefficient $ \widehat{O}_0 $ multiplying the representative operator.

Concrete examples are presented in the applications below. The rows of Eq. (58) yield T1A rules in the sunset family, while the coupled rows of Eq. (127) yield T2A rules in the planar double box. Furthermore, Eqs. (129) and (130) demonstrate how the latter are updated to T1A rules.

The T2-type reduction rules are less straightforward because they do not fully reduce the degree in a single step. Consequently, the solution strategy must be chosen with care: the algorithm requires a global ordering of operator indices such that higher-priority operators are expressed in terms of lower-priority operators of the same degree. Once this ordering is established, the remaining operators can subsequently be converted into T1-type reductions, ensuring that the degree is eventually reduced in a controlled manner.

These conditions ensure that our algorithm systematically reduces any integral with sufficiently high powers to simpler integrals. Importantly, in sectors with no master integrals (degenerate sectors), the algorithm naturally produces a reduction rule for the identity operator itself (i.e., $ \vec{b}=\vec{0} $), thereby proving that all integrals in that sector either vanish or reduce to subsector contributions.

An important application of the descendant concept arises when the highest-degree operator in a newly generated equation is itself a descendant of an operator previously resolved by one of the aforementioned reduction rules. In this case, the previously applied rule can be used to eliminate the descendant operator from the equation. We refer to this process as simplification; its algorithmic implementation is discussed in Section 3.2.3.

Consider an equation of the form $ 0=\widehat Q+\cdots $, where $ \widehat Q $ is the highest-degree operator. If $ \widehat Q $ is a descendant of a TA-type reduction rule, such as (31), then from (28) we immediately obtain

Substituting this result back into the original equation replaces the highest-degree operator $ \widehat Q $ with lower-degree operators. If, instead, $ \widehat Q $ is a descendant of a TB-type reduction rule, such as (32), then it follows immediately from (29) that

To eliminate $ \widehat Q $ from the equation, applying $ \widehat{O}_0 $ to both sides yields $ 0 = {\displaystyle \sum}_i c_i \widehat{O}_\delta \widehat{O}_i + \widehat{O}_\delta B + \widehat{O}_0(\cdots) $, where the highest-degree operator is again replaced by lower-degree structures. The same logic applies to the descendants of the T2A- and T2B-type reduction rules.

When solving differential equations, particularly in subsectors, one frequently encounters highest-degree operators whose index satisfies $ \vec{o}\,\, \not\geqslant \,\, 0 $ (see (5)). Although such operators are not suitable targets for reduction rules, they still affect the structure of the equations in which they appear. A practical approach to controlling these operators is to apply the corresponding relations only within restricted regions of the lattice where the problematic terms vanish identically. For example, if the operator $ \eta_2 \dfrac{ \partial^2} {\partial\eta_5^2} $ appears in the differential equation, the term $ (n_5+2)(n_5+1) F_{\vec{n}+2\vec{e}_5-\vec{e}_2} $ arises in the corresponding recurrence relation. Consequently, when applying this recurrence relation to reduce the lattice point $ \vec{n}=(n_1,0,n_3,\dots) $, the contribution of $ \eta_2 \dfrac{ \partial^2}{ \partial\eta_5^2} $ vanishes because $ F_{\vec{n}\not\geqslant 0}=0 $.

When designing the solution order for a T2-type equation, operators with index $ \vec{o}\not \geqslant 0 $ must therefore be handled carefully. To avoid infinite recursion, the reduction rule derived from such an equation should ideally be applied only to restricted lattice regions (for example, $ \vec{n}=(n_1,0,n_3,\dots) $ in the example above), unless the rule is specifically used to reduce the problematic index itself. In practice, this requires other reduction rules to first shrink the irreducible set to the appropriate boundary, after which the partially restricted rule can be applied safely and effectively.

Having established the generating-function formalism, we now present the algorithmic core of this work. The algorithm derives a complete set of symbolic reduction rules for a given integral family by operating directly on the differential equations satisfied by the generating function. Conceptually, the algorithm transforms these equations into a well-defined set of operator identities that act as recurrence relations on the coefficient lattice.

The overall workflow is iterative. Starting from the fundamental IBP-induced differential equations, we extract the immediately solvable reduction rules, apply them to simplify newly generated equations, and subsequently determine which regions of the lattice remain irreducible. If the resulting rule set is incomplete, the remaining irreducible structure dictates which descendant equations must be generated next. Equation generation, rule extraction, and completeness checking are therefore tightly coupled.

This section is organized as follows. In Section 3.1, we define the geometric interpretation of reduction rules and establish the completeness criterion in terms of reducible ($ \mathcal{S} $) and irreducible ($ \mathcal{U} $) lattice sets. We then describe the iterative procedure, which consists of three main modules:

● Module I (Section III.B): Generation of a system of differential equations.

● Module II (Section III.C): Solving these equations to derive valid reduction rules.

● Module III (Section III.E): Verifying the completeness of the derived rules with respect to the lattice criteria.

The overall workflow is illustrated in Fig. 1. If the completeness check in Module III fails, its output is fed back to Module I to guide the subsequent iteration of equation generation. The algorithm is therefore not a one-pass elimination procedure; rather, it is a structured loop in which the remaining irreducible lattice points determine the information required for the next step.

Figure 1.  (color online) The iterative structure of the reduction algorithm.

Before detailing the iterative steps, we first describe the geometric interpretation of a reduction rule on the coefficient lattice. A symbolic reduction rule associated with an operator of index $ \vec{o} $ specifies precisely which lattice points can be mapped to points of lower complexity. This perspective is valuable because completeness is ultimately a matter of lattice coverage, rather than merely the existence of formal operator identities.

Mathematically, the set of lattice points reducible by a rule of index $ \vec{o} $ is defined as:

where the inequality is understood componentwise. Conversely, the lattice points that cannot be reduced by this rule constitute the complementary set:

Our algorithm aims to identify a set of k rules $ \{ \vec{o}_1, \dots, \vec{o}_k \} $ such that the intersection of their irreducible sets,

is minimized. In particular, when the number of points in $\mathcal{U}_{\rm total}$ equals the number of master integrals (MIs), the reduction is complete.

The first module of each iteration generates the differential equations that feed the solver. These equations are derived from two sources: the fundamental IBP identities and the descendant equations ³ generated from rules identified in previous iterations. This module serves an exploratory role: it yields candidate operator structures that may lead to new reduction rules.

The algorithm begins with a set of seed equations derived by applying IBP identities to the generating function $ G(\vec{\eta}) $, yielding

where $ R=\ell_i $ for $ i=1,\ldots,L $ or $ R=p_j $ for $ j=1,\ldots,E $. These IBP identities can be transformed into differential equations in the variables $ \eta_i $ via the mappings

The explicit expressions for these IBP identities are as follows:

● For $ R=\ell_b $, where $ b=1,\ldots,L $, we have

where $ {\alpha}_{ab,k} $ and $ {\beta}_{ab,k} $ are defined by

● For $ R = p_b $, where $ b = 1, \ldots, E $, we have

where $ \widetilde {\alpha}_{ab,k} $ and $ \widetilde {\beta}_{ab,k} $ are defined by

Based on this derivation, several important observations can be made regarding the initial seed equations:

1. They are at most second-order differential equations.

2. The degree-two operators take the form $ \partial_i \partial_j $, where one index corresponds to a propagator and the other to an ISP.

3. The degree-one operators $ \partial_i $ exclusively involve propagator indices.

4. The boundary terms involve generating functions of the subsectors evaluated at $ \eta_k=0 $, reflecting a simpler topology.

These properties are specific to the form of the generating functions given in (41). When generalized to other contexts involving different generating functions, they must be modified accordingly.

After the initial round, subsequent iterations yield increasingly complex equations. This is accomplished by differentiating the equations or reduction rules derived from the previous solving module.

We apply the derivative operators $ \partial_i \equiv \partial/\partial\eta_i $ to the set of equations identified as "solvable" in the previous round (further details in Section III.C). In principle, one could apply all possible $ \partial_i $ to all available equations. However, this approach is computationally expensive and often redundant. We therefore employ a selection strategy to prune the generation process:

● Criterion 1 (Efficiency): If a reduction rule for a specific operator (e.g., $ \partial_a $) has been established, it follows that any power of the index $ n_a $ can be reduced. Consequently, the application of $ \partial_a $ to generate new equations can be excluded, as this case is already accounted for.

● Criterion 2 (Hierarchy): Equation generation can be prioritized from a specific subset of existing equations, such as exclusively from T1-type or degree-one rules. The selected subset may be updated during subsequent iterations of the algorithm.

● Criterion 3 (Cutoff): To manage computational complexity, we introduce a cutoff degree $ M_0 $ and solve equations only up to degree $ M \leqslant M_0 $ after full simplification (see Section III.B.3). Furthermore, equations of degree $ M<0 $ are not solved. As discussed in Section II.E, the target of a reduction rule must have an operator index $ \vec{o}\geqslant 0 $. For $ M<0 $, every highest-degree operator necessarily satisfies $ \vec{o}\,\,\not\geqslant\,\, 0 $; thus, such equations serve solely as sources for generating further equations rather than as direct reduction rules.

For equations of degree $ M=0 $, the situation is slightly more subtle because the only valid operator index is $ \vec{o}=\vec{0} $, corresponding to the zero-index operator. If the equation contains the identity operator as the highest-degree operator, we can solve it and conclude that the lattice point $ \vec{n}=\vec{0} $ reduces to subsectors, indicating that the sector contains no master integrals. Conversely, if the sector is known to contain nontrivial master integrals, no valid operator is available, and the equation should not be solved. Both situations arise in the examples below.

The IBP relations naturally suggest an initial value of $ M_0=2 $. In most practical examples, $ M_0=2 $ is sufficient; however, the algorithm permits $ M_0 $ to be increased if necessary.

In practice, Criterion 1 is highly effective. If the algorithm stalls (i.e., fails to identify a complete set of rules), these restrictions should be relaxed to encompass a broader range of source equations.

Each newly generated equation must be simplified using the complete set of currently known reduction rules. This is arguably the most critical step of the algorithm, as it reveals the novel operator structures that must be solved. If an equation simplifies to $ 0=0 $, it is classified as a trivial descendant and discarded, as it yields no new information. Otherwise, it is classified as a non-trivial descendant and added to the pool of equations to be solved in Module II.

Once an equation has been generated, whether through the action of $ \partial_i $ or by linear combination, it must be simplified using the complete set of known reduction rules. This step is essential because it eliminates operator content that has already been resolved, isolating only the genuinely novel structures that remain to be solved.

The simplification strategy relies on two ordered lists of reduction rules: $ \mathcal{R}_{A} $ (TA-type rules) and $ \mathcal{R}_{B} $ (TB-type rules) (see, for example, (31), (32), (34), and (35) in Subsection B.5). We prioritize simpler rules to maintain manageable expressions. Thus, we perform the simplification in two distinct rounds:

● First Round: Simplify using only $ \mathcal{R}_{A} $.

● Second Round: Simplify the first-round result using $ \mathcal{R}_{B} $.

Within each list, simpler reduction rules from earlier iterations and T1-type rules are assigned higher priority.

Given a rule list $ \mathcal{R} $, the procedure for a single round of simplification follows the workflow depicted in Fig. 2.

Figure 2.  (color online) Flowchart of the simplification.

● Step 1: Extract and Sort Highest-Degree Operators. Given an equation E, we identify all highest-degree differential operators and sort them to optimize the matching process. The sorting scheme is based on the finer index $ (\vec{a}, \vec{b}) $. Operators with a lower degree in $ \vec{b} $ take precedence. For operators of equal degree, we adopt a convention, such as prioritizing those involving ISPs over propagators. Let this sorted list of highest-degree operators be $ \mathcal{O}_{ops} = \{ \widehat{O}_1, \widehat{O}_2, \dots \} $.

● Step 2: Selection and Matching. For each operator $ \widehat{O}_i $ in $ \mathcal{O}_{ops} $, we attempt to match it against the rules in $ \mathcal{R} $ based on Section II.D.3. Given a specific operator $ \widehat{O}_i $ and a candidate reduction rule for the operator $ \widehat{R}_j $, we determine whether $ \widehat{O}_i $ is a descendant of $ \widehat{R}_j $ according to the criterion (27):

1. f it is not satisfied, we continue checking other rules.

2. If this condition is satisfied, $ \widehat{O}_i $ is a descendant of $ \widehat{R}_j $. Before finalizing this match, we verify whether substituting this expression into equation E yields a trivial identity ($ 0=0 $) ⁴:

(a) If the match is trivial, we discard it and proceed to evaluate the remaining rules.

(b) If it is non-trivial, we record the pair $ (\widehat{O}_i, \widehat{R}_j) $ as a valid simplification (to be applied later according to the prescription specified in (36) or (37)) and immediately terminate the search for $ \widehat{O}_i $.

● Step 3: Execute the simplification. After evaluating all operators in $ \mathcal{O}_{ops} $, we determine whether any valid simplification pairs have been identified.

1. If no pairs are found, the current set of operators cannot be reduced further using the list $ \mathcal{R} $, and the process terminates for the current iteration.

2. If pairs are found, we apply the substitutions for all recorded pairs as specified in (36) or (37). For T1-type reduction rules, the substitutions can be performed simultaneously. For T2-type reduction rules, it is advisable to apply the substitutions sequentially, verifying after each step whether the pair $ (\widehat{O}_i, \widehat{R}_j) $ still appears in the updated equation.

After completing all substitutions, we obtain the modified equation $ E' $.

● Step 4: Verification and Iteration To ensure the robustness of the results, we conduct a final verification of $ E' $. We determine whether the modified equation contains any remaining highest-degree operators that are candidates for reduction.

1. If any such operators remain, we update the current equation $ E \leftarrow E' $ and return to Step 1 to perform another iteration of extraction and matching.

2. If no such operators remain, the process terminates for this iteration. It should be emphasized that the resulting equation must be assigned to the lower-degree category for subsequent simplification.

In practice, the same high-degree descendant operator $ \widehat{O}_i $ frequently appears across multiple generated equations. To enhance computational efficiency, we implement a caching mechanism that records the mapping between a descendant operator and its corresponding parent rule, i.e., $ \widehat{O}_i \to \widehat{R}_j $. When the algorithm encounters an operator $ \widehat{O}_i $ previously matched in an earlier simplification instance, it retrieves the substitution directly from the cache, thereby bypassing the redundant search and derivative construction in Step 2.

For the second round of simplification using rules from $ \mathcal{R}_{B} $, additional notation beyond the basic logic is required. As noted in Section II.E, TB-type rules may require applying a zero-index operator $ \widehat{O}_{0} $ to the entire equation E to isolate the target term, transforming $ E \to \widehat{O}_{0} E $ (see (37)). Although rigorous, this procedure increases complexity. In practice, we frequently omit this global multiplication during the generation phase, reserving TB-type simplifications for the solving module, where they are applied strictly as needed.

Upon completion of the simplification rounds, the resulting system of equations is naturally filtered. Equations reduced to mathematical identities ($ 0=0 $) are discarded, as they provide no new constraints. The remaining set of non-trivial differential equations serves as the refined input for the subsequent Module II.

The primary input for this module is the simplified system of differential equations generated by Module I. The objective is to decouple this system and extract explicit reduction rules. Whenever a suitable high-degree operator is solved, the corresponding rule enlarges the reducible region $ \mathcal{S}_{\vec{o}} $ of the lattice. In this sense, Module II serves as the stage where algebraic information is translated into geometric progress.

To handle the complexity of the system systematically, we organize the equations by their total differential degree. The algorithm proceeds hierarchically, resolving the highest-degree sector (typically degree-2) before propagating the results to lower-degree sectors.

For a fixed degree M, the solution procedure comprises two distinct steps: Linearization (Sec. III.C.1), which classifies equations via Gaussian elimination, and Operator Selection (Sec. III.C.2), which applies logical criteria to determine reduction targets. Finally, an Update Mechanism (Sec. III.D) ensures that new reduction rules are propagated globally. These three operations must be repeated for each set of equations partitioned by degree and type.

Let $ \mathcal{E}_M $ denote the subset of equations with maximum degree M. We aim to isolate the linearly independent constraints on degree-M operators. Each equation in $ \mathcal{E}_M $ can be decomposed as:

where $ \widehat{O}_{\vec{o}_a}^{(M)} $ denotes the unique operator with a given operator index $ \vec{o}_a $ of degree M (i.e., $ \widehat{O}^{F}_{\vec{o}_a} $ in (24) in the preceding section), $ c_{i,a} $ are rational coefficients, $ \widehat{O}_{0,i,a} $ are zero-index operators, and $ \widehat{R}^{(<M)} $ encapsulates all terms of degree strictly less than M. The combinations $ c_{i, a} \widehat{O}_{0,i,a} $ are treated as scalars since zero-index operators commute with one another.

We construct a coefficient matrix C with entries $ c_{i,a}\widehat{O}_{0,i,a} $ and perform Gaussian elimination ⁶ to reduce it to row-reduced form. This linear algebra step arises naturally from the input structure and classifies the resulting rows into three categories:

1. Null Rows (Seeds for Lower Degrees): If a row yields zero coefficients for all degree-M operators, the equation reduces to the form $ \widehat{R}'^{(<M)} = 0 $. Such equations describe constraints entirely within lower-degree sectors. They are removed from the current degree-M processing and transferred to the corresponding lower-degree subsets only after being fully simplified by the known reduction rules detailed in III.B.3.

2. Type-1 Rows (Direct Reductions): A row containing exactly one non-zero coefficient for a degree-M operator $ \widehat{O}_k $ enables immediate reduction:

This unambiguously establishes a reduction rule for $ \widehat{O}_k $. We refer to these equations as T1-type equations.

3. Type-2 Rows (Coupled Relations): A row containing two or more non-zero coefficients defines a coupled constraint:

Solving such T2-type equations requires selecting one operator $ \widehat{O}_{p_i} $ to serve as the reduction target (LHS), expressing it in terms of the remaining operators. This selection is non-trivial and determines the efficiency of the reduction on the lattice $ \mathcal{S}_{\vec{o}} $, as detailed in the following subsection.

Importantly, the linear combinations (48) and (49) serve as "seeds" that generate non-trivial relations in subsequent steps, as discussed in Section III.B.2.

Solving T1-type equations is straightforward; they directly yield T1A-type or T1B-type reduction rules. In contrast, T2-type equations require more careful treatment. Specifically, a global priority ordering must be defined to determine which operator is solved for and placed on the left-hand side. The stability and convergence of the algorithm strongly depend on this choice.

The solution process is executed iteratively, prioritizing the simplest reductions. We begin by resolving all available T1-type equations. The resulting rules are then used to update the remaining T2-type equations (as described in Section 3.4), potentially simplifying some T2-type equations into new T1-type equations. We repeat this cycle until no further T1-type equations remain. Only then do we select targets within the remaining T2-type equations, according to the following ordered criteria:

● Criterion 1: Validity (Negative-Index Avoidance). Operators with an index $ \vec{o}\not\geqslant 0 $ should not be selected as solution targets. If an equation contains only such operators as its highest-degree terms, it should be retained as a source for subsequent equation generation rather than directly converted into a reduction rule.

● Criterion 2: Lattice reduction (the "good" operator). We prioritize solving for an operator $ \widehat{O}_i $ if it lowers a specific lattice index $ n_a $ such that this reduction is not offset by the competing operators in the same equation.

Formally, let $ \vec{v}_i $ and $ \vec{v}_j $ be the lattice shifts associated with operators $ \widehat{O}_i $ and $ \widehat{O}_j $ (i.e., the indices of the operators). We define $ \vec{v}_{ij} \equiv \vec{v}_i - \vec{v}_j $. If there exists a component a such that $ (\vec{v}_{ij})_a \gt 0 $ for all competing operators j in the equation, then solving for $ \widehat{O}_i $ ensures a net reduction in the index $ n_a $. In geometric terms, the recurrence moves lattice points toward the boundary of $ \mathcal{S}_{\vec{v}_i} $.

We emphasize that a single equation may contain multiple such "good" operators, each reducing a different lattice direction. Subsequent criteria are applied only when this ambiguity remains unresolved by the desired reduction flow.

● Crit. 3: Boundary-aware hierarchy. In certain cases, an equation contains an operator $ \widehat{O}_{\text{bad}} $ whose index has a negative component along the t-th direction. Such an operator vanishes identically on the boundary $ n_t=0 $. This observation can be exploited strategically: one first employs other equations to constrain the irreducible set to this boundary, and subsequently applies the current equation to derive an additional reduction rule once the problematic term has become ineffective.

A closely related situation arises when the same equation also contains a well-behaved operator that directly reduces $ n_t $, such as $ \partial/\partial\eta_t $. Under these conditions, the equation can still be applied safely, as the presence of $ \widehat{O}_{\text{bad}} $ does not obstruct the reduction of the index that controls its boundary behavior.

● Criterion 4: Empirical hierarchy for variables. If the preceding criteria remain insufficient to select a unique target, we employ an empirical ordering that has proven effective in practice. For degree-two operators, the suggested precedence is:

where the indices a and b denote ISP variables and i and j denote propagators. For degree-1 terms, we adopt the ordering $ \partial_a \succ \partial_i $.

Once a new TA-type reduction rule is derived, it must be propagated throughout the system via the simplification procedure (see Section III.B.3). This global update ensures the consistency of the iterative process: it eliminates solved operators from all remaining equations and frequently reveals new T1-type structures that were previously obscured.

During this update, we treat the unsolved equations and the library of previously established rules uniformly. The procedure is as follows:

1. Unification: We first convert all previously solved reduction rules back into equation form ($ \text{LHS} - \text{RHS} = 0 $) and combine them with the current set of unsolved equations.

2. Global Substitution: We apply the newly derived reduction rules to this unified set of equations to perform the simplification procedure (see III.B.3). Each occurrence of the reduced operator is replaced with its corresponding expression. This step eliminates the solved operator from the entire system, ensuring it does not reappear as a variable.

Following the global update, the system is effectively reset. The updated pool of equations is fed back into the solver, and the process is repeated until the set of reduction rules no longer changes. At this stage, the rules are maximally decoupled within the chosen strategy.

Note that, for simplicity, the update may be restricted to TA-type reduction rules; TB-type reduction rules are applied solely during the solving step.

Once a candidate set of reduction rules is obtained, we must determine its completeness. This is accomplished by analyzing the total irreducible set $\mathcal{U}_{\rm total}$ defined in Eq. (40). Thus, the completeness check serves a dual role as both a stopping criterion and a diagnostic tool. If the rule set is incomplete, the geometry of $\mathcal{U}_{\rm total}$ indicates the type of equations to be generated next.

We denote the number of lattice points in $ \mathcal{U}_{\rm total} $ by $ |\mathcal{U}_{\rm total}| $. Two cases must be distinguished:

If the number of master integrals, $N_{\rm MI}$, is known a priori (e.g., via geometric or algebraic methods), the stopping condition is simply:

If $|\mathcal{U}_{\rm total}| \gt N_{\rm MI}$, the rule set is incomplete, and the subsequent iteration must generate additional equations.

Another useful feature of the algorithm is its ability to determine $N_{\rm MI}$. Upon obtaining a finite set $\mathcal{U}_{\rm total}$ without prior knowledge of the master-integral count, it remains necessary to verify the completeness of this set. To this end, we perform a consistency test by reducing selected points outside $\mathcal{U}_{\rm total}$ along different reduction paths (i.e., using different combinations of rules).

● If all paths yield identical expressions, the reduction is consistent and $N_{\rm MI} = |\mathcal{U}_{\rm total}|$.

● If different paths yield inconsistent expressions, hidden relations (IBP identities) must exist among the points in $\mathcal{U}_{\rm total}$. These relations must be incorporated into the system, and $\mathcal{U}_{\rm total}$ must be further reduced until consistency is restored.

Having established the general algorithm in Section III, we now turn to explicit examples. We begin with the top sector of the sunset diagram because it is the minimal setting in which all the main ingredients appear. Although the topology is simple, it exhibits the essential interplay among generating functions, operator identities, descendant equations, and lattice-based completeness. Consequently, we treat the sunset example in greater detail than the subsequent case studies, using it as the primary pedagogical thread of the paper.

The sunset kinematics, the top-sector and subsector generating functions, and a representative fundamental IBP equation were introduced in the preview provided in Section II.C. In the remainder of this section, we focus on $G_{111}$ and present a complete iterative analysis: which equations are generated, how they are simplified, what remains irreducible after each round, and how the reduction closes in different kinematic regimes.

We begin with the fundamental IBP relations and derive the first layer of reduction rules. This initial iteration is crucial because it reveals the basic geometry of the reduction problem: it identifies which directions in the lattice are immediately controlled and which subsets remain unresolved, requiring subsequent treatment in later rounds.

Generating and organizing the IBP relations

Applying the general formulas in Eqs. (43) and (45) to the sunset topology yields six fundamental differential equations. We organize these equations according to the degree of the operators acting on the generating function. Because the present discussion focuses on the top sector, we suppress boundary terms involving subsector generating functions to maintain a transparent operator structure.

The first relation and the shorthand $f_{abc}$ were presented in Eqs. (20) and (19), respectively, as part of the preview. The remaining five differential equations are listed below:

Solving the system of equations

Having established the differential system, we isolate its highest-degree part and derive the associated reduction rules. In the present case, all six equations possess leading operators of degree two; therefore, we first analyze the corresponding coefficient matrix:

where $ \partial_i \equiv \dfrac{\partial}{\partial \eta_i} $. Applying Gaussian elimination, we transform the matrix into a form in which each row contains a leading term corresponding to a specific highest-degree operator:

Each row in (58) is of type T1A; consequently, each directly yields a degree-two reduction rule in terms of lower-degree operators. Thus, we obtain the following six primary reduction rules:

The six highest-degree operators on the left-hand side exhibit a distinct pattern:

which corresponds to the product of derivatives with respect to the propagator variables ($ \eta_1, \eta_2, \eta_3 $) and the ISP variables ($ \eta_4, \eta_5 $). This is the first concrete manifestation of the mixed propagator–ISP pattern anticipated in the general formalism.

Completeness check

We now interpret these rules on the coefficient lattice to determine what has and has not been reduced. Translating the differential equations into recurrence relations for the Taylor coefficients makes this interpretation explicit. For example, Eq. (59) implies:

where $ \vec{e}_i $ denotes the unit vector in the i-th direction. The term on the left-hand side, $ F_{\vec{n}+\vec{e}_3+\vec{e}_5} $, has the highest index sum. This relation allows us to reduce any integral with the index $ (n_1, n_2, n_3, n_4, n_5) $, provided that $ n_3 \geqslant 1 $ and $ n_5 \geqslant 1 $. Consequently, the integrals that cannot be reduced by this rule belong to the set $ (n_1, n_2, 0, n_4, n_5) \cup (n_1, n_2, n_3, n_4, 0) $.

By combining the constraints from the six reduction rules derived above, we identify the lattice points that remain irreducible at this stage. Following the completeness check specified in Eq. (40), we determine the irreducible sets to be:

Since both $ \mathcal{U}_1 $ and $ \mathcal{U}_2 $ still contain infinitely many lattice points, the first iteration is incomplete. Nevertheless, it has already partitioned the unreduced region into two well-defined components: one dominated by ISP powers and the other by propagator powers. This geometric decomposition provides precisely the information required to guide the subsequent iteration.

Before proceeding, it is instructive to separate the discussion into two kinematic branches. The degree structure of the descendant equations depends on the mass configuration: when $ m_i=0 $, several formally degree-two equations degenerate into degree-one equations, altering both the available reduction rules and the final master-integral count. We therefore discuss the massless and massive cases separately, beginning with the massless case where the mechanism of complete reduction is most transparent.

Generating the descendant equations

At this stage, no index has yet been reduced to a finite set; thus, the algorithm proceeds by generating descendant equations. Applying $ \partial_i $ for $ i=1,\dots,5 $ to the six first-round rules in Eqs. (59)–(64) initially yields $ 6\times 5=30 $ candidate degree-three equations.

However, not all generated equations are independent or relevant. We simplify them according to the descendant logic and integrability conditions discussed in Section III.B.2. We illustrate this simplification procedure with two representative examples:

1. Consider the action of $ \partial_5 $ on the reduction rule (60). This action generates the degree-three operator $ \partial_5 \partial_5 \partial_1 $, which is a descendant of the previously solved degree-two operator $ \partial_5 \partial_1 $. Substituting the known rule for $ \partial_5 \partial_1 $ reduces the equation identically to zero. Consequently, this equation is trivially generated and yields no new information.

2. Consider the action of $ \partial_4 $ on the reduction rule (60). This action generates the operator $ \partial_4 \partial_5 \partial_1 $, which has two parents: $ \partial_5 \partial_1 $ (from rule (60)) and $ \partial_4 \partial_1 $ (from rule (61)). Applying rule (60) yields zero. However, applying rule (61) reduces the degree-three operator to lower-degree terms. This procedure is equivalent to enforcing the integrability condition $ [\partial_4, \partial_5] = 0 $. After simplifying all resulting degree-two terms using the six fundamental rules, we obtain a non-trivial degree-one equation:

where the operators $ {\partial_i} $ appear as the highest-degree operators, potentially yielding a new reduction rule.

Of the 30 generated equations, 12 are trivially zero. Of the remaining 18, symmetry considerations (in particular, the integrability conditions) reveal that only 9 are distinct. A further reduction yields 8 independent descendant equations. The first is given by Eq. (68), and the remaining 7 are listed below:

1. Applying $ [\partial_4, \partial_5] $ to $ \partial_2 $: we have

2. Applying $ [\partial_4, \partial_5] $ to $ \partial_3 $: we have

3. Applying $ [\partial_1, \partial_2] $ to $ \partial_4 $: We obtain

4. From the action of $ [\partial_1, \partial_3] $ on $ \partial_4 $: we obtain

5. From $ [\partial_2, \partial_3] $ acting on $ \partial_4 $: we obtain

6. Applying $ [\partial_1, \partial_2] $ to $ \partial_5 $: we obtain

7. From the action of $ [\partial_1, \partial_3] $ on $ \partial_5 $: we obtain

8. From the action of $ [\partial_2, \partial_3] $ on $ \partial_5 $: we obtain

Solving the system in the massless case

A crucial feature of the new equations is that the coefficients of their highest-degree operators depend on the masses. Consequently, when $ m_i=0 $, several equations that are quadratic in the generic massive case become linear. This reduction is not merely technical; it directly reflects the smaller number of master integrals in the massless theory compared to the generic massive case ⁶. In the massless limit, only a single master integral remains.

We now set $ m_i = 0 $ in the equations derived above. Equations (73) and (75) degenerate into linear equations. Simplifying these using the rules from the previous iteration, we obtain:

We now solve the system, beginning with the remaining degree-two equations. The coefficient matrix corresponding to the operators $ \partial_1\partial_2 $, $ \partial_1\partial_3 $, and $ \partial_2\partial_3 $ is given by:

Thus, a combination yields a degree-one equation. Simplifying the new degree-one equation, $ (74)+ (71) $, using the descendant reduction rules yields:

There remain three degree-two T1A-type equations, which provide reduction rules for the operators $ { \partial_1} { \partial_2} $, $ { \partial_1} { \partial_3} $, and $ { \partial_2} { \partial_3} $ as

which enable the reduction of cross-terms between propagators.

Having derived three new reduction rules, we would, in principle, update the six previously solved reduction rules and the remaining degree-one equations as described in Section III.D. In the present case, however, no further updates are required.

Next, we consider the six degree-one equations. Combining the degenerate equations with the results of the degree-two simplification, we construct the coefficient matrix for the single derivatives as

where, for brevity, we define the following differential operators:

Applying Gaussian elimination yields the reduced matrix form:

where

We observe that the last row yields a trivial equation, confirming the redundancy anticipated from integrability conditions. The first row yields a degree-zero equation; however, it lacks a suitable operator to solve for and will therefore not be utilized further. We thus focus on the four independent equations represented by rows 2–5.

The second and third rows of the reduced matrix decouple from the propagator derivatives, allowing us to solve explicitly for $ \partial/\partial \eta_4 $ and $ \partial/\partial \eta_5 $, thereby yielding two T1B-type reduction rules:

which ensure the exact reduction of any integrals with nonzero indices $ (n_4, n_5) $ to $ (0, 0) $.

The fourth and fifth rows of the matrix relate the derivatives of the propagator variables. By rearranging Eq. (77), we can solve for the operator involving $ \eta_1 $ and $ \eta_2 $, yielding two T2B-type reduction rules:

These relations enable the reduction of any integral with indices $ (n_1, n_2, n_3) $ to a linear combination of integrals with indices $ (0, 0, N) $, where $ N \leqslant n_1+n_2+n_3 $, without increasing the complexity of the integrals.

Because all four new reduction rules are of the TB type, we do not perform the update step. The reason for this is evident from the present example: once the T1B-type reduction rule for $ \dfrac{ \partial}{ \partial\eta_4} $ in (89) is established, the previously derived degree-two operators $\dfrac { \partial}{ \partial\eta_4}\dfrac{ \partial}{ \partial\eta_i} $ with $ i=1,2,3 $ become descendants of that rule. To preserve the simpler and more convenient TA-type reduction rules, we therefore exclude TB-type rules from the global update.

Completeness check for the massless case

At this stage, we have obtained seven new reduction rules: four degree-one rules, namely (89), (90), (92), and (90), along with three degree-two rules, namely (71), (72), and (76). Applying the four degree-one rules to the remaining lattice points in (67), we find that the points in the set

cannot be reduced using these rules.

Applying the three degree-two reduction rules does not resolve this, as the set $ {\cal U}_3 $ remains unreduced. The reason is that these rules are descendants of the operators $ \dfrac{ \partial}{ \partial\eta_i} $ for $ i=1,2 $. Nevertheless, they remain useful for simplification in the final iteration.

The third iteration: reducing the final sector

Because only $ n_3 $ has not yet been fully reduced, we perform a third iteration by applying $ \partial_3 $ to the four degree-one equations derived in the previous step, namely (89), (90), (92), and (90).

First, applying $ \partial_3 $ to Eq. (89) yields an equation containing three degree-two operators: $ \partial_3 \partial_4 $, $ \eta_4 \partial_4\partial_3 \partial_4 $, and $ \eta_5 \partial_5\partial_3 \partial_4 $. All of these descend from the operator $ \partial_3 \partial_4 $. Upon simplification, the new equation involves five degree-one operators: $ \partial_3 $, $ \eta_3 \partial_3^2 $, $ \eta_1 \partial_1 \partial_3 $, $ \eta_2 \partial_2 \partial_3 $, and $ \eta_5 \partial_4 \partial_3 $. These operators can be simplified using (63) from the first iteration, together with (71), (72), and (76) from the second iteration. The equation then reduces to the following T1B-type equation:

with

Alternatively, applying $ \partial_3 $ to Eq. (90) yields another T1B-type equation involving the same operator $ \widehat{O}_{0;II} $. Indeed, when restricted to the irreducible set $ \mathcal{U}_3 $ (where $ n_1=n_2=n_4=n_5=0 $), this equation reduces to the constraints in Eq. (94), thereby confirming their consistency.

Finally, applying $ \partial_3 $ to Eqs. (92) and (91) initially yields second-degree equations involving the operators $ {\partial_3^2 } $ and $ \eta_3{\partial_3^3 } $. However, Gaussian elimination combines these two equations into a first-degree relation. Upon simplifying the first-degree terms using existing rules, the resulting relation assumes an instructive form:

which is a linear combination of Eqs. (92) and (91), and is therefore redundant for the reduction of $ \mathcal{U}_3 $.

Applying one of these T1B-type equations yields a T1B-type reduction rule for $ {\partial_3} $. Consequently, the remaining propagator index $ n_3 $ can be reduced to zero. Thus, the only integral that cannot be reduced by any of these rules is

This single point constitutes the master integral basis for the massless sunset topology. The massless branch thus serves as a complete example of how the iterative procedure closes: the first round isolates the unreduced subsets, the second collapses them to a single line, and the third removes the final propagator tower.

Having established the massless case, we now consider the fully massive sunset diagram. In this scenario, the three masses and the external invariant $ K^2 $ restore the generic degree structure of the equations and expand the final master-integral basis. Nevertheless, the logic of the algorithm remains unchanged: first reduce the ISP directions, then constrain the propagator directions, and finally close the remaining unreduced set.

From the six degree-two equations, the corresponding matrix is

where $ f_{\pm\pm\pm} $ was previously defined in (19). Using Gaussian elimination, the combination

is zero, leaving five independent T2-type equations for the six degree-two operators. To proceed, we adopt an ordering scheme in which five operators (e.g., $ \partial_1^2 $ and $ \partial_2^2 $) are expressed in terms of the sixth (e.g., $ \partial_3^2 $). No updates are required for these five derived reduction rules.

Next, we address the three degree-one equations. The corresponding matrix contains derivatives with respect to all five degree-one operators:

After Gaussian elimination, the system reduces to

where

We identify the last two rows as T2-type equations, which are suitable for solving for the ISP operators $ \partial_4 $ and $ \partial_5 $. Assuming the standard ordering, in which ISP operators take precedence over propagator operators, we derive the following explicit reduction rules:

Since these relations will be applied exclusively to the set $ (0,0,0,n_4,n_5) $, we retain only the terms that contribute nontrivially ⁷.

At this stage, we have seven new reduction rules, two of which are degree-one rules for $ n_4 $ and $ n_5 $. Applying these to the irreducible set (67) reduces $ n_4 $ and $ n_5 $ to zero, yielding the remaining irreducible set $ (n_1, n_2, n_3, 0, 0) $. Applying the five degree-two propagator equations further constrains the system, rendering this set reducible except for $ {\cal U}_3 $ in (93) along with the two isolated points $ (1,0,0,0,0) $ and $ (0,1,0,0,0) $. The existence of infinitely many irreducible points in the $ n_3 $ direction necessitates a third iteration.

The third iteration: reducing the final set

Because only $ n_3 $ remains unreduced, we proceed by applying $ \partial/\partial \eta_3 $ to the existing equations. We focus on the descendant of the degree-one equation $ (69) + (70) $. Applying $ \partial_3 $ to the full expression of this equation generates terms involving degree-two operators, such as $ \partial_3\partial_4 $, $ \partial_3\partial_2 $, $ \partial_3\partial_1 $, and $ \partial_3^2 $. Most of these can be simplified using the reduction rules obtained in the previous iterations. Crucially, the $ \partial_3^2 $ term does not cancel.

After simplification, the new equation yields a relation involving the following degree-two propagator structure:

where $ h_i $ and $ g_{ijk} $ are defined in (102). This new equation is linearly independent of the five T2-type equations from the system (98). Therefore, we can uniquely solve for all six degree-two operators $ {\partial_1^2} $, $ {\partial_2^2} $, $ {\partial_3^2} $, $ {\partial_1}{\partial_2} $, $ {\partial_1}{\partial_3} $, and $ {\partial_2}{\partial_3} $, yielding six T2A-type reduction rules.

Consequently, any integral with $ n_3 \geqslant 2 $ can be reduced. Combined with the previous steps, this ensures that the set of irreducible integrals is finite and consists of

These four points represent the master integrals for the generic massive sunset topology, consistent with standard expectations. The contrast with the massless branch highlights the role of kinematics: while the same algorithmic framework applies to both cases, the degree structure of the descendant equations determines the rate at which the irreducible lattice collapses and the number of surviving master integrals.

With both kinematic branches established, it is instructive to consider the practical application of the reduction rules. The algorithm typically yields more valid relations than are strictly necessary for completeness; consequently, an operational reduction strategy must be selected. For the sunset topology, two natural strategies emerge.

Strategy I: The minimal basis

The first approach prioritizes a minimal set of recurrence relations that remain valid across the entire kinematic lattice. A convenient choice consists of the following five general relations:

● Eqs. (89) and (90) for the ISP directions $ \eta_4 $ and $ \eta_5 $,

● Eqs. (92) and (91) for the propagator directions $ \eta_2 $ and $ \eta_1 $,

● Eq. (94) for the final direction $ \eta_3 $.

The reduction flow is then straightforward. First, the ISP relations are employed to reduce any integral with $ n_4>0 $ or $ n_5>0 $ to the sector with $ n_4=n_5=0 $. Next, the propagator relations are applied to reduce $ n_2 $ and $ n_1 $ to zero, at the expense of shifting weight to $ n_3 $. Finally, the last relation is applied to reduce the remaining $ n_3 $ tower to the master integral. While this strategy is compact, the individual equations retain their full algebraic complexity.

Strategy II: The hierarchical approach

The second approach exploits the fact that the boundary regions of the lattice satisfy significantly simpler equations. This strategy employs a larger set of rules but applies them in stages, ensuring that each step is algebraically simpler.

1. General reduction. Starting with the six fundamental degree-two rules in Eqs. (59)–(64), we reduce a generic point $ (n_1,\dots,n_5) $ to the boundary sets defined in Eq. (67).

2. ISP reduction on the boundary. For the subset $ (0,0,0,n_4,n_5) $, terms proportional to $ n_1 $, $ n_2 $, and $ n_3 $ vanish upon expansion; thus, the general ISP relations simplify considerably. Setting $ \eta_1=\eta_2=\eta_3=0 $ in Eq. (89) yields

and setting the same variables to zero in Eq. (90) yields

3. Propagator reduction along the constrained line. Along the line $ (n_1,n_2,n_3,0,0) $, setting $ \eta_4=\eta_5=0 $ simplifies Eqs. (92) and (91) to

Likewise, the specialized relation for the final subset $ (0,0,n_3,0,0) $, derived from Eq. (94), reduces to

Of the two strategies, Strategy II is often computationally superior. Although it necessitates a larger rule set, the individual algebraic manipulations on the boundary are substantially less demanding, thereby preventing the excessive growth of intermediate expressions. Because the detailed optimization of reduction paths is implementation-dependent, we defer this issue to future work. Nevertheless, the sunset example illustrates the primary point: the framework inherently generates both globally valid rules and specialized boundary rules, and either category can be organized into a consistent reduction strategy.

The rule derived thus far can be applied to reduce integrals with arbitrary propagator powers and tensor ranks. In practice, integrals may exhibit special properties, such as unit propagator powers. For the sunset integral, this corresponds to the lattice point where $ n_1=n_2=n_3=0 $. At such specific points, we can use Eqs. (89) and (90) to efficiently reduce the tensor structures. Explicitly, Eqs. (89) and (90) simplify to

and

In the literature, syzygy and module-intersection techniques achieve a comparable result, specifically concentrating the reduction on integrals with unit-power propagators. This approach is implemented in NeatIBP, a tool widely employed in state-of-the-art computations. Here, we show that our algorithm provides an alternative realization of this special case (see Eqs. (113) and (114)).

Next, we consider the two-loop planar massless double-box diagram shown in Fig. 3. Compared to the sunset example, this topology constitutes a more demanding case: it features more propagators, two independent kinematic scales, and a substantially larger reduction lattice. Consequently, it serves as the first rigorous test of whether the same algorithmic logic scales beyond the minimal two-loop family.

Figure 3.  The two-loop planar double-box diagram.

In this problem, four external momenta, $ k_1 $, $ k_2 $, $ k_3 $, and $ k_4 $, satisfy the on-shell conditions $ k_i^2=0 $ and momentum conservation, $ {\displaystyle \sum}_{i=1}^4 k_i =0 $, leaving two independent scales, $ s=2 k_1 \cdot k_2 $ and $ t=2 k_2 \cdot k_3 $. The complete set of Lorentz scalars is chosen as:

where $ \mathcal{D}_{1},\dots,\mathcal{D}_{7} $ are the standard propagators, and $ \mathcal{D}_8, \mathcal{D}_9 $ serve as ISPs. The generating function for the top sector is defined as:

For simplicity, we set $ G = G_{1111111} $ throughout this section. To keep the presentation focused, we explicitly display only the highest-degree and next-to-highest-degree operator structures, denoting lower-degree terms with ellipses. As our emphasis is on the top sector, we also omit subsector contributions from the displayed equations.

For the double-box topology defined above, the standard derivation via total derivatives $ \dfrac{\partial}{\partial \ell_i} \cdot \{k_j, \ell_k\} $ generates ten fundamental IBP relations. Applying these relations to the generating function G and extracting the corresponding differential operators, we obtain the primary system of equations:

Thus, the first iteration begins with two degree-one equations and eight degree-two equations. The immediate objective is to determine which lattice directions are already constrained by these relations and which irreducible subsets remain after the first pass.

We begin with the degree-two equations because they govern the mixed propagator–ISP directions that dominate the first iteration. Gaussian elimination yields:

From this reduced system, we derive four explicit T1A-type reduction rules for the composite operators $ \partial_2 \partial_9 $, $ \partial_5 \partial_9 $, $ \partial_5 \partial_8 $, and $ \partial_7 \partial_8 $. The two coupled rows yield two T2A-type reduction rules, such as those for $ \partial_2 \partial_8 $ and $ \partial_7 \partial_9 $. At this stage, the operator structure remains sufficiently simple that no additional update step is required.

Next, we combine the degree-one primary equations ((117) and (122)) with the degree-down combinations derived in the previous step ($ (118)-(119) $ and $ {(123)} - {(125)} $). Gaussian elimination yields:

Thus, four T1A-type reduction rules are derived for $ \partial_1 $, $ \partial_3 $, $ \partial_4 $, and $ \partial_6 $. These rules immediately eliminate four propagator directions from the irreducible lattice, substantially simplifying the geometry of the remaining problem.

Furthermore, these four degree-one T1A-type reduction rules can be used to construct the degree-two operators $ { \partial_4 \partial_9} $, $ { \partial_6 \partial_9} $, $ { \partial_1 \partial_8} $, and $ { \partial_3 \partial_8} $. Consequently, the two T2A-type reduction rules obtained in the previous step can be transformed into fully resolved T1A-type reduction rules:

Although the choices for intermediate T2 are not unique, the updated equations correspond to the same underlying operator relations. In this sense, the ambiguity is organizational rather than mathematical.

To summarize the first round, we derived ten reduction rules for the following differential operators:

We now interpret these first-iteration rules on the coefficient lattice. The degree-one relations in (131) reduce any integral with nonzero $ n_1 $, $ n_3 $, $ n_4 $, or $ n_6 $; consequently, all unreduced points are confined to the five-parameter subspace $ (0, n_2, 0, 0, n_5, 0, n_7, n_8, n_9) $. The degree-two rules in (132) then decouple the ISP indices from the remaining propagator indices whenever both appear together. In this sense, the first iteration reproduces, on a larger lattice, the same mixed-derivative reduction pattern previously observed in the sunset example. Applying these rules sequentially forces the indices to decouple into two distinct irreducible subsets ⁸:

The emergence of this decoupling pattern indicates that the reduction is already organizing into a small number of geometrically simple subsets. However, because the indices within $ \mathcal{U}_1 $ and $ \mathcal{U}_2 $ remain unbounded, the first iteration is incomplete. The second iteration must therefore probe the internal structure of these two residual subsets.

Before proceeding to the second iteration, we note that the remaining irreducible subsets $ \mathcal{U}_1 $ and $ \mathcal{U}_2 $ impose stringent kinematic restrictions on the operator content. Specifically, any term with a negative component in the i-th position of its operator index ($ i=1,3,4,6 $) can be omitted, as its contribution vanishes on the reduced points of $ \mathcal{U}_1 $ and $ \mathcal{U}_2 $.

Since the propagator indices $ n_1 $, $ n_3 $, $ n_4 $, and $ n_6 $ are already fully reduced, derivatives with respect to these variables yield only redundant equations. We therefore generate new equations by applying $ \partial_i $ for $ i \in \{2,5,7,8,9\} $ to the first-round reduction rules.

Applying these operators to the four previously derived degree-one reduction rules yields only trivial degree-two equations. Instead, the nontrivial information arises from applying them to the six degree-two rules and organizing the results via integrability conditions of the form $ [\partial_a, \partial_b]\mathcal{O} = 0 $. After simplification using the known descendant rules, we obtain nine linearly independent equations.

To express these compactly, we introduce two sets of zero-index operators. For the ISP parameters ($ i=8,9 $), we define:

For the remaining propagator parameters ($ i=2,5,7 $), we define:

On the subset $ \mathcal{U}_2 $, the descendant equations generated for the ISP are:

1. From $ [\partial_8,\partial_9]\partial_2 = 0 $:

2. from $ [\partial_8,\partial_9]\partial_5 = 0 $:

3. from $ [\partial_8,\partial_9]\partial_7 = 0 $:

On the subset $ \mathcal{U}_1 $, the descendant equations for the propagator are:

1. from $ [\partial_2,\partial_5]\partial_8 = 0 $:

2. From $ [\partial_2,\partial_7]\partial_8 = 0 $:

3. from $ [\partial_7,\partial_5]\partial_8 = 0 $:

4. from $ [\partial_2,\partial_5]\partial_9 = 0 $:

5. From $ [\partial_2,\partial_7]\partial_9 = 0 $:

6. from $ [\partial_7,\partial_5]\partial_9 = 0 $:

Of these nine equations, exactly six contain degree-two operators. After Gaussian elimination, we obtain

which yields three T1A-type reduction rules for $ { \partial_2 \partial_5} $, $ { \partial_2 \partial_7} $, and $ { \partial_5 \partial_7} $, one T1B-type reduction rule for $ { \partial_8 \partial_9} $, and two residual degree-one equations:

where the zero-index operators are:

We now combine the existing and newly derived degree-one equations. However, an important subtlety arises: a naive elimination of $ \partial_8 $ and $ \partial_9 $ using Eqs. (138) and (148) would yield a leading coefficient proportional to

Since this operator contains $ \eta_8\dfrac{\partial}{\partial\eta_8} $, it vanishes on the boundary $ n_8=0 $. Consequently, it cannot serve as a globally valid leading coefficient for the reduction. To circumvent this structural obstruction, we modify the elimination strategy and maintain the rules in a form that remains well-defined on the relevant lattice boundaries:

The rules in the third and fourth rows reduce $ (n_2, n_5, n_7) $ to $ (0, N_5, 0) $, where $ N_5 \leqslant n_2+n_5+n_7 $. The second row then reduces $ (n_8, n_9) $ to $ (N_8, 0) $, where $ N_8 \leqslant n_8+n_9 $. Finally, the first row reduces any point with $ n_8 \geqslant 2 $ to the cases $ n_8=0 $ or $ n_8=1 $.

Since all four new degree-one reduction rules are TB-type, we exclude them from the global update step. This preserves the simpler degree-two TA-type rules, which remain more convenient for practical reduction.

Substituting these newly derived rules into the irreducible subsets reduces $ n_2 $, $ n_7 $, and $ n_9 $ to zero. Additionally, the index $ n_8 $ is constrained to the values $ 0 $ or $ 1 $. The residual irreducible set resulting from $ \mathcal{U}_1 $ and $ \mathcal{U}_2 $ is therefore restricted to

The subset governed by $ n_5 $ retains an infinite tower of points. A third iteration is therefore required to complete the reduction.

At this stage, the residual set $ \mathcal{U}^{\text{2nd}} $ is already so constrained that most newly generated equations are either trivial or redundant. Consequently, it suffices to focus on the essential operators $ \partial_5 $ and $ \partial_8 $ applied to selected degree-one relations from the previous iteration. We present only the new relations required to eliminate the final irreducible subset.

First, we apply $ \partial_5 $ to the degree-one relation $ {(142)} $. After simplifying the result using previously established reduction rules, we obtain the following degree-two equation:

which is a degree-two T1B-type reduction rule for $ \partial_5^2 $. This allows us to rigorously reduce any integral with $ n_5 \geqslant 2 $ to $ n_5 = 0,\,1 $.

Next, we apply $ \partial_8 $ to the relation $ {(143)} $. Upon simplification, the resulting equation involves $ \partial_8 $, $ \partial_5 $, and $ \partial_2 $. By combining this with $ -\dfrac{s+t}{2 s_0}{(143)} $, the $ \partial_2 $ component cancels exactly, yielding a degree-one T1B-type equation relating $ \partial_8 $ and $ \partial_5 $:

which provides a rule to fully express $ \partial_8 $ in terms of $ \partial_5 $. Using this relation, we can systematically and unconditionally reduce $ n_8 $ to $ n_8 = 0 $.

At the conclusion of this third and final iteration, every index $ n_i $ with $ i \neq 5 $ has been reduced to zero, while the remaining propagator index $ n_5 $ is restricted to $ \{0,1\} $. Consequently, we obtain two master integrals:

which constitute the complete master-integral basis for the top sector. The planar double-box example thus terminates after three iterations, confirming that the same operator-based strategy remains effective well beyond the sunset family.

We next consider the two-loop non-planar massless double-box diagram shown in Fig. 4. This example is more complex than the planar double box, as the non-planar routing results in a less transparent organization of the irreducible lattice. Consequently, it provides a rigorous test of robustness: if the same iterative framework closes cleanly in this context, the method is not restricted to particularly simple routing structures.

Figure 4.  The two-loop non-planar double-box diagram.

As in the previous example, four external momenta $ k_1 $, $ k_2 $, $ k_3 $, and $ k_4 $ satisfy the on-shell conditions $ k_i^2=0 $ and momentum conservation $ {\displaystyle \sum}_{i=1}^4 k_i =0 $, leaving two independent kinematic scales: $ s=2 k_1 \cdot k_2 $ and $ t=2 k_2 \cdot k_3 $. The complete set of Lorentz scalars is explicitly chosen as:

where $ \mathcal{D}_1 \cdots \mathcal{D}_7 $ are the standard propagators defining the topology, and $ \mathcal{D}_8, \mathcal{D}_9 $ serve as ISPs. The generating function for the top sector is defined as:

For simplicity, we set $ G = G_{1111111} $ throughout the remainder of this section.

As the general workflow has been detailed in previous examples, we streamline the present discussion by displaying only the essential relations obtained at each iteration, along with the auxiliary equations required for the subsequent step. As before, we make only the leading and next-to-leading operator structures explicit, absorbing lower-degree terms into ellipses. Consequently, the emphasis of this section is not on reproducing every intermediate manipulation, but on demonstrating how the same iterative logic persists in a more entangled lattice geometry.

Following the same procedure as before, we begin with the ten fundamental IBP relations generated by acting with total derivatives on the generating function. The degree-two operator content of these relations is summarized by the matrix

After Gaussian elimination, we first extract the degree-one core of the system. This yields two T1A-type reduction rules for the operators.

which allow the reduction of the propagator indices $ n_1 $ and $ n_3 $ to zero.

Returning to the degree-two IBP relations, we update the operator content with the newly solved descendants of $ \partial_1 $ and $ \partial_3 $, yielding a set of T1A-type and T2A-type reduction rules. Six T1A-type rules are derived for the operators

and two T2A-type rules for the operators

Here, $ \rightarrow $ indicates that the reduction rule for $ \dfrac{ \partial}{ \partial\eta_6}\dfrac{ \partial}{ \partial\eta_9} $ is expressed in terms of $ \dfrac{ \partial}{ \partial\eta_7}\dfrac{ \partial}{ \partial\eta_9} $.

Using these first-round rules, the lattice reduction can be organized as follows:

1. Applying the degree-one relations for the operators (160) reduces all integrals with $ n_1 \geqslant 1 $ or $ n_3 \geqslant 1 $, yielding the irreducible subset $ (0,n_2, 0,n_4, n_5, n_6, n_7, n_8, n_9) $.

2. Using the degree-two operator relations (161) and (162), we completely decouple $ n_8 $ from $ n_2 $, $ n_4 $, $ n_5 $, $ n_6 $, and $ n_7 $, and decouple $ n_9 $ from $ n_2 $, $ n_4 $, and $ n_6 $, yielding $ 3 $ irreducible subsets:

The first iteration thus reduces the problem to three residual subsets with a significantly simpler structure than that of the full nine-dimensional lattice. Because each subset still contains infinitely many points, the reduction remains incomplete; however, the remaining geometry is now sufficiently simple to guide the subsequent iteration.

Since $ n_1 $ and $ n_3 $ are already fully reduced, the second iteration need only probe the remaining directions. Consequently, we apply $ \partial/\partial\eta_i $ for $ i\neq 1,3 $ to the degree-two rules established in the first iteration ⁹. After eliminating trivial descendants, the integrability conditions yield sixteen new equations, fourteen of which still contain degree-two operators. Their leading structure is summarized by

where we introduce the zero-index operator

Gaussian elimination yields ten degree-two T1A-type equations, two degree-two T2A-type equations, and two degree-one equations. From these, we derive ten T1A-type reduction rules for the operators $ { \partial_2}{ \partial_4} $, $ { \partial_2}{ \partial_5} $, $ { \partial_2}{ \partial_6} $, $ { \partial_2}{ \partial_7} $, $ { \partial_4}{ \partial_5} $, $ { \partial_4}{ \partial_7} $, $ { \partial_5}{ \partial_9} $, $ { \partial_5}{ \partial_6} $, $ { \partial_6}{ \partial_7} $, and $ { \partial_7}{ \partial_9} $ ¹⁰, along with one T2A-type reduction rule, $ \partial_5 \partial_7\rightarrow \partial_4 \partial_6 $, and one T2B-type reduction rule, $ \partial_8 \partial_9\rightarrow \partial_8 \partial_8 $. Applying these new reduction rules to simplify the two resulting degree-one equations alongside the two original ones yields the following coefficient matrix for the degree-one operators:

where the zero-index operators are defined as

One linear combination reduces to degree zero. The remaining structure yields three degree-one T2B-type rules for the operators.

These newly extracted rules further decompose the irreducible sets as follows:

1. Applying T2B-type rules to the operators $ \widehat{O}_{0;2}^{up}\dfrac{ \partial}{ \partial\eta_2} $, $ \widehat{O}_{0;6}^{up}\dfrac{ \partial}{ \partial\eta_6} $, and $ \widehat{O}_{0;7}^{up}\dfrac{ \partial}{ \partial\eta_7} $ reduces the set $ \mathcal{U}_{1}^{1st} $ to

Then, by applying the T1A-type rule, the operator $ \dfrac{ \partial}{ \partial\eta_4}\dfrac{ \partial}{ \partial\eta_5} $ can be decomposed into

2. Applying the T1A-type rule for the operator $ \dfrac{ \partial}{ \partial\eta_5}\dfrac{ \partial}{ \partial\eta_9} $ and the T2B-type rule for $ \widehat{O}_{0;7}^{up}\dfrac{ \partial}{ \partial\eta_7} $, the set $ \mathcal{U}_{2}^{1st} $ fragments into

3. Application of the T2B-type rule to reduce the operator $ 2\widehat{O}_{0;89} \dfrac{ \partial}{ \partial\eta_8} \dfrac{ \partial}{ \partial\eta_9} $ to $ \widehat{O}_{0;89} \dfrac{ \partial^2}{ \partial\eta_8^2} $ fragments the set $ \mathcal{U}_{3}^{1st} $ into

After the second iteration, every remaining irreducible subset is supported on a single nonzero index:

This already represents a substantial simplification: the coupled subsets from the first iteration have collapsed to one-dimensional towers. Nevertheless, a third iteration is required because each of these towers remains infinite.

Continuing the same procedure, and using the fact that $ n_2 $, $ n_6 $, and $ n_7 $ have been fully reduced, we apply only $ \dfrac{ \partial}{ \partial\eta_i} $ for $ i=4,5,8,9 $ to the three degree-one relations established in the second iteration. This yields four degree-two equations involving three degree-two operators:

and four valid degree-one equations involving four degree-one operators:

After eliminating trivially generated dependencies, simplifying via descendant reduction rules, and performing Gaussian elimination, we first obtain three degree-one T2B-type rules for the operators $ \dfrac{\partial}{\partial\eta_4} $, $ \dfrac{\partial}{\partial\eta_8} $, and $ \dfrac{\partial}{\partial\eta_9} $:

where the zero-index operators $ \widehat{O}_{0; 89,1} $, $ \widehat{O}_{0; 89,2} $, and $ \widehat{O}_{0; 89,3} $ are defined by

and $ \widehat{O}_{0;4}^{up}, \widehat{O}_{0;5}^{up} $ are defined in (169).

Additionally, we derive two crucial T1B-type reduction rules for the second-order derivatives associated with $ \eta_4 $ and $ \eta_5 $:

Under these final rules, the remaining irreducible sets collapse as follows:

1. Applying the T2B-type rules to the operators $ \dfrac{\partial}{\partial\eta_5} $, $ \dfrac{\partial}{\partial\eta_8} $, and $ \dfrac{\partial}{\partial\eta_9} $ unconditionally reduces the indices $ n_5, n_8, $ and $ n_9 $ to zero. Consequently, the infinite lines corresponding to the sets $ \mathcal{U}_{2}^{2nd} $, $ \mathcal{U}_{3}^{2nd} $, and $ \mathcal{U}_{4}^{2nd} $ are completely reduced to $ \mathcal{U}_{1}^{2nd} $.

2. Applying the T1B-type rule for the reduction operator $ \dfrac{ \partial^2}{ \partial\eta_4^2} $ ultimately restricts the remaining subset $ \mathcal{U}_{1}^{2nd} $ to $ n_4 \in \{0, 1\} $.

At this stage, the infinite towers have been fully eliminated, leaving only a finite set of master integrals:

This completes the reduction for the non-planar double box. Although the intermediate operator structure is more complex and the decoupling pattern less transparent than in the planar case, the same iterative logic drives the irreducible lattice to a finite master-integral basis.

The preceding examples focused on top sectors. To demonstrate that the formalism is not limited to this context, we now consider the subsector $ G_{110} $ of the sunset family. This example is conceptually important because it illustrates how the same operator-based strategy applies after part of the propagator structure has been removed. Using the kinematics introduced in Section IV, we again extract the reduction rules directly and track their action on the remaining lattice.

From the six fundamental IBP relations, we derive the coefficient matrix for four equations of degree two:

and the coefficient matrix of two degree-one equations:

Diagonalizing the degree-one system yields two explicit T1A-type rules for the operators $ \partial/\partial\eta_1 $ and $ \partial/\partial\eta_2 $:

Applying these T1-type rules unconditionally reduces $ n_1 $ and $ n_2 $ to zero, leaving only the irreducible set.

At this stage, it is more efficient to update the remaining degree-two relations of the original IBP system than to generate numerous new descendants. Substituting the newly solved degree-one relations, Eqs. (189) and (190), reduces the problem to a degree-one system with a coefficient matrix

Applying Gaussian elimination to this updated system yields four new degree-one equations. From these equations, we derive three key reduction rules for the operators $ \dfrac{\partial}{\partial\eta_3} $, $ \dfrac{\partial}{\partial\eta_4} $, and $ \dfrac{\partial}{\partial\eta_5} $:

where we introduce the respective zero-index coefficient operators:

Using these updated equations, the remaining irreducible set is resolved in a straightforward hierarchical manner:

1. We first apply the rule for $ \widehat{O}_{s;34,1} \dfrac{ \partial}{ \partial\eta_3} $ (Eq. (193)), which has the highest priority. Although the right-hand side contains an operator with a negative index component, namely $ \eta_3 \dfrac{ \partial^2}{ \partial\eta_5^2} $, this poses no difficulty for the reduction of $ n_3 $. The $ \eta_3 $ factor implies that the corresponding contribution begins only at $ n_3=1 $, so it does not impede the reduction from $ n_3=1 $ to $ n_3=0 $. The set $ \mathcal{U}_1 $ is therefore reduced to

2. When $ n_3 = 0 $, every term containing an explicit factor of $ \eta_3 $ vanishes identically. Consequently, the right-hand sides of Eq. (194) and Eq. (195) vanish, and these two relations reduce to T1A-type rules for $ \dfrac{\partial}{\partial\eta_4} $ and $ \dfrac{\partial}{\partial\eta_5} $. Applying these rules reduces the residual set to a single master integral with all indices equal to zero:

In summary, the sunset subsector closes without requiring a separate descendant-equation analysis beyond the updated IBP system. This example complements the top-sector discussion by demonstrating that the same framework naturally handles lower sectors and consistently yields a finite master-integral basis—in this case, a single master integral.

We conclude the sequence of worked examples with a nontrivial subsector of the two-loop massless non-planar box. Employing the kinematics introduced in Section 6, we consider the subsector in which $ \mathcal{D}_2 $ is treated as an ISP. This example is not redundant with respect to the top non-planar sector; rather, it demonstrates that the method remains applicable when the sector structure changes while the non-planar geometry is preserved. To maintain a focused presentation, we omit lengthy intermediate matrices that do not contribute conceptually and focus on the key reduction rules.

The original ten fundamental IBP relations yield the coefficient matrix of the degree-two equations as:

After Gaussian elimination, we derive four degree-two T1A-type reduction rules for the operators $ \partial_i \partial_8 $ ($ i=4,5,6,7 $), along with two explicit degree-one T1A-type scaling rules for the operators

Applying these T1-type rules reduces $ n_1 $ and $ n_3 $ to zero.

These initial rules yield an important strategic consequence: once $ \partial/\partial\eta_1 $ and $ \partial/\partial\eta_3 $ are resolved, the remaining degree-two relations can be updated to reveal a new degree-one reduction rule. In particular, updating the equations involving $ \partial/\partial\eta_1\,\partial/\partial\eta_8 $ and $ \partial/\partial\eta_3\,\partial/\partial\eta_8 $, and subsequently performing Gaussian elimination, yields the key T1B-type rule.

where the zero-index operator is given by:

Applying this rule reduces the index $ n_8 $ to zero.

Furthermore, we present two T2A-type reduction rules for the degree-two operators $ \dfrac{\partial}{\partial\eta_2}\dfrac{\partial}{\partial\eta_6} $ and $ \dfrac{\partial}{\partial\eta_4}\dfrac{\partial}{\partial\eta_9} $, namely

Applying these rules decouples the index $ n_2 $ from $ n_6 $ and $ n_4 $ from $ n_9 $, thereby fragmenting the irreducible set as follows:

1. Applying degree-one rules to the operators $ \dfrac{ \partial}{ \partial\eta_1} $, $ \dfrac{ \partial}{ \partial\eta_3} $, and $ \widehat{O}_{0;a,8} \dfrac{ \partial}{ \partial\eta_8} $, we obtain an irreducible set:

2. Using the degree-two rules for the operators $ \dfrac{\partial}{\partial\eta_2}\dfrac{\partial}{\partial\eta_6} $ and $ \dfrac{\partial}{\partial\eta_4}\dfrac{\partial}{\partial\eta_9} $, we decompose the irreducible set into four irreducible subsets:

Since infinitely many lattice points remain unreduced, a second iteration is required.

After the first iteration, the remaining unreduced set is spanned by $ n_2 $, $ n_4 $, $ n_5 $, $ n_6 $, $ n_7 $, and $ n_9 $. The second iteration therefore considers only these directions. Following the general algorithm, we apply the remaining derivatives $ \partial/\partial\eta_i $ for $ i=2,4,5,6,7,9 $ to the equations established in the first step. Particular attention is given to the degree-zero relation, as it can generate structurally novel degree-one constraints.

After updating the newly generated equations with the previously established reduction rules, we obtain twelve degree-one equations and seven degree-two equations. To achieve full reduction, the degree-one equations suffice; their coefficient matrix is:

where the zero-index operators are defined as

and the operators $ \widehat{O}_{0;i}^{up} $ are defined in (169).

Through Gaussian elimination and simplification using the existing rules, we derive four T1B-type rules for the degree-one operators:

Additionally, we extract two T2B-type rules, namely

where the previously defined operator $ \widehat{O}_{0;8}^{a} $ is applied because we employ the T1B-type reduction rule (201) to update the operator $ \eta_9 \dfrac{ \partial}{ \partial \eta_2} \dfrac{ \partial}{ \partial \eta_8} $.

These newly derived rules reduce the irreducible sets as follows:

1. By applying the four T1B-type rules, we reduce $ \mathcal{U}^{1st}_1 $ to $ (0, n_2, 0, 0, 0, 0, 0, 0, 0) $, $ \mathcal{U}^{1st}_2 $ to $ (0, 0, 0, 0, 0, 0, 0, 0, 0) $, and both $ \mathcal{U}^{1st}_3 $ and $ \mathcal{U}^{1st}_4 $ to $ (0, 0, 0, 0, 0, 0, 0, 0, n_9) $.

2. Applying the reduction rule for $ \widehat{O}_{0;8}^{a} \widehat{O}_{0;2}^{b} \dfrac{\partial}{\partial\eta_2} $ to $ (0, n_2, 0, 0, 0, 0, 0, 0, 0) $ maps this set to $ (0, 0, 0, 0, 0, 0, 0, n_8, n_9) $, which already contains the subset $ (0, 0, 0, 0, 0, 0, 0, 0, n_9) $.

3. At the bound $ n_2 = 0 $, the reduction rule for $ \widehat{O}_{0;9}^{b} \dfrac{\partial}{\partial\eta_9} $ simplifies to

Thus, together with the reduction rule (201), the set $ (0, 0, 0, 0, 0, 0, 0, n_8, n_9) $ can be reduced to $ (0, 0, 0, 0, 0, 0, 0, 0, 0) $.

Following the second iteration, we have derived a complete set of reduction rules for each $ \dfrac{\partial}{\partial\eta_i} $. Consequently, any integral in this subsector reduces to the unique master integral,

In this section, we examine a qualitatively different scenario: a family whose top sector is degenerate, meaning it contains no master integrals and reduces completely to subsectors. This example is significant because it demonstrates that the formalism extends beyond the construction of finite master-integral bases. Furthermore, it can identify cases where the top sector itself is absent as an independent component of the reduction problem. We consider the topology defined by the integral family

with the conventional two-loop configuration:

Since $ L=2 $ and $ E=1 $, the number of complete basis invariants is exactly $ N=5 $. Consequently, there are no Irreducible Numerator Parts (ISPs). The generating function for the top sector is explicitly given by:

To concisely describe the reduction to lower topologies, we introduce a compact notation for the subsector generating functions. For example, the boundary function corresponding to the contraction of the first propagator ($ \eta_1=0 $) is denoted by:

with analogous definitions for $ B_2 $, $ B_3 $, $ B_4 $, and $ B_5 $.

Applying the standard IBP differential operators, we obtain six fundamental relations. In the absence of ISPs, this initial system consists solely of degree-one operators acting on G, offset by subsector terms:

Consequently, there are six equations but only four independent degree-one operators governing the top-sector system, namely $ \{\dfrac{ \partial}{ \partial\eta_1}, \dfrac{ \partial}{ \partial\eta_2}, \dfrac{ \partial}{ \partial\eta_3}, \dfrac{ \partial}{ \partial\eta_4}\} $. Diagonalizing the corresponding coefficient matrix yields four T1A-type reduction rules, reducing any top-sector integral to the subset $ (0,0,0,0,n_5) $. More importantly, the remaining two independent linear combinations contain no degree-one terms, yielding pure degree-zero constraints. These degree-zero equations serve as a key indicator that the top sector may be degenerate.

Specifically, the combination of Eqs. (229) and (230) evaluates to:

while the combination of (233)–(232) yields:

These degree-zero constraints constitute the defining feature of this degenerate topology.

To exploit these constraints, we apply the four T1A-type reduction rules to the degree-zero equations. On the restricted subset $ (0, 0, 0, 0, n_5) $, Eq. (235) simplifies to

and Eq. (236) simplifies to

We now observe that both (237) and (238) are T1-type equations for the zero-index operator.

Using this operator, we reduce all remaining top-sector integrals to subsector contributions. This result is therefore stronger than ordinary completeness: the top sector itself contains no master integrals.

Furthermore, subtracting (237) from (238) yields a nontrivial relation among subsectors, indicating that certain subsectors are not independent. We do not pursue this point further here. This example highlights a further strength of the framework: the same operator-based analysis can certify the elimination of a top sector, rather than merely reducing it to a finite basis.

In this paper, we reformulated the reduction of multi-loop Feynman integrals in terms of generating functions and differential operators. Within this framework, integration-by-parts identities are recast as differential equations, while symbolic reduction rules emerge as operator identities acting on the lattice of expansion coefficients. This perspective enables the treatment of reduction not merely as a computational task, but also as a structured algebraic problem.

Building upon this formalism, we introduced an iterative algorithm that generates differential equations, simplifies them using previously derived rules, extracts new reduction rules by solving linear systems, and verifies completeness through the geometry of irreducible lattice subsets. The worked examples demonstrate that this logic applies across qualitatively distinct scenarios: simple top sectors, higher-point two-loop topologies, subsectors, and a degenerate case in which the top sector is entirely eliminated.

Several aspects of the present work are intentionally conservative. We have prioritized the structural formulation of the method and explicit examples over implementation-level optimization or exhaustive comparisons with alternative reduction strategies. These areas represent natural directions for future research. In particular, it would be valuable to investigate automated operator-selection strategies, a more systematic incorporation of symmetries, and the application of this framework to broader classes of multi-loop families.

The examples presented here do not establish a universal scaling advantage for higher-loop families over families with more external legs. In fact, both the sunset and double-box examples are two-loop systems. At a fixed loop order, increasing the number of independent external momenta expands the scalar product basis $N=L(L+1)/2+LE$; however, increasing the loop order also expands this basis. In practice, the computational cost also depends on the number of sectors, the rank and sparsity of the operator systems, the number of master integrals, coefficient growth, and the number of descendant rounds. A quantitative scaling comparison therefore necessitates a broader set of examples and a common implementation benchmark.

Overall, the generating function viewpoint provides a unified language in which symbolic reduction rules, descendant equations, and completeness conditions can be treated on an equal footing. We anticipate that this formulation will prove useful both for conceptual studies of IBP reduction and for the development of practical symbolic reduction algorithms.

We thank Mingxing Luo, Antonela Matijašić, and Stefan Weinzierl for valuable discussions.