Differential equations on unitarity cut surfaces

TL;DR

This paper introduces a novel approach to deriving differential equations for Feynman integrals that avoids doubled propagators and simplifies the process, especially for complex multi-loop integrals, using unitarity-inspired methods and computational techniques.

Contribution

It presents a new method to formulate differential equations on unitarity cut surfaces, eliminating the need for IBP reduction in certain cases, and extends the approach to two-loop nonplanar five-point integrals.

Findings

Direct derivation of DEs for one-loop box without IBP reduction.

Application to two-loop nonplanar five-point integrals at maximal cut.

Use of finite field techniques for computational efficiency.

Abstract

We reformulate differential equations (DEs) for Feynman integrals to avoid doubled propagators in intermediate steps. External momentum derivatives are dressed with loop momentum derivatives to form tangent vectors to unitarity cut surfaces, in a way inspired by unitarity-compatible IBP reduction. For the one-loop box, our method directly produces the final DEs without any integration-by-parts reduction. We further illustrate the method by deriving maximal-cut level differential equations for two-loop nonplanar five-point integrals, whose exact expressions are yet unknown. We speed up the computation using finite field techniques and rational function reconstruction.

Tables

Table 1. Table 1: Time required to obtain the results Eq. ( 84 ) with full dependence on the 8 8 8 master integrals including daughter topologies. The computation is performed on 1 1 1 CPU core on a laptop computer with an Intel Core i5 processor (clock frequency 2.5 GHz). Finite field techniques and rational function reconstruction are not used for the timing comparison here, but will be used for the more complicated nonplanar five-point integrals.

Software	Time taken by final run
FIRE5	141 seconds
Own code	37 seconds

Equations

$$\sum _i{\beta }_i^{\mu }\frac{\partial }{\partial p_i^{\mu }}\int d^{d}l\frac{\mathcal{N}}{{\prod }_j{z}_j},$$

$${\beta }_i^{\mu }={\beta }_{ij}{p}_j^{\mu }.$$

$$\int d^{d}l\left[\sum _i{\beta }_i^{\mu }\frac{\partial }{\partial p_i^{\mu }}\frac{\mathcal{N}}{{\prod }_j{z}_j}+\frac{\partial }{\partial l^{\mu }}\frac{\mathcal{N}{v}^{\mu }}{{\prod }_j{z}_j}\right]$$

$$=\int d^{d}l\sum _i\left({\beta }_i^{\mu }\frac{\partial }{\partial p_i^{\mu }}+v^{\mu }\frac{\partial }{\partial l^{\mu }}\right)\frac{\mathcal{N}}{{\prod }_j{z}_j}+\int d^{d}l\frac{\mathcal{N}}{{\prod }_j{z}_j}\frac{\partial v^{\mu }}{\partial l^{\mu }}.$$

$$\left({\beta }_i^{\mu }\frac{\partial }{\partial p_i^{\mu }}+v^{\mu }\frac{\partial }{\partial l^{\mu }}\right)z_j=f_j{z}_j,$$

$$z_i=0,\forall i\in \Delta .$$

$$\left({\beta }_i^{\mu }\frac{\partial }{\partial p_i^{\mu }}+v^{\mu }\frac{\partial }{\partial l^{\mu }}\right)z_k=f_k{z}_k=0.$$

$${\beta }_i^{\mu }\frac{\partial }{\partial p_i^{\mu }}+v^{\mu }\frac{\partial }{\partial l^{\mu }},$$

$$v^{\mu }\frac{\partial }{\partial l^{\mu }},$$

$${\mathcal{I}}_{\mathrm{box}}=\int d^{d}l\frac{\mathcal{N}}{{\prod }_{j=1}^4{z}_j},$$

$$z_1$$

$$z_2$$

$$z_3$$

$$z_4$$

$$(p_1+p_2{)}^2$$

$$(p_2+p_3{)}^2$$

$$(p_1+p_3{)}^2$$

$${\eta }^{\mu \nu }=\sum _{i=1,2,3}\sum _{j=1,2,3}{(G^{-1})}_{ij}{p}_i^{\mu }{p}_j^{\nu }+{\hat{\eta }}^{\mu \nu },$$

$$G_{ij}=p_i\cdot p_j=\frac{1}{2}\left(\begin{array}{ccc}0& s& -(s+t)\\ s& 0& t\\ -(s+t)& t& 0\end{array}\right).$$

$$l\cdot p_1$$

$$l\cdot p_2$$

$$l\cdot p_3$$

$$z_{ij}=z_i-z_j.$$

$$z_1={\eta }^{\mu \nu }{l}_{\mu }{l}_{\nu }={(G^{-1})}_{ij}(l\cdot p_i)(l\cdot p_j)-{\hat{\mu }}^2,$$

$${\hat{\mu }}^2=-{\hat{\eta }}^{\mu \nu }{l}_{\mu }{l}_{\nu }.$$

$${\hat{\mu }}^2=F(z)={(G^{-1})}_{ij}(l\cdot p_i)(l\cdot p_j)-z_1,$$

