Pentagon functions for massless planar scattering amplitudes

Loop amplitudes for massless five particle scattering processes contain Feynman integrals depending on the external momentum invariants: pentagon functions. We perform a detailed study of the analyticity properties and cut structure of these functions up to two loops in the planar case, where we classify and identify the minimal set of basis functions. They are computed from the canonical form of their differential equations and expressed in terms of generalized polylogarithms, or alternatively as one-dimensional integrals. We present analytical expressions and numerical evaluation routines for these pentagon functions, in all kinematical configurations relevant to five-particle scattering processes.


Introduction
Predictions for scattering process in elementary particle physics can be computed to high accuracy through a perturbation series expansion of the relevant scattering amplitudes. In this expansion, higher order corrections imply an increasing number of closed particle loops, leading to Feynman integrals over the loop momenta. While one-loop corrections are known for processes with arbitrary multiplicity [1], two-loop corrections to scattering amplitudes have up to now only been derived on a case-by-case basis, mostly for two-to-two scattering processes. In going to higher multiplicities for two-loop amplitudes, one faces two challenges: to express the large number of Feynman integrals in terms of a smaller set of basis integrals (for a minimal basis, these are often called master integrals) and to efficiently compute these basis integrals.
In the recent past, important progress has been made on integral reduction techniques for two-loop multi-particle processes, both in terms of semi-analytical approaches [2][3][4] as well as in optimizing algebraic reductions to master integrals [5][6][7][8][9][10] to cope with the high complexity of processes with five external particles [11,12] and beyond. As a result, expressions for two-loop five-gluon amplitudes in terms of a set of basis integrals were derived, first for specific helicity configurations [13][14][15][16], and most recently for the generalhelicity case [4,17].
These basis integrals can be expressed in terms of a set of master integrals [14,18]: massless two-loop five-point functions. A subset of these functions are two-loop fourpoint functions with one off-shell leg, which were computed in analytical form [19][20][21] already long ago in the context of lower multiplicity processes. The genuine five-point master integrals can be separated into planar and non-planar topologies, depending on the internal momentum routing. For the planar integrals, differential equations in the external momentum invariants [22,23] were derived and solved in [14] and [18]. In this paper, we build upon our work on the planar master integrals in [14] by taking a systematic approach aiming to combine the differential equations with knowledge on the kinematical analyticity structure of the master integrals to identify the minimal set of functions that can appear in them. These pentagon functions are the basic building blocks for two-loop five point master integrals, and our procedure used for their identification can be expected to generalise to non-planar integrals, to higher multiplicities and to higher loop orders. Various representations for these pentagon functions can be derived using the differential equation method.
Fully analytical expressions are found in terms of generalized polylogarithms [24,25] or Chen iterated integrals [26]. These functions have a long history in the mathematics of Feynman diagrams. In recent years, huge progress was made in understanding how to handle these functions systematically. This concerns in particular the multi-variable case. An important tool is the so-called 'symbol' of iterated integrals [27][28][29][30], which makes it easy to understand identities between different functions. Closely related to this is the idea of defining the special functions needed for Feynman integrals from certain canonical differential equations [23]. The latter encode the relevant data about the functions in a compact and unique way. In particular, they contain the 'symbol alphabet', denoting the integration kernels that are allowed to appear in the iterated integrals.
While the fully analytical expressions enable detailed studies of analyticity properties and of asymptotic properties, their numerical evaluation is rather inefficient. Instead, following [31], one can derive one-dimensional integral representations, that are optimised for numerical integration.
The paper is structured as follows. Following a brief description of the kinematics of five-particle scattering in Euclidean and Minkowskian space in Section 2 and of the notation used for the planar two-loop five-point master integrals in Section 3, we derive and analyse the differential equations for these master integrals in Section 4. Confronting the singularity structure of the differential equations with the physical requirements on the analyticity properties of the master integrals provides strong constraints on the pentagon functions that are allowed in the solutions of the differential equations, as discussed in Section 5. The pentagon functions are then computed in Section 6 by matching the generic solutions of the differential equations onto boundary conditions at specific kinematical points. In Sections 7 and 8, we describe a variety of consistency checks on these results and introduce a public numerical code which evaluates the pentagon functions and master integrals. We conclude with an outlook in Section 9.

Conventions and Lorentz invariants
The kinematics is described by five external momenta, p µ i , subject to the on-shell conditions p 2 i = 0, and momentum conservation 5 i=1 p µ i = 0. From the momenta, we can build ten scalar products s ij = 2p i · p j , of which five are independent.
We choose the following five Note that, in practice, we can always fix the overall scale of a quantity by dimensional arguments, so that effectively we need to deal with four-variable functions only. For most of the discussion, however, we prefer to keep working with the five variables v i , as they allow to see symmetries in an easier way. The above variables are parity even. There is also a parity odd invariant, We can also write this as (1234) = tr(γ 5/ p 1 / p 2 / p 3 / p 4 ). Note that µνρσ is a four-dimensional object. The calculations in this paper are done for general dimension D = 4 − 2 . However, since only four of the external momenta are independent, we can assume without loss of generality that they lie in some four-dimensional subspace.
We also introduce the following dimensionless Gram determinant of the four linearly independent vectors p i , i = 1 . . . 4, It is related to the parity-odd Lorentz invariant (1234) according to (2.4)

