Analytic result for the nonplanar hexa-box integrals

In this paper, we analytically compute all master integrals for one of the two non-planar integral families for five-particle massless scattering at two loops. We first derive an integral basis of 73 integrals with constant leading singularities. We then construct the system of differential equations satisfied by them, and find that it is in canonical form. The solution space is in agreement with a recent conjecture for the non-planar pentagon alphabet. We fix the boundary constants of the differential equations by exploiting constraints from the absence of unphysical singularities. The solution of the differential equations in the Euclidean region is expressed in terms of iterated integrals. We cross-check the latter against previously known results in the literature, as well as with independent Mellin-Barnes calculations.


Introduction
Scattering amplitudes for multi-particle processes start to play an increasingly important role in future collider physics analyses, as processes at higher multiplicity are being probed more and more accurately. Recently, rapid progress has been achieved for five-particle processes at next-to-next-to-leading order. This concerns several areas, such as the efficient computation of loop integrands [1][2][3][4], the analytic computation of the Feynman integrals [5][6][7] as well as advances in integral reduction techniques [8][9][10][11][12][13][14][15][16]. Most recently, two independent numerical determinations of all planar five-gluon scattering amplitudes [4,17] have been achieved.
Non-planar corrections are unfortunately considerably more difficult to handle, due to a variety of reasons. Owing to the richer cut-structure of non-planar amplitudes, they can contain a larger number of rational factors in the external invariants, leading to more complicated algebraic expressions, both in the integrand reduction and in the determination of the integrals. This article addresses the second challenge, specifically at the level of the computations of non-planar Feynman integrals. The first steps in this direction were taken in [18], where three of the present authors conjectured the function space describing the Feynman integrals, and proposed a bootstrap method for determining the functions. Furthermore, individual integrals were computed in ref. [19], using a method based on conformal symmetry. There are two non-planar integral families for five particles at two loops, namely the hexa-box integral family a) and the double pentagon integral family b), shown in figure 1. In this paper, we analytically compute all master integrals for this first one, namely the hexa-box integral family a).
We begin by deriving a basis of integrals with constant leading singularities, also known as d-log integrals [20][21][22]. This is done by adapting the algorithm described in [23] to the five-particle kinematics. We then use integral reduction programs to find a basis of 73 d-log integrals.
We follow this up by computing the differential equations for the basis integrals, and find that they obey the canonical form of [21], as expected by the conjecture made therein. We find that the differential equations can be expressed in terms of the non-planar pentagon alphabet of reference [18].
Having obtained a system of first-order differential equations, the solution is fully specified by providing a complete set of boundary constants. We do so by deriving constraints from the absence of unphysical singularities. In this way, we obtain analytical constraints for all boundary constants (up to an overall normalization, which is fixed by a trivial calculation). The constraints are written in terms of Goncharov polylogarithms. We evaluate the latter to high numerical precision, and give the solutions with 100-digit accuracy. This fully determines the solution of the differential equations, which can be expressed in terms of iterated integrals. These are straightforward to evaluate in the Euclidean region, as documented in detail for example in [7]. To obtain the results in the Minkowskian region requires either an analytic continuation of the results, or an independent determination of the boundary conditions in each Minkowskian scattering channel.
We validate our solution by comparing it to previously known results in the literature, for subtopologies that are planar or that correspond to four-point functions, as well as against an independent Mellin-Barnes calculation described in appendix A.
The paper is organized as follows. We begin in section 2 by describing our notation and the kinematics of the problem. We also discuss the integral reduction to master integrals, and the differential equation satisfied by the latter. Then, in section 3, we explain the determination of the d-log basis. Section 4 is dedicated to the canonical differential equations and their analytic solution. The appendix A discusses checks performed on the results. Finally, we conclude in section 5.