$$F$$

$$\displaystyle\quad+2st(-z_{1}z_{2}-z_{2}z_{3}-z_{3}z_{4}-z_{4}z_{1}+2z_{1}z_{3}+2z_{2}z_{4})\big{]}\,.$$

$$\Omega =\frac{d^{d}l}{{\prod }_j{z}_j}$$

$$=\frac{d{z}_{12}d{z}_{23}d{z}_{34}}{2^3\sqrt{\det G}\cdot {\prod }_j{z}_j}\cdot \frac{{\pi }^{(d-3)/2}}{\Gamma ((d-3)/2)}{({\hat{\mu }}^2)}^{(d-5)/2}d{\hat{\mu }}^2$$

$$=\frac{{\pi }^{(d-3)/2}}{2^3\Gamma ((d-3)/2)\sqrt{\det G}\cdot {\prod }_j{z}_j}d{z}_{12}d{z}_{23}d{z}_{34}d{\hat{\mu }}^2({\hat{\mu }}^2{)}^{(d-5)/2}$$

$$\times d{z}_1\delta \left(z_1-{(G^{-1})}_{ij}(l\cdot p_i)(l\cdot p_j)-{\hat{\mu }}^2\right)$$

$$=\frac{{\pi }^{(d-3)/2}}{2^3\Gamma ((d-3)/2)\sqrt{\det G}\cdot {\prod }_j{z}_j}d{z}_{12}d{z}_{23}d{z}_{34}d{z}_1{F(z)}^{(d-5)/2}$$

$$=\frac{{\pi }^{(d-3)/2}}{2^3\Gamma ((d-3)/2)\sqrt{\det G}}\frac{d{z}_1}{z_1}\frac{d{z}_2}{z_2}\frac{d{z}_3}{z_3}\frac{d{z}_4}{z_4}{F(z)}^{(d-5)/2}.$$

$$\int d^{d}l\frac{1}{{\prod }_j{z}_j}=\int \Omega ,$$

$$\int d^{d}l\frac{\partial }{\partial l^{\mu }}\frac{v^{\mu }}{{\prod }_j{z}_j}=\int d\left(v\circ \Omega \right).$$

$$\sum _{\mu }{v}^{\mu }{\partial }_{\mu }=\sum _i{v}_i\frac{\partial }{\partial z_i}.$$

Full text

††preprint: UCLA-17-TEP-102††arxiv: 1702.02355

1 Introduction

The generalized unitarity method Bern:1994cg ; Bern:1995db ; Bern:1996ja ; Bern:1997sc ; Britto:2004nc ; Ossola:2006us ; Forde:2007mi ; Ellis:2011cr ; Ita:2011hi has been very successfully applied to constructing loop integrands in quantum field theories, while applications to loop integration remain frontiers to explore. Integration-by-parts (IBP) reduction Chetyrkin:1981qh ; Laporta:2001dd ; Laporta:1996mq ; Smirnov:2014hma ; vonManteuffel:2012np was recently reformulated Gluza:2010ws ; Schabinger:2011dz ; Chen:2015lyz ; Sogaard:2014jla ; Georgoudis:2015hca ; Ita:2015tya ; Larsen:2015ped ; Georgoudis:2016wff ; Ita:2016oar ; Zhang:2016kfo in a unitarity-compatible manner without doubled propagators, relying on special combinations of loop momentum derivatives which form tangent vectors Ita:2015tya to unitarity cut surfaces.

After IBP reduction to a basis of master integrals, the master integrals still need to be evaluated, e.g. using the method of differential equations (DEs) Kotikov:1990kg ; Bern:1993kr ; Remiddi:1997ny ; Gehrmann:1999as ; Argeri:2007up ; Henn:2013nsa . A recent breakthrough was Henn’s canonical form of DEs Henn:2013pwa ; Henn:2014qga , allowing a large class of loop integrals to be expressed in terms of iterated integrals of uniform transcendentality. Various algorithms and software packages Lee:2014ioa ; Henn:2014qga ; Tancredi:2015pta ; Meyer:2016slj ; Meyer:2017joq ; Prausa:2017ltv ; Gituliar:2017vzm have appeared to find algebraic transformations of DEs to the canonical form, while a complementary approach is finding master integrals in the $d\log$ form with unit leading singularities ArkaniHamed:2010gh ; Henn:2014qga ; Bern:2015ple .

In constructing DEs, an important intermediate step is IBP reduction, which brings the RHS of the DEs into a linear combination of master integrals. External momentum derivatives increase the power of propagator denominators, so unitarity-compatible IBP reduction is not directly applicable. To solve this problem, one approach is to decrease the power of propagator denominators using dimension shifting Georgoudis:2016wff . However, we propose an alternative approach that completely avoids doubled propagators, even in intermediate steps. We promote unitarity cut surfaces to be objects embedded in the space of not only loop momenta, but also external momenta. By combining external and loop momentum derivatives to form tangent vectors to unitarity cut surfaces, doubled propagators cancel out, in direct analogy with unitarity-compatible IBP reduction.