Physical region
The physical regions of the five-point functions are given by 2 → 3 scattering kinematics (i + j → k + l + m). Any pair of two momenta (i, j) can be incoming, such that the physical region in Minkowski space corresponds to ten distinct regions (channels), which are commonly labelled by their initial-state invariant as s ij -channel. For each channel, there are constraints on the signs of the the kinematic invariants defined as the scalar products between two external momenta. These are summarised in Table 1.
Individual five-point integrals may be related among different channels through momentum permutations. The kinematical region in each channel is delimited by requiring all s-channel invariants positive and all t-channel invariants negative, plus negativity ∆ ≤ 0 of the Gram determinant (following from the real-valuedness of all momenta).
On the example of the s 12 -channel, these non-Gram-determinant constraints imply for the independent invariants: The last remaining independent invariant is then constrained by the positivity of the Gram determinant: For fixed values of s 12 , s 34 , s 45 , these constraints describe an ellipse in the (s 23 , s 15 )-plane, as shown in an example in Figure 1.  3 Two-loop five-point planar master integrals The family of penta-box integrals is defined as with p 2 i = 0, i = 1, . . . 5, and 5 i=1 p µ i = 0, and where a 1 , . . . a 8 ≥ 0 are propagators and a 9 , a 10 , a 11 ≤ 0 numerator factors. See Figure 2. Integral reduction [32] (using for example FIRE [33] or Reduze [34]) shows that there are 61 master integrals for this family of integrals. The master integrals can be organized in terms of integral sectors, corresponding to the 8-propagator sector, and subsectors with fewer propagators. In total, one needs 46 distinct sectors. One can further organize the latter by grouping together sectors that are related by relabelling the external legs. In this way, one can group all integrals in terms of 17 sectors, which are classified in the following. First, there are a number of integrals that are products of one-loop integrals. They are shown in Fig. 3. In this figure, the I j indicate at which position these integrals appear in the basis that we chose (the latter will be discussed below). Next, integrals corresponding to four-point functions with one off-shell leg are shown in Fig. 4. They are known from ref. [19,21], where they are expressed in terms of generalized harmonic polylogarithms [24,29], for which efficient numerical evaluations are available [35][36][37]. The genuine five-point sectors are shown in Fig. 5.
For what follows, we need to choose a basis of 61 integrals. Traditionally, integrals were typically chosen relatively randomly, as implied by the lexicographic ordering [32] that was used in the integral reduction. In ref. [23], it was proposed to choose the basis such that the integrals have simple properties. This considerably simplifies their computation, and allows to obtain the result in a form that is as compact as possible. As we will see in the next section, in the expansion, all integrals evaluate to multiple polylogarithms. In general, such Feynman integrals will involve linear combinations of multiple polylogarithms of varying transcendental weight, and with various prefactors that depend rationally or algebraically on the kinematics. It is desirable to disentangle the latter, so that such factors are moved into overall normalizations of the integrals, and such that only functions of homogeneous weight appear.
Understanding for predicting which integrals have this property came initially from studies in N = 4 super Yang-Mills. Conjecturally, integrals whose integrands can be written as a 'd-log'-form, and hence have constant leading singularities [38], have this property. The initial examples satisfying this conjecture were massless, planar, finite, dual conformal integrals. It has since been generalized to more generic integrals within dimensional regularization [23].
What is important to emphasize is that the basis choice can be done a priori, by analyzing the loop integrand. In principle, one could classify all integrands having the  desired properties, and then select a linearly independent (under integral reduction) subset. This can be done algorithmically, see e.g. [39].
In practice, it may not be necessary to classify all such integrals, but just to construct a sufficient number of them. It is possible to construct many 'd-log' integrals directly, for example by iteratively using lower-loop building blocks. See refs. [40,41] for examples. Our choice of 61 basis integrals is given in the ancillary file pentabox basis2.txt.

Differential equations for the master integrals
To compute the master integrals, the explicit and complicated integration over the loop momenta can often be avoided for multi-scale integrals by using differential equations in kinematical invariants, as first demonstrated for the two-loop four-point functions in [22].
We use the integral basis I j , with j = 1 . . . 61, discussed in the previous section, and compute the differential in all variables v j , j = 1, . . . 5. We find the following canonical form of the differential equations [23] d I(v i ; ) = dÃ I(v i ; ) , (4.1) with the matrixÃ The ensemble of the letters W i is called the alphabet of the problem under consideration. We derived the planar pentagon alphabet in our previous work [14], which was subsequently extended to account also for letters relevant to non-planar pentagon functions in ref. [42], whose notation we adopt here: 3) with W 1+i , W 6+i , W 11+i , W 16+i , W 21+i , W 26+i , with i = 1 . . . 5, defined by cyclic symmetry. Note that the W i , with i = 26, . . . 30, are parity-odd, in the sense that they go to their inverse under ∆ → −∆, while all other letters are parity-even under that transformation. Finally, the a i are constant 61 × 61 matrices. The W i notation covers that case of planar [14] as well as non-planar pentagon functions [42]. Following the notation of that reference, the full 31-letter alphabet is called A NP . Here we will only need the planar case, which consists of the 26  The planar 26-letter alphabet A P covers the pentabox integrals (and all their subintegrals), as well as cyclic rotations thereof. If one focuses only on the orientation of the pentabox shown in Fig. 2, only 24 of those letters occur; the letters that are absent in that case are W 8 and W 10 .