Kinematics and notation
We denote the momenta of the on-shell particles in pentagon kinematics by p µ i , i = 1, . . . 5, with p 2 i = 0. Momentum conservation reads 5 i=1 p µ i = 0. We introduce the following five independent Mandelstam variables, Here µ is an arbitrary scale, e.g. the scale appearing in the dimensional regularization, making the v i dimensionless. In the following, we will set µ 2 = 1 GeV without loss of generality, as the dependence on µ can always be restored by dimensional analysis. We parametrize the integrals of the integral family shown in figure 1a) using the following notation for the individual integrals. In the above, a 1 , . . . , a 8 ≥ 0 are propagators and a 9 , a 10 and a 11 ≤ 0 numerator factors.
To perform the integral reduction [24] for this hexa-box family, we use the program Re-duze2 [25], which yields a basis of 73 master integrals. Aiming for the differential equations for the hexa-box family, we need to go beyond the reduction of integrals with unit powers on all propagators (which was accomplished previously, [15,16]), which are sufficient for scattering amplitudes, and include the reduction of integrals with single squared propagators. To limit the size of intermediate algebraic expressions in this reduction, we perform independent reductions on spanning cuts, i.e. by projecting onto subspaces of integrals that are required to contain a specific combination of propagators. The hexa-box family has in total 11 spanning cuts (single-scale three-point or four-point subtopologies that each do not contain any further subtopologies). A sufficient practical mitigation of the complexity of the integral reduction can be achieved by combining them into four spanning topologies, identified by two-particle cuts, i.e. requiring the non-vanishing of (a 5 , a 7 ), (a 5 , a 8 ), (a 6 , a 7 ) or (a 6 , a 8 ). The full integral reduction is then assembled by adding the cut-reductions, with individual integrals weighted by appropriate inverse multiplicity factors, which correct for their appearance in more than one cut-reduction tree.
The master integrals in the hexa-box family can be classified as follows. There are 54 planar integrals, 9 are non-planar with up to four external legs (four-point functions with one off-shell leg, which were computed in [26][27][28] in terms of generalized harmonic polylogarithms [29][30][31][32][33]) and 10 that are non-planar with five external legs. The latter type of genuine non-planar five-point integrals in the hexa-box integral family are depicted in figure 2. The second one, (h), can also be flipped upside down, and hence there are four such integral sectors. Together they have 3+3+3+1 = 10 master integrals. The reduction selects a basis of master integrals in each topology according to lexicographic ordering, typically JHEP03(2019)042 Figure 2. Non-planar integral sectors with genuine five-particle kinematics. The labelling follows that of [22].
containing the scalar integral and integrals with simple numerator factors. Differential equations for the hexa-box integrals in an alternative basis in terms of pure integrals (containing higher propagator powers) were derived most recently in [34]. We will be interested in a different basis, in which the integrals have a d-log form. Such d-log integrals have properties that significantly simplify their computation. In particular, in the ǫ expansion all such integrals evaluate to multiple polylogarithms of homogeneous weight. Determining this basis is the subject of section 3. We note already here that this basis choice can be done algorithmically [23] by analyzing just the loop integrand.

The alphabet
The 73 integrals that we shall compute can be expressed through iterated integrals of the type d log W i 1 · · · d log W i L , where the algebraic functions W i of the kinematic variables are called letters. The ensemble of the letters {W i } is called the alphabet of the problem under consideration. We recall the notation of [18], where the 31 letters of the non-planar pentagon alphabet were introduced. They fall into six classes of five letters W 1+i , W 6+i , W 11+i , W 16+i , W 21+i , W 26+i , with i = 0 . . . 4, that are mapped into each other by cyclic permutations together with one lonely letter W 31 . Explicitly, the first twenty letters are 3) while the next ten are

JHEP03(2019)042
and the last one is Here, ∆ is the Gram determinant that can be written explicitly as Note that the letters W i , with i = 26, . . . 30, are parity-odd, in the sense that they go to their inverse under √ ∆ → − √ ∆, while all other letters are parity-even under that transformation. Furthermore, only the letters can appear as first entries (see section 4.3 for information regarding the symbols). There is also a hypothetical second-entry condition that the integrals should obey that forbids some combinations of pairs of letters from appearing. We refer to [18] for more details.

