Pentagon Functions for One-Mass Planar Scattering Amplitudes

We present analytic results for all planar two-loop Feynman integrals contributing to five-particle scattering amplitudes with one external massive leg. We express the integrals in terms of a basis of algebraically-independent transcendental functions, which we call one-mass pentagon functions. We construct them by using the properties of iterated integrals with logarithmic kernels. The pentagon functions are manifestly free of unphysical branch cuts, do not require analytic continuation, and can be readily evaluated over the whole physical phase space of the massive particle production channel. We develop an efficient algorithm for their numerical evaluation and present a public implementation suitable for direct phenomenological applications.


Introduction
Electroweak boson production in association with photons and jets plays a key role in the physics program of hadron colliders such as the Large Hadron Collider (LHC). It provides a window into understanding the electroweak symmetry breaking mechanism and constitutes a prominent irreducible background in new-physics searches. Our ability to correctly interpret the experimental data relies to a large extent on the availability of precise theoretical predictions for this class of processes. Generally the knowledge of the next-to-next-to-leading order (NNLO) in the perturbative expansion in the strong coupling constant is highly desirable in order to harness the wealth of data to be collected. To access the important W/Zbb and W/Zγγ production processes, as well as W/Z + multijet associated production, the calculation of 2 → 3 NNLO corrections is required. The latter represents a major challenge and remains on the forefront of modern computational methods' capabilities. The remarkable progress in the computation of massless two-loop five-particle scattering amplitudes in the last few years [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20] has recently culminated in the first NNLO theoretical predictions for 2 → 3 processes [21][22][23][24][25]. This success was enabled by the availability of analytic results for the required two-loop Feynman integrals [1,3,10,[26][27][28][29][30][31]. In particular, it was crucial to express the latter in terms of judicious bases of transcendental functions [27,31]. On the one hand, these bases make the analytic structure of scattering amplitudes manifest, enable the analytic cancellation of ultraviolet and infrared divergences, and facilitate the derivation of compact representations. On the other hand, they are amenable for efficient and stable numerical evaluation suitable for phenomenological applications [31]. In this paper, we construct a similar basis of transcendental functions for the planar two-loop five-particle scattering amplitudes with a single external massive particle.
The planar one-mass two-loop five-point Feynman integrals are currently the subject of extensive study. The four-point integrals with two external masses were computed in refs. [32][33][34]. As for the genuine five-point Feynman integrals, bases of independent integrals -alias master integrals -satisfying systems of canonical differential equations (DEs) [35] have been constructed in ref. [36]. Their numerical evaluation was achieved by solving the DEs in terms of generalized power series expansions [36][37][38]. Expressions for the master integrals in terms of multiple polylogarithms (MPLs) [39][40][41] have been derived in refs. [26,42,43]. In ref. [44], a basis for the subset of transcendental functions contributing to completely color ordered amplitudes has been constructed. Similarly to ref. [36], their numerical evaluation was performed through series expansion of their DEs. These developments have led to the first analytic calculations of two-loop five-particle amplitudes with an external massive leg, for the production of a W and a Higgs bosons associated with a pair of bottom quarks [44,45], as well as of the two-loop four-point form factor of a length-3 half-BPS operator in planar N = 4 sYM [46]. First results for the non-planar integrals have also been reported recently [47][48][49].
In this paper we construct a basis of transcendental functions that is sufficient to express any planar two-loop corrections for scattering cross sections involving four massless and one massive particle. We call the functions in this basis one-mass pentagon functions. We follow the approach of ref. [31]. We start from the canonical DEs of ref. [36] and generalize them to all permutations of the external massless legs. We compute the initial values for the DEs using the expressions of the master integrals in terms of MPLs from refs. [42,43]. We write the solutions order by order in the dimensional regularization parameter in terms of iterated integrals [50] with logarithmic kernels. The solutions are pure functions of uniform transcendental weight [35,51] and satisfy a shuffle algebra, which we use to linearize and eliminate all algebraic relations in the space of functions generated by the DEs. The algebraic independence and closure over momenta permutations of the one-mass pentagon functions is very advantageous for the modern workflow to compute scattering amplitudes analytically, which is based on the application of finite field arithmetic and functional reconstruction [52,53]. Indeed, writing all permutations of the Feynman integrals in terms of the same basis of transcendental functions gives access to simplifications of the rational coefficients that would otherwise be missed. We find explicit representations of the basis functions up to weight two in terms of logarithms and dilogarithms, and up to weight four in terms of one-dimensional integrals. All the expressions are branch-cut free within the physical scattering region of the massive-particle production channels. We thus demonstrate that the constructive algorithm of ref. [31] can be straightforwardly generalized beyond the purely massless kinematics. The numerical evaluation of the weight three and four pentagon functions on the other hand is significantly more challenging compared to their massless counterpart. This is due to the non-trivial geometry of the phase space, which we study in detail, and to the presence of non-linear spurious singularities in the one-dimensional integral representations. Nevertheless, we design an algorithm for the numerical evaluation that is fast and stable across the whole phase space of the 2 → 3 massive-particle production channels. We achieve performance comparable to that of the massless pentagon functions from ref. [31], and provide a public implementation which meets the demanding requirements of phenomenological applications.
Compared to a more traditional approach of attempting to express the master integrals in terms of MPLs, our method is rather insensitive to the presence of square roots in the symbol alphabet. Even in cases where the canonical " d log"-form of DEs can be obtained, it can be difficult, if not impossible, to find MPL solutions if the symbol alphabet contains non-rationalizable square roots [42,[54][55][56][57]. Very importantly, even when these attempts did succeed, it was frequently observed that the numerical evaluation in the physical scattering regions is very challenging [26,27,42,58]. In addition, contrary to naïve expectations, the physical properties may even be more obscured due to the presence of spurious branch cuts, which are instead absent in the iterated integral representation. Our approach elegantly bypasses these issues, and therefore can be a useful alternative for cases when the integration of the DEs in terms of MPLs is impractical or impossible.
The rest of the paper is structured as follows. We begin in section 2 by presenting our notation, discussing the kinematics, and defining the 2 → 3 physical scattering regions. In section 3 we define the planar families of five-particle integrals with one external massive leg up to two loops, and discuss how to express the master integrals in terms of Chen's iterated integrals through the canonical DEs they satisfy. In section 4 we discuss the construction of our function basis. First, in section 4.1 we use the properties of Chen's iterated integrals to construct systematically a basis of algebraically independent functions to express all permutations of the master integrals up to the order in required for two-loop corrections to cross sections. Then, in section 4.2 we design explicit expressions for the basis functions which are well suited for their numerical evaluation, and show that they are well defined in the chosen physical scattering region. We discuss the implementation in a C++ library of routines for the numerical evaluation of our function basis in section 5. We draw our conclusions in section 6. We also include a number of appendices, where we discuss in detail the physical phase-space geometry and the positivity properties of the symbol alphabet.

Scattering kinematics
In this section we discuss the scattering kinematics and introduce certain quantities which will play an important role in the rest of the paper. We consider the scattering of five particles, four massless and one massive. We take their momenta p i to be outgoing, and the momentum conservation reads (2.1) Following the notation of refs. [36,48], we assume that momentum p 1 is massive, while the others are massless, where ε µνρσ is the fully antisymmetric Levi-Civita symbol. The square of tr 5 is algebraically related to the Mandelstam invariants in eq. (2.3) through the Gram determinant of four external momenta 1 as We emphasize that, despite the relation in eq. (2.5), the sign of Im tr 5 is required to fully specify a point in the scattering phase space. The sign of Im tr 5 is changed under parity conjugation or permutations of external momenta, whereas √ ∆ 5 stays unchanged. Therefore one should be careful to distinguish these objects.
The other three Gram determinants that give rise to square roots which are relevant for the analytic structure of the five-particle one-mass scattering amplitudes are where λ denotes the Källen function, 3 } is closed under the permutations of external massless momenta S 4 , namely the ∆ (i) 3 's transform into each other under the action of σ ∈ S 4 . In this paper we denote by the permutation σ ∈ S 4 such that σ(i) = σ i . The action of σ on the external massless momenta is given by 12) and the permutation σ of a generic function f of the external kinematics -e.g. a Feynman integral or a scattering amplitude -is given by From the definition (2.4) it follows that σ • tr 5 = sign(σ) tr 5 , (2.14) where sign(σ) is the signature of the permutation σ. The Feynman integrals are multi-valued functions of the kinematics with a very intricate analytic structure. It is therefore fundamental to specify the domain of the kinematic variables X. We consider the physical scattering regions associated with the production of one massive and two massless particles, where i, j, k, l take distinct values in {2, 3, 4, 5}. These regions are relevant for instance for the production of a massive vector boson in association with two jets at a hadron collider. Without loss of generality we assume that the incoming momenta are p 4 and p 5 , i.e. we focus on the s 45 channel. The s 45 channel is defined by requiring that the momenta correspond to physical configurations of the scattering process 45 → 123, i.e. they are real and are associated with real scattering angles and positive energies. These constraints translate into the following inequalities for the scalar products of the momenta, as well as the following Gram-determinant inequalities [59][60][61], where i, j, k, l take distinct values in {1, 2, 3, 4, 5}. We give a thorough discussion of the role of the Gram determinants in appendix A. Here we content ourselves with noting that the only non-vanishing Gram determinant involving only one momentum is G(p 1 ) = p 2 1 , which we assume to be positive, and that all the Gram determinants involving four momenta are proportional to ∆ 5 . The negativity of ∆ 5 can be understood as a consequence of the reality of the momenta, which implies that tr 5 is purely imaginary (see eq. (2.4)) and hence, through eq. (2.5), that ∆ 5 < 0. Putting together all the inequalities above and simplifying them (see appendix A) gives the following definition of the s 45 scattering channel 45 → 123, which we label by P 45 : 2 In other words, given a set of Mandelstam invariants (2.3) satisfying the inequalities (2.18), there is a corresponding configuration of particle momenta describing the scattering process 45 → 123. 3 Moreover, since tr 5 is purely imaginary for real momenta, the scattering region P 45 is separated into two halves, corresponding to Im tr 5 > 0 and Im tr 5 < 0, Each of the halves is a connected region in the momentum space [59]. The two halves are mapped onto each other by a parity transformation. Apart from the massive-particle production channels (2.15), other possible physical channels are the decay into four particles (1 → 2345) and three-particle production with a single massive particle in the initial state (1j → klm where j, k, l, m take distinct values in {2, 3, 4, 5}). Here we focus on the case of massive-particle production, which is the most important for the hadron-collider phenomenology.