Chen iterated integrals for solving the differential equation
Given the first-order system of differential equations (4.1), the solution is completely determined upon giving a boundary condition. It turns out that the latter can be obtained by simple physical considerations. This will be discussed in section 6.
What functions will appear in the solution of eq. (4.1)? In practice, we need to solve this equation in an expansion in , i.e. (4.10) Note that our integral basis I is normalized, without loss of generality, such that this expansion starts at 0 . Inserting (4.10) into (4.1), we see that the equation decouples order by order in , This means that we have where C is an integration contour in the space of kinematic variables v i , and I We see that the solution, to all orders in the expansion, is given in terms of iterated integrals. The integration kernels are given by logarithmic differential forms, see (4.2). The possible arguments of the logarithms (the letters) are given by the set A P . We call the functions arising from this alphabet planar pentagon functions. Along the lines of this terminology, let us observe that the matrixÃ of eq. (4.2) dictates, via eq. (4.12), how words are build up from the letters. In this sense we can think of the constant matrices a i in eq. (4.2) as the 'grammar' needed to form words.
It is convenient to introduce a shorthand notation for the iterated integrals. For simplicity, let us first give the definition for iterated integrals of a single variable x, and choose (here) as boundary point of the integration x 0 = 1. We denote iterated integrals by brackets [. . .]. Integrals are defined iteratively, namely (4.13) Here the integration goes along the path The iteration starts with the empty bracket [ ] ≡ 1. The number of entries n is called weight. For example, we have the weight-one function 14) and the weight two function An important difference of the pentagon functions w.r.t. the single-variable case is that they are defined for the five-dimensional kinematic space {s i,i+1 }. Just as in the above examples, the integrals are defined along a path, parametrized by a variable t ∈ [0, 1]. Let us start with a single integral. We denote where integration path C starts at the boundary point (−1, −1, −1, −1, −1), and goes to the function argument (which we assume to be in the Euclidean region for now), without picking up monodromies. To be completely explicit, we could choose a straight path Other kinematic regions can be obtained via analytic continuation. Just as above, we can define iterated integrals. There is one new feature of the multi-variable case with respect to the single-variable one. In order for the integral to be well-defined, it is important that it is independent of the choice of contour, as long as singularities are not crossed. In other words, the integral has to be homotopy invariant. These conditions are called integrability conditions [27]. For a generic term of the form This equation puts a constraint on the consecutive entries in the iterated integrals. In section 5, we will classify all iterated integrals built from the alphabets A P and A P;1 satisfying the integrability conditions. Letter v notation momentum notation cylic Table 2. Interpretation of pentagon alphabet in terms of particle momenta.

Properties of the pentagon alphabet
The iterated integrals that can appear in our case are characterized by the alphabet A P . We wrote the alphabet in terms of the Mandelstam invariants s i,i+1 . As we will see, it is instructive to rewrite it in terms of different variables. This will allow us to see an underlying simplicity of this alphabet. Let us discuss the different types of letters, and reveal their simple dependence on the external momenta p i . This will give us insights into simple parametrizations of the alphabet, and into the singularity structure of the associated functions. The first 25 letters W 1 to W 25 (out of which 20 belong to A P ) are simple scalar products of the loop momenta. Their interpretation is straightforward: these are possible singularities of Feynman integrals. One could have discovered these, for example, by an analysis of the Landau equations. Feynman integrals are multivalued functions, so it is expected that they can have branch cuts. However, for planar integrals, only the first five singularities correspond to branch cuts on the first sheet of the functions. All these letters appear already in four-point integrals with one leg off-shell [19,20], with the planar ones shown in Fig. 4 (and cyclic permutations thereof).
Next, we have the odd letters W 26 , . . . W 30 , that are genuine to five-particle kinematics. They appear already at one loop, in particular in the six-dimensional pentagon integral.
Remarkably, as pointed out in ref. [14], they can be written as ratios of traces with a very simple dependence on the momenta, e.g.
Note that writing the letters as ratios is a choice. Using the property of the logarithm d log(αβ) = d log α + d log β, we could have equally well chosen their numerators as independent letters. This is possible since, e.g.
which does not contain new letters. We prefer to use the ratios, as they have simple transformation properties under parity. Finally, we have the Gram determinant as an independent letter. Since log(∆) = 2 log √ ∆ we could equally use (1234) = tr[γ 5/ p 1 / p 2 / p 3 / p 4 ]. We summarize the letters of the alphabet, and the equivalent ways of expressing them in Table 2. From the equivalent representations of the alphabet letters, we can deduce three remarkable properties of the alphabet.
The first property we comment on is that the alphabet is to a large part determined by the singularity structure of Feynman integrals, i.e. the locations in kinematic space where singularities can occur. This can be analyzed, in principle, via Landau equations. We remark that the alphabet knows not only about singularities on the first sheet of the multi-valued functions, but about all singularities. In principle, the latter can be rather complicated. What we observe looking at Table 2 is that in our case the allowed singularities follow a simple pattern.
Remarkable property #1: the singularities of the pentagon alphabet (4.3) all correspond to exceptional configurations of the external momenta p i · p j = 0 or p i · (p j + p k ) = 0, or to restricting them to a lower-dimensional subspace D < 4, where ∆ = 0.
It is an interesting open question whether these properties will continue to hold at higher loop orders, or whether the alphabet may have to be enlarged in that case. This is relevant for understanding the function space at higher loops, and can have bootstrap applications, see [43] for a recent discussion.
The next property we wish to emphasize is somewhat related to the first one, and concerns the linearity of the alphabet in the momenta. This has an immediate application. Sometimes, it can be useful to parametrize the alphabet in a way that is rational in a given variable, and as simple as possible. We remark that the above linearity suggests a way of introducing a variable parametrizing one particular direction of the kinematic space. The idea is to deform the external momenta, while preserving momentum conservation and the on-shell conditions. This concept is well-known and important in the context of BCFW recursion relations for scattering amplitudes [44]. Writing the momenta in terms of spinors, p i = λ iλi , one defines for example the deformation In this way, one can study the differential equation along the direction z. It turns out that the alphabet depends linearly on the shift parameter z. The simplicity of our alphabet is also related to a similar parametrization considered in ref. [18]. Summarizing, we have Remarkable property #2: the pentagon alphabet (4.3) can be written in a way that is linear in the external momenta. This allows to introduce BCFW shifts, with a linear dependence on the parameter.
The same comment applies to the extension to non-planar pentagon functions of ref. [42].
In fact, it is possible to describe not just a single direction of the alphabet in a rational way, but in fact define a change of variables that completely rationalizes them. There are various ways of doing this.
One way is using a parametrization suggested by the spinor helicity variables. In four dimensions, we can write the on-shell momenta in terms of spinors, and use the standard bra-ket notation for spinor invariants. In that language, it is clear that letters Note that while the spinors manifestly solve the on-shell conditions, they are not independent, due to momentum conservation. In ref. [45], an independent set of spinor products was used to parametrise the five-particle kinematics. Given the fact that all letters are simple rational functions in terms of the spinor brackets, it is not surprising that this leads to a rational version of the alphabet. We collect the variables, as well as the alphabet in this parametrization in Appendix A.2.
We remark in passing that eq. (4.22) shows that the parity-odd letters W 26 , . . . W 30 can be interpreted as phases for real momenta (in Minkwoski space), i.e. when λ i andλ i are related by complex conjugation. Indeed, in that case we have |W 26 | = 1.
A parametrization closely related to spinors is that of momentum twistors. The latter variables have the advantage of solving both the on-shell and momentum conservation constraints, and are hence free variables. Additionally, they offer a geometric interpretation of the singularities of the functions. Like the spinor parametrization, they rationalize the alphabet. We review this in Appendix A.1.
In summary, we see a third remarkable feature of the pentagon alphabet: Remarkable property #3: the pentagon alphabet can be naturally parametrized in a rational way using spinor or momentum twistor variables.
Both property #2 and property #3 imply that, if desired, the pentagon functions can be represented by Goncharov polylogarithms. In summary, the solutions of eq. (4.1) can be expressed as iterated integrals in the space of pentagon functions built from the alphabet A P . In principle, we could proceed with solving the differential equations for the master integrals directly. We find it instructive, however, to first study the function space from a slightly more general point of view. This will then allow us to express the solutions of the differential equations in terms of a minimal number of functions.