The canonical differential equations
Using integration by parts identities (IBP), one can reduce the general integral (2.2), to a linear combination of a basis set of integrals I(v i ; ǫ), called master integrals. The next step is then to find a way to compute those master integrals. We can accomplish this by using the method of differential equations [21,[35][36][37], which works as follows. We first differentiate the set of master integrals I with respect to the variables (2.1). This can be done at the level of the integrals (2.2) by using appropriate derivatives in the external momenta, respecting the on-shell conditions. The derivatives obtained in this way can then also be expressed as a linear combination of the master integrals I, which means that we obtain a set of first order linear differential equations with five different matrices {A j } 5 j=1 that depend in a non-trivial way on the v i as well as on ǫ = 2 − D/2. Now, if the set of master integrals is chosen to be of a d-log form, as is discussed in section 3, then the differential equations simplify significantly. For such a basis, after combining the five derivatives in a 1-form, we obtain the following canonical form of the differential equations [21] with the matrix being independent of ǫ. We note that once (2.8) and the value of I at some boundary point are known, then the problem of computing the master integrals I at any kinematic point in an ǫ expansion is solved [21]. The value of the integrals at the boundary point will be derived in section 4. We wish to emphasize here that the construction of

JHEP03(2019)042
the canonical basis is done at the integrand level and as such does not require the a priori knowledge of the differential equation.
Finally, let us anticipate that one can write theÃ(v i ) matrix in a nice way by using the algebraic functions W i of section 2 as where theã i are constant 73 × 73 matrices (with rational entries). We remark thatÃ is independent of seven of the letters, namely of the letters 8, 9, 10, 21, 22, 23 and 24. The correspondingã i matrices are zero.

Construction of a basis of d-log integrals
In this section we describe how we obtained a d-log basis with constant leading singularities. An algorithm for doing this is provided in [23]. Let us briefly summarize the method. We start from a given propagator structure, in the present case that of part a) of figure 1.
Then, an ansatz for all possible numerator structures is made. The degree of the latter is constrained by the requirement of the absence of double poles. Computing all leading singularities of a general linear combination of such numerators, we obtain a complete solution of all d-log integrands for the corresponding propagator structure. We perform the analysis in four dimensions, expressing the loop momenta in a basis built from the spinor helicity variables of the external momenta. Furthermore, we find it convenient to parametrize the kinematics as in eq. (3) of ref. [38], as the latter rationalizes the Gram determinant √ ∆ (2.6) that can be built from four of the five external momenta. The computation of the leading singularities can be combined nicely with the computation of cuts. In the case of the maximal cut of the full topology there are no integration variables left, so we obtain the leading singularities in this case without further calculations. Computing the leading singularities on cuts has several advantages. First, it drastically reduces the amount of different leading singularities that have to be computed. Second, we can split the calculation in several smaller parts that can be computed in parallel. Third, we can choose for each cut an optimized parametrization of the loop momenta and this way minimize the appearance of square roots in intermediate steps.
In order to find a d-log solution in a given sector we proceed as follows: first, we compute the leading singularities on the maximal cut of that sector in order to get all solutions projected on that sector. Then, for each solution we add a linear combination of integrals of the subsectors and fix their coefficients by computing the leading singularities on that subsector. In this way, we can iteratively construct a list of d-log integrals.
As a check that the integrals obtained with this procedure are correct, we verified them all by computing the leading singularities for each solution individually without taking cuts. Along this way we also checked that the solutions for the hexa-box family provided in [22] are d-log integrals with constant leading singularities. For the verification we used a seminumerical approach, setting all but one of the external kinematical variables to numerical JHEP03(2019)042 values, thus proving that the leading singularity does not depend on the one kinematical variable that was not replaced by a numerical value. Repeating this for all external variables ascertains that a given possible solution has constant leading singularities. This seminumerical approach simplifies the calculation substantially.
We remark on a subtlety in this approach. As the above analysis is done in four dimensions, it cannot detect certain Gram determinants that vanish in four dimensions. Therefore the latter represent an ambiguity. While we expect a refined version of the leading singularity analysis to also fix this ambiguity, here we chose a pragmatic solution. We aimed for finding 'simple' solutions without the admixtures of Gram determinants (that necessarily involve many numerator terms, and hence typically integrals of several topologies). Unwanted and complicated solutions of this type can in most cases be easily identified and removed.
In this way, we obtained 157 d-log integrals for the hexa-box family. The integrals obtained are all expected to be pure functions of uniform transcendental weight [20,21]. We perform the following consistency check on this assertion. The number of d-log integrals (in our case 157) is much bigger than the number of master integrals (in our case 73). We first choose a set of linearly independent d-log integrals as a basis of master integrals. Then, we express the remaining integrals in terms of this basis, using the reduction obtained in section 2 above. If all integrals have uniform transcendental weight, then the basis coefficients must be numerical constants (for general Feynman integrals these coefficients would be functions of the external variables s ij = 2p i · p j and of D, the parameter of dimensional regularization). Indeed, we explicitly found that all relations involved constants only.