It has been proposed that the maximal cut can provide valuable information about differential equations and the function space of the integrals CaronHuot:2012ab ; Primo:2016ebd ; Frellesvig:2017aai . The latter two references rely on consistent definitions of unitarity cuts in the presence doubled propagators (see also Lee:2012te ; Sogaard:2014ila ). This issue is bypassed in our approach by construction.

Section 2 sketches the basics of our formalism. Section 3 uses inverse propagator coordinates, also known as the Baikov representation in the $d$-dimensional case, to present the detailed formalism through the one-loop box example. Section 4 applies the formalism to the nonplanar pentabox at the maximal cut level, and finds a system proportional to the dimensional regularization parameter $\epsilon$, for tensor integrals with unit leading singularities. Section 5 discusses the use of finite field techniques and rational function construction in speeding up the nonplanar pentabox computation. Some concluding remarks are given in Section 6.

2 Basic formalism

2.1 Avoiding doubled propagators

Generalized unitarity cuts replace propagators by delta functions. However, when the propagator is doubled (i.e. squared), it is no longer straightforward to impose unitarity cuts Primo:2016ebd ; Frellesvig:2017aai ; Lee:2012te ; Sogaard:2014ila . Therefore, a unitarity-compatible approach to integration-by-parts reduction uses special IBP relations that do not involve doubled propagators Gluza:2010ws ; Schabinger:2011dz ; Chen:2015lyz ; Sogaard:2014jla ; Georgoudis:2015hca ; Ita:2015tya ; Larsen:2015ped ; Georgoudis:2016wff ; Ita:2016oar ; Zhang:2016kfo . This allows IBP relations to be put on unitarity cuts, which can be exploited to construct multi-loop generalizations of the OPP parameterization Ossola:2006us of one-loop integrands, as well as allowing analytic IBP reduction to be achieved by merging results from a spanning set of unitarity cuts Larsen:2015ped .

We will explore a similar unitarity-compatible approach to differential equations, with no Feynman integrals involving doubled propagators appearing on the LHS or RHS of the differential equations, even before IBP reduction is performed to simplify the RHS. The advantage is two-fold. First, the differential equations may be put on unitarity cuts, and therefore can be constructed by merging incomplete results on a spanning set of unitarity cuts. Second, unitarity-compatible IBP reduction can be used to reduce the RHS into the original set of master integrals, since no integrals with doubled propagators are ever generated.

We wish to compute the derivative of a Feynman integral with propagators $1/z_{j}$ and a tensor numerator $\mathcal{N}$,

[TABLE]

where $\beta^{\mu}_{i}$ is a linear combination of external momenta $p_{j}$,

[TABLE]

We are free to add total divergences to Eq. (1) without changing its value after integration, obtaining

[TABLE]

In the above expressions, $v^{\mu}=v^{\mu}(l,p)$ has polynomial dependence on internal and external momenta, with one free Lorentz index. The final expression Eq. (4) has no doubled propagators $\sim 1/z_{j}^{2}$, if the following condition is satisfied,

[TABLE]

where $f_{j}$ has polynomial dependence on internal and external momenta.

2.2 Relation to unitarity cut surfaces

As we will see, the condition Eq. (5) has a nice geometric interpretation in terms of unitarity cut surfaces, similar to what was observed Ita:2015tya in unitarity-compatible IBP reduction.

Consider a $L$-loop Feynman diagram topology with $N$ internal propagators. For a subset $\Delta$ of all inverse propagators ${1,2,\dots,N}$, the unitarity cut surface labeled by $\Delta$ is defined as the hypersurface of all points $(l_{1},l_{2},\dots,l_{L})$ in the complex loop momentum space which solves the generalized unitarity cut condition,

[TABLE]

Notice that generalized unitarity cuts differ from traditional unitarity cuts in QFT textbooks, treated by e.g. the optical theorem, Cutkosky rules Cutkosky:1960sp , and the largest time equation Veltman:1963th ; Remiddi:1981hn , in several respects,

Complex rather than real loop momentum space is considered. In particular, the energy component of the cut loop momentum is not required to be (real) positive. The algebraic closure of the complex field guarantees that the unitarity cut condition has solutions for generic external momenta, as long as not too many propagators are cut (which will be assumed to be the case for the rest of the paper).111For example, for one-loop topologies in $4$ spacetime dimensions, the unitarity cut condition can always be solved when $4$ or fewer propagators are cut. 2. 2.

There is no direct connection with the discontinuity of Feynman diagrams across branch cuts.

When $\Delta$ contains all inverse propagators, i.e. when the unitarity cut condition sets all the propagators on-shell, the unitarity cut is called a maximal cut. Many problems, such as integrand construction and integration by parts, are simplest at the level of the maximal cut. When a proper subset of propagators are set on-shell, we have a non-maximal cut.