Classification of planar pentagon functions
The differential equations we provided can in principle be solved to any desired order in the expansion. For physical applications to two-loop amplitudes, one typically wishes Although there is no general proof to our knowledge, it is generally believed that twoloop Feynman integrals in four dimensions involve weight four functions at most. For this reason we focus on this case in particular, and discuss in detail the functions that make an appearance.

Planar pentagon functions to weight four
Physical amplitudes are expected to have branch cuts only originating at v i = 0. Therefore we will focus on functions where only v i appear in the first letter. We proceed by writing all possible words in the alphabet of a given weight, and imposing this first entry condition, as well as the integrability conditions. Table 3 summarizes the classification of pentagon functions up to weight four. In the second line, the number of integrable symbols at a given weight is given.
When working with iterated integrals, products of lower-weight integrals can be expanded (via the shuffle algebra) into sums of higher-weight functions. Inverting these relations, we can remove all product terms from a given expression, so that only 'irreducible' functions remain. This organisation is very useful, as the product terms are faster to evaluate (numerically). Let n 1 , n 2 , n 3 , n 4 be the number of irreducible functions at weight 1, 2, 3, 4. Then, the number of product functions at weights 2, 3, 4 is #prod(2) = 1 2 (1 + n 1 )n 1 , #prod(4) = 1 24 n 1 (1 + n 1 )(2 + n 1 )(3 + n 1 ) + 1 2 n 1 (1 + n 1 )n 2 + 1 2 respectively. The number of products functions and of irreducible (new) functions are given in lines 3 and 4 of Table 3. The functions can be further classified by the number of entries containing the parityodd letters letters W 26 to W 30 . This is done in the last three lines of Table 3. Finally, integrals can be organised according to cyclic symmetry, into either quintets or singlets. For example, at weight one, the five functions all are obtained by cyclic symmetry from one basic function. At this stage we could proceed and select specific functions representing the various entries on the last three lines of Table 3. These would be guaranteed to cover the solution space of the differential equations. However, we anticipate a further simplification. As we will see, most functions that are needed actually depend on fewer than the 26 letters. This is closely related to a conjectured second entry condition of ref. [42]. It turns out that up to weight four, the letters v 1 + v 2 and cyclic only appear in the slashed box function shown in Fig. 4(h).
This fact motivates to repeat the above classification of pentagon functions for a smaller 21-letter alphabet, with letters W 6 , . . . W 10 removed. As mentioned earlier, this is the oneloop alphabet A P;1 . We see that in this case, the number of functions needed is reduced considerably, see Table 4. In the following subsections, we will explicitly construct this basis of iterated integrals. We introduce all functions needed to describe planar five-particle amplitudes at two loops, up to weight four. The notation we use is f (k) i,j , where i refers to the weight, j is a label, and k = 1 . . . 5 refers to different cyclic orderings of the same function. For functions that are singlets of the cyclic group, we simply use the notation f i,j .

Weight one functions
At weight one, the allowed symbols are dictated by the first entry condition, i.e. the branch cut structure of the integrals, These integrals evaluate to logarithms. Using the definition of the iterated integral, with the boundary point v i = −1, we have This formula is manifestly well-defined in the entire Euclidean region v i < 0. To define the function in other regions, one adds the usual Feynman +i0 prescription, i.e. −v i → −v i −i0, and analytically continues [21].