Determination of the boundary conditions
In this section, we will determine the boundary conditions of the differential equations (2.8) for the hexa-box integrals, such that their complete solution becomes uniquely specified. The method for computing the boundary conditions starts by picking a convenient reference point where the integrals are finite. Then, one integrates the differential equation along a path joining the boundary point with kinematic points where letters of the alphabet vanish and where singularities can thus appear. By demanding the absence of spurious singularities, we obtain constraints on the values of the integrals at the reference point. This turns out to be sufficient to determine the boundary conditions (up to a trivial overall normalization).

The origin point and spurious singularities
The Euclidean region is given by the conditions One may verify that in this case the Feynman denominator of any integral of the integral family under consideration is positive definite. This implies that the solution of the differential equation (2.8) is real within that region. The latter observation is useful, as the equations we will obtain are in general complex.

JHEP03(2019)042
In order to provide an explicit solution to the differential equation, we need to determine the value of I(v i ; ǫ) at one point. A suitable candidate is the point p E in the Euclidean region corresponding to setting v 1 = −3, v 2 = −3, v 3 = v 4 = v 5 = −1 and choosing the positive sign for the root, √ ∆ = 3 √ 5. Let us denote the value of I(v i ; ǫ) at p E as I E (ǫ). We can now impose conditions on the value of I E (ǫ) by demanding that the integrals stay finite when taking certain suitable limits. This is justified as follows: all integrals in the d-log basis are ultraviolet finite, by construction. We take ǫ < 0 in order to regulate the integrals in the infrared and consider limits in which some of the letters of the alphabet vanish. Taking such limits does not change the UV structure of the integrals, and so we require that the integrals remain finite in the limit, provided that ǫ < 0. In other words, we constrain the boundary condition I E (ǫ) by demanding that these spurious singularities are absent.