In general, for a diagram topology with $N$ propagators, there are $2^{N}$ different unitarity cut surfaces, since we may choose each propagator to be either cut or not cut. Since this paper is concerned with differential equations w.r.t. external momenta, we define “extended unitarity cut surfaces” which is embedded in the space of not only loop momenta but also external momenta. Here the loop momenta are still allowed to be complex but external momenta are required to be real. As usual, the surface is defined by the cut condition, Eq. (6).

For any of the $2^{N}$ extended unitarity cut surfaces $\Delta$ with $z_{k}\in\Delta$, the cut condition sets $z_{k}=0$, while the condition for the absence of doubled propagators, Eq. (5), becomes

[TABLE]

Therefore, the expression

[TABLE]

which we refer to as a “DE vector”, is a tangent vector to every extended unitarity cut surface embedded in the space of both internal and external momenta. The loop part,

[TABLE]

is called an “IBP vector”. In the special case that the DE vector has no external momentum derivative, the IBP vector itself is a tangent vector to unitarity cut surfaces, and is used in unitarity-compatible IBP reduction. Computational algebraic geometry is used to find IBP vectors in the literature, and will also be used in this paper to find DE vectors.

3 Detailed formalism in inverse propagator coordinates

3.1 Inverse propagator coordinates

We re-examine the simple example of the one-loop box Bern:1993kr using our method. We assume that all internal and external lines are massless. The scalar box integral with some numerator $\mathcal{N}$, shown in the leftmost diagram of Fig. 1, is

[TABLE]

with

[TABLE]

The kinematic invariants are

[TABLE]

It is well known that after IBP reduction, there are $3$ master integrals. For our purposes, they are conveniently chosen as the scalar box $\mathcal{I}_{\rm box}$, the $s$-channel triangle $\mathcal{I}_{\rm triangle}^{(s)}$, and the $t$-channel triangle integral $\mathcal{I}_{\rm triangle}^{(t)}$, shown in Fig. 1.222Although there are $4$ possible daughter triangles of the box, the two $s$-channel triangles are identical after integration, and the same is true for the two $t$-channel triangles.

The loop integral can be parameterized in the inverse propagator coordinates, in either $4$ or $d$ dimensions. This parameterization goes back to Cutkosky’s proof of the cutting rules Cutkosky:1960sp , and has been systematically studied by Baikov for the $d$-dimensional case Baikov:1996rk ; Baikov:1996iu ; Grozin:2011mt . More recently, this parameterization was used in Ita:2015tya ; Larsen:2015ped for unitarity-compatible IBP reduction. A detailed explanation of the Baikov representation recently appeared in Frellesvig:2017aai , where a public code was made available, and applications to differential equations were discussed. We now derive this parameterization for the one-loop box, and refer the readers to the literature for the multi-loop case.

Using the Van Neerven-Vermaseren basis vanNeerven:1983vr , the metric tensor $\eta^{\mu\nu}$ is written as the sum of a “physical” component in the 3-dimensional space spanned by external momenta, and a “transverse” component $\hat{\eta}^{\mu\nu}$ in the remaining $(d-3)$-dimensional space which is orthogonal to every external momentum,

[TABLE]

where we used the inverse of the Gram matrix $G$,

[TABLE]

We will need the relations,

[TABLE]

In the above equations we have defined

[TABLE]

From Eq. (11),

[TABLE]

where $\hat{\mu}^{2}$ denotes the Euclidean norm of the $(d-3)$-dimensional transverse components of $l$,

[TABLE]

Substituting Eqs. (20)-(22) into the above equation, $\hat{\mu}^{2}$ is expressed in terms of the inverse propagators $z_{i}$,

[TABLE]

where $F$ is called the Baikov polynomial. For the one-loop box, it evaluates to

[TABLE]

The loop integration measure, multiplied by the propagators, is re-written as

[TABLE]

In the second-to-last line of the above equations, we integrated out $\hat{\mu}^{2}$ against the delta function, and in the last line we made a linear transformation of integration variables with unit Jacobian. Eq. (28) accomplishes the transformation of the loop momentum from the Lorentzian / Euclidean component coordinates to the inverse propagator coordinates.

We will adopt the differential form notation for IBP relations Ita:2015tya ; Larsen:2015ped , which helps to transform between different coordinate systems for the loop momenta. We write

[TABLE]

where $\Omega$ is a maximal differential form in the loop momentum space. For convenience, we define $\Omega$ to include both the integration measure and the propagators. The total divergence integral of any IBP vector, $v=v^{\mu}\partial_{\mu}$, is

[TABLE]

3.2 Vectors for differential equations

We first discuss IBP vectors, i.e. the loop component of DE vectors. The IBP vector Eq. (9) can be transformed into the new parameterization,

[TABLE]

Since $v^{\mu}$ is required to have polynomial dependence on external and internal momenta, it is easy to show that in the new coordinates, $v_{i}$ must have polynomial dependence on the inverse propagators $z$. Furthermore, $v^{\mu}\partial_{\mu}$ is required to satisfy rotational invariance in the $(d-3)$-dimensional transverse space. An example of such an expression is

[TABLE]

with $\eta^{\mu\nu}$ defined in Eq. (18). This vector, and in fact all rotational invariant vectors satisfy