Weight two functions
It turns out that the combination of first entry conditions and integrability of the iterated integrals is rather restrictive. Apart from products of weight-one functions, only one new type of function is needed. It is given by Note that when generalizing to the non-planar pentagon alphabet A NP , a new type of weight two parity-odd function appears [42]. At this stage it is crucial to note that the weight one and two functions are very important. Following [31], all functions up to weight four can be written in terms of one-fold integral representations. We will discuss this in section 5.6. For this reason, in general at higher weight it will be sufficient to define the necessary functions, and give the representations we use for their numeric evaluation. For simple cases, we can of course use explicit representations in terms of polylogarithms similar to the ones given above.

Weight three functions
At weight three, according to Table 4, we have three even functions, as well as one odd function.
Here, we need for the first time functions that depend on three of the v i at the same time. This is typical of four-point kinematics, with one off-shell leg. For example, for I 39 shown in Fig. 4(i), the momenta are p 1 , p 2 , p 3 , p 4 + p 5 , such that the relevant invariants are v 1 , v 2 , v 4 .
Two of the functions we define are very simple, These iterated integrals can be done explicitly in terms of polylogarithms. For this, we use the definition of the iterated integrals. By definition, the first two integrations can be done in terms of the weight one and weight two functions given above. We have In the second equality, we have used eqs. (5.6) and Table 2. Here the integration path C goes from the boundary point (−1, −1, −1, −1, −1) to the function argument (which we assume to be in the Euclidean region for now), without picking up monodromies. To be completely explicit, we could choose a straight path v i (t) = t − 1 + t v i , parametrized by t ∈ [0, 1]. Carrying out the integration, we find These expressions are valid in the full Euclidean region. For other regions, one analytically continues, taking into account the Feynman i0 prescription. The third function is slightly more complicated. We define 3,3 for i = 2, 3, 4, 5 obtained from cyclic symmetry. Here the attentive reader might wonder about the role of the ζ 2 W 16 W 4 term, especially since W 16 is not part of the expected branch cuts. The reason for this term becomes transparent upon explicitly carrying out the first two integrations. We have 1,1 f 1,1 − f 1,1 f Recall that the integrals should be well-defined, and be free of branch cuts, in the Euclidean region v i < 0. One can check that the ζ 2 term in the second line of eq. (5.13) is necessary for the integral to be well-defined. Specifically, it ensures that the function multiplying d log W 16 vanishes at v 4 = v 1 + v 2 , i.e. when W 16 vanishes. Carrying out the integration, we find the following representation This formula for f 3,3 is valid for the region v 4 < v 1 < 0, v 4 < v 2 < 0, a subset of the Euclidean region. The corresponding formulas for other regions are obtained via analytic continuation. We do not print them here, but they are encoded in an ancillary file. This, and other practical questions about numerical evaluation, are discussed in more detail in sections 6 and 8.
Finally, there is one parity-odd function at weight three that depends on the full fivepoint kinematics. This function can be identified with the (normalized) six-dimensional one-loop pentagon integral Φ 5 , We denote it by f 3,4 = −Φ 5 . In our two-loop integral family, it appears as the weight three part of integral I 37 , see Fig. 5

(a). In terms of iterated integrals, it is given by
Here −d 37,3 is the value of Φ 5 at the symmetric point. Its analytic expression is given in an ancillary file. It follows from the previous subsections that the two inner integrations in the iterated-integral definition of Φ 5 can be expressed in terms of logarithms and dilogarithms.
In this way, we obtain the one-dimensional integral representation where the Mandelstam variables v i are given an implicit t dependence via the parametrization v i (t) = 1 + t (v i − 1).

Weight four functions
We have a total of 9 parity-even functions without odd letters, one parity-odd function, and one parity-even function with up to two odd letters. The definitions can be found in an ancillary file. Here we summarize the main properties of the functions. Three functions just depend on letters W 1 , W 4 , W 11 , similarly to functions f (i) 3,1 and f (i) 3,2 at weight three. Five further functions are expressed in terms of the kinematics of the box integrals with one off-shell leg, and depend on W 1 , W 2 , W 4 , W 11 , W 14 , W 16 letters only. As was mentioned earlier, one function corresponds to the slashed box integral I 22 shown in Fig. 4(h). It is the only function containing the letter W 9 = v 1 + v 2 . These functions are sufficient to describe all integrals shown in Figs. 3 and 4, thereby providing a minimal functional basis for the results of ref. [19].
Finally, in order to describe the integrals of Fig. 5, only three additional functions are needed. Let us present their main features (the precise definition is provided in an ancillary file.) One parity even function is defined as All functions are defined in an ancillary file. The total number of 'irreducible' functions at weights 1, 2, 3, 4 is 1, 1, 4, 12, respectively. The total number of 'irreducible' functions is 18. The final results for the integrals will be given in terms of these functions.
Summary of pentagon function classification: All planar pentagon functions up to weight four are expressed in terms of a basis of 18 irreducible functions, and permutations thereof. Only four of these functions depend on the genuine pentagon kinematics. When expressed in our basis, all identities between functions are manifest.

