Pentagon Functions for Scattering of Five Massless Particles

We complete the analytic calculation of the full set of two-loop Feynman integrals required for computation of massless five-particle scattering amplitudes. We employ the method of canonical differential equations to construct a minimal basis set of transcendental functions, pentagon functions, which is sufficient to express all planar and nonplanar massless five-point two-loop Feynman integrals in the whole physical phase space. We find analytic expressions for pentagon functions which are manifestly free of unphysical branch cuts. We present a public library for numerical evaluation of pentagon functions suitable for immediate phenomenological applications.


Introduction
Scattering amplitudes are among the central objects of interest in quantum field theories (QFT). On the one hand, they are the building blocks for scattering cross sections, which are the crucial theoretical input for phenomenological studies of high-energy particle collisions. On the other hand, they exhibit intriguing mathematical properties which provide us an opportunity to understand fundamental structure of QFTs. The coupling constants of the Standard Model are small in the high-energy regime, which implies that the scattering amplitudes can be consistently approximated by their perturbative expansion. Beyond the leading order in the expansion, the amplitudes are represented as sums of Feynman integrals with increasing number of loops. Only one-loop integrals are known for arbitrary scattering processes [1][2][3]. Evaluation of Feynman integrals with two or more loops is an open problem and an active area of studies in theoretical physics and mathematics. While two-loop integrals for many 2 → 2 scattering processes have been already obtained (for a recent review see [4,5]), 2 → 3 processes are on the current frontier of research. Massless Feynman integrals play a special role in QFT. The most abundantly produced particles in hadron collisions are the partons of quantum chromodynamics (QCD): gluons and quarks. Both can be treated as massless at sufficiently high energies. On a formal side, mathematical structure of QFT is more transparent in the absence of (spontaneously) broken symmetries.
A large number of Feynman integrals contributing to scattering amplitudes can be reduced to a smaller set of master integrals with the help of integration-by-parts identities [6]. It is a formidable challenge for multi-scale processes, and a number of novel ideas and algorithms has been developed to tackle integral reduction of five-particle processes [7][8][9][10][11][12][13][14][15][16][17][18]. Thanks to the advances in integral reduction and functional reconstruction techniques [8,9,19], we have witnessed a tremendous progress in calculation of two-loop five-particle amplitudes. All planar five-point QCD helicity amplitudes have been obtained in [17,[20][21][22][23][24][25]. The first results for non-planar five-point amplitudes were obtained in N = 4 super-Yang-Mills theory [26,27] and in N = 8 supergravity [28,29], followed by the fullcolor five-gluon amplitude with all positive helicities [30]. Also the full-color six-gluon all-plus helicity amplitude was obtained in [31]. Important progress has been made in evaluation of five-point amplitudes and integrals with one massive leg [32][33][34]. The first cross section computation of a 2 → 3 process was carried out in [35], where the planar two-loop amplitudes for the qq → γγγ process were evaluated on a small set of phase space points to construct an interpolating function.
The master integrals for massless five-point scattering processes have been a subject of extensive studies in recent years. The method of differential equations (DE) [36][37][38][39][40] in their canonical form [41][42][43][44][45], and systematic understanding of the transcendental functions appearing in calculations of multi-scale Feynman integrals [46][47][48][49][50] proved to be indispensable to obtain analytic results for five-point massless master integrals for planar [51][52][53] and non-planar [27,30,[54][55][56] topologies. Differential equations in canonical form provide a natural framework for expressing master integrals in terms of functions of uniform transcendental (UT) weight order by order in the dimensional regulator. It is advantageous, both for analyzing analytic structure of scattering amplitudes and for their efficient nu-merical evaluations, to have a good grasp on the analytic understanding of the relevant space of transcendental function. Finding a minimal set of transcendental functions that is sufficient to express all master integrals is essential for deriving compact analytic representations of scattering amplitudes and studying their asymptotic behavior in singular limits (soft, collinear, high-energy, etc.). Successful applications of modern semi-numerical approaches to analytic reconstruction of amplitudes [8,9,17,21,23,25] rely to a large extent on the knowledge of this set. At the same time, the representation of amplitudes in terms of a minimal set of transcendental functions achieves the efficiency required in phenomenological applications, where they have to be numerically evaluated on huge samples of phase space points. In the context of five-particle scattering we refer to this basis set as pentagon functions.
For planar massless five-point integrals, a set of pentagon functions was constructed in [53] and a reference implementation of their numerical evaluation was provided. For the scattering processes involving solely QCD partons, only planar Feynman integrals contribute in the leading-color approximation. Obtaining complete NNLO predictions and assessing the accuracy of this approximation requires calculation of amplitudes involving non-planar Feynman integrals. The scattering processes with photons in the final state also involve non-planar contributions originating from closed fermion loops, which in general cannot be considered subleading. The knowledge of the full set of pentagon functions is vital for obtaining scattering amplitudes for these classes of processes. In this work, we construct a basis set of pentagon functions that is sufficient to express all planar and nonplanar massless five-point two-loop master integrals in any physical scattering channel. We consider series expansion of the master integrals in the dimensional regulator up to the order sufficient for calculation of two-loop hard functions and next-to-next-to-leading order (NNLO) cross sections. The pentagon functions manifestly posses only physical branch cuts. In particular, they are branch-cut-free in the whole physical phase space and do not require analytic continuation. We find explicit representations of pentagon functions which admit efficient and stable numerical evaluations, and we implement the latter in a C++ library (section 7.2). In addition to extending the set of pentagon functions to the non-planar sector, at the same time we reconsider the analysis of the planar sector carried out in [53]. The planar subset of pentagon functions presented in the current work is explicitly closed under permutations of external momenta and involves a much smaller set of transcendental constants. Furthermore, the numerical evaluation of the pentagon functions is significantly improved both in speed and precision. Thus, for the first time, we provide an implementation that is immediately applicable to computations of NNLO cross sections of any scattering process involving five massless particles.
To find a minimal set of pentagon functions, we follow a constructive approach, which relies almost entirely on the information contained in the canonical differential equations [27,[53][54][55][56][57][58]. We consider the DEs for the planar pentagon-box, hexagon-box, double pentagon, and the one-loop pentagon integral topologies (see section 3) in all 5! permutations of external momenta. We solve each DE in terms of iterated integrals [49] with an initial point X 0 in the physical scattering region (section 4). To completely fix the solutions of DEs one needs to provide the initial values -values of the master integrals at X 0 . Building on the results of [30,59], we obtain a complete set of the initial values of all DEs at X 0 from the requirement of absence of unphysical singularities and identify a generating set of 19 algebraically-independent transcendental constants. We then employ the shuffle algebra of iterated integrals (see e.g. [48]) to find a set of linear-independent irreducible iterated integrals up to transcendental weight four in section 5. We evaluate the iterated integrations up to weight two in terms of logarithms and dilogarithms, and we derive onefold integral representations [53,60] for the iterated integrals of weight three and weight four. In this way, we find expressions for all master integrals in any scattering channel sidestepping a difficult problem of analytic continuation. The obtained analytic expressions allow us to perform a detailed analysis of their behavior in singular limits. As an example, in section 6 we investigate the behavior of pentagon functions on boundaries of the physical phase space where all five momenta belong to a three-dimensional subspace, but none of the external momenta are soft or collinear. Confirming the observation of [59], we find that certain weight three and weight four pentagon functions contributing to non-planar master integrals are divergent on these boundaries.
All results of the paper are made available through data files and can be explored with the Mathematica package presented in section 7. We elaborate on the implementation details of pentagon function numerical evaluation by the C++ library and demonstrate its performance. In section 8, we discuss validation of our results, and we conclude in section 9.

Kinematics
We study the scattering of five massless particles in four-dimensional Minkowski space-time. The particles momenta p i are subject to momentum conservation 5 i=1 p i = 0, and on-shell conditions p 2 i = 0. We parametrize points X of the physical phase space as where s ij := (p i + p j ) 2 are the Mandelstam invariants, and ε(·, ·, ·, ·) is the fully antisymmetric Levi-Civita symbol. The parity-odd invariant 5 is related to the determinant of the Gram matrix G(p 1 , p 2 , p 3 , p 4 ) as In the physical phase space ∆ < 0 [61], so it is convenient to define such that δ ∈ R. It is worth noting that although |δ| is not algebraically independent from v i , the sign of δ is necessary to fully specify a point in the physical phase space. Depending on the problem at hand, one can choose different parametrizations of the scattering kinematics (see e.g. appendix of [53]). Our choice of the parametrization in eq. (2.1) is motivated by the fact that X transforms linearly upon permutations of momenta p i .

Integral Topologies
We consider all Feynman integral topologies required in computation of two-loop scattering amplitudes of five massless particles. There are four topologies, see fig. 1. Their legs are decorated with particle labels what we call the topology permutation. To regularize the divergences of loop integrals we employ dimensional regularization and extend the integration measure to D = 4 − 2 dimensions. We define the integral families G τ,σ for each topology τ in permutation σ as where γ E is the Euler-Mascheroni constant 1 , L τ is a number of loops in topology τ , D τ,σ is an ordered set of inverse propagators of integral topology τ in permutation σ, the exponents a i ∈ Z for i ∈ [1,8] and a i ∈ Z ≤0 for i ∈ [9,11], and µ 2 is an arbitrary regularization scale, which preserves the integer dimensions of the integrals. In this paper, we choose the units of energy such that µ = 1. The explicit dependence on the regularization scale can then be restored by the dimensional analysis. For each of the four integral topologies we choose the standard permutation σ 0 = (1,2,3,4,5), and define the sets D τ,σ 0 as This choice is illustrated in fig. 1, where we show the top topology (the topology with maximal number of denominators) for each integral family in the standard permutation σ 0 . The denominator-variable sets in other permutations σ = (σ 1 , . . . , σ 5 ) ∈ S 5 are generated by the action of the symmetric group S 5 on the set of external momenta p i , The Feynman integrals from each family G τ,σ in eq. (3.1) form a linear vector space. For each G τ,σ (separately) we choose a set of basis elements, which are independent under the linear relations generated from integration-by-parts identities [6]. We refer to these sets as master integrals. One can further decrease the number of master integrals in ∪ τ,σ G τ,σ by identifying integrals among different topologies and permutations. However, as we explain in section 5, we find that it is more convenient to resolve these relations together with the functional relations while constructing a basis of transcendental functions.
A choice of a basis in the vector space G τ,σ is in general arbitrary. Frequently it is specified by an ordering relation on the set of exponents a i in eq. (3.1) [62] (see also [63,64]). However, it was observed in [41] that certain integrals have particularly nice properties. Following the approach of [54,56], we choose the master integrals with constant leading singularities in D dimensions. The differential equation for such integrals can be cast into the canonical form [41], and the integrals can be expressed as Q-linear combinations of pure functions with uniform transcendentality (UT). In the following we call them UT master integrals. We discuss the construction and the solution of the differential equations in the next section.

Differential Equations
To obtain analytic expressions for all master integrals from integral families defined in eq. (3.1), we construct the corresponding differential equations [36][37][38][39][40] in the canonical form [41]. The differential equations (DE) for all integral topologies in fig. 1 have been extensively studied in literature. The canonical form of differential equations has been obtained in [27,[53][54][55][56]. The sub-topologies of fig. 1 with less than five external momenta were also studied in [57,58].
For the double-pentagon topology, we directly use the canonical DE of [56]. For hexagon-box, planar pentagon-box, and one-loop pentagon topologies we repeat the analysis of [56] to find master integrals with unit leading singularities in D-dimensions. Their four-dimensional integrands have d log form.
In this section we provide details of the construction and integration of DEs, which are necessary for the construction of a basis of transcendental functions in section 5.1.

Construction of canonical differential equations
We would like to find the analytic expressions for all integral topologies in fig. 1 in all 5! permutations. One can think of two different approaches to constructing DE solutions for all permutations. In the first approach, one would consider a single permutation of each topology, e.g. the one depicted in fig. 1, solve it analytically by the method of differential equations, and then obtain analytic expressions for all other permutations of the topology by means of analytic continuation. The latter is a highly nontrivial task for integrals depending on many scales. In particular, some of the non-planar master integrals develop discontinuities and even divergences inside the physical region on subvarities of ∆ = 0 without collinear or soft momenta [59]; we discuss this in section 6. In this paper we follow an alternative approach which was advocated in [53]. We work simultaneously with all 5! permutations of each topology in fig. 1 and consider canonical differential equations (DE) for each of them where f τ,σ is a vector of UT master integrals of topology τ taken in permutation σ. Entries of the matrices a (i) τ,σ are rational constants, and {W i } 31 i=1 are letters of the pentagon alphabet [65]. The letters are algebraic functions of the Mandelstam variables. We review the pentagon alphabet in appendix A. The matrices of the DE in permutation σ are related to the DE in the standard permutation σ 0 as follows, where σ permutes the external momenta according to (3.3). The pentagon alphabet is closed under S 5 permutations of the external momenta, and the letters of the alphabet have where the action of σ ∈ S 5 on the kinematic point X (see eq. (2.1)) is induced from the action of σ on momenta p i by eq. (3.3). In addition, for arbitrary η ∈ S 5 , the following relation holds, The UT master integrals in the standard permutation are related to the integrals from (3.1) by linear transformations (4.5) The UT master integrals f τ of our basis are split in the parity-even f The number of parity-even and parity-odd master integrals (i.e. dimensions of f (±) τ ) are given in table 1. We provide explicit transformations (4.5) in ancillary files (see section 7.1).

Integrating DE and iterated integrals
To integrate DEs (4.1) order-by-order in , we define -expansions of the UT master integrals as τ,σ are of uniform transcendental weight w. The -expansion of the two-loop master integrals starts with 1 4 pole, and with 1 2 pole for one-loop pentagons, i.e. the soft-collinear pole 1 2 per loop order. We omit in the following the topology and permutation labels τ, σ to avoid bulky notations.
Let us denote a point of the kinematic space by X (2.1). We specify it by the set of five adjacent Mandelstam invariants and the sign of 5 . We choose an initial point X 0 and integrate the DE along a path γ connecting X 0 and X. Thus we express weight-w solutions at an arbitrary kinematic point X as iterated integrals of the initial values f (w ) (X 0 ), i.e. solutions with w ≤ w evaluated at X = X 0 , At weight 0 the previous equation simplifies to f (0) (X) = f (0) (X 0 ), i.e. it is a constant vector of rational numbers. Moreover, the vector is the same for all permutations of the given topology, Equation (4.8) can be rewritten explicitly as a linear combination of iterated integrals built upon the pentagon alphabet The coefficients κ (w−w ) of the linear combination are transcendental constants of weight w − w , Chen iterated integrals [49] of weight w along the path γ are defined recursively as with [] X 0 = 1. The iterated integrals vanish at the initial point X = X 0 by construction, The DE guaranties that only homotopy invariant, i.e. invariant under small deformations of the integration contour γ, linear combinations of the iterated integrals are present in the solution (4.10). The iterated integrals are in general multi-valued functions since they pick up a nontrivial monodromy upon integrating around a pole, i.e. zero locus of an alphabet letter. Thus we have to specify an analyticity region P 0 within which the iterated integrals are single-valued and real-analytic functions. Then the result of the iterated integration depends only on the end points X and X 0 of the integration path γ ⊂ P 0 .
The iterated integral representation is a powerful tool which enables us to classify the solutions of the DE within the analyticity domain by doing simple algebraic calculations. In particular, the iterated integrals satisfy the shuffle algebra relations (see e.g. [48]), which specify how to rewrite a product of several iterated integrals as a sum of iterated integrals, where we sum over all {k 1 , . . . , k w 1 +w 2 } in the shuffle product {i 1 , . . . , i w 1 } ¡ {j 1 , . . . , j w 2 }.
After applying the shuffle algebra relations, all polynomial identities among the functions represented by iterated integrals become linear, and only the trivial combination of the iterated integrals vanishes. In this way we take into account the functional relations among the DE solutions.

Physical region
As we have already mentioned, the iterated integrals (4.10) are not single-valued and realanalytic at any kinematic point X. We choose the following analyticity domain  It is a half of the physical s 12 -channel scattering region, i.e. 12 → 345 scattering process. Fixing signs of the Mandelstam invariants implies that the particle energies are positive and scattering angles are real. ∆ < 0 implies the reality of momenta (see also eqs. (2.1) and (2.2)). In addition, we also fix the branch of the square root √ ∆ by the condition δ > 0. The boundaries of P 0 corresponding to vanishing of one or several s ij describe the soft/collinear limits. The boundary ∆ = 0 of P 0 lies inside the physical s 12 -channel and splits it into two halves. It corresponds to the kinematics with all five momenta lying in a three-dimensional hyperplane. Crossing of the ∆ = 0 variety separating Im( 5 ) > 0 and Im( 5 ) < 0 regions is not innocuous since the master integrals could diverge there [59].
In the following we work strictly inside P 0 (4.15) and classify the solutions of the DE (4.10) only for X ∈ P 0 . Since we consider all 5! permutations of each topology, we can immediately translate our results to any physical scattering region. This, we provide analytic expressions for the master integrals taken in arbitrary permutation (3.1) through the whole phase space, see section 5.8.
Given an arbitrary point X ∈ P 0 , we evaluate the iterated integrals by choosing the path γ to be a line segment connecting X 0 and X, eq. (2.1). We parametrize the segment If we are to avoid the problem of analytic continuation, the integration path γ must never leave the analyticity domain P 0 eq. (4.15), i.e. for any X ∈ P 0 , {X(t)} 0≤t≤1 ⊂ P 0 must be satisfied. To this end, we note that P 0 is not a convex 2 . Nevertheless, a weaker statement holds: a line segment connecting X 0 (4.16) with an arbitrary point X ∈ P 0 lies entirely inside P 0 . We outline the proof in appendix B. Thus, integrating in (4.12) along the straight lines connecting X 0 with any X ∈ P 0 we obtain real-analytic single-valued solutions (4.10) throughout P 0 .

Initial values
In order to be able to integrate DEs (4.1), we need to know -expansion of all UT master integrals at the reference point X 0 (4.16) -the initial values of the DEs. As we pointed out in section 4.1, we would like to trade the problem of analytic continuation of the UT master integrals in the standard permutation σ 0 to all possible permutations, for solving the DEs with initial values at X 0 in all 5! permutations: { f τ,σ 0 (X 0 )} + analytic continuation of f τ,σ 0 (X) to any physical s ij -channel We restrict our consideration to weight w ≤ 4 initial values, i.e. we truncate theexpansion (4.7a) of the two-loop topologies at the finite part, and at O( 2 ) for the one-loop pentagon (4.7b).
Initial values of the DE for two-loop five-point topologies have been extensively studied previously. Weight-0 initial values f (0) are rational numbers. They are enough to construct 2 For example, let us consider the following line segment γ parametrized by 0 ≤ t ≤ 1, We find ∆(X(t = 0)) = − 87 16 < 0 and ∆(X(t = 1)) = − 231 100 < 0 so the end points of γ belong to P0 (4.15). However, at the intermediate point ∆(X(t = 1 2 )) = 1323 50 > 0, and the segment does not lie inside the physical s12-channel. symbols [46][47][48] of the UT master integrals, and are relatively easy to obtain. Calculation of higher weight initial values is much more tedious. The planar pentagon-box topology ( fig. 1b) has been solved in any physical region in [53]. In [54] the initial values for one permutation of the hexagon-box topology ( fig. 1c) were evaluated in the Euclidean region. In [56] the initial values in a physical region for one permutation of the double pentagon topology ( fig. 1d) were presented with 50 digits precision. In calculation of a five-point nonplanar amplitude in [30], the initial values at X 0 were computed for all permutations of all four topologies in fig. 1 with 200 digit precision, but they were not explicitly reported. The initial values have been fixed by requiring the absence of the unphysical singularities in the DE solutions (4.8), see [53,54,66,67], and special care have been taken owing to the singular behaviour of the nonplanar Feynman integrals at ∆ = 0. In the present paper, we publish for the first time the complete set of initial values at X 0 (all permutations σ of all four topologies τ ).
In [59], by integrating the DEs, the initial values at X 0 were transported to a point X Regge ∈ P 0 in the Regge asymptotic regime. In this regime the pentagon alphabet enormously simplifies, which leads to a more simple form of the initial values at X Regge as compared to X 0 . The available numerical precision is enough for fitting to a basis of transcendental constants. A small generating set S Regge (see table 2 in [59]) of algebraically independent over Q transcendental constants was identified, and all initial values at X Regge were written as polynomials Q[S Regge ] graded by the transcendentality degree. These analytic expressions for the initial values at X Regge were transported back to X 0 (4.16) by integrating the DEs (4.1) in terms of multiple polylogarithms (MPLs) [46,68] and evaluating the latter with GiNaC [69] with 10 4 digit precision. In this way the initial values at X 0 have been found with at least 9 · 10 3 digit precision. This precision is enough to identify the generating set S 0 of algebraically independent over Q transcendental constants, see table 2, and to fit the initial values at X 0 to graded polynomials Q[S 0 ]. The analytic form of the initial values can be found in the data files supplied with the Mathematica package (see section 7.1).
The set S 0 consists of only 19 transcendental constants, which we classify as follows. We assign the Z ≥0 × Z 2 -charge (or grading) to the constants where the first factor refers to the transcendentality degree while the second factor Z 2 = {+, −} = {even, odd} counts parity. Then the weight-w initial values f (w,+) τ,σ (X 0 ) of the parity-even master integrals are homogeneous polynomials Q[S 0 ] (w,+) , while the initial values f (w,−) τ,σ (X 0 ) of the parity-odd master integrals are homogeneous polynomials Q[S 0 ] (w,−) . In order to be able to consistently assign the parity to the initial values we have to introduce two copies of i π of opposite parity, i.e. parity-odd i π and parity-even i • π. For example, log(3) + i • π carries charge (1, +) and it is an admissible weight-one initial value of a parity-even UT master integral; π • π + i Im Li 2 e i π 3 carries charge (2, −) and it is an admissible weight-two initial value of a parity-odd UT master integral. Obviously, i π and i • π are numerically identical, and we are allowed to identify π 2 = ( • π) 2 . We notice that all parity-odd constants are pure imaginary, and all parity-even constants are real (except for i • π). The reality properties of the transcendental constants imply that the initial-values of the parity-even master integrals Weight even (+) Table 2: The basis S 0 of algebraically independent over Q transcendental constants specifiying the initial values { f If there existed an Euclidean region for all master integrals G τ [ a] (3.1) where they take real values, then both parity-even and parity-odd UT master integrals would also be real in that region. Then, by analytic continuation from the Euclidean region to a physical scattering region the parity-odd integrals f (−) τ would become imaginary up to the contributions from discontinuities, which are proportional to i π. This would be in agreement with eq. (4.19). However, non-planar integrals from the double-pentagon topology ( fig. 1d) do not have an Euclidean region and they are complex-valued everywhere. Nevertheless, we find the observed correspondence between the reality properties of the initial values and the parity of the UT master integrals very intriguing.

Parity of the UT master integrals
As we discussed in the previous section, the initial values at X 0 , or equivalently the transcendental numbers κ in eq. (4.11), obey Z ≥0 × Z 2 -grading. The same is true for the iterated integrals (4.12). Indeed, counting of the parity-odd letters (see appendix A) in the iterated integral [W i 1 , . . . , W in ] corresponds to the Z 2 -grading, and the number of iterated integrations (weight) corresponds to Z ≥0 -grading. Both gradings are compatible with the shuffle algebra (4.14). It then follows that the UT master integrals f (w) τ,σ (X) inherit the Z 2 -grading of the initial values and iterated integrals. This fact allows us to establish an equivalence between the Z 2 -grading and the parity of the UT master integrals f (±) τ,σ implied by the parity-conjugation properties of the coefficients T τ in definition (4.5).
We would like to emphasize that this compatibility is not trivial. Indeed, had we not introduced two copies of i π to represent the initial values, the equivalence would not hold. Moreover, if taken literally, the parity conjugation maps the initial point X 0 (4.16) and the path γ of an iterated integral (4.12) out of the chosen analyticity region P 0 into a region with δ < 0 in addition to parity conjugation of the d log-kernels (see appendix A). The relations between parity-conjugated iterated integrals could then be established through analytic continuation. As discussed in the previous section, our strategy is to avoid analytic continuation in favor of considering momenta permutations of the master integrals, so we do not pursue this approach in what follows.

Classification of Functions
Upon integration of DEs (4.1) for all four topologies τ , each in 5! permutations, we obtain a number of iterated integrals (4.10) which are not independent. We would like to reduce them to a minimal set of functions which are sufficient to express all solutions f (w) τ,σ up to weight w ≤ 4. Before we delve into the classification procedure, it is worth noting that from (11+61+73+108)×5! UT master integrals involved in DEs (4.1) only 1865 two-loop and 52 one-loop UT master integrals are linear-independent under the topological identifications among different topologies and permutations of integrals G τ,σ [ a] (see eq. (3.1)) [30]. These relations in the set of UT master integrals are trivialized by solving their DEs in terms of iterated integrals. For this reason, we do not explicitly implement them in our classification. In this section, we find the linear-independent solutions at weights w ≤ 4 and show that their number is smaller than the one obtained from the topological analysis. 3 We then further reduce the set of linear-independent iterated integrals to a smaller set of irreducible iterated integrals, i.e. the ones which cannot be represented as products of lower-weight iterated integrals. We claim that the latter set is minimal, and we denote it as pentagon functions.

Classification strategy
Let us briefly outline our classification strategy. We will assume that the solutions vanish at X = X 0 , or in other words we consider f τ,σ (X 0 ). We proceed recursively in weight. At weights 1 ≤ w 1 < w we have already identified {I . We rewrite them schematically in the following form splitting out the term I (w) with the maximal number of iterated integrations, (5.1) In other words, are transcendental constants of weight w − w and R (w ) represents Q-linear combination of weight-w iterated integrals. We mod out lower-weight iterated integrals, i.e. we apply the symbol map [46][47][48] and choose the subset of M w ≤ N w linear independent 3 It is expected that this redundancy should be lifted by considering higher orders in -expansion (4.7a). Nw }. Then we need to eliminate reducible iterated integrals from (5.2), i.e. iterated integrals which are products of lower weight iterated integrals. In order to achieve it, we consider symbols of weight-w products of lower weight iterated integrals, {I with w 1 < w, already classified at the previous steps of our procedure, e.g.
which are linear independent by induction. Using the shuffle algebra (4.14) we rewrite them as sums and complement them with the symbols (5.2). The resulting set of symbols is overcomplete. Then we choose a basis in their Q-linear span. We include the maximal possible number of the products (5.3) in the basis and complement them by a subset of (5.2), {I The linear span of the subset does not contain products of lower weights, i.e. it is irreducible. Let us now relabel the iterated integrals by 1, . . . , L w . Complementing the symbols to the complete solutions of the DEs by means of (5.1) we would like to choose as an irreducible set of iterated integrals at weight w. However, it could happen that not all weight-w solutions of the DE are expressible in terms of (5.4) and products of already classified lower-weight iterated integrals, {I We could encounter a solution J of the DE which is expressible in the constructed basis at the symbol level, but it also contains "beyond-the-symbol" terms,

products of lower weights
Here b k , b n,m , . . . are rational numbers, andR (w ) a is a Q-linear span of weight-w with w < w iterated integrals which have not been included in the weight-w basis of iterated integrals at the previous steps of the classification. Then to classify weight-w iterated integrals we would need to reconsider the classification of all lower weights. Fortunately, this complication can be very easily resolved. We find that it is sufficient to only extend the set of weight-1 iterated integrals {I . These extra integrals do not appear in the weight-1 solutions f (1) τ,σ (see eq. (5.13)). But their powers and their products with {I terms in eq. (5.5). As a result of this classification, we find that we need 24 weight-2 functions, 111 weight-3 functions and 472 weight-4 functions to express all UT master integrals for all four topologies in all orientations (up to weight 4), see table 3. This counting of pentagon functions should be compared with the counting of integrable symbols summarized in table 1 from [65]. Predictably, integrating the DEs we find considerably fewer solutions than one obtains by imposing the second-entry restriction, inspired by the Steinmann relations [70,71], on the space of integrable symbols. In the rest of this section we elaborate the details of the pentagon function classification.  Table 3: Number of the parity even|odd independent iterated-integral solutions of DEs (4.1) for all four topologies and all 5! orientations combined.

Parity-even letters of the alphabet in the analyticity region
The iterated integrals (4.12) involve d log kernels d log(W i ), i = 1, . . . , 31, so while implementing integrations in the region P 0 (4.15) we need to keep track of possible singularities of d log(W i ). Let us consider first the parity-even letters of the alphabet (A.2). The parity-odd ones are discussed in section 5.3.2.
Most of the parity-even letters have definitive sign within P 0 , and W 31 = 5 ≡ i δ with δ > 0 inside P 0 . The corresponding d log(W i ) kernels are realanalytic inside P 0 and integration along a path γ, γ ⊂ P 0 , is well defined. Missing in eq. (5.6) are the parity-even letters which all vanish at the initial point X 0 (4.16). Since they are linear in the Mandelstam invariants, along any line segment γ = [X 0 ; X] parametrized by 0 ≤ t ≤ 1 (see eq. (4.17)), In the latter case, the d log kernel has a simple pole at t = 0, and we should verify that the corresponding integration in eq. (4.12) is well-defined.

Weight-one solutions
Weight-1 solutions of DEs (4.1) have a very simple form. According to (4.8) The previous equation involves only letters (4.9)) is annihilated by the components of dÃ τ,σ (4.1) corresponding to the remaining letters, The one-fold integrals [W i ] X 0 (X) from (5.9) are calculated straightforwardly, see (4.12) and (4.16). For example, Tab. 2). In fact, the transcendental constant log (3), which comes from the one-fold integrations and from weight-1 initial values, cancels out in (5.9).
The spurious transcendental constants related to the specific choice of X 0 (4.16) is one of the reasons why we prefer not to evaluate iterated integrals as in (5.11). Instead, we introduce a set of the following ten functions {g 1,1 = log(s 12 ) , g 1,2 = log(−s 23 ) , g 1,3 = log(s 34 ) , g 1,4 = log(s 45 ) , g 1,5 = log(−s 15 ) , g 1,6 = log(−s 13 ) , g 1,9 = log(−s 14 ) , g 1,10 = log(−s 25 ) . (5.12) They are well-defined in the analyticity region P 0 (4.15). The arguments of logarithms in (5.12) are equal up to a sign to the letters , and they are listed in the first and the third columns of (5.6). We can assign even parity to {g Then we represent functions (5.12) as one-fold integrals with the initial point X 0 , resolve the one-fold integrals in terms of functions (5.12), and substitute the former in (5.9). In this way, we find weight-1 solutions f (1,±) τ,σ (5.9) in the parity-even and parity-odd sectors where b, c are vectors of rational numbers. Of course, at weight 1 the two approaches are completely equivalent, but at higher weights the second one is more practical.

Extra weight-one functions
As we have noted at the end of section 5.1, we need to supplement functions (5.12), appearing in the weight-1 solutions (5.13), by some extra weight-1 functions that are needed to describe higher weight solutions of the DEs. Let us start with the parity-even functions. We define ten functions which are logarithms of the remaining arguments from (5.6) 2,2 = log(−s 13 − s 14 ) , g 2,3 = log(s 45 − s 13 ) , g 2,5 = log(s 35 + s 45 ) , g 2,10 = log(W 31 ) = log( i δ) . (5.14) They are well-defined everywhere inside P 0 (4.15). Let us note that we do not introduce logarithms of letters (5.7) which do not have definitive signs in P 0 . We assign positive parity to {g . We will also need parity-odd weight-1 functions. They are one-fold integrals (4.12) of the parity-odd letters {W i } 30 i=26 along a path γ connecting X 0 and X ∈ P 0 , In the next section, we explain that they are well-defined and single-valued within P 0 (4.16) and provide explicit expressions (5.21), (5.22) for them.
Thus, −2π < ϕ k < 0 at k = 1, 2, 3, and 0 < ϕ k < 2π at k = 4, 5. We need to match them with the principal branch of the logarithm. We cross the branch cut and go to another Riemann sheet of the logarithm only at a k = 0, see (5.20). If we go from the region a k > 0 to the region a k < 0, then we decrease the phase ϕ k , and we should add −2 i π to the principal value of the logarithm. If we go from the region a k < 0 to the region a k > 0, then we increase the phase ϕ k , so we should add +2 i π to the principal value of the logarithm. We illustrate this in fig. 2 Thus we obtain the following following continuous and real-analytic expressions for integrals (5.15) inside P 0 (4.15): for k = 1, 2, 3, and g 3,k (X) = θ(−a k ) log W 25+k + θ(a k ) (log W 25+k + 2 i π) + i π δ 0a k − 2 i π 3 (5.22) for k = 4, 5. Here δ ij is the Kronecker delta, and the function θ(x) is defined as

Weight-two solutions
According to (4.8), weight-2 solutions of DEs (4.1) can be represented as The first term in the previous equation corresponds to the symbol of the solution. Instead of evaluating one-fold and two-fold iterated integrals from (5.24), and then looking for cancellation of the spurious transcendental constants among the three terms in (5.24), we prefer to start with a set of 15 parity-even and 9 parity-odd weight-2 functions and to express (5.24) in terms of them.

All master UT integrals at weight two
Now we have all necessary ingredients to express weight-2 solutions f τ,σ . We transform 15 parity-even (5.25) and 9 parity-odd (5.28) functions into the iterated-integral representation which involves one-fold [W i ] X 0 (X) and two-fold [W i , W j ] X 0 (X) iterated integrations. We also transform all weight-1 functions in eqs. (5.12), (5.14) and (5.15) into the iterated integral representation, i.e. one-fold integrations. Then we resolve the iterated integrals [W i ] X 0 and [W i , W j ] X 0 for the weight-1 and weight-2 functions, and, by substituting these relations into (5.24), we express all weight-2 solutions in terms of the functions g.
The parity-even and parity-odd solutions, defined in section 4.5, involve different subsets of functions. We find that the parity-even solutions have the following form 1,i + where b, c, and d are vectors of rational numbers. These expressions respect Z ≥0 × Z 2grading. As we can see upon identifying two copies of i π the only transcendental constants in the solution are i π and π 2 . In fact, not all extra weight-1 functions (5.14), (5.15) are present. The weight-2 solutions contain only parity even g 2,k with k = 2, 3, 6, 7, 8, 9 and parity odd g 3,k with k = 1, 3, 4, 5.
As we can see, the weight-1 functions in eqs. (5.14) and (5.15) that are not needed to express the weight-1 solutions (5.13) start appearing in the weight-2 solutions eqs. (5.32a) and (5.32b) in the "beyond-the-symbol" part. This phenomenon illustrates the statements from section 5.1.

Weight-three pentagon functions
We follow the classification procedure from section 5.1 and find 90 parity-even and 21 parity-odd weight-3 irreducible iterated integrals. Let us recall that they have the form of eq. (5.1) with w = 3. At weights one and two we started with the nice choices of logarithmic and dilogarithmic pentagon functions and expressed all iterated integrals in terms of these functions. It was argued in [60] that at higher weights explicit representations of iterated integrals in terms of MPLs are not always beneficial, especially for the purpose of their efficient numerical evaluation. We find that this also applies to the iterated integrals studied in this paper (see section 5.7). We then choose a set {g i=1 of irreducible iterated integrals as our weight-3 pentagon functions. Taken together with already classified iterated integrals at weights one and two, they are sufficient to express any weight-3 solution f i } 111 i=1 take the form of a one fold integral along a path γ connecting X 0 (4.16) and X ∈ P 0 , where h (2,±) i,j are weight-2 polynomials of definite parity, 2,i } i=2,3,6,7,8,9 (2,+) , (5.35b) We must verify that the integrations in eq. (5.33) are well-defined. As in eq. (4.17), we choose the path γ as the line segment γ = [X 0 ; X]. The functions h (2) i,j are well-defined on γ, since they are polynomials of weight-1 and weight-2 pentagon functions. The only possible source of potential problems is d log W k with k ∈ α := {7, 10, 12, 21, 22, 23}, i.e. d log-forms of the letters that do not have a definite sign in the analyticity region P 0 (see eq. (5.7)). As we showed in eq. (5.8), on the path X(t) either they are d log(t) with a pole at t = 0 or they are identically zero. Fortunately, the pole is compensated by the accompanying h (2) i,k which vanishes at t = 0. Indeed, inspecting h (2) i,k for k ∈ α we find that they involve very simple pentagon functions: {g (1) 1,i } 10 i=1 (eq. (5.12)) and g 1,i for i = 4, 7, 8, 11, 14, 15 (eq. (5.25)). Resolving W k = 0 (A.2) as a constraint on {v i } 5 i=1 and evaluating h (2) i,k on this subspace we find that it vanishes, Thus the integrations in the definition of the weight-3 pentagon functions (5.33) are welldefined for any point of the analyticity domain P 0 (see eq. (4.15)).

Weight-four pentagon functions
At weight 4 the classification procedure from section 5.1 results in 316 parity-even and 156 parity-odd irreducible iterated integrals having the form of eq. (5.1) with w = 4. As at weight 3, we choose a set of 472 irreducible iterated integrals {g (4) i } 472 i=1 as our weight 4 pentagon functions. To bring them into a more explicit form, we use the definition in eq. (4.12) and perform the two innermost integrations. The functions g (4) i are then expressed as two-fold iterated integrals over functions g (1) i,j and g (2) i,j introduced above, where the weight-2 functions h (2) i,j,k have definite parity, and they are of the same form as in eq. (5.34). The transcendental constants κ (3) ∈ Q[S 0 ] (3,±) are given in table 2. Equation (5.38) respects Z ≥0 × Z 2 -grading, i.e. the transcendental weight and parity counting. In order to render the weight-4 pentagon functions (5.38) to a form better adapted for numerical evaluations, we rewrite them as one-fold integrals in the next section.

One-fold integral representation of weight-four pentagon functions
We apply the technique from [60] to rewrite the two-fold iterated integrals in (5.38) into one-fold integrals. We introduce parametrization (4.17) of the path γ which we choose as the line segment, γ = [X 0 ; X], such that X(1) = X and X(0) = X 0 , and interchange the order of integrations where W i (t) := W i (X(t)). Thus, one of the integrations in eq. (5.39) becomes trivial, and naively we obtain 1 u dt ∂ t log(W j (t)) = log(W j (u = 1)) − log(W j (u)) .

(5.40)
However, the right-hand-side of (5.40) should be well-defined for 0 ≤ u ≤ 1. It is straightforward to express it in terms of the weight-1 pentagon functions (eqs. (5.12), (5.14) and (5.15)). Indeed, for the parity-odd letters we use eqs. (5.21) and (5.22), 3,j−25 (X(1)) − g 3,j−25 (X(u)) , j = 26, . . . , 30 . Let us note that here we have to consider letters (5.7), i.e. j ∈ α := {7, 10, 12, 21, 22, 23}, which do not have definitive sign in P 0 . Still, they do have definitive sign on the line segment γ, see (5.8), and they vanish at u = 0. Unlike other letters, which do not vanish inside P 0 , they introduce logarithmic singularity at u = 0, and thus they are not among pentagon functions (5.12) and (5.14). Finally, for the remaining parity-even letter W 31 , 1 u dt ∂ t log(W 31 (t)) = log(δ(1)) − log(δ(u)) . (5.43) In this away, we arrive at the one-fold integral representation of the weight-4 pentagon functions (5.38), We need to verify that integrations in (5.44) are well-defined. The analysis is similar to the one from section 5.5.1. The weight-2 functions h (2) are polynomials in the pentagon functions as in (5.34), and they are real-analytic. A potential source of pole singularities in the integrand is ∂ u log(W k (u)) at k ∈ α. We find that the pole is suppressed since Summarizing, we find that the integrations in eq. (5.44) are well-defined. All terms of the integrands are real-analytic at 0 ≤ u ≤ 1, except for the terms with j ∈ α, which are real-analytic at 0 < u ≤ 1. The only singularity of the integrands is the logarithmic singularity log(u) in the terms with j ∈ α. This singularity is integrable. Nevertheless, some care should be taken in an algorithm for numerical evaluations, see section 7.2 for details.

All master UT integrals at weight four
We express all weight-4 solutions (4.8) of the DEs as homogeneous polynomials of definite parity in the pentagon functions of weight one (5.12), (5.14), (5.15), weight two (5.25), (5.28), weight three (5.33), and weight four (5.38), and in transcendental constants from table 2, and we find that the parity-even solutions have the following form f (4,+) τ,σ ∈ Q g (4,+) , g (3,+) , {g 2,10 , {g Thus, we have classified all pentagon functions up to transcendental weight four, and we identified the minimal generating set in the pentagon function space. All constructed pentagon functions are well-defined within the physical region P 0 (4.15). We also provided -expansion of all two-loop UT master integrals that describe the massless five-particle scattering in terms of the pentagon functions.

Alternative representation of the pentagon functions
We provided expressions for the pentagon functions in terms of the familiar polylogarithmic functions only at weights one and two, i.e. logarithms and dilogarithms, respectively. We preferred to express the pentagon functions of higher weights as one-fold integrations, see (5.33) and (5.44). This approach provides a convenient setup for numerical evaluations of the pentagon functions, and thus it is completely sufficient for all imaginable phenomenological applications of the pentagon functions. We implemented this approach in a public C++ library and a public Mathematica package, which we describe in sections 7.1 and 7.2.
Nevertheless, one could ask a question how to express weight three and four pentagon functions in terms of polylogarithmic functions. We found expressions for all 90|21 weight-3 pentagon functions (5.33) in terms of logarithms, dilogarithms and trilogarithms with arguments built from the letters of the pentagon alphabet. Using this weight-3 result it is straightforward to obtain and alternative integral representation for all 315|156 weight-4 pentagon functions (5.38). Indeed, we explicitly implement all inner three-fold iterated integrations and obtain g (4) where h (3) are weight-3 polylogarithmic functions, which involve logarithms, dilogarithms and trilogarithms. Thus, we have at hand two alternative ways to evaluate weight-3 and weight-4 pentagon functions -the one extensively described in the previous subsections and the one briefly outlined in this subsection. We stick to the first approach and use our private Mathematica implementation of the second approach as a highly-nontrivial test for the public library and the public package.

Master integrals in arbitrary channel
All previous considerations were restricted to the subset P 0 (4.15) of the physical region, the s 12 -scattering channel. We classified pentagon functions up to weight four, which are welldefined inside P 0 , and we provided -expansion of all UT master integrals in all orientations in terms of the pentagon functions. Thus, we are able to evaluate all master integrals within P 0 . In this section we demonstrate that having at hand results for all 5! orientations of the master integrals is equivalent to knowing them in any physical region. Let X be a kinematic point in an arbitrary scattering channel of the physical region. We can always find an elementσ in S 5 which maps the point X into a pointX from the s 12 -scattering channel,X =σX ∈ P 0 . (5.48) The previous equation implies that X belongs to the scattering channelσ 1σ2 →σ 3σ4σ5 whereσ := (σ) −1 and the following inequalities which specify the scattering channel hold sσ 1σ2 , sσ 3σ4 , sσ 3σ5 , sσ 4σ5 > 0 , X : sσ 1σ3 , sσ 1σ4 , sσ 1σ5 , sσ 2σ3 , sσ 2σ4 , sσ 2σ4 < 0 , Then we can use eqs. (4.3) and (4.4) to evaluate master integrals f τ,σ (in arbitrary orientation σ) at the point X from the arbitrary scattering channel as follows, Let us note that the UT master integrals by definition in eq. (4.6) are eigenvectors of the parity conjugation. Hence, if we evaluated the UT master integrals at a point X = (v 1 , . . . , v 5 ; 5 ), we can obtain their values at the parity-conjugated point X P = (v 1 , . . . , v 5 ; − 5 ) by inverting the signs of the parity-odd integrals, The Z 2 -grading discussed in sections 4.5 and 5 guaranties that also the pentagon functions and the transcendental constants are the eigenvectors of parity conjugation. Consequently, one can parity-conjugate each function and constant individually in the way that is compatible with eq. (5.51). Finally, we note that we could restrict our attention to even smaller portion of the physical phase space than P 0 . Indeed, the s 12 -channel is invariant under the S 2 × S 3permutations, which preserve signs of the Mandelstam invariants in (4.15). Then eq. (5.50) atσ ∈ S 2 × S 3 and sgn(σ) = +1 relates all UT master integrals evaluated at a pair of points of P 0 . If sgn(σ) = −1, then eq. (5.50) should be supplemented with eq. (5.51). Thus, knowing the values of the master integrals in the region P 0 /S 2 ×S 3 , which is six times smaller than P 0 , we can reconstruct values of the master in the whole P 0 , and consequently, in the whole physical phase space.
In conclusion, we reduced the problem of evaluating the master integrals in arbitrary physical channel to evaluating their permutations in the s 12 -channel. Our classification of the pentagon functions in the s 12 -channel is thus sufficient to evaluate the master integrals in an arbitrary physical channel.

Behavior near the boundary ∆ = 0
The obtained analytic expressions for the pentagon functions, enable us to evaluate master integrals in the physical channels and also to study their asymptotic regimes. In [59], the multi-Regge asymptotics in the non-planar sector of five-particle massless amplitudes has been studied. One could also study soft or collinear asymptotics by approaching s ij = 0 boundaries of the physical scattering channels. We are not going to plunge here into the detailed study of all possible singular regimes of the master integrals. Instead, we consider asymptotic behaviour of the master integrals (pentagon functions) when approaching ∆ = 0 boundary of P 0 . The surface ∆ = 0 separates any physical channel into to halves: δ > 0 and δ < 0. It was demonstrated in [59] that the five-particle non-planar Feynman integrals have discontinuities on subvarieties of ∆ = 0 and can even be divergent there. This is a peculiar feature of the non-planar five-particle scattering, which does not manifest itself in simpler planar master integrals studied in the past. To gain more experience with the non-planar master integrals we consider ∆ → 0 asymptotics of the pentagon functions.
We should stress that discontinuities and divergences at ∆ = 0 appear in the Feynman integrals, but the scattering amplitudes are expected to be free of these singularities in the physical region. In other words, only certain combinations of the Feynman integrals are allowed to contribute to the physical amplitude. The superamplitudes presented in [59] have been tested to satisfy this property. One could try to reverse the argumentation and apply the bootstrap approach to amplitudes [65,[72][73][74]. In the spirit of the Steinmann relations [70,71] for the planar hexagon scattering in N = 4 super-Yang-Mills theory, one could exploit the absence of discontinuities at ∆ = 0 as a nontrivial dynamical input on an amplitude ansatz consisting of the pentagon functions. It would be interesting to see how strong is this restriction, and to which extent it fixes nonplanar five-point two-loop amplitudes, in particular the QCD helicity amplitudes.
We choose a generic kinematic point X b on the boundary ∆ = 0 of P 0 (4.15). It describes a configuration of momenta {p µ i } 5 i=1 lying in a 3-dimensional hyperplane, but none of the momenta are soft or collinear, i.e. none of the Mandelstam invariants s ij vanish, One can easily check that all parity-odd letters are equal to 1 at X = X b , a k (X b ) = 0 , W 25+k (X b ) = 1 for k = 1, . . . , 5 , (6.2) and the statement does not depend on the path 5 we choose to approach X b . Let us inspect the pentagon functions in the asymptotic regime X → X b . The results of this subsection are implemented in the Mathematica package (see section 7.1).

Weights one and two
We start with weights one and two. In view of eq. (6.2), the parity-odd functions in eqs. (5.21), (5.22) and (5.25) take the following form for k = 4, 5 , (6.3c) while among the parity-even pentagon functions in eqs. (5.12), (5.14) and (5.25) only the form of g (1) 2,10 (X) = log W 31 (X) changes at X → X b : it is divergent on ∆ = 0. Let us note that g (1) 2,10 is absent in the weight-1 (5.13) and weight-2 (5.32a), (5.32b) solutions of the DE. Thus, they are finite at ∆ = 0. Approaching the surface ∆ = 0 from the opposite sides -δ > 0 and δ < 0 -inside any physical channel, we find a discontinuity in the parity-odd UT master integrals if they do not vanish at ∆ = 0. Inspecting the parityodd weight-2 solutions given in eq. (5.32b) at ∆ = 0, which involve only weight-1 pentagon functions g (1) 3,k (X b ), we find that they are not identically zero [59]. More precisely, they vanish on some parts of ∆ = 0 carved out by {a k (X) = 0} k=1,3,4,5 , while are constant and nonzero on the remaining parts. This is illustrated in fig. 3.

Weight three
Let us find asymptotics of the weight-3 pentagon functions at ∆ → 0 using their integral representation in eq. (5.33). The integration path γ is a line segment (4.17) with the end point X(t = 1) = X b . The functions h (2) (X(t)) from eq. (5.33) are finite at t = 1. On the line segment γ, we find (see eqs. (6.2) and (A.2)) that Then the integration kernels in eq. (5.33) involve the following singularities at t → 1: The latter are integrable singularities, and only d log W 31 (t) introduces a divergence. We regularize it with ε → 0 as follows, where f (t) is a regular function, and the last term contains logarithmic divergence if f (1) = 1, Inspecting all weight-3 pentagon functions {g (3) i } 111 i=1 (5.33), we find that only parity-odd g (3) 98 and g (3) 103 are singular at ∆ = 0. Parity-odd weight-3 solutions (5.37b) of the DEs involve logarithmically divergent at ∆ → 0 weight-1 and weight-3 pentagon functions g (1) 2,10 , g 98 , g 103 , and indeed we observe divergences of the UT master integrals at ∆ = 0.

Weight four
In a similar spirit, we study asymptotics of the weight-4 pentagon functions given in eq. (5.44). Only terms with k = 31 and/or j = 31 produce singularities at X = X b . Thus, we need to regularize three types of terms: where f (u) is a regular function. We have already regularized the first term (6.8a) in (6.6). Taking into account eq. (5.43), we regularize the second term, eq. (6.8b), as follows: Let us note that log ( 5 (u)) is an integrable singularity at u → 1. Finally, we regularize the third term (6.8c) as follows: The integrations in eq. (6.10) are convergent, and singularities are revealed in the form of divergent log ( 5 (1 − ε)) as ε → 0. The divergent terms in eqs. (6.8a) to (6.8c) are present only in the parity-odd pentagon functions of weight four. Thus, only they contain logarithmic divergences log 5 and log 2 5 in the limit ∆ → 0. These divergences of the pentagon functions are also inherited by the parity-odd weight-4 solutions of the DEs given by eq. (5.46b).

Numerical Evaluation
In section 5, we constructed a complete set of pentagon functions required to analytically represent (up to weight 4) all UT master integrals of the topologies shown in fig. 1. In this section, we describe our implementation of numerical evaluation of the pentagon functions and the UT master integrals.
All results of sections 5.8 and 6 are implemented in a Mathematica package, discussed in section 7.1. The package is provided mainly for the purpose of demonstration, and it is not intended for the use cases where high throughput and/or numerical robustness is required. For the later, we provide a C++ library, which we present in section 7.2. The library is optimized for performance, and its numerical efficiency makes it well-suited for evaluation of phase-space integrals with the Monte-Carlo method -the key ingredient for obtaining theoretical predictions for any observable cross section.

Mathematica package PentagonMI
The Mathematica package PentagonMI implements numerical evaluation of all UT master integrals f τ,σ , defined in eqs. (4.5), (4.7a) and (4.7b), at any point of the physical phase space. The master integrals are expressed in the basis of the pentagon functions, constructed in section 5. The package consists of three main components. The first component is data files containing the definitions of the objects employed in this paper. The data files are in Mathematica format. However, they are made to be self-consistent and understandable as plain text, such that they can be used outside of this package. 6 The second component uses P pentagon PB planar pentagon-box HB hexagon-box DP double pentagon  To install the package, follow the instructions in the "Installation" section of the README.md file in the root directory of the distribution.
Let us first describe the data files in the directory datafiles/ provided with the package. Our choice of UT master integrals, i.e. eq. (4.5), is specified for each Feynman integral topology by a corresponding file in the directory datafiles/UT-MI/ (see abbreviations of topologies τ in table 4). To reduce the size of files, we write all UT master integrals as

Q-linear combinations of a smaller subset of 1917 UT integrals
which we found by identifying Feynman integrals among different topologies and permutations [30]. Each UT master integral is rewritten as a linear combination of UT integrals G in the file datafiles/MI_in_G.m. UT integrals G expressed in terms of pentagon functions, as given by eqs. The parity grading of the alphabet, master integrals, pentagon functions, and transcendental constants plays an important role in our classification. We list all parity-odd objects in the file datafiles/parity-grading.m. Further details can be found in the datafiles/README.md file supplied along with the distribution. The main interface of PentagonMI is given by the function EvaluateMI, which accepts a list of master integrals to be evaluated, and a kinematical point in the physical region given EvaluateMI uses eq. (5.50) to find a permutation that maps it to P 0 . By default, δ > 0 is assumed, and evaluation of master integrals at a parity-conjugated point (δ < 0) can be requested with the option "ParityConjugation" -> True. Several other options that can be used to modify certain aspects of evaluation are available; we refer to the documentation provided in the file PentagonMI.m. For an example, see the program test/all_master_integrals.m, which evaluates all master integrals in all permutations at a single phase-space point.
EvaluateMI is only responsible for constructing a representation of each UT master integral in terms of pentagon functions. In order to obtain numerical values of the master integrals, numerical values of the pentagon functions at the given kinematical point are required. Numerical evaluation of the pentagon functions is carried out either with the (sub-)package PentagonFunctionsM, described in the next section, or through a Mathematica interface of the C++ library (see section 7.2). By default, the latter is chosen if available, and the option "UseCppLib" -> False can be used to choose the Mathematica implementation instead.

Numerical evaluation of pentagon functions
We implemented numerical evaluation of pentagon functions in a Mathematica package PentagonFunctionsM. The package can be used independently or as a part of PentagonMI.
All functions are evaluated in the analyticity region P 0 (4.15). We evaluate weight-1 and weight-2 functions, explicitly given by eqs.
17 , g (4) 113 , g (4) 470 . The requested numerical integration error of weight-3 and weight-4 functions can be set with the option "IntegrationPrecisionGoal". Its default value is 10, which means that the integration is terminated when the (negative) log 10 of an estimate of either relative or absolute error reaches 10. The requested functions are evaluated in parallel by default, using all available CPUs. Parallelization can be disabled by setting the option "Parallel" -> False.
Kinematical points are allowed to lie on the surfaces of spurious singularities exactly. In this case, the corresponding d log-kernels are set to zero, see discussion around eqs. (5.8), (5.35b) and (5.42).
As a special case, EvaluatePentagonFunctions also evaluates asymptotics of the pentagon functions at ∆ → 0 as discussed in section 6. The special case is activated automatically whenever evaluation at a kinematical point sitting on the boundary ∆ = 0 is requested. Concretely, for weight-3 and weight-4 pentagon functions, we implemented the regularization of the divergent one-fold integrals as introduced in eqs. (6.6), (6.9) and (6.10). The asymptotics, which is divergent in the limit ∆ → 0, is (at most quadratic) polynomial in log( 5 ) with numerical coefficients resulting from the regularized one-fold integrations. An example can be found in the test program test/functions_delta_singular.m. In this program all pentagon functions are evaluated at a kinematical point X b with ∆(X b ) = 0. Also, as a consistency check, asymptotics of the divergent at the point X b weight-4 functions are compared to their values at a point that is slightly deformed away from ∆ = 0.

C++ library PentagonFunctions++
One of the main goals of this paper is to take advantage of analytic understanding of the five particles massless scattering to derive a representation of the corresponding two-loop master integrals that is suitable for phenomenological applications. In particular, this representation should lend itself to a numerically efficient and stable implementation. We believe that the classification of pentagon functions, which we carried out in section 5, indeed provides such a representation. Nonetheless, we find that the Mathematica implementation described in the previous section does not realize the full potential of our method. To this end, we implement numerical evaluation of pentagon functions in a C++ library PentagonFunctions++, which we present in this section.

Features
PentagonFunctions++ is a C++14 library, which implements numerical evaluation of the pentagon functions, classified in section 5, in their analyticity region P 0 (4.15).
For numerical evaluation of weight-1 and weight-2 functions, we use their explicit representation in eqs. (5.12), (5.14), (5.15), (5.25) and (5.28) in terms of the log, Li 2 , and Cl 2 functions. The latter are evaluated numerically with a custom C++ implementation [76] based on the algorithms of [77]. For weight-3 and weight-4 functions, we use the one-fold integral representations in eqs. (5.35a), (5.35b) and (5.44) and evaluate the integrals numerically. The integrands of certain weight-3 and weight-4 functions are somewhat lengthy. Thus, to speed up their numerical evaluation, we optimize the integrands with respect to the number of floating-point operations with the code-optimization facilities [78] of the computer algebra system FORM [79]. 7 The choice of a numerical integration algorithm for the evaluation of the one-fold integrals (quadrature) can significantly impact evaluation times. Thus, it is essential to choose an algorithm that is suitable for the problem at hand. We employ the double exponential tanh-sinh quadrature [80]. The quadrature exploits a change of an integration variable t ∈ (0, 1), which maps the endpoints of the integration region to infinities, x ∈ (−∞, +∞), and the transformed integrand decays double exponentially, i.e. as exp − π 2 exp(|x|) with x → ±∞. It can then be shown [80] that the integral can be approximated remarkably well by a simple trapezoidal rule. In fact, it was proven in [80,81] that the tanh-sinh quadrature is the optimal choice for integrands that are analytic inside the integration domain (excluding, perhaps, the endpoints) in a sense that it requires the least number of evaluations of the integrand to reach a given integration error. For this class of integrands, the tanhsinh quadrature converges exponentially, i.e. the number of correct digits in the numerical approximation is proportional to the number of evaluations of the integrand. Integrable singularities at the endpoints of the integration domain, such as the logarithmic singularity in the integrands of the weight-4 functions (see section 5.6.2), do not introduce any complications for this quadrature, hence no special handling is required. PentagonFunctions++ uses an adapted implementation of the tanh-sinh quadrature from Boost C++ [82].
We pointed out in eqs. (5.7) and (5.8) that several letters W k∈α of the pentagon alphabet do not have a definite sign inside the analyticity region P 0 and vanish at the endpoint t = 0 of the integration interval. Thus, their d log-forms have a simple pole at t = 0. As we discussed around eqs. (5.36) and (5.45), these poles are compensated by vanishing combinations of weight-2 functions. The quadrature algorithm discussed above might require evaluation of the integrands very close to the endpoints. It is thus important to ensure that the cancellation of the poles is numerically stable. To this end, in the neighborhood of t = 0, t <t we evaluate the kernels d log (W k (t)) together with their coefficients h k (t) through their generalized series expansion around t = 0 as such that no numerical cancellation of the pole has to occur. The thresholdt is chosen such thatt 2 ε T , where ε T is the roundoff error (or machine epsilon) of the floating-point number type T (e.g. ε double 10 −16 ).
On certain subvarieties of the physical phase space (spurious singularities) the letters W k∈α might be identically zero. As we mentioned around eq. (5.8), on these subvarieties the corresponding integration kernels also vanish. However, in small neighborhoods of these subvarieties the letters W k∈α (γ(t)) = b k (X) t almost vanish along the whole line segment γ = [X 0 ; X], i.e. 0 < |b k (X)| 1. Then the contribution from h k d log W k to the integral is rendered small by potentially large cancellations in h k (t). In principle, this situation can be avoided by using an appropriate representation of h k (t) and/or path deformations. However, we find that it is sufficient to simply set the integration kernels d log (W k ) to zero exactly whenever |b k (X)| is below a certain threshold. We note that the pentagon functions are analytic on the surfaces of spurious singularities, and only the functions that vanish on a particular spurious-singularity surface can be significantly impacted by this procedure. The threshold can thus be adjusted in such a way that only insignificant neighborhoods of the spurious-singularity surfaces are potentially affected. We demonstrate this a posteriori in section 7.2.3. We leave a more refined analysis of the spurious singularities for future study.
PentagonFunctions++ is able to perform all evaluations in three fixed-precision floatingpoint types: double, quadruple and octuple precision, which respectively represent significands of approximately 16, 32, and 64 decimal digits. We use a C++ implementation of quadruple and octuple numerical types from the qd library [83]. Numerical evaluation in multiple fixed-precision types is indispensable for understanding numerical stability of the implementation, as well as for adaptively balancing precision against performance.

Usage
The library can be obtained from the git repository [84] with git clone https://gitlab.com/pentagon-functions/PentagonFunctions-cpp.git To install the library, follow the instructions in the "Installation" section of the README.md file in the root directory of the distribution [84].
The intended way to use PentagonFunctions++ is to write a C++ program, which links to the provided static or shared library. Further details can be found in the "Usage" section of the README.md file found in the root directory of the distribution [84].
The main interface of the library is provided by a struct FunctionID, which is declared in the header file src/FunctionID.h. An instance of FunctionID, constructed with integer arguments (w,i,j) or (w,i), represents the pentagon function of weight w and indices i,j, according to the definitions in eqs. (5.12), (5.14), (5.15), (5.25), (5.28), (5.35a), (5.35b) and (5.44). The instances of FunctionID can be used to obtain callable function objects of numerical type T (e.g. double, dd_real, qd_real) with the method get_evaluator<T>(). For example, with FunctionID fid{4,471}; auto f = fid.get_evaluator<double>(); one creates a function object f, which can be used to numerically evaluate the pentagon function g (4) 471 at any number of kinematical points in double precision. The method get_evaluator<T> performs an initialization stage of the integration framework that need not be repeated for each subsequent evaluation. Several example programs can be found in examples/ directory of the distribution.
The termination condition of the numerical integration is controlled by the global variable template <typename T> extern T IntegrationTolerance; which specifies the tolerance for each numerical type T independently. It is declared in the header file src/Constants.h. The numerical integration is terminated when the difference of two subsequent estimates of the integral have absolute value less than the tolerance multiplied by an estimate of the L 1 norm of the integral. The default value is chosen such that the integration error is close to the rounding error of the numerical type T. Finer control over tolerance might be exploited for improving either integration speed or precision of the results.
For convenience, we also provide a Mathematica interface. It is realized as a Mathematica package PentagonFunctions, which interacts with the program mathematica_interface/evaluator_math.cpp. The interface is similar to the one of the package PentagonFunctionsM, described in the previous section. An example of the interface usage can be found in examples/math_interface.m.

Performance
In this section, we demonstrate performance of our implementation with respect to evaluation speed and numerical stability, which are the most important properties of a numerical algorithm.
To characterize evaluation speed of our implementation, we evaluate all pentagon functions with PentagonFunctions++ (with the standard settings) at a random generic point from the physical phase space. We perform the evaluation on a single core of Intel(R) Core(TM) i7-7700 CPU @ 3.60GHz. We show the evaluation times as well as the (minimal) number of correct digits 8 for the three supported floating-point types in table 5. Evaluation of the planar subset of the pentagon functions 9 takes approximately 40% of the total evaluation time. Comparing to the evaluation times of the planar pentagon functions of [53] reported in [35], we observe that the planar subset of our implementation evaluates approximately 100 times faster.
Further, we demonstrate the numerical stability of our implementation by evaluating all pentagon functions on a sample of 90000 phase-space points, drawn from a typical distribution employed in computations of differential cross sections for processes with five Precision Correct digits Timing (s)   double  13  2.5  quadruple  29  180  octuple 60 3900 Table 5: Evaluation times of all pentagon functions on a typical phase-space point. Evaluation is performed in a single thread. massless particles 10 . We evaluate all pentagon functions in double and quadruple precision at each phase space point, and we use the latter to compute the accuracy of the former. We characterize the accuracy of the evaluationĝ i (X) of the i-th pentagon function on a kinematical point X by the logarithmic relative error ("correct digits") r i (X) which we define as i (X) is the numerical evaluation of the same function in quadruple precision. We define the minimal logarithmic relative error among all pentagon functions at the kinematical point X as We display the distribution of R(X) over the phase space in fig. 4. We observe very good numerical stability in the bulk of the phase space: only 0.1% of the phase-space points evaluate with less than 8 correct digits. All 12 kinematical points X (R<6) with R < 6 are from the region of the phase space with 0 < δ 1 . As we discussed in section 6, some pentagon functions diverge in the limit δ → 0. But the divergence is only logarithmic. So, with min X (R<6) [δ] 10 −7 , the absolute values of the divergent pentagon functions still remain relatively small. However, the condition number κ of e.g. the function g (1) 2,10 = log( i δ) diverges much faster, In other words, numerical evaluation of the function g 2,10 in the regime δ 1 becomes dominated by the rounding error of the input data (Mandelstam invariants) much earlier than the function itself becomes large. Let us mention that it should be possible to circumvent this issue in applications to numerical evaluation of complete two-loop amplitudes expressed in terms of the pentagon functions. We leave these considerations for future studies. We would like to note that the desired precision of pentagon function evaluation heavily depends on the intended application. For this reason, we did not attempt to reach any given accuracy threshold. Instead, we have studied the numerical stability of doubleprecision evaluation in the default running mode of PentagonFunctions++. If a certain accuracy threshold is to be guaranteed, the following strategy can be attempted. Only the kinematical points close to either spurious singularities or δ = 0 are expected to be problematic, which can be detected before any integrations take place. Then the potentially problematic functions can be evaluated in higher precision with the help of multi-precision facilities of PentagonFunctions++.
We conclude that the performance of our numerical implementation of the pentagon functions will almost certainly be sufficient for the main anticipated application in phenomenological studies.

Validation
The presented analytic expressions for the master integrals as well as their numerical implementations have passed a number of cross checks which we report in this section.
The classification of all pentagon functions (planar and nonplanar) elaborated in this paper does not directly rely on the planar classification from [53], which brings forward the cyclic symmetry. In other words, the analytic expressions for the master integrals of the planar pentagon-box topology provided here are not literally the same as in [53], but of course both should be in agreement. Indeed, using our code, we reproduce numerical values of the planar integrals calculated by the code from [53]. Furthermore, we employed an implementation of numerical unitarity framework [16,17,22] in Caravel [87] to successfully reproduce the numerical evaluations of the leading-color two-loop five-gluon helicity amplitudes in the physical region presented in [20].
One of the crucial ingredients in the analytic solution of the DEs for the master integrals are initial values, see section 4.4. For the double-pentagon family we use the initial values from [56], which have been already validated independently with pySecDec [88,89] for several permutations of this topology 11 . Multi-digit numerical values of the full set of initial values presented in this paper have been already employed in [26,59], where expected physical properties have been observed. In particular, the weight-3 and weight-4 initial values have been probed in [26], where the weight-drop of the five-gluon two-loop hard function has been explicitly verified. This provides an indirect consistency check.
We check our analytic solutions of the DEs comparing their numerical evaluations with the numerical integration of the DEs via generalized series expansions [90], implemented in DiffExp [91]. We employed DiffExp to transport the initial values at the reference point eq. (4.16) to an arbitrary point in the physical region that can be reached without leaving the analyticity region and compared the result with our evaluations. We found perfect agreement within numerical precision. We should mention that we integrate each permutation of the DE separately and never cross the physical region boundaries. It would be interesting to analytically continue nonplanar topologies from one scattering channel to another with DiffExp.
Also we numerically evaluated a number of the nonplanar master integrals with pySecDec [88,89] at arbitrary points of the physical phase space and compared the result with our evaluations of the pentagon functions finding a perfect agreement. The latter check is sensitive to the initial values and the DE integration.
In order to verify the one-fold representations of the pentagon functions, see (5.33) and (5.44), and their numerical implementation, we constructed alternative representations of the weight-three and four pentagon functions. The weight-three functions are written explicitly as polynomials in the polylogarithmic functions, and weight-four functions are one-fold d log integrals with the weight-three integrands (see section 5.7). We find that both implementations of the pentagon functions agree numerically. Obviously, the analytic expressions for the pentagon functions -the one implemented in the public library and the alternative one -are completely different. They correspond to different rewritings of the iterated integral form of the pentagon functions into a more tangible form. The numerical agreement of the two representations is a strong check for both of them as well as for their numerical implementations.
Last but not least, the two public numerical implementations of the pentagon functions -the Mathematica package and the C++ library -enables us to control the quadrature accuracy.

Conclusions
In this paper, we constructed a complete set of transcendental functions, pentagon functions, which are sufficient to represent all planar and non-planar master integrals for massless five-point two-loop and one-loop scattering amplitudes up to transcendental weight four. Our analytic results are valid over the whole physical phase space. This was achieved by expressing all external momenta permutations of master integrals in the same set of pentagon functions. We expressed the pentagon functions of weights one and two in terms of (di-)logarithms, and we constructed one-fold integral representations for the pentagon functions of weights three and four. We presented an implementation of numerical evaluation of pentagon functions in a public C++ library. The latter was designed to satisfy demands of phenomenological applications. Thus, for the first time, all massless five-point two-loop Feynman integrals are available for immediate application in computations of fully differential cross sections. Together with the ongoing advances in reduction of five-point two-loop amplitudes, our results pave the way for computation of NNLO predictions for a number of key scattering processes at hadron colliders. The latter include production of three hard jets, two-photon and three-photon production in association with jets.
The pentagon functions are crucial for finding compact analytic representations of scattering amplitudes as well as for studying their asymptotic limits. In fact, the analytic form of amplitudes can be directly reconstructed from their numerical evaluations over finite fields [8]. This approach has lead to remarkable progress in calculation of planar five-point amplitudes [17,[20][21][22][23][24][25][26][27][28][29]35]. We expect that our results will open a possibility of extending these methods to non-planar sector.
Our strategy of constructing bases of transcendental functions can be generalized to scattering processes with even larger number of scales, such as five-point processes with massive particles. It might be also possible to apply our approach for the cases where Feynman integrals evaluate to iterated integrals over more complicated differential forms, such as modular forms (see e.g. [92][93][94]). It would be interesting to explore this in the future. Finally, an interesting question is to compare efficiency of the numerical evaluation of our analytic results to the approach of solving DEs numerically via generalized series expansions [90,91].

A Pentagon alphabet
We recall the definition of the nonplanar pentagon alphabet [51,65]. We use shorthand notations where we assume cyclicity of v-variable labels, see eq. (2.1), and for the sake of presentation we split 31 letters {W i } 31 i=1 of the alphabet in the orbits of the cyclic group Z 5 , Thus, {W i } 25 i=1 and W 31 are parity-even and {W i } 30 i=26 are parity-odd. Note that, since 5 is imaginary in the physical region, the parity-conjugation is equivalent to the complex conjugation in this region.
Working in the non-planar sector we have to consider all S 5 permutations of the external momenta p µ 1 , . . . , p µ 5 . The alphabet is closed under this action, which induces representation of S 5 in the space of the letters. The decomposition in the irreducible representations of S 5 can be found in [65]. Here we prefer to work with reducible representations, and we summarize the action of S 5 on the letters: where all indices run cyclically over 26, . . . , 30, so only the five parity-odd letters appear in the transformation rules.
• W 31 is mapped to itself up to sign, i.e. W 31 → W 31 or W 31 → −W 31 depending on the signature of S 5 permutation.
• Considering d log-forms of the letters, the previous transformations simplify. For the parity-even letters, the action of i=21 and {d log W 31 } is by permutations, and for the parity-odd letters, {d log W i } 30 i=26 transform linearly.

B Physical region geometry
The analyticity region P 0 (4.15) is not convex, but any line segment from X 0 (4.16) to any X ∈ P 0 lies inside P 0 . Let us prove this statement. We are going to show that any ray from X 0 crosses the boundary of P 0 only once, i.e. it cannot enter back inside P 0 . We need to consider only ∆ = 0 boundary of P 0 . For other boundaries of P 0 , i.e. s ij = 0 with i, j = 1, 2 or i, j = 3, 4, 5 or i = 1, 2 and j = 3, 4, 5, the statement is obvious. Let us choose an arbitrary X 1 on the boundary of P 0 , ∆(X 1 ) = 0 , s 12 , s 34 , s 35 , s 45 | X 1 ≥ 0 , s 13 , s 14 , s 15 , s 23 , s 24 , s 25 | X 1 ≤ 0 , (B.1) and connect it with X 0 by the line segment γ = [X 0 ; X 1 ] parametrized by 0 ≤ t ≤ 1. We need to prove that ∆(t) ≡ ∆(γ(t)) < 0 at 1 − ε < t < 1 for small positive ε. Since ∆ is a degree-4 polynomial in the Mandelstam variables, it is enough to show that none of the following four inequalities is compatible with (B. This can be verified using a computer algebra system. Extra simplification of the inequalities is achieved by fixing value of one of s ij (since all expressions are homogeneous polynomials in the Mandelstam invariants) and by using ∆ (k−1) (1) = 0 to lower the degree of ∆ (k) (1).