[TABLE]

where the proportionality constant has polynomial dependence on $z_{i}$. In the inverse propagator coordinates, this means

[TABLE]

for some $v_{F}$ which is a polynomial in the $z$ variable.333We thank Harald Ita for giving a two-loop version of this argument in private communications. The above equation is a crucial criterion for a valid IBP vector, derived in a different way in Larsen:2015ped .

The IBP relation from an IBP vector $v=v_{i}\partial_{z_{i}}$ is obtained by substituting the last line of Eq. (28) into the exterior derivative expression Eq. (30),

[TABLE]

with $v_{F}$ defined in Eq. (34).

Setting the tensor numerator $\mathcal{N}$ to $1$ in Eq. (3), the 2nd term in the square bracket is given by Eq. (35), while the first term is simply, suppressing the $i$ index,

[TABLE]

Adding (36) and the vanishing integral Eq. (35), we obtain,

[TABLE]

The above expression has no doubled propagators when

[TABLE]

for some polynomial $f_{k}=f_{k}(z)$ for every inverse propagator $z_{k}$. This is equivalent to Eq. (7).

In summary, in inverse propagator coordinates, a unitarity-compatible DE vector that does not lead to doubled propagators is required to satisfy two requirements, Eq. (34) and (40). Combining these two equations gives

[TABLE]

Given the derivative against external momenta, $\beta^{\mu}\partial_{\mu}$, Eq. (42) is a polynomial equation in the unknown polynomials $f_{j}$ and $v_{F}$. This is an inhomogeneous version of the “syzygy equation” of Ref. Larsen:2015ped .444We note that Ref. Frellesvig:2017aai used a different, but related, polynomial equation to construct differential equations within the Baikov representation, while our hybrid approach involves objects such as $\partial z_{j}/\partial p^{\mu}$, which are in turn expressed as polynomials in the $z$ variables. We use the computational algebraic geometry package SINGULAR DGPS to find a particular solution for $f_{j}$ and $v_{F}$ in Eq. (42), which in turn fixes $v_{k}$ in Eq. (40). This allows us to compute the derivative of the loop integral w.r.t. external momenta using Eq. (38).

3.3 Evaluating DEs: one-loop box warm-up

Consider the derivative with respect to $t=(p_{2}+p_{3})^{2}$,

[TABLE]

which annihilates $s=(p_{1}+p_{2})^{2}$ and keeps external momenta on-shell. Combined with Eqs. (11)-(14), we find the following expressions for the RHS of Eq. (42),

[TABLE]

Solving Eq. (42) for the unknown polynomials $f_{i}$ and $v_{F}$, we find the following solution, which can be easily checked,

[TABLE]

which in turn fixes $v_{k}$ in Eq. (40). Substituting the results into Eq. (38) gives the $t$-derivative of the box integral in a form that has no doubled propagators,

[TABLE]

where we set $d=4-2\epsilon$. The triangle integrals are one-scale integrals whose derivatives against $t$ can be fixed by simple dimensional analysis, so we omit the calculation. Transforming to the basis suggested in Refs. Henn:2014qga ; Henn:2013pwa , we obtain

[TABLE]

which agrees with the result obtained from the standard approach involving IBP reduction of integrals with doubled propagators. In this simple example, we directly obtain the DEs without any IBP reduction, but unitarity-compatible IBP reduction of tensor integrals is needed in more complicated cases.

3.4 4D leading singularities

Although this paper mainly concerns unitarity cuts in $d$ dimensions, the inverse propagator coordinates also make it easy to compute cuts in $4$ dimensions. Setting $d=4$ in Eq. (28), the maximal cut residue of the box integral $\mathcal{I}_{\rm box}$ is equal to

[TABLE]

where $F(0)=F(z_{i}=0)$. Therefore, $st\,\mathcal{I}_{\rm box}$ is an integral whose leading singularity is a $\mathbb{Q}$ number with no dependence on $s$ or $t$. This is an important criterion for selecting candidate integrals with uniform transcendentality.

We now translate the pentagon integrals with unit leading singularities ArkaniHamed:2010gh to the inverse propagator coordinates, which serves as a building block for the master integral choice for the nonplanar pentabox. This loop topology is shown in Fig. 2. The five kinematic invariants can be chosen, in a cyclic invariant fashion, to be $s_{12}$, $s_{23}$, $s_{34}$, $s_{45}$, $s_{51}$, where $s_{ij}=(p_{i}+p_{j})^{2}$. The Gram matrix is again defined as $G_{ij}=p_{i}\cdot p_{j},1\leq i,j\leq 4$. In $D$ dimensions, the pentagon integral with a numerator $\mathcal{N}$ can be expressed in terms of the five inverse propagator coordinates, $z_{i}=(l-\sum_{j=1}^{i-1}p_{j})^{2}$,

[TABLE]

whereas in $4$ dimensions, the Baikov polynomial resides in a Dirac delta function since the $z_{i}$’s become linearly dependent Ita:2015tya ,

[TABLE]