One-fold integral representations for basis functions
In the previous subsections, we defined a complete basis of functions. As we already mentioned, many functions appeared previously in the context of two-loop massive box integrals with one off-shell leg, and reliable numerical codes exist for their evaluation, for all relevant kinematic regions [36].
Therefore here we only need to discuss the new functions, f 3,4 = −Φ 5 at weight three, and the three functions f There are different choices of representations suitable for numerical evaluation. These include expressing the answer in terms of a minimal function basis (consisting of polylogarithms and Li 2,2 functions), in terms of Goncharov polylogarithms, or in terms of an iterated integral representation [31]. For a more complete discussion, including examples, we refer the reader to [31], where a similarly complicated alphabet was discussed. The upshot of that discussion is that while the first two types of representations offer advantages, such as a fast and reliable implementation of the basic functions, there are also disadvantages, such as a proliferation of terms. On the other hand, the iterated integrals one starts with are relatively compact.
What was suggested in [31] is a hybrid approach, where one carries out the first two integrations explicitly, in terms the logarithms and dilogarithms discussed above. Naively, it would then seem that one has two integrations remaining to get to weight four. However, by changing the order of integration, another integration can be done explicitly in terms of a logarithm. In this way one gets a one-fold representation that reads, schematically As a concrete example on how the parametrization proposed in ref. [31] allowed us to express integrals of a function of given weight as integrals of combinations of lower-weight functions let us revisit eq. (5.16) as follows: where the dependence on the variable t assigned to the letters W i and to the iteratedintegral functions indicates that all kinematic variables v i are parametrized as below eq. (5.17). If we were to carry out this integration explicitly we would obtain the result presented in eq. (5.17), where dilogarithms appear in the integrand. Alternatively, we can exchange the order with which we integrate over the variables t and t , obtaining Just as in [31], we coded up these one-fold integral representations and evaluated them numerically. These routines have been implemented in a public code, as described in section 8.

Integral basis in terms of pentagon functions
As a result of the classification of the previous sections, we can express the integrals appearing in the differential equations in terms of that function basis. Here, we give examples for some of the integrals that have genuine five-particle kinematics. For I 37 shown in Fig. 5(a) 1,1 + 2f 2,1 + 2f 1,1 + 3f 1,1 + 3f where the dots in the expression of I 60 indicate that, at each order in , only the functions with the highest weight have been included, whereas products of lower weight functions and constants have been omitted; notice that the expressions of I 37 , I 59 and I 61 are complete.
In the above formula, we have already used the explicit boundary values for the integrals obtained in the following section 6. Some of the boundary constants appearing at weight four are rather lengthy, so that we have abbreviated them as bc i . Their values are provided in pentagox basis2 bdry weight4.txt. The results for all integrals expressed in terms of pentagon functions can be found in masters-Walphabet f.txt.

Boundary conditions
In section 4, we derived a system of first-order differential equations that the master integrals satisfy, and subsequently identified the basis of functions that can appear in the solution to these differential equations. To fully specify the solution, we need to find the boundary conditions to the differential equations, i.e. the values of all integrals at specific kinematical points.
One might think that this requires a separate calculation of Feynman integrals. However, experience shows that boundary conditions for Feynman integrals can usually be obtained by the differential equations themselves, together with some physically motivated constraints. We also find this to be the case here. All except one boundary constant are obtained from an analysis of theÃ matrix in eq. (4.2). The final boundary constant represents an overall normalization, and is fixed by evaluating one of the trivial bubble type integrals. The latter can be given, to all orders in , in terms of Γ functions.
The conditions we will use come from the expectation that the integrals should be non-singular at several hypersurfaces. We have • no branch cuts within the Euclidean region, i.e. at v 1 − v 3 = 0; • no branch cuts in the u-channels: p i · p i+2 = 0; • we also expect the integrals to be finite at ∆ = 0, which corresponds to the external momenta lying in a three-dimensional subspace.
The predictive power of these conditions comes from the fact that the differential equations, and hence the matrixÃ, contain singularities at these locations. This means that the general solution to the equations has singularities at these locations, while these singularities should be spurious for the actual Feynman integrals (which may still contain discontinuities there).

Boundary conditions in the Euclidean region
The conditions described above give relations between values of the integrals evaluated at different hypersurfaces. It is desirable to have an explicit boundary value for all integrals at the same point. In the Euclidean region, where s i,i+1 < 0, the symmetric point s i,i+1 = −1 is a particularly convenient choice. 1 In order to compare the boundary values at the symmetric point, we transport them there, using the differential equation. In other words, we integrate the differential equation from some boundary point back to the symmetric point. As the problem is homotopy invariant, one is free to choose a convenient path. Often, this can be done in such a way that the answer is relatively simple. Let us give an example of this. We choose the parametrization Along this path, we find the following reduced alphabet By this we mean that eq. (4.1) becomes with the b i being constant matrixes, and β i ∈ B.
The boundary point ∆ = 0 corresponds to x = −1. The symmetric point at s 12 = −1 corresponds to 1 − 3x + x 2 = 0. We take the smaller solution x 0 = 1 2 (3 − √ 5). So, we can use eq. (6.3) to compute the connection between these two points. Some care is required due to the singular point at x = 0. For definiteness, we start in the Euclidean region, at x = x 0 . When reaching x = 0, we analytically continue, using the i0 prescription of the Feynman integral. Finally, we integrate from x = 0 to x = −1. The result of these integrations is expressed analytically in terms of Goncharov polylogarithms.
We proceed in a similar way for all boundary points. Finally, we fix the overall constant of integration (recall that (4.1) is a homogeneous equation) from a trivial bubble integral (8.1). In this way, all constants of integration are fixed.
Having fixed the boundary condition for all master integrals, in principle we could proceed and obtain a boundary value for each of the physical regions described in section 2 by analytic continuation. See section 6 of [46] for an example of this in a similar context. We find it most convenient, however, to repeat the above analysis of consistency equations, directly in the physical regions. This is discussed in the following section.