The paths to the spurious singularities
Let us now explain how this is implemented in detail. We begin by choosing an index j ∈ {1, . . . , 25} and considering the limit W j ≡ y → 0. Without loss of generality, let us take j = 11 to illustrate the situation. We decompose the matrixÃ as A =Ã sing log(y) + non-singular for y → 0 . (4.2) In the neighborhood of y = 0, the solution of (2.8) has the form where J(ǫ) is a constant boundary vector. Computing explicitly the matrix exponential in the above equation, we obtain many terms proportional to y aǫ , where a is an integer. Since we demand that the integrals are finite at y = 0 for negative values of ǫ, the coefficients of y aǫ in (4.3) have to vanish for a > 0. This imposes conditions on the constant vector J(ǫ), which we now have to translate to conditions on the value of the integral at the Euclidean point I E (ǫ), see for example [39]. In order to do this translation, we need to consider a path γ(x) that starts at p E and continues to a point p in which W 11 vanishes. It is advantageous to choose the parametrization of the path in such a way as to resolve the square root in √ ∆. Explicitly, the path reads where x parametrizes the path. This path reduces to p E for the beginning point x = x 0 = √ 5 and leads to the vanishing of W 11 (and also some other letters) for the end point x = x 1 = √ 13. Along the path from x 0 to x 1 , we have to go around the (spurious) singularity atx = 3 where the letters 18, 19, 27, 28 vanish. Since such a singularity can introduce a branch cut, we need to go around it by adding a small imaginary part. We can do in two ways, namely above or below the cut. Thus, in general, we obtain in principle two solutions for the value of the integrals in the vicinity of the end point p and two corresponding path parameterizations. In practice, we do not need to worry about this and can take just JHEP03(2019)042 one of the two, say the one going over the cut. We will then in general obtain complex equations for the unknown real vector I E (ǫ), but we can simply decompose them into real and imaginary parts. We illustrate the path γ in figure 3.
We now expand the boundary values of the integral as On the one hand, integrating the differential equation, we get the iterated integrals expression while on the other hand we get from (4.3) the expansion (4.6) In the above, we have to first impose on J(ǫ) the vanishing of the terms proportional to y aǫ with a > 0 in (4.3). Furthermore, the parameter y needs to be matched to the parametrization of the path as y = x 1 − x = √ 13 − x. The matching of (4.5) with (4.6) imposes conditions on the I E,n .
We now need to evaluate explicitly the iterated integrals like γ dÃ I (0) E in (4.5). This task is performed explicitly in terms of Goncharov polylogarithms G(a 1 , . . . , a k ; z) in three steps. First, we perform the iterated integrations along the path γ(x) around the beginning point x 0 using the definitions of the Goncharov polylogarithms: G(a 2 , . . . , a k ; t) . (4.7) In a second step, approaching the end point x 1 , we need to use the shuffle algebra for the Goncharov polylogarithms in order to make the terms containing log(y) = log(x 1 − x) explicit so that we can match (4.5) to (4.6). This means that we obtain an explicit expression for the integrals like γ dÃ I E in the vicinity of y = 0 that is of the type m c m log(y) m where the coefficients c m are y-independent and explicitly given in terms of the values of the Goncharov polylogarithms that are finite for y = 0. The specific value of these constants depends in principle on the path chosen to avoid the spurious singularity atx. It suffices for our purposes to choose one of the two.
We can now perform the matching (4.5) to (4.6) and obtain analytic conditions on the I (n) E vectors that contain many different Goncharov polylogarithms. Repeating the JHEP03(2019)042 same procedure that we did for the path γ for many other paths going to other spurious singularities, we obtain many constraints on the boundary conditions.
Finally, we used one more constraint, which comes from the analysis of the leading singularities. One may classify the integral basis according to parity. Then, the parity odd integrals are expressed in terms of certain F 's of eq. (2.2), and normalized by a factor proportional to √ ∆ to make them pure integrals. Since ∆ → 0 is not a physical singularity of the Feynman integrals, the odd pure integrals have to vanish at ∆ → 0. Similarly to the previous analysis we consider a path γ(x), which rationalizes √ ∆, joining the Euclidean point p E where ∆(p E ) = 3 √ 5 and a singular point where ∆ = 0, and we integrate the differential equation along this path in terms of Goncharov polylogarithms. The analytic conditions on the vector I (n) E come from vanishing of I(v i ; ǫ) at the boundary point. More specifically, we use the shuffle algebra for the Goncharov polylogarithms to extract logarithmic singularities and we demand vanishing of the coefficients in front of all powers of the logarithms. Taking the union of all the constraints discussed above we find that they are sufficient to fix the boundary conditions analytically. We note that performing the matching at weight L imposes additional conditions needed to fix the coefficients at weight L − 1. Thus, we need to go to weight 5 for some of the paths, in order to obtain enough conditions to fix all the coefficients.
In fact, we obtain an overdetermined system of equations and solving it requires using many identities for the Goncharov polylogarithms. Thus, we choose to solve the equations numerically, which leads us to the third step, namely the numerical evaluation of those Goncharov polylogarithms that are finite at the end point x 1 of the path. In order do that, we use the GiNaC implementation of [33]. While in principle we can solve the equations to arbitrary numerical precision using GiNaC, to limit computing time, we chose 100 digits precision.
The consistency conditions from matching (4.5) to (4.6) for all the possible paths γ going from p E to points at which some even letters W j vanish is enough to fix all the unknown coefficients in I (n) E , up to an overall normalization condition. The latter reflects the fact that eq. (2.8) is homogeneous in I(v i ; ǫ). To fix the normalization, it is sufficient to compute one of the trivial integrals in the d-log basis I(v i ; ǫ) analytically. Factoring out the overall divergence and common factors from dimensional regularisation, the first component of I(v i ; ǫ) is defined and expressed as follows: The result (4.8) provides the normalization fixing all the remaining coefficients. In particular, we have for the first vector Furthermore, for illustration we show explicitly the complete solution for I(v i ; ǫ) up to linear order in ǫ, Inserting the values of v i for the Euclidean point p E in the above, one obtains our analytic expression for I (1) E , which is proportional to log(3). As was already mentioned, the other boundary vectors I Having fixed the boundary values up to weight 4 one can easily find an analytic solution of the differential equation (2.8) up to the same order in ǫ-expansion. Given a point in the Euclidean region of the kinematical space one connects it with the point p E by a path and integrates the differential equation along the path according to (4.5). Choosing a piecewise linear path one can rationalize the integration kernels dÃ and reduce all integrations to Goncharov polylogarithms (4.7). For more details on this procedure see e.g. [19,40].