The Baikov polynomial $F(z)$ is invariant under a cyclic permutation of the indices, $z_{i}\rightarrow z_{i+1},s_{ij}\rightarrow s_{i+1\,j+1}$, modulo $5$. Let us consider the $4$-dimensional cut $z_{1}=z_{2}=z_{3}=z_{4}=0$, which localizes Eq. (51) to, omitting constant factors,

[TABLE]

where $F(0)$ is again a shorthand for

[TABLE]

and $z_{5}^{\pm}$ are the two solutions for $z_{5}$ on the cut,

[TABLE]

Choosing

[TABLE]

then Eq. (52) produces the leading singularity $1$ on the cut solution $z_{5}=z_{5}^{-}$, and [math] on the other cut solution $z_{5}=z_{5}^{+}$. One important feature which makes Eq. (55) a valid ansatz is that the $z_{5}$-independent term, $-F(0)/\sqrt{\det G}$, is a cyclic invariant expression which is the same on all the five possible $4D$ maximal cuts. So Eq. (55) can be promoted to the full uncut expression for $\mathcal{N}$ by combining the $z$-dependent terms obtained on the individual cuts, and the result corresponds to the chiral pentagon integral of Ref. ArkaniHamed:2010gh up to terms that vanish on $4$-dimensional box cuts.

4 A non-planar five-point topology at two-loops

To illustrate our method for computing differential equations on unitarity cuts, we consider two-loop five-point scattering. While planar master integrals have been computed Gehrmann:2015bfy ; Papadopoulos:2015jft , the nonplanar counterpart remains unknown. Ref. Frellesvig:2017aai suggested probing the properties of such integrals via the maximal cut. We focus on the nonplanar pentabox in Fig. 3, and compute differential equations at the maximal cut level.

The inverse propagators are

[TABLE]

while the irreducible numerators are initially chosen as

[TABLE]

Using inverse propagator coordinates, the loop integral can be written as, omitting overall factors,

[TABLE]

where $\det G$ is the Gram determinant for the external momenta, defined in the same way as for the one-loop pentagon, and $F(z)$ is defined by the $(d-4)$-dimensional components of $l_{1}$ and $l_{2}$, which are $\mu_{1}$ and $\mu_{2}$, by

[TABLE]

The five independent kinematic invariants are defined in the same way as for the one-loop pentagon, as

[TABLE]

The loop integrals, up to an overall power of $s_{12}$, depend on the four dimensionless ratios,

[TABLE]

We perform IBP reduction at the maximal level for tensor numerators with up to $6$ total powers of $z_{9}$, $z_{10}$ and $z_{11}$, using an in-house implementation of the algorithm of Larsen:2015ped . This involves solving a syzygy equation, which is the same as our Eq. (42) with the RHS set to zero, to find vectors whose total divergences generate IBP relations without doubled propagators. Three master integrals, $\mathcal{I}_{1}$, $\mathcal{I}_{2}$, and $\mathcal{I}_{3}$, are found, with tensor numerators

[TABLE]

respectively.

Next, the derivatives w.r.t. the four scaleless kinematic invariants, $\partial/\partial\chi_{ij}$, are each written in terms of external momentum derivatives $\sim\beta^{\mu}(\partial/\partial p^{\mu})$, in the fashion of Eq. (43). This fixes the RHS of Eq. (42). On a unitarity cut, Eq. (42) is simplified, allowing a solution to be found by the computer more quickly.555If a particular $z_{k}$ is set to zero, the $f_{k}$ term can simply be dropped from the LHS of the equation, while the RHS of the equation needs to be evaluated before $z_{k}$ (and any other inverse propagators in the cut) is set to zero, since these two operations do not commute. Finding one particular solution for $f_{j}$ and $v_{F}$ using SINGULAR gives us the needed ingredients to compute the derivative of the loop integral in Eq. (4). We obtain a combination of tensor integrals without doubled propagators, which are then reduced to the three master integrals using the IBP reduction procedure described above. The end results are the differential equations relating the three master integrals with their partial derivatives against the four scaleless kinematic invariants,

[TABLE]

The full results for the matrix $M_{AB}^{ij}$ are attached in an ancillary file pentaCrossBox.m in the Mathematica format. We will elaborate on the computation techniques in Subsection 5.2. A sample component is

[TABLE]

which have poles that can be identified with kinematic singularities, such as $\chi_{45}\propto(p_{4}+p_{5})^{2}$, and $(\chi_{23}-\chi_{45}+1)\propto(p_{1}+p_{3})^{2}$. Since Eq. (63) allows second derivatives to be computed, we are able to perform another check using the consistency condition Meyer:2016slj ,

[TABLE]

As is the case for the component shown in Eq. (64), the complete matrix is linear in $\epsilon$ but not proportional to $\epsilon$, because we have not yet transformed the result into a “good” basis of master integrals with unit leading singularities. Such a “good” basis can be found easily using the 4D cut method in Henn:2013pwa . Roughly speaking, cutting the $l_{2}$ box sub-loop produces the Jacobian