Integral families and differential equations
The modern approach to computing scattering amplitudes analytically relies on the fact that they can be expressed in terms of linear combinations of scalar Feynman integrals. The latter are then grouped into families based on their propagator structure. The advantages are that within each family only a finite number of integrals are linearly independenttypically much fewer than the Feynman diagrams contributing to the amplitudes -and that it is possible to express the amplitudes in terms of them algorithmically. These linearly independent integrals, called master integrals in the literature, constitute the basis of the linear span of all the integrals in the family. Crucially, the master integrals satisfy a system of first-order differential equations (DEs) [35,[62][63][64][65], which constitutes a powerful tool to compute them analytically. In this section we discuss the planar Feynman integral families relevant for the scattering of one massive and four massless particles at one and two loops with internal massless propagators. We closely follow ref. [36] and lift their results to all S 4 permutations of the external massless particles. We begin by defining the integral families. Then we discuss the DEs satisfied by the master integral bases proposed in ref. [36], and show how they can be solved in terms of Chen's iterated integrals [50]. We finish the section by discussing how we computed the initial values necessary to fix uniquely the solution of the DEs. The resulting expression of the master integrals in terms of iterated integrals is then the starting point of the construction of the function basis, which we discuss in section 4.

Integral families
The Feynman integrals appearing in planar five-particle amplitudes with a single massive external leg up to two-loop order belong to the families shown in figure 1 and permutations thereof. We adopt the definitions of ref. [36]. There is a single one-loop pentagon family, which we label by 1L, and three two-loop pentabox families, which we label by mzz, zmz and zzz depending on the location of the massive leg, as shown in figure 1. We consider each family in the S 4 permutations of the external massless particles. The divergences in the integrals are regulated in dimensional regularization with D = 4 − 2 . The one-loop one-mass pentagon integrals are defined as while the integrals of the two-loop one-mass pentabox families are given by Here γ E is the Euler-Mascheroni constant, µ is the regularization scale, σ is a permutation label, and a list of integers a encodes the propagator exponents. In this paper we set µ = 1,  which is equivalent to considering dimensionless ratios of the Mandelstam invariants in eq. (2.3) to µ 2 . The definitions of the inverse propagators D τ,σ are gathered in table 1.
The last three propagators of the two-loop families are irreducible scalar products, i.e. auxiliary propagators introduced in order to express the numerators. Therefore we assume that {a j } 11 j=9 are non-positive integers. The two-loop scattering amplitudes also contain Feynman integrals which factorize into the product of two one-loop integrals of the pentagon family 1L. Some of these "one-loop squared" integrals are linearly independent from the two-loop pentaboxes. They constitute an additional family shown in figure 2, which we dub 1L 2 . One can express the integrals from this family in terms of transcendental functions by multiplying the expressions of the one-loop integrals which constitute them. Since we are interested in constructing a set of algebraically independent functions, we can omit the one-loop squared integrals from our considerations in section 4.
The linear span of all scalar integrals G τ,σ [ a] in a given family τ and orientation of the external legs σ, for any (allowed) set of indices a, has a finite-dimensional basis. Any scalar integral G τ,σ [ a] can thus be expressed as a Q ( , X)-linear combination of the master integrals, namely a linear combination with coefficients which are rational functions of the Mandelstam invariants and the dimensional regulator . This is typically achieved by solving linear systems of integration-by-parts identities (IBPs) [66][67][68]. The choice of the master integrals is to a large extent arbitrary and can be exploited to simplify selected Table 1: Inverse propagators of the one-loop pentagon (1L) and of the two-loop pentabox families (mzz, zmz and zzz) for all permutations σ = (σ 2 σ 3 σ 4 σ 5 ) of the external massless legs. The representative diagrams are shown in figure 1 for the standard ordering of the external legs, σ id = (2345). aspects of their computation. In this paper, we take advantage of the integrals which satisfy the mathematical property known as transcendental purity [35,51]. We denote by g τ,σ the vector of pure master integrals for the family τ in the permutation σ. The pure master integrals g τ,σ are expressed as linear combinations of scalar Feynman integrals, where the coefficients T τ,σ [ a] are rational functions of the X and . The pure master integral bases for the families in figure 1 have been constructed in ref. [36] for the orientation σ id = (2345). We obtain the other orientations by permutations according to eq. (2.13).
There are relations also among the integrals belonging to different families and in different permutations. We identified the relations in the Q-linear spans of the one-and two-loop integrals respectively by analyzing their Symanzik polynomials (see e.g. ref. [69] Table 2: Number of independent one-and two-loop master integrals in the unions of integral families {1L} and mzz, zmz, zzz, 1L 2 respectively. The second column corresponds to the integrals in the standard permutation σ id , and the third column corresponds to the complete permutation orbit. The last row shows the number of non-factorizable two-loop master integrals in each case.

Permutations and physical channels
The master integrals of the families shown in figs. 1 and 2 appear in the scattering amplitudes in various permutations of the external massless legs. Following refs. [10,31,44], we aim to express all permutations of the master integrals in terms of a common basis of special functions which are well-defined and can be evaluated numerically in a whole physical scattering region. Without loss of generality we chose the latter to be the s 45 channel. We show that considering all S 4 permutations of the master-integral families allows us to relate any massive-particle production channel (2.15) to 45 → 123. Let us consider a phase-space point X of the physical channel σ 4 σ 5 → 1σ 2 σ 3 for some permutation σ = (σ 2 σ 3 σ 4 σ 5 ) ∈ S 4 , In order to evaluate the master integrals g τ,σ of the family τ in the orientationσ at a point X ∈ P σ 4 σ 5 , we first map X to a phase-space point X ∈ P 45 using the inverse permutation, Then, we evaluate the same master integrals in the composed permutation σσ at the corresponding point X in the s 45 channel, as shown by the following chain of equalities: This is to be contrasted with an alternative approach often adopted in the literature, which consists in considering the master integrals for a fixed ordering of the external legs and analytically continuing them to the other scattering regions. The permutations of the master integrals in this case are performed only at the evaluation stage, This however means that the integrals have to be evaluated a number of times equal to the number of required permutations for each desired phase-space point. Moreover, all the analytic simplifications which stem from cancellations among integrals in different permutations of the external legs are missed. This not only makes the expressions for the amplitudes more obscure, but also prevents from subtracting the ultraviolet and infrared poles analytically, and introduces numerical cancellations which may affect the numerical stability of the evaluations.
If the master integrals are pure, their DEs take a particularly simple canonical form [35] where only logarithmic differential forms contribute. In this work we employ the canonical DEs for g τ,σ id that were constructed in ref. [36], where a (i) τ,σ are matrices of rational numbers, and W i are algebraic functions of the kinematic variables called letters of the alphabet A D 4 . Here the integrals are normalized such that their expansion around = 0 starts with 0 , so that their fourth order in is the highest required for the computation of two-loop scattering amplitudes up to order 0 . The alphabet A D 4 is closed under the cyclic permutations of the external massless legs, We are interested in a basis of transcendental functions sufficient to represent any planar two-loop five-point one-mass scattering amplitude through transcendental weight four. We thus consider the master integrals in all orientations and construct an S 4 -closure A S 4 of the cyclic alphabet A D 4 by choosing a basis in the linear span of the permutation orbits σ • d log W i , for any i ∈ A D 4 and σ ∈ S 4 . We then straightforwardly generalize the canonical DEs (3.8) to all permutations σ ∈ S 4 of the integral families, d g τ,σ (X) = dÃ τ,σ (X) g τ,σ (X) , The alphabet A S 4 is a subset of the non-planar alphabet recently established in ref. [48], and for convenience we adopt the notation of the latter. We summarize here the embedding of A D 4 and A S 4 into the non-planar alphabet of ref. [48]. We explicitly spell out the expressions of the letters involved in our work in appendix C and in the ancillary file alphabet.m [70]. The cyclic alphabet contains 58 letters, which can be separated into two sets, (3.14) The permutation closure A S 4 contains 156 letters, which again can be split into two sets,

Solution of the DEs in terms of iterated integrals
The canonical DEs eq. (3.8) can be solved straightforwardly order-by-order in theexpansion (3.9) in terms of Chen's iterated integrals [50] (see also [71]) as where g (w) τ,σ (X 0 ) are the values of the master integrals at an arbitrary initial point X 0 , which we discuss in section 3.5. The iterated integrals [W i 1 , . . . , W in ] over d log kernels are defined as [] X 0 (X) := 1, (3.19) where γ is a path in the space of the kinematic variables connecting the initial point γ(0) = X 0 and γ(1) = X. The number of iterations n corresponds to the transcendental weight. We consider the initial point X 0 as fixed and treat the iterated integrals as functions of the endpoint X. 5 The iterated integrals form a shuffle algebra with the product where a = {a 1 , . . . , a m } and similarly for b and c, and a ¡ b denotes all possible ways to shuffle a and b. The shuffle product allows us to linearize the relations between products of transcendental functions of different weights, and therefore to systematically identify all functional relations through linear algebra. This equips us with a powerful tool in the construction of a function basis which we discuss in section 4.1.
The components g τ,σ of the -expansion of the pure master integrals are graded by the transcendental weight w. This grading is manifest in eq. (3.18) and is compatible with the shuffle product in eq. (3.20). In other words, g (w) τ,σ and the iterated integrals can be assigned a Z ≥0 charge. Further gradings can be associated with changing signs of the square roots of the problem: √ ∆ 5 (or equivalently tr 5 ) and ∆ (i) 3 with i = 1, 2, 3. The pure master integrals which are normalized by one of these square roots are odd with respect to that square root, while the others are even. We can thus assign a Z 2 = {+, −} = {even, odd} charge to each pure master integral for each of the four square roots. One can easily verify that the d log forms in eq. (3.18) can also be assigned the same charges (see appendix C). Thus, it follows from eqs. (3.18) and (3.20) that the shuffle algebra of the iterated integrals with logarithmic integral kernels drawn from the alphabet A S 4 is graded by Z ≥0 × (Z 2 ) 4 . This grading constitutes a useful organizing principle in the construction of the function basis discussed in section 4.1.
We stress that the choice of signs of the square roots is arbitrary and does not change within P + 45 . In this paper we choose The results for the opposite sign of each of the square roots can be obtained by flipping the sign of corresponding negatively-charged integrals. 5 The (Q-linear combinations of) iterated integrals considered in this work are homotopy functionals, i.e. they are invariant under continuous deformations of the integration path γ. While this is in general not true for the separate iterated integrals shown in eq. (3.19), it does hold for the combinations given in eq. (3.18), since the latter arise as the solution of DEs which satisfy integrability conditions. We assume that the same path γ is used for all iterated integrals, and that it lies entirely in the analyticity region of the integrals P + 45 . The considered iterated integrals therefore do not depend on the specific choice of γ.

Initial values for the differential equations
In order to construct the DE solutions (3.18) up to transcendental weight w we need to determine the initial values g (w) τ,σ (X 0 ). Since our aim is to find a function basis in the space of solutions, we have to take the algebraic relations between the initial values into account. Our approach to this end is the following. We evaluate the master integrals at an initial point with sufficiently large accuracy, and use the PSLQ algorithm [72] to identify integer relations among them. This allows us to extract a generating set κ of algebraically independent transcendental constants over Q, and to express the initial values g We choose an initial point X 0 ∈ P + 45 , which satisfies the following requirements: 1. X 0 introduces a minimal number of distinct prime factors.
2. X 0 is invariant under the automorphisms of P 45 , In other words, it is invariant under the exchanges of momenta 2 ↔ 3 and 4 ↔ 5.
3. The four letters which have indefinite sign inside P 45 and depend linearly on the Mandelstam invariants vanish at X 0 (see appendix D).
In order to understand the first two requirements we recall that the master integrals considered in this work can be expressed in terms of multiple polylogarithms whose arguments are rational/algebraic functions of the external kinematics [42,43]. It is natural to expect that many of these arguments coincide if the phase-space point is degenerate. Choosing an initial point which is as symmetric as possible therefore reduces the number of distinct multiple polylogarithms in the initial values. Minimising the number of distinct prime factors in the initial point -and hence in the rational arguments -further simplifies the arguments of the polylogarithms evaluated there. These requirements therefore decrease the number of algebraically-independent initial values and thus facilitate the discovery of the relations through the PSLQ algorithm. We will make use of the last two properties in sections 4.2 and 5, where we construct an explicit representation of our function basis and design an algorithm for their numerical evaluation. The next step consists in obtaining high-precision values of the master integrals at the initial point X 0 (3.22). We employ the analytic expressions in terms of Goncharov polylogarithms (GPLs) [39][40][41] constructed in refs. [42,43]. We use GiNaC [73,74] to evaluate numerically the GPLs with 3000-digit accuracy. We evaluate the permuted master integrals g τ,σ at X = X 0 by evaluating the expressions of refs. [42,43] for the standard ordering σ id of the external legs at permuted points, X = σ • X 0 , as shown in eq. (3.7).
The explicit GPL expressions of refs. [42,43] are given in an Euclidean region, while we are interested in the values of the Feynman integrals in the physical scattering regions.
We perform the analytic continuation by adding a small positive imaginary part to the values of the independent Mandelstam invariants at the initial point X 0 as where the superscript (0) denotes the values at the initial point X 0 (3.22), η is an infinitesimal positive parameter, and the c i 's are positive rational numbers. We choose the c i 's randomly so that all scalar products s ij (including the non-adjacent ones) have a positive imaginary part and all square roots in the alphabet are evaluated on the chosen branch (see eq. (3.21)). An additional complication we need to take care of is that some of the GPLs are divergent at X 0 (and permutations thereof). Since the master integrals have no singularity in P + 45 , these divergences are spurious and must cancel out. We handle this by using the parameter η in eq. (3.24) as regulator and compute the limit η → 0 + of the divergent GPLs analytically using PolyLogTools [75]. The cancellation of divergences at all permutations of X 0 is a consistency check of our GPL-evaluation routine. We further checked our results by evaluating all permutations of the master integrals also with the generalized power series expansion method [37]. We used the initial values provided in ref. [36] and integrated the canonical DEs with the Mathematica package DiffExp [38] at 15-digit accuracy. We evaluated the permuted master integrals by permuting the canonical DEs as discussed in section 3.3, rather than permuting the initial point as done for the GPLs. This way we verified also the consistency of the two approaches. We found full agreement.

Linear span (⊕ products) Irreducible
Weight Re Im Re Im  Table 3: Number of independent transcendental constants at each weight. The second column shows the dimensions of the linear spans of constants together with (weight-graded) products of lower-weight constants. The third columns shows the number of irreducible constants. * Only three real constants appear at weight one, but we add an extra constant to resolve reducible higherweight constants.
We construct the generating set κ (w) at weight w recursively starting from weight 0 where κ (0) ≡ {1}. We consider the Q-linear span C τ,σ (X 0 ), and weight-graded products of the lower-weight constants κ (w ) with w < w.
We employ an adapted version of the program PSLQM3 from the MPFUN2020 package [76], which implements the three-level multipair PSLQ algorithm [77,78] to search for integer relations using the initial values evaluated to 2000 digits. We then confirm the identified relations at the 3000-digit precision level. We then use these relations to construct a basis in C (w) 0 , preferring the lower-weight products. The remaining basis elements of C (w) 0 are then deemed irreducible and assigned to κ (w) . Note that we require the generating set κ to be real, while the initial values C We conclude this section with some remarks. First, we should stress that our numerical analysis cannot guarantee that the enumerated transcendental constants do not satisfy any additional algebraic relations over Q. However, as we explain in section 4.1, this would not jeopardize our function basis. Second, we emphasize that, while finding the relations among the transcendental constants is essential for the construction of the function basis, we do not need to know the constants in the generating set analytically. Furthermore, as we discuss in section 4.1, it is possible to absorb the transcendental constants into the definition of the basis functions such that only i π and ζ 3 explicitly appear in the expressions for the master integrals in term of these functions.

Function basis
In this section we present a basis of special functions in terms of which we can express all the planar two-loop five-particle amplitudes with one external off-shell leg up to order 0 . We dub these functions one-mass pentagon functions, and we denote them by {f i }, where w is the transcendental weight and i is a label. We begin by discussing our strategy for constructing function bases using their iterated-integral representations in section 4.1. Then, we present the expressions for the functions of the basis which are well suited for efficient numerical evaluation in section 4.2.

Construction of the basis using Chen's iterated integrals
The starting point are the canonical DEs (3.11) for the master integrals of the planar families shown in figure 1, each considered in all the S 4 permutations of the external massless legs. We solve the DEs in terms of iterated integrals up to order 4 , or equivalently up to transcendental weight four, with the base point given by eq. (3.22). The components of the expansion of the master integrals are not linearly independent, and the properties of the iterated integrals expose in a manifest way all the functional identities. Our aim is to identify a minimal set of linearly independent components, and to express any solution of the DEs as a polynomial in them. This minimal set constitutes the basis of one-mass pentagon functions, {f (w) i }. We refer the interested reader to ref. [31] for a thorough discussion of this procedure, and give here just a brief outline. Since the functional relations are uniform in the transcendental weight, we proceed recursively weight by weight. At each weight w, we find a basis of the linear space spanned by the weight-w master integral components over Q. The dimensions of these bases are shown in the second column of table 4. Comparing to the third column in table 2, we observe that at each weight the number of linearly independent integrals is less than the one implied by the topological identifications. Next, we use the shuffle algebra (3.20) obeyed by iterated integrals to remove the reducible functions, i.e. those functions which can be expressed as products of lower-weight functions. Assuming that the lower-weight pentagon functions are already found, {f such that w 1 +w 2 +. . .+w k = w, with 1 ≤ w i < w for i = 1, . . . , k and 1 ≤ k ≤ w. This way we can mod out the products from the weight-w basis. Repeating this procedure weight by weight up to weight four gives a basis algebraically independent functions {f (w) i }, which we call the one-mass pentagon functions. Their number is shown in the third column of table 4.
The iterated-integral representation of the master integrals in eq. (3.18) involves the initial values g (w) τ,σ (X 0 ). Thus, for our strategy to succeed, the relations between these transcendental constants have to be taken into account. We identified these relations in section 3.5 with a PSLQ-based analysis. As a result, it is in principle possible that certain algebraic relations were overlooked. Even if it were the case, the pentagon function classification of this section would not be affected. Indeed, the linear and algebraic independence of the functions persists upon the symbol projection [79][80][81], which is obtained by putting to zero all the non-rational constants in the iterated integral representation. In other words, the number of pentagon functions coincides with the number of linearly independent and irreducible symbols, and is therefore minimal.

Weight
Linearly independent  Irreducible  Irreducible,  cyclic   1  11  11  6  2  86  25  8  3  483  145  31  4 1187 675 113 Table 4: Counting of the master integral components at each transcendental weight. The second column corresponds to the number of the Q-linearly independent components. The third column shows the number of irreducible components, i.e. the number of functions in the basis. The fourth column shows the number of functions sufficient to express the master integrals in the cyclic permutations D 4 (3.10) only.
The one-mass pentagon functions inherit the Z ≥0 ×(Z 2 ) 4 grading from the pure master integrals (see section 3.4). Since the pure master integral normalizations involve at most one square root, the sets of functions with negative charge with respect to each square root are disjoint. The number of basis functions with negative Z 2 charge with respect to each of the square roots is shown in table 5. The charges of all one-mass pentagon functions are provided in the ancillary file pfuncs charges.m [70]. Weight 87 (8) 18 (9) 18 (0) 18 (1) Table 5: Number of one-mass pentagon-functions with negative Z 2 charge with respect to each of the square roots. In parentheses are the corresponding numbers for the cyclic subset specified in eq. (4.3).
The procedure described above to construct a basis out of the pure master integral components leaves a lot of freedom in choosing the basis elements {f (w) i }. We profit from this freedom to construct a basis which fulfills several additional constraints.
First, it is convenient to reduce as much as possible the amount of transcendental constants in the -components of the pure master integrals written as polynomials over Q of the pentagon functions. We absorb the majority of those constants in the definition of the one-mass pentagon functions. The only leftover transcendental constants are powers of i π and ζ 3 . For example, at weight two we have where by A Q we denote the linear span of the set of pentagon functions and constants A over Q. As a consequence, and in contrast with the basis of functions defined for the massless case in ref. [31], no constant in our basis has a non-trivial charge with respect to changing the signs of the square roots. We observe that, while the weight-four components of the master integrals contain products of all weight-one and two pentagon functions, only a subset of the weight-three ones -{f i } 120 i=4 -appear. The second useful constraint is the separation of the minimal subset of functions relevant for the cyclic permutations of the integral families. If one is interested in fully color-ordered scattering amplitudes, it suffices to consider only the permutations of the master integrals which preserve the cyclic ordering of the particles, namely σ ∈ D 4 (3.10). It is therefore meaningful to separate a minimal subset of the one-mass pentagon functions which contribute in the two permutations D 4 . We call them cyclic pentagon functions. They are counted in the fourth column of table 4. Only the 49 relevant letters of the cyclic alphabet A rel D 4 (3.13) contribute in the iterated integral expressions of these master integrals. Explicitly, the cyclic pentagon functions are The expressions of the master integrals g τ,σ in the cyclic permutations σ ∈ D 4 thus involve only these functions. The third property we made manifest stems from the observations made in the computations of the two-loop amplitudes in refs. [44,45,82]. Of the 49 letters of the cyclic alphabet A rel D 4 (3.13), the following six, drop out from the amplitudes truncated at order 0 , despite being present in the expressions of the two-loop master integrals. In the same works it was also observed that the letter W 198 = √ ∆ 5 , which is present in the amplitudes, drops out from the properly defined finite remainders. This is reminiscent of the behavior of the letter √ ∆ 5 of the massless pentagon alphabet [83], which has been observed to drop out from several massless two-loop fiveparticle finite remainders. This phenomenon is expected to be a general feature, and the first steps towards an explanation have been taken in the context of cluster algebras [84]. Based on the experience with the massless scattering, it is thus natural to conjecture that W 198 drops out of the finite remainders for any one-mass five-particle scattering process up to two loops. To make this explicit, we further refine the cyclic pentagon functions (4.3) by confining the letters Z ∪ { √ ∆ 5 } to a minimal subset of them. This way, expressing the known two-loop five-particle amplitudes with one external massive leg [44,45,82] in terms of our cyclic pentagon function basis (4.3), we avoid the spurious appearance of the letters Z altogether. Similarly, the finite remainders are free from the pentagon functions involving the letter √ ∆ 5 . In addition to analytic simplifications, this has the benefit of making the numerical evaluation more efficient, as we avoid evaluating functions involving spurious integration kernels which would otherwise lead to numerical cancellations.
We provide the expressions of the independent master integrals in terms of one-mass pentagon functions in the ancillary file mi2pfuncs.m [70]. The other master integrals can be obtained through the mappings in master integral mappings.m (see section 3.1). We spell out the one-mass pentagon functions in terms of Chen's iterated integrals and transcendental constants in pfuncs iterated integrals.m.
The one-mass pentagon functions are by construction closed under permutations of the external massless legs. In other words, any permutation σ ∈ S 4 of a pentagon function can be expressed in terms of (graded) polynomials in the pentagon functions, (4.5) This allows to evaluate the pentagon functions in any physical 2 → 3 scattering region with production of a massive particle following the same procedure discussed in section 3.2 for the master integrals. We provide all S 4 permutations of the one-mass pentagon functions in the folder pfuncs permutations/ of the ancillary files [70].

Explicit representations
In the previous section we discussed how we identified a minimal set of irreducible iterated integrals which is sufficient to write down the -expansion of any pure master integral up to transcendental weight four. For this purpose it was crucial to express the master integrals, and thus the functions in the basis, in terms of Chen's iterated integrals. In this section we present expressions for the basis functions which are well suited for their numerical evaluation. We give them in the ancillary file pfuncs expressions.m [70]. Following refs. [27,31] we provide expressions in terms of logarithms and dilogarithms up to transcendental weight two, and one-fold integral representations for the weight-three and four functions, as first suggested in ref. [85]. We constructed the expressions in terms of logarithms and dilogarithms using the strategy proposed in ref. [80] (see also ref. [31]).

Weight-one pentagon functions
We choose the 11 weight-one pentagon functions as follows,  such that only i π appears explicitly at weight one in the expressions of the pure master integrals in terms of pentagon functions.

Weight-two pentagon functions
From the analysis of the iterated integral solution of the canonical DEs for the master integrals we have identified 25 irreducible weight-2 functions. In order to present them in a more compact way, we introduce the following short-hand notations,  Then we have 6 weight-two functions whose iterated integral representation also involves alphabet letters which are quadratic in the Mandelstam invariants, (4.10) The first 5 of these functions are well-defined in the physical region P 45 , since the dilogarithm arguments cannot cross the branch cut within P 45 . In order to see that f 22 is well-defined as well we need to verify that the argument of the logarithm is positive inside P 45 . Indeed we find that s 14 s 15 −p 2 1 s 23 = W 33 , which in appendix D we prove to be positive in P 45 .
We choose the remaining 3 weight-two functions to be given by the so-called threemass triangles, i.e. they correspond (up to an algebraic prefactor) to the finite part of the one-loop triangle integrals with three-external off-shell legs (see e.g. ref. [86]). In order to write them down explicitly, we introduce the three-mass-triangle function, where ρ = ±1 and the Källen function λ is defined in eq. (2.10). In terms of this auxiliary function the last 3 weight-two pentagon functions are given by (4.12) Each of these three functions has a negative charge with respect to changing the sign of the three-mass triangle square root ∆ (i) 3 with the same arguments in eq. (2.9). Indeed, these are the three functions with non-trivial charge at transcendental weight two in table 5. In appendix D we show that the functions in eq. (4.12) are well-defined in P 45 .
Finally, let us recall that 8 weight-two pentagon functions, {f

Weight-three pentagon functions
At transcendental weights three and four we construct a representation of the pentagon functions in terms of one-dimensional integrals, which can be efficiently evaluated numerically [31]. The starting point is their expression in terms of Chen's iterated integrals, which for the 145 weight-three functions is given schematically by where τ . The criteria discussed in section 4.1 therefore do not forbid any weight-three pentagon function to appear in the amplitudes and finite remainders. Invoking the weight-one and two pentagon functions defined previously we can implement the two inner integrations in the three-fold integral in the first term on the right-hand side of eq. (4.13), and the inner integration in the second term. As a result, the weight-three pentagon functions take the following one-fold integral form, where c i,j,k are rational numbers, τ  k (X) are weight-two functions. The latter are polynomial in i π and the weight-one and two pentagon functions, The integration in eq. (4.14) is pulled back to the interval [0, 1] by an arbitrary contour γ, which connects the base point X 0 = γ(0) to an arbitrary point X = γ(1) ∈ P + 45 . To simplify the expression we used the short-hand notations The construction of the explicit contour γ is discussed in detail in section 5. For the purposes of this section it suffices to know that it lies entirely within the physical region P + 45 . For the sake of completeness, let us note that not all 108 letters (3.16) appear in the logarithmic integration kernel in eq. (4.14): the letter W 198 = √ ∆ 5 and 24 of the 54 letters which are quadratic in the Mandelstam invariants are absent.
Prior to attempting the numerical evaluation of the one-dimensional integrations in eq. (4.14) we need to verify that the integrations are well-defined for a path inside the physical region P + 45 . The functions h (2) k (X) are polynomials in the pentagon functions of weights one and two, and thus they do not have singularities nor branch cuts anywhere in P + 45 . We only need to inquire about possible singularities of the algebraic functions ∂ t log (W j (t)), which are located at W j (X) = 0. In appendix D we show that 20 letters can vanish inside P + 45 . We call them spurious letters. There are 4 spurious letters which are linear in the Mandelstam invariants,  While the linear spurious letters S vanish at the base point X 0 = γ(0), none of the quadratic spurious letters S q vanishes at the base point. Moreover, the weighttwo functions h k (X) which in eq. (4.14) accompany the logarithms of the spurious letters, log W j (t) with W j ∈ S, vanish on the locus where the corresponding letter vanishes, {X : W j (X) = 0}. The latter is straightforward to verify since the h (2) k (X) which accompany the d logs of spurious letters in the expressions of the weight-three functions do not involve the complicated three-mass-triangle pentagon functions. We can thus conclude that the logarithmic integration kernels in the weight-three one-fold integrals (4.14) may have at most a simple pole singularity, in the neighbourhood of t = t 0 such that W j (t 0 ) = 0, and that such simple poles are suppressed by the vanishing of the corresponding weight-two functions h For the linear spurious letters W j ∈ S the spurious singularity is at the base point t 0 = 0, while the quadratic spurious letters W j ∈ S q produce spurious singularities for t 0 ∈ (0, 1]. It is convenient to confine the weight-three functions with spurious singularities in a minimal set of 2 cyclic and 14 non-cyclic pentagon functions. This classification can be summarised schematically as We address the problem of the numerical evaluation in section 5.

Weight-four pentagon functions
The 675 weight-four pentagon functions are given in terms of iterated integrals by where τ 113 alone. Due to this choice, we expect the finite remainders of color ordered amplitudes to require only {f (4) i } 106 i=1 out of the 675 weight-four pentagon functions, as is the case for the results already available in the literature [44,45,82]. Beyond the cyclic sector, the letter W 198 = √ ∆ 5 is also present in the pentagon functions {f (4) i } 675 i=665 . We therefore expect that these 11 weight-four functions, together with the cyclic function f (4) 113 , are not required in the finite remainders. Of course, since we want to be able to evaluate all the master integrals and not only the finite remainders, we consider all the 675 weight-four pentagon functions in what follows.
We now proceed to constructing a one-fold integral representation for the weight-four pentagon functions. Instead of expanding the weight-four pentagon functions in terms of iterated integrals as in eq. (4.24), we represent them as one-fold d log integrals of the weight-three pentagon functions, where the weight-three pentagon functions f with their one-fold integral representations (4.14) in eq. (4.25) we rewrite the weight-four functions as two-dimensional integrals over a triangular domain, In order to obtain a one-fold integral representation, we swap the order of the integrations in eq. (4.27) [85]. This way one of the integrations becomes trivial, and we obtain the desired expression in terms of one-fold integrals, We rely on this representation of the weight-four pentagon functions for their numerical evaluation. In this view, we need to verify that the integrations in eq. (4.29) are welldefined. As argued for the analogous expression of the weight-three pentagon functions given by eq. (4.14), the only possible pole singularities originate from the algebraic factors ∂ s log(W k (s)). The spurious letters, W k ∈ S, may in fact vanish on the integration path γ. The weight-four pentagon functions however inherit the pole cancellation mechanism from the weight-three ones, see eqs. (4.21) and (4.22). There is an additional complication with respect to the weight-three case. Despite the pole singularity suppression, there are integrable logarithmic singularities, which arise from those terms in eq. (4.29) with W j ∈ S, so that L j (s) provides a logarithmic singularity. 6 We discuss the treatment of these integrable singularities in view of the numerical integration in section 5. The last piece of the weight-four pentagon functions which needs to be discussed are the logarithmic factors L j defined by eq. (4.28). In appendix D we study the behavior of the alphabet letters in the physical region P + 45 . If the letter W j is real-valued and has constant sign along the path γ, then L j (s) = log W j (X) W j (γ(s)) . (4.30) The complex-valued letters, {W j } 137 j=130 ∪ {W j } 188 j=186 , also do not pose any problem. They take values on the unit circle for any X ∈ P + 45 (see appendix D.4). We can define their phases ϕ j (X), to be continuous functions ranging in the interval (0, 2π), so that the integration in eq. (4.28) is also trivial, Finally, we argue that for the spurious letters W j ∈ S, which are all real-valued, we can choose Indeed, if a spurious letter W j (t) takes opposite signs at t = 1 and t = s, L j (s) gets a branch contribution which depends on the prescription for the integration contour deformation. However, the branch contribution cancels out among the first and second term on the right hand side of eq. (4.29). The weight-three pentagon function f k (t) accompanying ∂ t log(W j (t)) in eq. (4.25) with spurious W j ∈ S vanishes at t = t 0 for each t 0 such that W j (t 0 ) = 0, i.e.  thus suppressing the pole singularity originating from ∂ t log(W j (t)). This ensures the cancellation of the branch cuts of L j (s) in eq. (4.29), which can then be completely ignored. Note that the factors of L j (0) in eq. (4.29) corresponding to the linear spurious letters W j ∈ S , which vanish at the base point X 0 , are absent in the weight-four pentagon functions. The remaining letters do not vanish at the base point, so that the corresponding L j (0)'s are well-defined.
So far we tacitly assumed that the weight-four pentagon functions f (4) i (X) are evaluated at a phase-space point X = γ(1) which is not in a zero locus of the spurious letters. This special case however deserves some care, as L j (s) is divergent if W j (X) = 0. In fact, this singularity cancels out from the pentagon functions, which are not singular on any of the spurious surfaces {X : W j (X) = 0, W j ∈ S}. In order to show this, let us assume that the spurious letter W j vanishes on the path γ(t) at t = 1 + ε with |ε| 1, so that |W j (1)| 1. In other words, we assume that the integration along γ ends in the proximity of a spurious singularity. Using eq. (4.35) with t 0 = 1 + ε and eq. (4.34), we can rewrite the contributions to the right-hand side of eq. (4.29) with j : W j ∈ S as where each term is manifestly finite as ε → 0. Indeed, the second term contains an integrable logarithmic singularity, and the last term vanishes at ε → 0. The remaining contributions to the right-hand side of eq. (4.29) with j : W j / ∈ S are clearly finite as well. As a result, the weight-four pentagon functions f (4) i (X) are finite also when evaluated at a phase-space point X where any of the spurious letters vanishes.
From the considerations above it follows that all the weight-four pentagon functions are well-defined in the entire physical scattering region P + 45 . Some of them have bounded and explicitly finite integrands. Others, which involve d log's of spurious letters S in their one-fold integral representation (4.29), demand additional work in view of the numerical integration. For this reason we isolate the latter in a minimal subset, as done already at transcendental weight three. This separation, on top of splitting the functions into cyclic and non-cyclic sectors, and isolating the letters Z (4.4) and W 198 = √ ∆ 5 , gives the following classification of the weight-four pentagon functions {f . Now that we have expressions for all the one-mass pentagon functions that are welldefined throughout the physical region P + 45 we can move on to discussing their implementation in a C++ library which allows for their fast and stable numerical evaluation. This is the topic of the next section.

Numerical evaluation
We implement the numerical evaluation of the one-mass pentagon functions in the C++ library PentagonFunctions++ [87], which already provides the numerical evaluation of the massless pentagon functions from ref. [31].
The main idea of our approach follows the one from ref. [31]. We implement the evaluation of the weight-one and two functions through their explicit representations given in sections 4.2.1 and 4.2.2. For the weight-three and four functions we employ the one-fold integral representations derived in sections 4.2.3 and 4.2.4. The one-fold integrals are computed numerically with the tanh-sinh quadrature [88], which guarantees fast convergence for the integrands that are analytic and bounded everywhere inside the integration domain, with the exception of the endpoints where integrable divergences can occur [88,89].
In PentagonFunctions++ we use an adapted implementation of the tanh-sinh quadrature from Boost C++ [90], and we generate optimized code for the integrands with FORM [91,92]. We implement numerical evaluation in three fixed-precision floating-point types: double, quadruple, and octuple, approximately corresponding to significands of 16, 32, and 64 decimal digits respectively. The last two floating-point types are available in qdlib [93]. We refer for more details to ref. [31].
In this section, we discuss the practical implications for our implementation of the two features of the one-mass five-particle phase space which are new with respect to the fully massless case. The first new feature is the fact that the phase space P 45 , as an open subset of the affine space with coordinates (2.3), is not starshaped, i.e. there is no phasespace point from which all other points are reachable by straight lines fully within P 45 (see appendix B for details). The second new feature is the existence of spurious singularities associated with the quadratic letters S q that can occur anywhere along the integration path, as opposed to only at the initial point.

Integration path
As we show in appendix B, some points X ∈ P + 45 cannot be reached from X 0 by a line segment γ ∈ P + 45 . Therefore, instead of line segments, we consider integration paths γ which are polygonal chains, i.e. connected series of line segments. Given a set of points {A 1 , A 2 , . . . , A n } in the affine space, we denote by [A 1 , A 2 , . . . , A n ] the polygonal chain with vertices A 1 , A 2 , ..., A n . To avoid explicit analytic continuation and crossing of physical singularities we require that γ lies entirely within the physical phase space P + 45 . We begin by noticing that a line segment γ j (t) : [0, 1] → [X 0 , X] lies within P + 45 if and only if it does not intersect the variety ∆ 5 = 0, namely if the pull back ∆ 5 | γ j (t) = 0 for t ∈ [0, 1]. Therefore, using Sturm's theorem to count the number of roots of the quartic equation ∆ 5 | γ j (t) = 0, we can determine if any given line is inside P + 45 or not. We can then construct a two-chain γ = [X 0 , X , X] with the following algorithm. (b) X is close to the region boundary or to subvarieties of vanishing S ; 7 3. If both line segments [X 0 , X ] and [X , X] lie within P + 45 , accept X and terminate, else go to 1.
Here step 2 is designed to avoid potential numerical instabilities arising from the path unnecessarily entering near-singular configurations, as well as to avoid unnecessary crossing 7 The closeness is determined by a technical cutoff on the absolute values of the relevant letters. of spurious singularities. Some care should be taken when the coefficients of the ∆ 5 (t) polynomial become small and the number of roots cannot be reliably calculated from Sturm's theorem in step 3. In these cases we reject the point and try another one.
This simple strategy allows us to quickly check thousands of points and find a suitable polygonal two-chain. In principle, more line segments may be required and the presented algorithm can be straightforwardly extended. However, in practice we observe that two connected line segments are sufficient provided the base point X 0 is given by (3.22), and rarely more than 100 attempts are needed to find a suitable point X . We thus conjecture that any point of P + 45 can be reached from X 0 by at most two connected line segments. If PentagonFunctions++ ever encounters a point for which this is not the case, an error message is generated and the evaluation is aborted. We observe that, from a sample of O 10 6 randomly generated points, only a small fraction (few percents) is not visible from X 0 . We can view this as a statistical measure of non-starshapedness of the physical region.

Spurious singularities
In sections 4.2.3 and 4.2.4 we demonstrated that all integrands of the one-fold integrals in eqs. (4.14) and (4.29) are either finite or have integrable divergences on subvarieties where any of the spurious letters W i ∈ S is vanishing. Hence we showed that the integrals involving spurious divergences are well-defined from the mathematical point of view. Nevertheless, their numerical evaluation does require special treatment to ensure adequate numerical stability.
Whenever one of the quadratic spurious letters W i ∈ S q has opposite signs at the base point X 0 and at the evaluation point X, the integration path γ(t) has to cross the spurious singularity W i = 0 for some t 0 ∈ [0, 1]. The corresponding one-form d log W i then has a simple pole at t 0 , which is compensated by the corresponding multiplier h(X) which vanishes at t = t 0 (see section 4.2.3). This commonly causes significant loss of precision due to numerical cancellations. We address this problem by evaluating the contribution proportional to d log W i in the neighborhood where W i is small as where h 1 , h 2 are the first non-vanishing coefficients of the series expansion of h(X) in W i . Therefore the spurious poles are always canceled analytically. We choose the threshold for when to switch into the series expansion such that the O W 2 i corrections are estimated to be smaller than the roundoff error at given precision.
At weight three the coefficients h 1 , h 2 are explicitly finite by construction. At weight four they can contain logarithmic singularities of the form log |W i | originating from eq. (4.34). These singularities are integrable. Nevertheless, they violate the assumptions of the tanhsinh quadrature if they are encountered anywhere except at the endpoints. If left unattended, they spoil the expected double-exponential convergence rate and the associated error estimates. To avoid this issue, we further subdivide the integration chain at the locations of the spurious singularities such that they stand at the endpoints by construction. We do this on-demand for each pentagon function individually. As in the case of linear spurious singularities at X 0 , we must additionally ensure that, close to the spurious singularity, the vanishing letters are evaluated in a numerically stable manner, and that the rounding errors do not cause the integrand to be evaluated exactly at the singularity. We accomplish this by evaluating the letters that can vanish along a line segment γ j as where t 0 is a solution of W i | γ j (t) = 0, and the coefficients c (k) i;γ j are computed explicitly from the letter's derivatives. If a letter W i ∈ S q vanishes more than once along the same line segment, we construct multiple representations of the form (5.2) for each of the singularity locations t k . Then for any t ∈ γ j we select the representation which corresponds to the singularity that is the closest, i.e. such that |t − t k | is minimized.
The same considerations apply also to the case of linear spurious singularities W i ∈ S except that their location is always fixed to t 0 = 0 by construction. The representation in eq. (5.1) remains unchanged, while in eq. (5.2) the coefficients c (2) i;γ j = 0. Finally, let us note that the same considerations apply to the case when the point X itself is close to the zero locus W i = 0. As far as floating-point evaluations are concerned, no extra precautions are needed in this case, as long as W i (X) = 0 exactly. Indeed, the remaining spurious divergence is only logarithmic, hence it does not become large even for W i as small as machine epsilon. If a specifically fine tuned phase-space point is of interest, we refer to the representation in eq. (4.36).

Performance
To characterize the performance of our implementation we sample all one mass pentagon functions on points from a physical phase-space distribution that is representative of the expected phenomenological applications. For the purpose of demonstration, we construct a phase-space distribution by requiring optimal Monte Carlo integration of the leading order pp → eν e jj production cross section at √ s = 13TeV with Sherpa 2.2 [95]. We define the phase space by requiring two anti-k T jets with R = 0.4 [96], and we apply rather loose cuts, where p l T and p j T are the transverse momenta of the leptons and jets respectively, M W is the dilepton invariant mass, and η e , η j are the pseudo-rapidities of the electron and jets respectively. The renormalization scale is chosen to be half the scalar sum of the transverse momenta of all final-state particles.
We evaluate all pentagon functions in double and quadruple precision on 10 5 phasespace points from this distribution. We characterize the accuracy of the double-precision evaluationf (w) i;double (X) of a pentagon function on a kinematical point X by the logarithmic relative error r w,i , i;quad (X) is the numerical evaluation of the same function in quadruple precision. For convenience we will refer to this quantity as correct digits. We define the smallest number of correct digits (i.e. worst precision) among all pentagon functions at the kinematical point X as w, i ∈ {all functions up to weight four}. (5.5) We display the distribution of R(X) over the phase space as well as the average evaluation time in fig. 3. We observe excellent precision in the bulk of the phase space. A few of the phase-space points in the tail with low number of correct digits can be associated to the fact that the evaluation of tr 5 on those points becomes ill-conditioned. Comparing to the case of the massless pentagon functions in [31], the evaluation time is even lower despite the more complex phase space and the larger number of functions. We attribute this to the implementation details of the C++ library. This fact does however demonstrate that the overhead from integrating over multiple line segments is indeed negligible.  To characterize the scaling of our method with precision, we show in table 6 the timings of evaluation of all one-mass pentagon functions in double, quadruple, and octuple precision on a typical phase-space point. The number of correct digits is calculated by comparing with the results obtained from DiffExp as described in section 3.5.
The termination condition for the numerical one-fold integration for the benchmarks in this section is chosen as specified by the default settings of PentagonFunctions++: it is picked in such a way that the precision of the integration result is only limited by rounding Precision Correct digits Timing (s)   double  12  0.19  quadruple  28  159  octuple 60 1695 Table 6: Evaluation times of all one-mass pentagon functions on a typical phase-space point. The evaluations are performed on a single thread of Intel(R) Xeon(R) Silver 4216 CPU @ 2.10GHz.
errors in the integrand and the condition number which characterizes the sensitivity of the evaluation result to small perturbations in the input. One can therefore adjust this threshold to improve average evaluation time if lower precision in the bulk is sufficient for the desired applications. For example, it might be useful if quadruple precision is necessary only to overcome rounding errors in the integrand, while the target precision is still that of double precision. Finally, let us emphasize that in this section we considered the complete set of onemass pentagon functions sufficient to express any two-loop five-point scattering amplitude with four massless and one massive external momenta. It is however expected that not all of the functions contribute in particular applications. For instance, it has been observed that the functions involving certain subsets of letters drop out from suitably defined finite quantities [44,45,84]. Moreover, some scattering processes require only the subset of pentagon functions corresponding to the cyclic permutations (see e.g. [44]). Both of these properties are manifest in our function basis. Hence, the results presented in this section can be considered as the worst-case scenario analysis.

Validation
We validated our results against three sources: the function basis of ref. [44], the benchmark values for the pure master integrals in ref. [36], and the numerical evaluation of the master integrals using DiffExp.
We verified that the basis functions of ref. [44], which cover only the cyclic permutations D 4 of the planar integral families, can be expressed as (graded) polynomials in our cyclic one-mass pentagon functions, meaning that there is full compatibility between the two bases. We made use of this relation to cross-check the numerical values of the cyclic one-mass pentagon functions obtained with PentagonFunctions++ against those of the functions of ref. [44] in several kinematic points, including some in 2 → 3 physical scattering regions different from P 45 . We evaluated the one-mass pentagon functions at the latter points in double precision as discussed in section 3.2, using the analytic permutations of the pentagon functions provided in the ancillary files [70]. We found full agreement.
Reference [36] provides the numerical values of the pure master integrals of the planar families in the standard orientation of the external legs in six kinematic points, one for each distinct 2 → 3 physical scattering region with the massive particle in the final state. We reproduced these values through the one-mass pentagon functions evaluated with octuple precision and crossed to the scattering regions different from P 45 as discussed in section 3.2.
Finally, we used DiffExp to evaluate numerically all master integrals at several kinematic points in P 45 , including some close to spurious singularities. We used the benchmark values of ref. [36] as initial values. We evaluated the pure master integrals of the permuted families by permuting the evaluation point as in eq. (3.7) and flipping the sign of the parityodd integrals for the odd-signature permutations. We then matched the expressions of the pure master integrals in terms of the one-mass pentagon functions with their numerical values, and solved the ensuing system of equations to get the values of the pentagon functions. The fact that this system of equations admits a unique solution confirms the correctness of the expressions of the integrals in terms of pentagon functions. We then cross-checked the resulting values against the evaluations done with PentagonFunctions++ in octuple precision, finding full agreement to the expected accuracy.

Usage
The library PentagonFunctions++ can be used either as a C++ library or through its Mathematica interface implemented by the package PentagonFunctions`. We give here only a brief demonstration and refer to the documentation supplied together with the library distribution [87].
The library accepts as input the kinematical point given by the six Mandelstam invariants in eq. (2.3). The latter are assumed to be in ratios to the regularization scale, see eq. (3.1). The pentagon functions are evaluated in their analyticity region P + 45 defined in eq. (2.19), which means that Im tr 5 > 0 by construction. Function values for parity-conjugated points with Im tr 5 < 0 can be obtained by manually changing the sign of the parity-odd functions, which we list explicitly in the supplementary file pfuncs charges.m [70].
The interface of the library is designed in such a way that for each needed pentagon function the user must first obtain a callable function object or evaluator. In the C++ interface the evaluator of numerical type T is obtained by creating an instance of FunID<KinType::m1> type, and invoking its template method get_evaluator<T>(). For example, with the code 14 . In line 4 the function is evaluated on a point X specified as in eq. (2.3). It is worth noting that the created evaluators are completely independent of each other and will each evaluate only the requested pentagon functions. This means, for instance, that the user can choose to evaluate only a needed subset of pentagon functions.
Similarly, in the Mathematica interface this can be achieved with Here the second argument is a list of required functions, and the last argument can be used to select the required precision. In the Mathematica interface the functions are evaluated in parallel using all available threads. In line 3 the evaluation result is returned as replacement rules.
In both interfaces m1 is used as an identification for the set of one-mass pentagon functions introduced in this paper, 8 and the evaluator objects can be used for evaluations on any number of subsequent phase-space points. For more details we refer to the examples provided with the library distribution.

Conclusions
We constructed a basis of transcendental functions which allows to express any planar one-and two-loop five-particle Feynman integral with a single external massive leg up to transcendental weight four. This is sufficient for the computation of the planar twoloop corrections to cross sections. We call the functions in our basis one-mass pentagon functions. Our results are valid throughout the entire physical phase space for any 2 → 3 scattering process with production of a massive particle. We achieved this by constructing expressions for the one-mass pentagon functions which are well-defined in a particular 2 → 3 channel, and by providing the expressions of the master integrals in terms of our function basis for all permutations of the external massless legs. We expressed the one-mass pentagon functions in terms of logarithms and dilogarithms up to transcendental weight two. For the weight three and four functions, instead, we constructed a one-fold integral representation which is well suited for the numerical integration. We presented a public C++ library [87] which allows to evaluate numerically the one-mass pentagon functions in a stable and fast way. The supplemental material can be obtained from the git repository [70].
The results presented in this paper open the door for broad phenomenological applications, and constitute a crucial step forward towards the computation of NNLO QCD predictions for several processes which feature prominently in the LHC physics program, such as the production of a Higgs or an electroweak vector boson in association with two jets. Indeed, the one-mass pentagon functions have already been instrumental in the recent calculation of the planar two-loop amplitudes for four partons and an electroweak vector boson [82], and we envisage that they will facilitate calculations of other processes with similar scattering kinematics in the near future.
More broadly, we expect that the approach we adopted in this paper will help to tackle many further problems in the future. In primis, we envisage that it will be possible to classify the non-planar one-mass pentagon functions in a completely analogous way, once the canonical DEs for all the non-planar integral families are available (see refs. [47,48] for progress in this direction). Similarly, the expressions for the one-mass pentagon functions may be generalized to the 1 → 4 decay region, which would for instance enable a two-loop description of e + e − → 4j at lepton colliders. More generally, our approach can be useful in other multi-scale applications, in particular those where rationalizing all square roots is very difficult or even impossible.

A Gram determinants and the physical region
In this section we discuss in detail the properties of the Gram determinants which are relevant to define the physical scattering region and to determine the signs of the letters of the alphabet. Our discussion essentially follows refs. [59][60][61], specializing to the case of massless particles. We begin by defining the Gram determinant of two sets of momenta, G(q i 1 , . . . , q in ; q j 1 , . . . , q jn ) = det G(q i 1 , . . . , q in ; q j 1 , . . . , q jn ) , where G is the Gram matrix, We call the number of momenta in the two sets, n in eq. (A.1), the order of the Gram determinant. The Gram determinant is invariant under the exchange of the two sets of momenta, which corresponds to a transposition of the Gram matrix. We further introduce a short-hand notation for the Gram determinants of the external momenta The Gram determinant vanishes if the momenta in any of the two sets are linearly dependent, and is invariant under shifts of any of the momenta by any linear combination of the other momenta in its argument. As a result, when considering four-dimensional external momenta, the Gram determinants with order n ≥ 5 vanish. Moreover, the order-n Gram determinant is invariant under any permutation σ ∈ S n of the momenta in its argument, Finally, Sylvester's determinant identity implies a number of relations among Gram determinants with different orders which will be useful in the following discussion, Now that we have stated the main properties of the Gram determinants we can move on to discussing how they constrain the physical region for five-particle scattering [59][60][61]. In particular, we prove the constraints given by eq. (2.17). Let us consider the Gram determinants G i 1 ...in , with n = 1, . . . , 4 and distinct indices i 1 , . . . , i n . The corresponding Gram matrix G(p i 1 , . . . , p in ) is symmetric and thus has n real eigenvalues. Clearly if any of the eigenvalues is equal to zero, the determinant G i 1 ...in also vanishes. We can prove that the Gram matrix G(p i 1 , . . . , p in ) is either identically zero -which is of no interest -or it has exactly one positive and n − 1 non-positive eigenvalues. Since the Gram determinant G i 1 ,...,in is the product of all eigenvalues of the corresponding Gram matrix, this property straightforwardly implies the constraints given by eq. (2.17), which we rewrite here for the convenience of the readers, We begin by proving that the Gram matrix G(p i 1 , . . . , p in ) cannot have more than one positive eigenvalue. If it had two positive eigenvalues, the corresponding eigenvectors would be real linear combinations of p i 1 , . . . , p in with the following properties: linear independence, orthogonality, and time-likeness. This is impossible in the (1, 3) Minkowski space-time metric. Therefore either all n eigenvalues of G(p i 1 , . . . , p in ) are non-positive, or one eigenvalue is positive and n − 1 are non-positive.
Secondly, we prove that the Gram matrix G(p i 1 , . . . , p in ) has at least one positive eigenvalue or it vanishes identically. This follows from physical constraints on the kinematics. In particular, we assume that the momentum p 1 is time-like, p 2 1 > 0. The sum of the eigenvalues of G(p i 1 , . . . , p in ) equals its trace, i.e. n k=1 p 2 i k . Then either i k = 1 for some k, in which case the trace of the Gram matrix evaluates to p 2 1 and is thus positive, or no index is equal to 1, in which case the trace vanishes. The sum of the eigenvalues is thus non-negative. As a result, either all the eigenvalues are zero or at least one eigenvalue is positive. Since the Gram matrix G(p i 1 , . . . , p in ) can have at most one positive eigenvalue, this implies that it has exactly one positive and n − 1 non-positive eigenvalues, or that it vanishes. q.e.d.
The study of the positivity of the alphabet letters in appendix D makes use of another useful identity, with distinct indices i, j, k taken from {2, 3, 4, 5}. This relation holds in the physical scattering regions, i.e. where the Gram-determinant constraints (A.9) are satisfied. We prove this by showing that G(p 1 , p i + αp j + βp k ) is a non-positive polynomial in α and β, Using the identity (A.6) we see that the determinant of the Hessian matrix of this polynomial is non-negative, This, together with the physical constraint that G 1j < 0, implies that the polynomial has a global maximum where the first derivatives in α and β vanish. With the help of (A.6) and (A.8), we rewrite the maximum in terms of the Gram determinants so that it is explicitly non-positive, which implies the inequality (A.10). q.e.d.
In order to understand the implications of the Gram-determinant constraints (A.9) in terms of the scalar invariants it is useful to spell out the various Gram determinants. We begin with the order-one Gram determinants. Only one of them is non-vanishing, whereas the other four vanish identically, (A.14) Higher order Gram determinants vanish only in degenerate configurations of the momenta corresponding to the boundaries of the physical region defined by eq. (2.16) (see also footnote 2). The order-two Gram determinants are given by Clearly the G ij 's, for distinct i and j, can vanish only where p i · p j = 0, which indeed defines the boundary of the physical regions. The order-three Gram determinants have different expressions depending on whether one of the indices is equal to 1, where i, j, k take distinct values in {2, 3, 4, 5}. In order to show that the order-three Gram determinants vanish only on the boundaries of the physical regions we make use of eq. (A.7). Setting G ikl = 0 in eq. (A.7) in fact implies that Because of the Gram-determinant constraints (A.9) in the physical scattering regions, the left-hand side is non-negative, while the right-hand side is non-positive. Hence, both sides of the equation vanish. The vanishing of the order-three Gram determinant G ikl = 0 thus implies that either G kl = 0 or G ijkl = 0. Both cases correspond to boundaries of the physical regions. We have already proven this above for the order-two Gram determinants and we discuss it presently for the order-four ones. The latter are all proportional to ∆ 5 = tr 2 5 , which clearly can only vanish at the boundary ∆ 5 = 0 of the physical region (see also footnote 2). Finally, let us establish the inequalities in eq. (2.18), which define the s 45 scattering region in terms of the Mandelstam invariants s ij and p 2 1 . We begin by rewriting the physical constraints on the scalar products (2.16) in terms of the Mandelstam variables, (A. 20) We see that the latter are less strict than the constraints given in eq. Given two points x, y in some affine space A, we say that y is visible from x (or vice versa) if the line segment connecting x and y lies in A. First of all, we note that one can easily find two points within P 45 which are not visible from each other, i.e. P 45 is not convex. A more refined question about the geometry of P 45 , which has important practical implications, is whether it is starshaped. A set A is called starshaped if it contains a point x ∈ A such that all points in A are visible from x. The point x is then called a star center of A. We refer to ref. [98] for a detailed overview of properties of starshaped sets.
We are going to prove that P 45 is not starshaped. We will assume that P 45 is starshaped and show that this assumption leads to a contradiction. The domain P 45 owes its nontrivial geometry to the quartic constraint ∆ 5 < 0. The remaining constraints in eq. (2.18) are in fact linear, so if they are satisfied at points X 1 and X 2 then they are also satisfied along the line segment connecting X 1 and X 2 . For this reason we will focus on the quartic constraint ∆ 5 < 0 hereafter.
Let us assume that P 45 is starshaped and let the point X 0 ∈ P 45 be a star center. Since the inequalities defining P 45 are homogeneous, we parametrize X 0 without loss of generality as One can easily check that the inequalities (2.18) are satisfied for z 1. Physically, the asymptotics z 1 is the Regge limit with high energies s 13 ∼ s 45 ∼ z. We connect X 0 and X (z) Regge by a line segment γ z (t) parameterized with 0 ≤ t ≤ 1. Along the line segment the leading term of ∆ 5 for z 1 is The requirement that X (z) Regge is visible from X 0 implies the negativity of ∆ 5 (γ z (t)) for all t ∈ [0, 1]. Consequently, the leading term in the z-expansion (B.3) has to vanish, giving the following constraint on X 0 , In order to obtain more constraints on X 0 we consider those permutations of the one-parameter family (B.2) which are automorphisms of the domain P 45 , namely the permutations in the set Σ (P 45 ) (3.23). Since ∆ 5 is invariant under all S 4 permutations, we automatically obtain three more constraints on X 0 : where we tacitly imply that s Combining together the four constraints on X 0 we restrict it to the following oneparameter family, This asymptotics corresponds physically to the collinear limit p 2 ||p 3 . We connect X 0 and X col . One can then verify that the quadratic constraint ∆ 5 < 0 cannot be satisfied for all t ∈ [0, 1] for z 1. We do this by showing that, for z 1 and for any a > 0, the sign of ∆ 5 (γ z (t)) changes at some value of t ∈ (0, 1). For this purpose we find it convenient to pull back ∆ 5 to [0, +∞) by the map t = u 1+u with u ≥ 0, Clearly, the sign of the leading term of this expansion changes as u varies along the positive axis for any value of a. In other words, ∆ 5 cannot be negative along the entire line for any choice of a > 0, and X (z) col is not visible from any X 0 of the form (B.6). Therefore we can conclude that P 45 is not starshaped.

C Permutation closure of the planar alphabet
In this appendix we spell out the expressions of the relevant letters A rel S 4 listed in eq. (3.16). Following ref. [48], we present them grouped into orbits of the permutation group S 4 starting from a generating set of letters. The generating letters may be invariant under certain permutations, so that only a subset of S 4 is required to generate the full orbit. We give in the ancillary file alphabet permutation orbits.m [70] the subsets of permutations which, starting from the generating letters shown here, give the letters in the corresponding orbits in the correct ordering.
We also highlight the behavior of each letter with respect to changing the sign of one of the square roots. We recall that a letter W i is called odd with respect to a square root √ δ if d log W i changes sign when changing the sign of √ δ, namely if W i goes to its inverse 1/W i . The letters with non-trivial behavior with respect to changing the sign of the square roots are listed in table 7.
The purely rational relevant letters are given by 9 The letters with non-trivial behaviour with respect to changing the sign of one of the three-mass-triangle square roots but independent of tr 5 are 3 and ∆ 3 , respectively. The parity-odd letters -namely those whose d log changes sign when changing the sign of tr 5 10 -which are free from the three-mass-triangle square roots are The letters {W19, . . . , W24} require also a factor of −1 on top of the permutation of the generating letter. 10 In ref. [48] the letters are defined in terms of √ ∆5, rather than tr5, and are thus even under parity. We adopt their definition of the letters but replace √ ∆5 with tr5, to match the conventions used in the computation of the planar integral families [36]. Since √ ∆5 is invariant under any permutation of the massless legs whereas tr5 changes sign with an odd-signature permutation, the sets of permutations used here may differ from those of ref. [48].  Table 7: Relevant letters which depend on the sign of the square roots of the problem. On each row the letters in the right column are odd with respect to changing the sign of the square root in the left column. We recall that tr 5 is related to √ ∆ 5 through eq. (2.5).
There are three parity-odd letters which also depend on one of the three-mass-triangle square roots, 3 and ∆ 3 , respectively. The definition of W 188 in eq. (C.5) requires an inverse because we are using tr 5 rather than √ ∆ 5 (see also footnote 10). The last four relevant letters are given by the four square roots of the problem, We stress that these letters are even with respect to changing the sign of any square root. D Positivity of the alphabet in the physical scattering region 45 → 123 The one-mass pentagon functions are expressed as iterated integrals with integration kernels given by the logarithms of the alphabet letters (see section 3.4). The representation in terms of iterated integrals highlights that the pentagon functions may have singularities only on the loci where any of the involved letters vanishes or diverges. Therefore we need to inquire about the positivity of the alphabet letters in the region P 45 , i.e. whether they have definite sign in the region P 45 . This has important implications for the construction of the explicit representation of the pentagon functions in section 4.2. In this appendix we study systematically the positivity of the 108 relevant letters of the alphabet listed in eq. (3.16 We denote by S the set of all spurious letters, We prove this by rewriting the letters in such a way that their positivity properties in P 45 become obvious consequences of the inequalities (2.18). For the sake of presentation we split the alphabet into several subsets of letters having similar structure, to each of which we devote a subsection. Finally, in subsection D.5 we show that the three-mass triangle functions defined in eq. (4.12) have no branch cuts within P 45 .  i.e. both terms in each of the previous sums are simultaneously either positive or negative.

D.1 Linear letters
In order to establish the positivity of the remaining quadratic letters, the linear inequalities in eq. (2.18) are not sufficient. We have to appeal to the Gram determinant inequalities. We thus express the letters in terms of order-3 Gram determinants, G 1ij , and profit from their positivity in P 45 , namely G 1ij > 0 for i, j = 2, . . . , 5 and i = j. In this way we can straightforwardly establish that the following quadratic letters do not change sign in P 45 ,  Among the square-root letters, the letters {W i } 129 i=118 contain only the square roots of three-mass-triangle type. We observe that they are neither singular nor vanish in the physical region P 45 , and their signs are summarized as follows: 3 . The latter constraint is equivalent to p 2 1 s 45 = 0, which is impossible inside the region P 45 . Therefore, W 118 cannot change sign in P 45 .
The second example is given by the letter W 126 , W 126 = s 12 − s 15 + ∆ . (D.14) The latter vanishes or becomes singular only at (s 12 − s 15 ) 2 = ∆ 3 , which is equivalent to s 12 s 15 = p 2 1 s 34 . Resolving the previous relation in one of the Mandelstam variables, e.g. p 2 1 = s 12 s 15 /s 34 , and eliminating this variable from ∆ 5 , we observe that ∆ 5 turns into a perfect square, which cannot be negative, in contradiction with the quadratic constraint ∆ 5 < 0 in eq. (2.18). For the remaining algebraic letters involving only the three-mass triangle square roots, one of the two arguments above is applicable.

D.4 Parity-odd letters and
√ ∆ 5 The remaining relevant letters that we need to consider, Moreover, these letters cannot be equal to 1, as that would imply ∆ 5 = 0. In summary, for the parity-odd letters we have that

D.5 Three-mass triangles
In eq. (4.12) we defined three weight-two pentagon functions corresponding, up to an overall algebraic factor, to the finite part of the three-mass triangle Feynman integrals. Their explicit expressions (4.12) rely on the representation (4.11) of the three-mass triangle function. Here we show that they are well defined, namely real analytic, in the physical channel P 45 .