The symbol of the solution
Thanks to the boundary vector (4.9), we can easily derive an explicit expression for the symbol of all of the 73 integrals. We refer to [41]

JHEP03(2019)042
to [18] for specific information concerning the integrable symbols relevant for this article. It follows directly from the differential equation (2.8) and the definition (2.9), that the symbol of the integrals we have computed are given by (4.11) It should be noted that the standard ordering in symbols is the opposite to that of the Goncharov polylogarithms in (4.7). For symbols [W i 1 , . . . , W in ], derivatives act on the last entry.

Checks on the solution
Several independent checks were performed on the hexa-box integrals derived above. The full set of integrals can be compared to the purely numerical evaluation, obtained using sector-decomposition with the FIESTA [42] code. The comparison is performed in the Euclidean point p E , and yields good agreement within the available numerical precision. However, the error margins on the FIESTA results increase with increasing weight, and agreement can be established for I E only to 1% and for I E only to 2%. The hexa-box integral family contains subtopologies corresponding to planar and nonplanar four-point functions with one off-shell leg [26][27][28] and to planar five-point functions [5][6][7]. Analytical expressions for all these integrals were derived previously. Working again in the Euclidean point p E , we performed a detailed numerical comparison for all previously available integrals (63 of the 73 integrals from the hexa-box family), using the routines described in [31,32] for the four-point functions and [7] for the five-point functions, observing full agreement of the results. It is worth noting that the Euclidean five-particle kinematics translates for some of the subsector integrals into (space-like) Minkowskian four-particle kinematics [28], where the integrals nevertheless remain real.
Finally, in appendix A we perform a direct check of the symbols of some of the components of I(v i ; ǫ) by deriving their Mellin-Barnes representation, which can then be used to bootstrap their symbol using the methods explained in [18]. We obtain a perfect agreement with the expression in (4.11).

Conclusion and discussion
In this paper, we computed an analytical expression for all massless non-planar two-loop five-point integrals belonging to the hexa-box integral family. Our computation is based on the identification of a basis of integrals with constant leading singularities, which fulfil a system of differential equations in canonical form. Inspection of this system verifies a conjecture made in [18] about the function space governing these integrals. By construction, the d-log form of the differential equation system is solved trivially in terms of iterated integrals.