[TABLE]

So we can recycle the one-loop chiral pentagon expression Eq. (55) to write down three tensor integrals, $\tilde{\mathcal{I}}_{1}$, $\tilde{\mathcal{I}}_{2}$, and $\tilde{\mathcal{I}}_{3}$ with unit leading singularities. Their numerators are (ignoring the dependence of the chiral pentagon numerator $\mathcal{N}^{\rm chiral}$ on $z_{1},z_{2},z_{3},z_{4}$ which vanish on the maximal cut of the nonplanar pentabox),

[TABLE]

respectively. The first two of these integrals are among the nonplanar $\mathcal{N}=4$ SYM integrands for this particular topology given in Bern:2015ple .

The new basis is related to the old one via a matrix $T$,

[TABLE]

and the differential equations in the new basis,

[TABLE]

are related to the old system, Eq. (63), by the transformation formula,

[TABLE]

After the transformation, the matrices are proportional to $\epsilon=-(d-4)/2$, and now involve not only polynomials but also square roots of the Gram determinant. The square roots are eliminated by switching to momentum-twistor variables Hodges:2009hk using the parameterization given in Appendix (A.2) of Ref. Badger:2013gxa . The momentum-twistor variables, $x_{i}$ with $i=1,2,3,4,5$, are related to the usual kinematic invariants by

[TABLE]

where $\operatorname{tr}_{5}$ is defined via the Gram determinant,

[TABLE]

Since there are only $4$ dimensionless ratios of kinematic invariants, we will fix $x_{1}=s_{12}=1$, effectively only looking at the dependence of the integrals on $x_{2},x_{3},x_{4},x_{5}$. After re-writing the differential equations in terms of the above momentum-twistor variables, we used the CANONICA software package Meyer:2017joq to transform the differential equation into a form that contains “dlogs” Henn:2013pwa (a few seconds of computation time is used),

[TABLE]

In the above equation, $\vec{I}$ is a column vector consisting of the $3$ maximal-cut master integrals of unit leading singularities, $\tilde{\mathcal{I}}_{1}$, $\tilde{\mathcal{I}}_{2}$, and $\tilde{\mathcal{I}}_{3}$, defined via the tensor numerators in Eq. (67). The $s_{i}$ variables are the so called symbol letters, with polynomial dependence on the momentum-twistor variables. The $3\times 3$ matrices $\mathbb{M}_{i}$ are purely numerical, with no dependence on the dimension $d$ or the kinematic / momentum-twistor variables. The explicit expressions for the symbol letters are,

[TABLE]

and the explicit expressions for the $3\times 3$ matrices are,

[TABLE]

This form of the maximal-cut differential equations indicate that the solutions are multiple polylogarithms involving the $11$ symbol letters in Eq. (74), with uniform transcendentality in the $\epsilon$ expansion, as explained in Refs. Henn:2013pwa ; Henn:2014qga .

5 Computation techniques and timing comparisons

5.1 Double box and comparison with FIRE5

Having illustrated our method at the one-loop level in Section 3 and at the two-loop level in Section 4, we test the method for the double box and compare with existing methods. Being not too simple or too complicated, the double box topology can be handled by both traditional and unitarity-based methods, allowing for a meaningful comparison.

The double box with massless internal and external lines is shown in Fig. 4, with kinematic variables defined in the caption. It is well known Smirnov:1999wz that there are $2$ top-level master integrals, which may be chosen as the scalar integral $I_{\rm dbox}^{s}$ and the tensor integral $I_{\rm dbox}^{t}$,

[TABLE]

With daughter topologies included, there are $12$ master integrals in total. Among them, $8$ master integrals are independent after accounting for the discrete symmetry,

[TABLE]

which leaves the $7$ propagators invariant up to a permutation. We re-compute the most non-trivial part of the system of differential equations, i.e. the $t$-derivative of the $2$ top-level master integrals, expressed in terms of the $8$ master integrals.

Our method is applied to a spanning set of $6$ different unitarity cuts,

[TABLE]

The discrete symmetry Eq. (77) is used to speed up the calculation by relating the cut $\Gamma_{2a}$ with $\Gamma_{2b}$, and $\Gamma_{3a}$ with $\Gamma_{3b}$. In particular, the IBP relations generated on one cut automatically become valid IBP relations on the related cut after the symmetry transformation. Not surprisingly, these are essentially the same cuts used in Ref. Larsen:2015ped for IBP reduction of double box integrals.

By merging the results on the $6$ cuts, we reproduce the following differential equations,

[TABLE]

where daughter topology integrals are fully computed but omitted in the above sample results.

Both FIRE5 and our own unitarity-based code begin with “preparation runs” purely in Mathematica to process the diagram topology information supplied by the user. For the “final runs”, FIRE5 is used in the C++ mode, while our own code is written in Mathematica and SINGULAR. Only the final runs, which reflect the true computational complexities, are included in our timing comparison. The time required to obtain the results in Eq. (84) is shown in Table 1.