Boundary conditions in the physical regions
To evaluate the Feynman integrals for physical kinematics, we integrated the differentialequation system separately in each of the ten kinematical regions of Table 1, choosing as boundary point the value of the integral in points with a high degree of symmetry located inside each region (see Table 5).
The strategy we follow to evaluate such values is to exploit the cancellations of spurious singularities occurring on the boundary of the physical phase space to set constraints on value of the integrals. In particular, as the Gram determinant is negative definite for physical kinematics, each physical region is delimited by a hypersurface defined by the equation ∆ = 0 . On the other hand, all s i,j invariants, including non-adjacent ones, have a definite sign in the physical region, therefore also the hyperplanes defined by the equations s i,i+2 = 0 delimit the physical region. Points of tangency between the ∆ = 0 hypersurface and the s i,i+2 = 0 hyperplanes can therefore be reached from the boundary points through  Table 5. List of boundary points for each of the ten physical regions, together with the s i,i+2 = ∆ = 0 points used to obtain contraints from spurious-singularity cancellations. The ( * ) indicates the physical region where spurious-singularity cancellations in the four points reported leave the boundary condition on integral I 38 unconstrained. The value of the integral in this region can be found by analytically continuation from the Euclidean as described above eq. (6.6).
paths entirely contained in the physical regions. In these points of tangency, the differentialequation system features non-physical singularities. By not specifying the boundary values of integrals and by leaving them as unknown parameters, such divergences will also affect the integrated results, such that constraints on the boundary values themselves can be obtained by imposing the cancellation of these divergences. Considering only one point of tangency for each physical region does not provide enough constraints to be able to determine the boundary values of all integrals. The path to go from the boundary point (B) to the point in which spurious singularities appear (S) is parametrised as follows , For all points s B and s S given in the table above, the resulting reduced alphabet contains the square root √ −1 + t 2 . Such alphabet can be linearised with a change of variable , (6.5) which has been chosen so that the integration path is still parametrized by t ∈ [0, 1]. For the regions in which the incoming momenta are adjacent, it is necessary to consider three different s i,i+2 = ∆ = 0 points. For some of the regions corresponding to non-adjacent incoming momenta three points are not enough to constrain the system and a fourth point needs to be considered. For one region (s 13 channel) we obtained one boundary value (for integral I 38 ) by analytically continuing the solution from the Eulidean region. This was done by using the one-fold-integral representation of the Chen iterated integral expression, and by assigning a small imaginary part to the integration variable as follows where α represents any letter given in Table 2. The numerical evaluation of this Chen iterated integral is far more time-consuming than the evaluation of the pentagon basis functions described in Section 5, which is why this method is not used as a default for the implementation. Table 5 summarises, for each of the physical regions, the chosen boundary point and the auxiliary u-channel-zero points used to constrain the boundary value.

Checks on the results
The two-loop pentagon functions derived in the previous section can be combined to yield analytical expressions for the full set of planar massless five-point master integrals. Previous results on these master integrals concern the four-point sub-topologies (two-loop four-point functions with one off-shell leg [19,36]) and a numerical representation of the five-point integrals [18]. We checked numerically that our pentagon functions reproduce these known master integrals, finding machine precision agreement in the Euclidean region as well as in the Minkowskian region in all ten physical scattering channels. In the Euclidean region, we have moreover performed a numerical in-depth comparison between our expression in terms of Goncharov polylogarithms, and the one-dimensional integral representation, described in section 5.6. Again, machine precision agreement is found, with the observation that the one-dimensional integral representation evaluates considerably faster.
These master integrals appear in the expressions for all-massless two-loop five-point scattering amplitudes. Results for these amplitudes in the all-gluon case were first obtained for specific helicity configurations [13,15,16], and more recently for all general-helicity cases [4,17]. These results expressed the respective amplitudes as linear combinations of (numerically or analytically determined) integral coefficients, multiplying basis integrals. In these previous works, the five-point basis integrals were then evaluated numerically (restricted up to now to Euclidean kinematics), using the FIESTA sector decomposition code [47]. By re-expressing these basis integrals in terms of our master integrals, analytical expressions for two-loop five-gluon amplitudes can be obtained, such as in [14] for the allplus helicity amplitude of [13]. An early check of this result was obtained by an independent derivation of this specific amplitude [16].
Working with the authors of [4], we performed [48] a detailed comparison of all individual integrals and amplitudes in the Euclidean region, obtaining agreement with [4] on at least six significant digits in each case. Our analytical expressions were then also used for the evaluation of the five-gluon amplitudes of [4] in several physical kinematic points [48], which provides further opportunities for the validation of our pentagon functions. The correctness of the two-loop pole structure [49,50] for all helicity configurations provides a check of the integrals up to weight three; moreover, the cancellation of terms proportional to (D s − 2) 0 in the all-but-one-minus helicity configurations at order 0 provides a non-trivial consistency check of the weight-four contributions to the integrals.
Finally, point-wise checks of selected master integrals against numerical results using sector decomposition were performed with the FIESTA code [47] in the Euclidean region and with the SecDec [51] code in different Minkowskian regions.
Taken together, these checks provide strong support for the correctness of the analytical expressions that we obtained for the pentagon functions, and for their numerical implementation, whose usage is described below.