JHEP03(2019)042
To uniquely determine the integrals from their differential equations requires knowledge on their boundary values in one specific kinematical point. Using physical insights on the singularity structure of the integrals, we infer boundary conditions from their behaviour in spurious kinematical points where the differential equations become singular, but the integrals themselves should remain regular. These boundary conditions are combined into a boundary value for all integrals in one specific point in the Euclidean region, from where the integrals can be evaluated straightforwardly, for example in terms of Goncharov polylogarithms.
The integrals are real and single-valued in the Euclidean region. For practical applications in scattering amplitude calculations, their analytical continuation to the Minkowskian regions corresponding to all kinematical crossings is required. This can in principle be performed on the iterated integrals with an appropriate deformation of the integration contours. Aiming for an efficient numerical representation in all regions, a more systematic approach is in order, analogous to the work on the planar two-loop five-point integrals [7]. In there, the minimal basis of planar pentagon functions was identified from their required analyticity properties, expressed in entry conditions on their symbol. All these functions were written in terms of one-dimensional integrals containing simple logarithmic and polylogarithmic integrands, with boundary values determined separately in each Minkowskian region. A similar procedure should equally be feasible for the non-planar five point integrals from the hexa-box family. It will be subject of future work, aiming for an efficient numerical representation for arbitrary physical kinematics.
The master integrals considered in this paper are relevant for two-to-three scattering processes in arbitrary theories with massless particles. Any integral of the hexa-box family can be expressed in terms of the master integrals computed here. In some cases, this may involve additional integral reduction identities beyond the ones used here for deriving the differential equations. There are by now several approaches [15,16] for finding such integral reductions. On the other hand, in the case of five-particle scattering in N = 4 super Yang-Mills theory, no further integral reductions are necessary. This is due to the fact that all hexa-box integrals appearing there are directly part of our integral basis. Therefore, with the work presented here, all integrals of one of the two non-planar integral families contributing to the amplitude [22] are known. The second non-planar integral family is beyond the scope of the present paper.
Note added. After completion of this paper, important related progress was made by the analytic reconstruction of two-loop five-gluon QCD amplitudes [43,44] at leading color level in terms of master integrals, and the analytic determination of the full-color fivepoint two-loop amplitude in N = 4 SYM theory [45,46] at symbol level. Moreover, the calculation of the integrals [45,47]

A Comparison with Mellin-Barnes calculation
In this appendix, we compute the symbols of several integrals using the Mellin-Barnes technique, which provides a useful check of the results of the main text. The integrals we discuss are also of direct importance for amplitudes in N = 4 super Yang-Mills theory.
We are interested in checking the symbols (4.11). For this, we shall compute the symbols of a few integrals using the Mellin-Barnes bootstrap technique of [18]. The integrals that we shall consider are the four members of the hexa-box integral family that can be found in [22], see figure 2. In the notation of [22], these are the integral (i), the integrals (h) with numerators N The symbol of the integral (i) was computed in [18] up to and including the finite part in the ǫ-expansion by means of the Mellin-Barnes technique. Here we outline how to obtain the symbols of integrals (h) and (c) by this method as well. In order to be more specific, let us from now on concentrate on the integral of type (h) with the numerator N (h) 1 , which we dub I (h 1 ) . The steps that we perform can be done, with minimal modifications for the other (h) integral as well as for the integral (c).
We start by deriving a neat Feynman representation for I (h 1 ) , which will then allow us to obtain its Mellin-Barnes representation. This is done by using the fact that I (h 1 ) contains a box sub-integral (see figure 2) which can be rewritten as a two-fold integral of a propagator raised to the power 2 + ǫ as:
interested in the limits for which the simplified Mellin-Barnes integrals can be evaluated explicitly by means of the Cauchy theorem. Furthermore, the same limits considerably simplify the 31-letter alphabet. We specialize to those limits leading to 2dHPL and HPL alphabets. Then, by considering the computed asymptotics of the Mellin-Barnes integrals and by comparing them to the symbol ansatz, we can fix the unknown coefficients in the ansatz. In this way, we obtain the symbol of the integral I (h 1 ) up to and including the finite part. but we stress that we computed all the terms up to and including the ǫ 4 terms.
To summarize, we obtain the symbol of I (h 1 ) by first deriving a Feynman representation by getting rid of the box sub-integral, then trading that Feynman representation for a Mellin-Barnes one which is very convenient for taking suitable kinematical limits for which the integral can be evaluated explicitly such that finally one obtains constraints for an inspired symbol ansatz. The symbol (A.11) can now be compared directly to the results we have obtained in the main text. Specifically, I (h 1 ) = − I(v i ; ǫ) 56 and we obtained the symbol of the right hand side in (4.11). We find that both sides are in complete agreement.
Identical calculations for the symbols of the other integral of type (h) as well as of the integral (c) have also been done. They are given in terms of the integrals as I (h 3 ) = I(v i ; ǫ) 57 and I (c) = I(v i ; ǫ) 71 , and the symbol results agree with the computations of the main text.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.