Though there could be more improvements by switching to statically compiled computer languages such as C++, our code already offers significant improvement in speed over a calculation based on FIRE5, for the following possible reasons,

The lack of doubled propagators reduces the number of different integrals that appear in IBP relations, because for combinatorial reasons, the number of integrals grow rapidly when propagators are allowed to be raised to higher powers. 2. 2.

The use of unitarity cuts reduces the computational complexity. 3. 3.

Besides solving syzygy equations, SINGULAR is also used to solve the sparse linear system formed by the IBP relations. Unfortunately, due to lack of a controlled comparison, we do not know how this affects performance compared to Fermat Smirnov:2008iw used in the C++ version of FIRE5.

5.2 Finite field techniques and rational function reconstruction

The computation of the differential equations for the nonplanar pentabox posed challenges in terms of CPU time and memory consumption. Given that the results, e.g. Eq. (64), are rational functions in $\chi_{ij}$, we use the rational function reconstruction technique of Peraro:2016wsq to fit the analytic result from numerical inputs of $\chi_{ij}$.

The algorithm of Peraro:2016wsq reconstructs multivariate rational functions in two steps, (i) fitting univariate rational functions, and (ii) fitting multivariate polynomials, using the input from many iterations of step (i). We use a simple private implementation which performs most of the work by exploiting the built-in capabilities of Wolfram Mathematica.666In Mathematica 10, step (i) is accomplished by the command FindSequenceFunction with the option FunctionSpace -> "RationalFunction". Step (ii) is accomplished by the command InterpolatingPolynomial. The latter command allows the option of computing in a finite field $\mathbb{Z}_{p}$. For the nonplanar pentabox computation, step (i) requires $18$ kinematic points for each iteration, and $495$ iterations are performed to produce the input for step (ii). Step (ii) fits polynomials in $3$ variables, with degrees up to $8$. Finite field techniques are used to accelerate step (ii): using a large prime $p$, the full result can be constructed from its image in $\mathbb{Z}_{p}$ probabilistically using a minor modification of the extended Euclid algorithm Peraro:2016wsq ; vonManteuffel:2014ixa . After completing steps (i) and (ii), the fitted results are validated against new computations with additional random rational values of $\chi_{ij}$.

Here we give more information to quantify the performance gains from rational function reconstruction and finite field techniques. The computation of the differential equations for the nonplanar pentabox, with analytic dependence on the kinematic invariants, is very time consuming and does not finish after $48$ hours.777The computation gets stuck at the first stage, namely finding IBP-generating vectors that do not cause doubled propagators. However, with (rational) numerical kinematic invariants, the computation finishes in a few seconds per kinematic point on a modern computer. A total of $28$ hours is used in evaluating differential equations on $8910$ kinematic points and reconstructing the full analytic results. The last step of the calculation, i.e. multivariate polynomial fitting, dramatically benefits from finite field techniques and takes about $52$ seconds to reconstruct all the $36$ entries of the four $3\times 3$ matrices. In contrast, when finite field techniques are turned off, about $395$ seconds are needed to reconstruct only one of these $36$ entries (with the computation aborted afterwards), which is slower by more than $2$ orders of magnitude.

6 Conclusions

We have proposed a new method for constructing differential equations for Feynman integrals, which avoids generating integrals with doubled propagators, instead producing tensor integrals to be reduced by unitarity-compatible IBP reduction. In fact, for the simplest cases such as the one-loop box, no IBP reduction is needed at all. Our method allows constructing differential equations from a spanning set of unitarity cuts in $d$ dimensions, with IBP reduction also performed on the cuts.

Applying our method to the nonplanar pentabox, we obtained the homogeneous differential equations on the maximal cut in Henn’s canonical form. This allows us to confirm that the master integrals, at least when evaluated on the maximal cut, are multiple polylogarithms with uniform transcendentality in the $\epsilon$ expansion. We have extracted the $11$ symbol letters, which are polynomials of momentum-twistor variables.

We also demonstrated that finite field techniques and rational function reconstruction, which are emerging as new tools in studying scattering amplitudes Peraro:2016wsq ; vonManteuffel:2014ixa ; vonManteuffel:2016xki , are useful in computing differential equations for Feynman integrals.

There are several possible directions for follow-up studies. One direction is extending the calculation to other nonplanar five-point topologies, which can be done straightforwardly. It would be desirable to construct an automated implementation of our method, perhaps as an extension to unitarity-compatible IBP reduction software packages, such as Azurite Georgoudis:2016wff , since many computation steps can be shared. Eventually, we would like to construct full DEs for nonplanar two-loop five-point integrals, which are relevant for NNLO QCD corrections for $2\rightarrow 3$ scattering processes at the LHC Gehrmann:2015bfy ; Papadopoulos:2015jft .

7 Acknowledgment

We thank Zvi Bern and Harald Ita for enlightening discussions and comments on the manuscript, and Yang Zhang for enlightening discussions and sharing ideas for efficient implementations of unitarity-compatible IBP reduction. The work of MZ is supported by the Department of Energy under Award Number DE-SC0009937.