Numerical implementation
The two-loop planar pentagon functions and the master integrals constructed from them have been implemented as C++ functions, using their one-dimensional integral representation for the numerical evaluation. The code requires as prerequisites the installation of GiNaC [52] (and its dependencies) and of the GNU Scientific Library GSL [53]. The computation of the master integrals is performed in two stages, which correspond to two distinct functions.
The C++ function evaluating the linearly independent pentagon functions is defined in the file evaluate pentagonfunctions.cpp: evaluate pentagonfunctions(double v [5], double error) The arguments of this function are a vector of double precision floating-point numbers containing the five kinematic invariants v 1...5 and the relative size of the target error for the numerical integration. It is recommended to keep this target error well below machine precision. The function returns a 82-dimensional vector of double-precision complex numbers containing the values of all pentagon functions (including the cyclic permutations of the kinematic invariants).
The C++ function assembling the 61 independent master integrals (in the UT basis) and their cyclic permutations from these pentagon functions is defined in the file evaluate pentagonintegrals.cpp: evaluate pentagonintegrals(vector<complex<double> > eval pen) Here, the argument eval pen of the function is the result of a preceding function call to the above evaluate pentagonfunctions: the 82-dimensional complex double precision vector of the pentagon functions. The output is a nested vector of dimension 61 × 5 × 5 of complex double-precision numbers containing the values of the five cyclic permutations (third index) of the different coefficient of the the expansion in (second index) of the 61 master integrals in the UT basis basis2 (first index).
To illustrate the usage of these functions, and example program example-code.cpp is provided. This program evaluates all pentagon functions and all master integrals for a given kinematic point (either selected from predefined 16 example points, or entered manually as list of the v 1...5 ), and writes them to an output file. The target integration error is fixed to 10 −6 in this example code.
After potentially adjusting the pathnames for GiNaC and GSL in the makefile, the code can be complied by make and run by ./example-code.
The code is available as an ancillary file to the arXiv submission of this article. It will be part of a new www.hepforge.org repository PentagonFunctions, where also the GPL implementation and computer algebra expressions for all results from this article will be made available.
The master-integral basis basis2 is normalised such that its first element (two-loop massless sunrise with scale s 23 ) is (8.1)

Conclusion and discussion
In this paper, we gave details of the computation of all master integrals of massless, planar five-particle integrals at two loops [14]. The integral basis was chosen following the ideas of [23], and as a result, the analytic answers are in a simple and transparent form. In particular, all master integrals are pure functions of uniform transcendental weight in their expansion. They are expressed in terms of iterated integrals, with integration kernels given by the pentagon alphabet (4.3) that we identified in [14].
We classified all planar pentagon functions up to weight four. They are expressed in terms of a basis of 18 irreducible functions, and permutations thereof. Only four of these functions depend on the genuine pentagon kinematics. Next, we expressed all master integrals in terms of this basis of functions.
We used properties of the iterated integrals to provide efficient one-fold integral representations for all of the basis functions. For convenience of the user, we coded them in a publicly available program. We validated this code as described in section 7.
The main focus of the present paper was to provide the full analytic solution for the master integrals, and to provide a code for their fast and reliable numerical evaluation. Beyond this, we remark that the information provided here is very flexible, and can be used, for example, for obtaining asymptotic expansions in any desired limit. We refer the interested reader to ref. [54] for more details in a similar setting.
As a result of our analysis, all integrals needed for planar massless two-loop five-particle scattering amplitudes are available analytically, and can be evaluated numerically using the computer code provided with this publication.
The calculation of planar two-loop five-particle scattering amplitudes is already well underway [4,[13][14][15][16][17]. Our results have already been used in [4] for obtaining numerical results. We expect that in the future, they can be used for obtaining full analytic answers for all helicity configurations, extending the initial result of [14] for the all-plus helicity configuration.
A very interesting open problem are the non-planar massless five-particle integrals. First results [42,55,56] suggest a natural extension of the pentagon alphabet obtained in ref. [14]. It would be interesting to prove this conjecture, and to compute all master integrals, using the methods employed here.  Figure 6. Momentum twistor geometry of (a) five-particle non dual-conformal scattering amplitudes (5 + 2 twistors) and (b) seven-particle scattering with dual conformal symmetry.
To measure distances we also need an "infinity twistor", which we write as an auxiliary bitwistor I = Z 6 Z 7 . With these ingredients at hand, we can write The light-like conditions are encoded in momentum twistor space geometrically: the corresponding momentum twistor lines intersect, and we have, e.g. x 2 23 ∼ (2334) = 0. See Fig. 6.
It is interesting that there is a closely connected case that can also be described by seven twistors, namely seven-particle scattering in N = 4 super Yang-Mills. The latter theory has an additional symmetry, dual conformal symmetry, that implies that the infinity twistor does not appear. So, dual conformal seven-point functions in that theory are described by the data shown in Fig. 6 (b). It is clear that there are some differences between the two cases, e.g. due to different cyclic symmetry groups. On the other hand, it will be interesting to compare the functions found here to observations about seven-point functions in the N = 4 sYM literature.
Of course, the formulation above is a redundant formulation for just five independent variables v i . We can choose specific Z's to reduce or remove this redundancy. In the literature this is sometimes referred as 'gauge fixing'. The only constraint is that the parametrization we choose must be invertible. We choose 3 This leads to the following change of variables, The variable x 1 is dimensionful, whereas the four remaining x i , i = 2, . . . 5 are dimensionless. In other words, x 1 is a trivial overall scale. The inverse of this transformation has two branches. We chose the inversion for which the boundary point (B.P.) This parametrization has the advantage of making the Gram determinant a perfect square, leading to √ ∆ = − x 2 1 x 2 x 4 (x 5 − 1) + x 3 1 + x 2 x 5 + x 4 (−2 − x 2 + x 5 ) . (A.8) Here, the overall sign was chosen in such a way that √ ∆ is positive in the symmetric point (A.7).
The inverse relations admit two solutions, one of which is The second solution is obtained from the former by flipping the sign in front of the square roots.

A.3 Scattering equations parametrization
At five particles, the scattering equations [59] have two complex solutions [60] that one can use to parametrize the four scale-invariant ratios appearing at five particles. This leads to the following parametrization 4 Like in the other cases, these variables allow one to rationalize the square root √ ∆.