Hexagon functions and the three-loop remainder function

We present the three-loop remainder function, which describes the scattering of six gluons in the maximally-helicity-violating configuration in planar N=4 super-Yang-Mills theory, as a function of the three dual conformal cross ratios. The result can be expressed in terms of multiple Goncharov polylogarithms. We also employ a more restricted class of"hexagon functions"which have the correct branch cuts and certain other restrictions on their symbols. We classify all the hexagon functions through transcendental weight five, using the coproduct for their Hopf algebra iteratively, which amounts to a set of first-order differential equations. The three-loop remainder function is a particular weight-six hexagon function, whose symbol was determined previously. The differential equations can be integrated numerically for generic values of the cross ratios, or analytically in certain kinematics limits, including the near-collinear and multi-Regge limits. These limits allow us to impose constraints from the operator product expansion and multi-Regge factorization directly at the function level, and thereby to fix uniquely a set of Riemann-zeta-valued constants that could not be fixed at the level of the symbol. The near-collinear limits agree precisely with recent predictions by Basso, Sever and Vieira based on integrability. The multi-Regge limits agree with the factorization formula of Fadin and Lipatov, and determine three constants entering the impact factor at this order. We plot the three-loop remainder function for various slices of the Euclidean region of positive cross ratios, and compare it to the two-loop one. For large ranges of the cross ratios, the ratio of the three-loop to the two-loop remainder function is relatively constant, and close to -7.


Contents 1 Introduction 4
2 Extra-pure functions and the symbol of R For roughly half a century we have known that many physical properties of scattering amplitudes in quantum field theories are encoded in different kinds of analytic behavior in various regions of the kinematical phase space. The idea that the amplitudes of a theory can be reconstructed (or 'bootstrapped') from basic physical principles such as unitarity, by exploiting the link to the analytic behavior, became known as the "Analytic S-Matrix program" (see e.g. ref. [1]). In the narrow resonance approximation, crossing symmetry duality led to the Veneziano formula [2] for tree-level scattering amplitudes in string theory.
In conformal field theories, there exists a different kind of bootstrap program, whereby correlation functions can be determined by imposing consistency with the operator product expansion (OPE), crossing symmetry, and unitarity [3,4]. This program was most successful in two-dimensional conformal field theories, for which conformal symmetry actually extends to an infinite-dimensional Virasoro symmetry [5]. However, the basic idea can be applied in any dimension and recent progress has been made in applying the program to conformal field theories in three and four dimensions [6,7,8].
In recent years, the scattering amplitudes of the planar N = 4 super-Yang-Mills theory have been seen to exhibit remarkable properties. In particular, the amplitudes exhibit dual conformal symmetry and a duality to light-like polygonal Wilson loops [9]. The dual description and its associated conformal symmetry mean that CFT techniques can be applied to calculating scattering amplitudes. In particular, the idea of imposing consistency with the OPE applies. However, since the dual observables are non-local Wilson loop operators, a different OPE, involving the near-collinear limit of two sides of the light-like polygon, has to be employed [10,11,12,13].
Dual conformal symmetry implies that the amplitudes involving four or five particles are fixed, because there are no invariant cross ratios that can be formed from a five-sided light-like polygon [14]. The four-and five-point amplitudes are governed by the BDS ansatz [15]. The amplitudes not determined by dual conformal symmetry begin at six points. When the external gluons are in the maximally-helicity-violating (MHV) configuration, such amplitudes can be expressed in terms of the BDS ansatz, which contains all of the infrared divergences and transforms anomalously under dual conformal invariance, and a so-called "remainder function" [16], which only depends on dual-conformally-invariant cross ratios. In the case of non-MHV amplitudes, one can define the "ratio function" [17], which depends on the cross ratios as well as dual superconformal invariants. For six external gluons, the remainder and ratio functions are described in terms of functions of three dual conformal cross ratios.
At low orders in perturbation theory, these latter functions can be expressed in terms of multiple polylogarithms. In general, multiple polylogarithms are functions of many variables that can be defined as iterated integrals over rational kernels. A particularly useful feature of such functions is that they can be classified according to their symbols [18,19,20], elements of the n-fold tensor product of the algebra of rational functions. The integer n is referred to as the transcendental weight or degree. The symbol can be defined iteratively in terms of the total derivative of the function, or alternatively, in terms of the maximally iterated coproduct by using the Hopf structure conjecturally satisfied by multiple polylogarithms [21,22]. Complicated functional identities among polylogarithms become simple algebraic relations satisfied by their symbols, making the symbol a very useful tool in the study of polylogarithmic functions. The symbol can miss terms in the function that are proportional to transcendental constants (which in the present case are all multiple zeta values), so special care must be given to account for these terms. The symbol and coproduct have been particularly useful in recent field theory applications [23,12,24,25,26]. In the case of N = 4 super-Yang-Mills theory, all amplitudes computed to date have exhibited a uniform maximal transcendentality, in which the finite terms (such as the remainder or ratio functions) always have weight n = 2L at the L loop order in perturbation theory.
Based on the simplified form of the two-loop six-point remainder function obtained in ref. [23] (which was first constructed analytically in terms of multiple polylogarithms [27]), it was conjectured [24,25] that for multi-loop six-point amplitudes, both the MHV remainder function and the next-to-MHV (NMHV) ratio function are described in terms of polylogarithmic functions whose symbols are made from an alphabet of nine letters. The nine letters are related to the nine projectively-inequivalent differences z ij of projective variables z i [23], which can also be represented in terms of momentum twistors [28]. Using this conjecture, the symbol for the three-loop six-point remainder function was obtained up to two undetermined parameters [24], which were later fixed [29] using a dual supersymmetry "anomaly" equation [29,30]. The idea of ref. [24] was to start with an ansatz for the symbol, based on the above nine-letter conjecture, and then impose various mathematical and physical consistency conditions. For example, imposing a simple integrability condition [19,20] guarantees that the ansatz is actually the symbol of some function, and demanding that the amplitude has physical branch cuts leads to a condition on the initial entries of the symbol.
Because of the duality between scattering amplitudes and Wilson loops, one can also impose conditions on the amplitude that are more naturally expressed in terms of the Wilson loop, such as those based on the OPE satisfied by its near-collinear limit. In refs. [10,11,12,13], the leading-discontinuity terms in the OPE were computed. In terms of the cross ratio variable that vanishes in the near-collinear limit, the leading-discontinuity terms correspond to just the maximum powers of logarithms of this variable (L−1 at L loops), although they can be arbitrarily power suppressed. These terms require only the one-loop anomalous dimensions of the operators corresponding to excitations of the Wilson line, or flux tube. That is, higher-loop corrections to the anomalous dimensions and to the OPE coefficients can only generate subleading logarithmic terms. While the leading-discontinuity information is sufficient to determine all terms in the symbol at two loops, more information is necessary starting at three loops [24].
Very recently, a new approach to polygonal Wilson loops has been set forth [31,32], which is fully nonperturbative and based on integrability. The Wilson loop is partitioned into a number of "pentagon transitions", which are labeled by flux tube excitation states on either side of the transition. (If one edge of the pentagon coincides with an edge of the Wilson loop, then the corresponding state is the flux tube vacuum.) The pentagon transitions obey a set of bootstrap consistency conditions. Remarkably, they can be solved in terms of factorizable S matrices for two-dimensional scattering of the flux tube excitations [31,32].
In principle, the pentagon transitions can be solved for arbitrary excitations, but it is simplest to first work out the low-lying excitations, which correspond to the leading power-suppressed terms in the near-collinear limit in the six-point case (and similar terms in multi-near-collinear limits for more than six particles). Compared with the earlier leading-discontinuity data, now all terms at a given power-suppressed order can be determined (to all loop orders), not just the leading logarithms. This information is very powerful. The first power-suppressed order in the six-point near-collinear limit is enough to fix the two terms in the ansatz for the symbol of the three-loop remainder function that could not be fixed using the leading discontinuity [31]. At four loops, the first power-suppressed order [31] and part of the second power-suppressed order [33] are sufficient to fix all terms in the symbol [34]. At these orders, the symbol becomes heavily over-constrained, providing strong cross checks on the assumptions about the letters of the symbol, as well as on the solutions to the pentagon transition bootstrap equations. In short, the application of integrability to the pentagon-transition decomposition of Wilson loops provides, through the OPE, all-loop-order boundary-value information for the problem of determining Wilson loops (or scattering amplitudes) at generic nonzero (interior) values of the cross ratios. We will use this information in the six-point case to uniquely determine the three-loop remainder function, not just at symbol level, but at function level as well.
A second limit we study is the limit of multi-Regge kinematics (MRK), which has provided another important guide to the perturbative structure of the six-point remainder function [35,36,37,38,39,40,24,41], as well as higher-point remainder functions [42,43] and NMHV amplitudes [44]. The six-point remainder function and, more generally, the hexagon functions that we define shortly have simple behavior in the multi-Regge limit. These functions depend on three dual-conformally-invariant cross ratios, but in the multi-Regge limit they collapse [45] into single-valued harmonic polylogarithms [46], which are functions of two surviving real variables, or of a complex variable and its conjugate. The multi-Regge limit factorizes [41] after taking the Fourier-Mellin transform of this complex variable. This factorization imposes strong constraints on the remainder function at high loop order [41,45,47].
Conversely, determining the multi-loop remainder function, or just its multi-Regge limit, allows the perturbative extraction of the two functions that enter the factorized form of the amplitude, the BFKL eigenvalue (in the adjoint representation) and a corresponding impact factor. This approach makes use of a map between the single-valued harmonic polylogarithms and their Fourier-Mellin transforms, which can be constructed from harmonic sums [45]. Using the three-and four-loop remainder-function symbols, the BFKL eigenvalue has been determined to next-to-next-to-leading-logarithmic accuracy (NNLLA), and the impact factor at NNLLA and N 3 LLA [45]. However, the coefficients of certain transcendental constants in these three quantities could not be fixed, due to the limitation of the symbol. Here we will use the MRK limit at three loops to fix the three undetermined constants in the NNLLA impact factor. Once the four-loop remainder function is determined, a similar analysis will fix the undetermined constants in the NNLLA BFKL eigenvalue and in the N 3 LLA impact factor.
In general, polylogarithmic functions are not sufficient to describe scattering amplitudes. For example, an elliptic integral, in which the kernel is not rational but contains a square root, enters the two-loop equal-mass sunrise graph [48], and it has been shown that a very similar type of integral enters a particular N 3 MHV 10-point scattering amplitude in planar N = 4 super-Yang-Mills theory [49]. However, it has been argued [50], based on a novel form of the planar loop integrand, that MHV and NMHV amplitudes can all be described in terms of multiple polylogarithms alone. Similar "dlog" representations have appeared in a recent twistor-space formulation [51]. Because six-particle amplitudes are either MHV (or the parity conjugate MHV) or NMHV, we expect that multiple polylogarithms and their associated symbols should suffice in this case. The nine letters that we assume for the symbol then follow naturally from the fact that the kinematics can be described in terms of dual conformally invariant combinations of six momentum twistors [28].
Having the symbol of an amplitude is not the same thing as having the function. In order to reconstruct the function one first needs a representative, well-defined function in the class of multiple polylogarithms which has the correct symbol. Before enough physical constraints are imposed, there will generally be multiple functions matching the symbol, because of the symbollevel ambiguity associated with transcendental constants multiplying well-defined functions of lower weight. Here we will develop techniques for building up the relevant class of functions for hexagon kinematics, which we call hexagon functions, whose symbols are as described above, but which are well-defined and have the proper branch cuts at the function level as well. We will argue that the hexagon functions form the basis for a perturbative solution to the MHV and NMHV six-point problem.
We will pursue two complementary routes toward the construction of hexagon functions. The first route is to express them explicitly in terms of multiple polylogarithms. This route has the advantage of being completely explicit in terms of functions with well-known mathematical properties, which can be evaluated numerically quite quickly, or expanded analytically in various regions. However, it also has the disadvantages that the representations are rather lengthy, and they are specific to particular regions of the full space of cross ratios.
The second route we pursue is to define each weight-n hexagon function iteratively in the weight, using the three first-order differential equations they satisfy. This information can also be codified by the {n − 1, 1} component of the coproduct of the function, whose elements contain weight-(n−1) hexagon functions (the source terms for the differential equations). The differential equations can be integrated numerically along specific contours in the space of cross ratios. In some cases, they can be integrated analytically, at least up to the determination of certain integration constants.
We can carry out numerical comparisons of the two approaches in regions of overlapping validity. We have also been able to determine the near-collinear and multi-Regge limits of the functions analytically using both routes. As mentioned above, these limits are how we fix all undetermined constants in the function-level ansatz, and how we extract additional predictions for both regimes.
We have performed a complete classification of hexagon functions through weight five. Although the three-loop remainder function is a hexagon function of weight six, its construction is possible given the weight-five basis. There are other potential applications of our classification, beyond the three-loop remainder function. One example is the three-loop six-point NMHV ratio function, whose components are expected [25] to be hexagon functions of weight six. Therefore, it should be possible to construct the ratio function in an identical fashion to the remainder function.
Once we have fixed all undetermined constants in the three-loop remainder function, we can study its behavior in various regions, and compare it with the two-loop function. On several lines passing through the space of cross ratios, the remainder function collapses to simple combinations of harmonic polylogarithms of a single variable. Remarkably, over vast swathes of the space of positive cross ratios, the two-and three-loop remainder functions are strikingly similar, up to an overall constant rescaling. This similarity is in spite of the fact that they have quite different analytic behavior along various edges of this region. We can also compare the perturbative remainder function with the result for strong coupling, computed using the AdS/CFT correspondence, along the line where all three cross ratios are equal [52]. We find that the two-loop, three-loop and strong-coupling results all have a remarkably similar shape when the common cross ratio is less than unity. Although we have not attempted any kind of interpolation formula from weak to strong coupling, it seems likely from the comparison that the nature of the interpolation will depend very weakly on the common cross ratio in this region.
The similarity of the weak and strong coupling limits of remainder functions has been noticed before. For the eight-point case in two-dimensional (AdS 3 ) kinematics, there are two real kinematical parameters. Refs. [53,54] found impressive numerical agreement, to within 3% or better, between rescaled versions of the two-loop [55] and strong-coupling [56] octagon remainder function, as a function of both parameters. The octagon, decagon, and general 2n-point remainder functions were evaluated at strong coupling, analytically in an expansion around the regular polygon limit, and the rescaled functions were found to be very similar to the two-loop result [57]. The six-point case we consider in this paper has also been studied at strong coupling using a Z 4 symmetric integrable model [58] and using the homogeneous sine-Gordon model and conformal perturbation theory [59]. In the latter work, the strong-coupling result was compared with the two-loop one along a one-dimensional curve in the space of the three cross ratios, corresponding to the trajectory of an integrable renormalization group flow. Again, good numerical agreement was found between the two rescaled functions.
The remainder of this paper is organized as follows. In section 2 we recall some properties of pure functions (iterated integrals) and their symbols, as well as a representation of the twoloop remainder function (and its symbol) in terms of an "extra pure" function and its cyclic images. We use this representation as motivation for an analogous decomposition of the threeloop symbol. In section 3 we describe the first route to constructing hexagon functions, via multiple polylogarithms. In section 4 we describe the second route to constructing the same set of functions, via the differential equations they satisfy. In section 5 we discuss how to extract the near-collinear limits, and give results for some of the basis functions and for the remainder function in this limit. In section 6 we carry out the analogous discussion for the Minkowski multi-Regge limit. In section 7 we give the final result for the three-loop remainder function, in terms of a specific integral, as well as defining it through the {5, 1} components of its coproduct. We also present the specialization of the remainder function onto various lines in the threedimensional space of cross ratios; along these lines its form simplifies dramatically. Finally, we plot the function on several lines and two-dimensional slices. We compare it numerically to the two-loop function in some of these regions, and to the strong-coupling result evaluated for equal cross ratios. In section 8 we present our conclusions and outline avenues for future research. We include three appendices. Appendix A provides some background material on multiple polylogarithms. Appendix B gives the complete set of independent hexagon functions through weight five in terms of the {n − 1, 1} components of their coproducts, and in appendix C we provide the same description of the extra pure weight six function R ep entering the remainder function.
In attached, computer-readable files we give the basis of hexagon functions through weight five, as well as the three-loop remainder function, expressed in terms of multiple polylogarithms in two different kinematic regions. We also provide the near-collinear and multi-Regge limits of these functions.
2 Extra-pure functions and the symbol of R (3)  6 In this section, we describe the symbol of the three-loop remainder function as obtained in ref. [24], which is the starting point for our reconstruction of the full function. Motivated by an alternate representation [25] of the two-loop remainder function, we will rearrange the threeloop symbol. In the new representation, part of the answer will involve products of lower-weight (hence simpler) functions, and the rest of the answer will be expressible as the sum of an extrapure function, called R ep , plus its two images under cyclic permutations of the cross ratios. An extra-pure function of m variables, by definition, has a symbol with only m different final entries. For the case of hexagon kinematics, where there are three cross ratios, the symbol of an extrapure function has only three final entries, instead of the potential nine. Related to this, the three derivatives of the full function can be written in a particularly simple form, which helps somewhat in its construction.
All the functions we consider in this paper will be pure functions. The definition of a pure function f (n) of transcendental weight (or degree) n is that its first derivative obeys, where φ r are rational functions and the sum over r is finite. The only weight-zero functions are assumed to be rational constants. The f (n−1) r and φ r are not all independent of each other because the integrability condition d 2 f (n) = 0 imposes relations among them, Functions defined by the above conditions are iterated integrals of polylogarithmic type. Such functions have a symbol, defined recursively as an element of the n-fold tensor product of the algebra of rational functions, following eq. (2.1), In the case of the six-particle amplitudes of planar N = 4 super Yang-Mills theory, we are interested in pure functions depending on the three dual conformally invariant cross ratios, (2.4) The six particle momenta k µ i are differences of the dual coordinates x µ i : x µ i − x µ i+1 = k µ i , with indices taken mod 6.
Having specified the class of functions we are interested in, we impose further [24,25] that the entries of the symbol are drawn from the following set of nine letters, and relations obtained by cyclically rotating the six points. The variables y u , y v and y w can be expressed locally in terms of the cross ratios, where Note that under the cyclic permutation z i → z i+1 we have u → v → w → u, while the y i variables transform as y u → 1/y v → y w → 1/y u . A three-fold cyclic rotation amounts to a space-time parity transformation, under which the parity-even cross ratios are invariant, while the parity-odd y variables invert. Consistent with the inversion of the y variables under parity, and with eq. (2.7), the quantity √ ∆ must flip sign under parity, so we have altogether, The transformation of √ ∆ can also be seen from its representation in terms of the z ij variables, upon letting z i → z i+3 . It will prove very useful to classify hexagon functions by their parity. The remainder function is a parity-even function, but some of its derivatives (or more precisely coproduct components) are parity odd, so we need to understand both the even and odd sectors. Since the y variables invert under parity, y u → 1/y u , etc., it is often better to think of the y variables as fundamental and the cross ratios as parity-even functions of them. The cross ratios can be expressed in terms of the y variables without any square roots, where we have picked a particular branch of √ ∆. Following the strategy of ref. [24], we construct all integrable symbols of the required weight, using the letters (2.5), subject to certain additional physical constraints. In the case of the six-point MHV remainder function at L loops, we require the symbol to be that of a weight-2L parity-even function with full S 3 permutation symmetry among the cross ratios. The initial entries in the symbol can only be the cross ratios themselves, in order to have physical branch cuts [12]: first entry ∈ {u, v, w} . (2.12) In addition we require that the final entries of the symbol are taken from the following restricted set of six letters [60,24]: Next one can apply constraints from the collinear OPE of Wilson loops. The leading-discontinuity constraints [10,11,12] can be expressed in terms of differential operators with a simple action on the symbol [24]. At two loops, the leading (single) discontinuity is the only discontinuity, and it is sufficient to determine the full remainder function R 6 (u, v, w) [11]. At three loops, the constraint on the leading (double) discontinuity leaves two free parameters in the symbol, α 1 and α 2 [24]. These parameters were determined in refs. [29,31], but we will leave them arbitrary here to see what other information can fix them.
The two-loop remainder function R (2) 6 can be expressed simply in terms of classical polylogarithms [23]. However, here we wish to recall the form found in ref. [25] in terms of the infrared-finite double pentagon integral Ω (2) , which was introduced in ref. [61] and studied further in refs. [62,25]: 6,rat (u, v, w) . (2.14) The function R 6,rat can be expressed in terms of single-variable classical polylogarithms, R 6,rat = − 1 2 We see that R 6,rat decomposes into a product of simpler, lower-weight functions Li 2 (1 − 1/u i ), plus the cyclic images of the function r(u), whose symbol can be written as, The symbol of Ω (2) can be deduced [25] from the differential equations it satisfies [62,63]. There are only three distinct final entries of the symbol of Ω (2) (u, v, w), namely Note that three is the minimum possible number of distinct final entries we could hope for, since Ω (2) is genuinely dependent on all three variables. As mentioned above, we define extrapure functions, such as Ω (2) , to be those functions for which the number of final entries in the symbol equals the number of variables on which they depend. Another way to state the property (which also extends it from a property of symbols to a property of functions) is that p-variable pure functions f of weight n are extra-pure if there exist p independent commuting first-order differential operators O i , such that O i f are themselves all pure of weight (n − 1). More explicitly, the symbol of Ω (2) can be written as [25], The functions Q φ and Q r will be defined below. The functionΦ 6 is the weight-three, parity-odd one-loop six-dimensional hexagon function [63,64], whose symbol is given by [63], where Ω (1) is a finite, four-dimensional one-loop hexagon integral [61,62], Although we have written eq. (2.21) as an equation for the symbol ofΦ 6 , secretly it contains more information, because we have written the symbol of a full function, Ω (1) (u, v, w), in the first two slots. Later we will codify this extra information as corresponding to the {2, 1} component of the coproduct ofΦ 6 . Another way of saying it is that all three derivatives of the functionΦ 6 , with respect to the logarithms of the y variables, are given by −Ω (1) (u, v, w) or its permutations, including the ζ 2 term in eq. (2.22). Any other derivative can be obtained by the chain rule. For example, to get the derivative with respect to u, we just need, which leads to the differential equation found in ref. [63], (2.24) HenceΦ 6 can be fully specified, up to a possible integration constant, by promoting the first two slots of its symbol to a function in an appropriate way. In fact, the ambiguity of adding a constant of integration is actually fixed in this case, by imposing the property that the functioñ Φ 6 is parity odd.
Note that for the solution to the differential equation (2.24) and its cyclic images to have physical branch cuts, the correct coefficients of the ζ 2 terms in eq. (2.22) are crucial. Changing the coefficients of these terms in any of the cyclic images of Ω (1) would correspond to adding a logarithm of the y variables toΦ 6 , which would have branch cuts in unphysical regions.
The other weight-three symbols in eq. (2.19) can similarly be promoted to full functions. To do this we employ the harmonic polylogarithms (HPLs) in one variable [65], H w (u). In our case, the weight vector w contains only 0's and 1's. If the weight vector is a string of n 0's, w = 0 n , then we have H 0n (u) = 1 n! log n u. The remaining functions are defined recursively by Such functions have symbols with only two letters, {u, 1 − u}. We would like the point u = 1 to be a regular point for the HPLs. This can be enforced by choosing the argument to be 1 − u, and restricting to weight vectors whose last entry is 1. The symbol and HPL definitions have a reversed ordering, so to find an HPL with argument 1−u corresponding to a symbol in {u, 1−u}, one reverses the string, replaces u → 1 and 1 − u → 0, and multiplies by (−1) for each 1 in the weight vector. We also use a compressed notation where (k − 1) 0's followed by a 1 is replaced by k in the weight vector, and the argument (1 − u) is replaced by the superscript u. For example, ignoring ζ-value ambiguities we have, The combination occurs frequently, because it is the lowest-weight extra-pure function, with symbol u ⊗ u/(1 − u).
In terms of HPLs, the functions corresponding to the weight-three, parity-even symbols appearing in eq. (2.19) are given by, (2.28) Here we have added some ζ 2 terms with respect to ref. [25], in order to match the {3, 1} component of the coproduct of Ω (2) that we determine later. Note that the simple form of the symbol of R 6,rat in eq. (2.15) means that it can be absorbed into the three cyclic images of Ω (2) (u, v, w) without ruining the extra-purity of the latter functions. Hence R (2) 6 is the cyclic sum of an extra-pure function. With the decomposition (2.14) in mind, we searched for an analogous decomposition of the symbol of the three-loop remainder function [24] into extra-pure components. In other words, we looked for a representation of S(R (3) 6 ) in terms a function whose symbol has the same final entries (2.18) as Ω (2) (u, v, w), plus its cyclic rotations. After removing some products of lowerweight functions we find that this is indeed possible. Specifically, we find that, (2.29) Here P 6 is the piece constructed from products of lower-weight functions, The function R ep is very analogous to Ω (2) in that it has the same (u ↔ v) symmetry, and its symbol has the same final entries, (2.31) In the following we will describe a systematic construction of the function R ep and hence the three-loop remainder function. As in the case just described forΦ 6 , and implicitly for Ω (2) , the construction will involve promoting the quantities S(R u ep ) and S(R yu ep ) to full functions, with the aid of the coproduct formalism. In fact, we will perform a complete classification of all well-defined functions corresponding to symbols with nine letters and obeying the first entry condition (2.12) (but not the final entry condition (2.13)), iteratively in the weight through weight five. Knowing all such pure functions at weight five will then enable us to promote the weight-five quantities S(R u ep ) and S(R yu ep ) to well-defined functions, subject to ζ-valued ambiguities that we will fix using physical criteria.

Hexagon functions as multiple polylogarithms
The task of the next two sections is to build up an understanding of the space of hexagon functions, using two complementary routes. In this section, we follow the route of expressing the hexagon functions explicitly in terms of multiple polylogarithms. In the next section, we will take a slightly more abstract route of defining the functions solely through the differential equations they satisfy, which leads to relatively compact integral representations for them.

Symbols
Our first task is to classify all integrable symbols at weight n with entries drawn from the set S u in eq. (2.5) that also satisfy the first entry condition (2.12). We do not impose the final entry condition (2.13) because we need to construct quantities at intermediate weight, from which the final results will be obtained by further integration; their final entries correspond to intermediate entries of R ep .
The integrability of a symbol may be imposed iteratively, first as a condition on the first n − 1 slots, and then as a separate condition on the {n − 1, n} pair of slots, as in eq. (2.2). Therefore, if B n−1 is the basis of integrable symbols at weight n − 1, then a minimal ansatz for the basis at weight n takes the form, and B n can be obtained simply by enforcing integrability in the last two slots. This method for recycling lower-weight information will also guide us toward an iterative construction of full functions, which we perform in the remainder of this section. Integrability and the first entry condition together require the second entry to be free of the y i . Hence the maximum number of y entries that can appear in a term in the symbol is n − 2. In fact, the maximum number of y's that appear in any term in the symbol defines a natural grading for the space of functions. In table 1, we use this grading to tabulate the number of irreducible functions (i.e. those functions that cannot be written as products of lower-weight functions) through weight six. The majority of the functions at low weight contain no y entries.
The y entries couple together u, v, w. In their absence, the symbols with letters {u, v, w, 1 − u, 1 − v, 1 − w} can be factorized, so that the irreducible ones just have the letters {u, 1 − u}, plus cyclic permutations of them. The corresponding functions are the ordinary HPLs in one variable [65] introduced in the previous section, H u w , with weight vectors w consisting only of 0's and 1's. These functions are not all independent, owing to the existence of shuffle identities [65]. On the other hand, we may exploit Radford's theorem [66] to solve these identities in terms of a Lyndon basis, where H u lw ≡ H lw (1 − u), and Lyndon(0, 1) is the set of Lyndon words in the letters 0 and 1. The Lyndon words are those words w such that for every decomposition into two words w = {u, v}, the left word u is smaller 1 than the right word v, i.e. u < v. Notice that we exclude the case Weight y 0 y 1 y 2 y 3 y 4 18 4 13 6 -6 27 4 27 29 18 Table 1: The dimension of the irreducible basis of hexagon functions, graded by the maximum number of y entries in their symbols.
l w = 0 because it corresponds to ln(1 − u), which has an unphysical branch cut. Further cuts of this type occur whenever l w has a trailing zero, but such words are excluded from the Lyndon basis by construction. The Lyndon basis of HPLs with proper branch cuts through weight six can be written explicitly as, In order to parametrize the full space of functions whose symbols can be written in terms of the elements in the set S u , it is useful to reexpress those elements in terms of three independent variables. The cross ratios themselves are not a convenient choice of variables because rewriting the y i in terms of the u i produces explicit square roots. A better choice is to consider the y i as independent variables, in terms of which the u i are given by eq. (2.11). In this representation, the symbol has letters drawn from the ten-element set, We appear to have taken a step backward since there is an extra letter in S y relative to S u . Indeed, writing the symbol of a typical function in this way greatly increases the length of its expression. Also, the first entry condition becomes more complicated in the y variables. On the other hand, S y contains purely rational functions of the y i , and as such it is easy to construct the space of functions that give rise to symbol entries of this type. We will discuss these functions in the next subsection.
which they appear in the argument of "Lyndon(0, 1)", i.e. 0 < 1. Later we will encounter words with more letters for which this specification is less trivial.

Multiple polylogarithms
Multiple polylogarithms are a general class of multi-variable iterated integrals, of which logarithms, polylogarithms, harmonic polylogarithms, and various other iterated integrals are special cases. They are defined recursively by G(z) = 1, and, G(a 1 , . . . , a n ; z) = z 0 dt t − a 1 G(a 2 , . . . , a n ; t) , G(0, . . . , 0 Many of their properties are reviewed in appendix A, including an expression for their symbol, which is also defined recursively [67], S(G(a n−1 , . . . , a 1 ; a n )) = n−1 i=1 S(G(a n−1 , . . . ,â i , . . . , a 1 ; a n )) ⊗ (a i − a i+1 ) −S(G(a n−1 , . . . ,â i , . . . , a 1 ; a n )) ⊗ (a i − a i−1 ) , where a 0 = 0 and the hat on a i on the right-hand side indicates that this index should be omitted. Using eq. (3.6), it is straightforward to write down a set of multiple polylogarithms whose symbol entries span S y , (3.7) The set G also emerges naturally from a simple procedure by which symbols are directly promoted to polylogarithmic functions. For each letter φ i (y u , y v , y w ) ∈ S y we write ω i = d log φ i (t u , t v , t w ). Then following refs. [19,68], which are in turn based on ref. [18], we use the integration map, The integration is performed iteratively along the contour γ which we choose to take from the origin t i = 0 to the point t i = y i . The precise choice of path is irrelevant, provided the symbol we start from is integrable [18,19]. So we may choose to take a path which goes sequentially along the t u , t v , t w directions. Near the axes we may find some divergent integrations of the form y 0 dt/t•. . .•dt/t. We regularize these divergences in the same way as in the one-dimensional HPL case (see the text before eq. (2.25)) by replacing them with 1 n! log n y. In this way we immediately obtain an expression in terms of the functions in G, with the three subsets corresponding to the three segments of the contour.
The set G is larger than what is required to construct the basis of hexagon functions. One reason for this is that G generates unwanted symbol entries outside of the set S u , such as the differences y i − y j , as is easy to see from eq. (3.6); the cancellation of such terms is an additional constraint that any valid hexagon function must satisfy. Another reason is that multiple polylogarithms satisfy many identities, such as the shuffle and stuffle identities (see Appendix A or refs. [67,69] for a review). While there are no relevant stuffle relations among the functions in G, there are many relations resulting from shuffle identities. Just as for the single-variable case of HPLs, these shuffle relations may be resolved by constructing a Lyndon basis, G L I ⊂ G, (3.9) A multiple polylogarithm G(w 1 , . . . , w n ; z) admits a convergent series expansion if |z|≤ |w i | for all nonzero w i , and it is manifestly real-valued if the nonzero w i and z are real and positive. Therefore, the set G L I is ideally suited for describing configurations for which 0 < y i < 1. In terms of the original cross ratios, this region is characterized by, We will construct the space of hexagon functions in Region I as a subspace of G L I with good branch-cut properties.
What about other regions? As we will discuss in the next subsection, multiple polylogarithms in the y variables are poorly suited to regions where ∆ < 0; in these regions the y i are complex. For such cases, we turn to certain integral representations that we will describe in section 4. In this section, we restrict ourselves to the subspace of the unit cube for which ∆ > 0. As shown in fig. 1, there are four disconnected regions with ∆ > 0, which we refer to as Regions I, II, III, and IV. They are the regions that extend respectively from the four points (0, 0, 0), (1, 1, 0), (0, 1, 1), and (1, 0, 1) to the intersection with the ∆ = 0 surface. Three of the regions (II, III and IV) are related to one another by permutations of the u i , so it suffices to consider only one of them, Region II : (3.11) In Region II, the set G L I includes functions G(w 1 , . . . , w n ; z) for which |w i |< |z| for some i. As mentioned above, such functions require an analytic continuation and are not manifestly realvalued. On the other hand, it is straightforward to design an alternative basis set that does not suffer from this issue, (3.12) Like G L I , G L II also generates symbols with the desired entries. It is therefore a good starting point for constructing a basis of hexagon functions in Region II.

The coproduct bootstrap
The space of multiple polylogarithms enjoys various nice properties, many of which are reviewed in appendix A. For example, it can be endowed with the additional structure necessary to promote it to a Hopf algebra. For the current discussion, we make use of one element of this structure, namely the coproduct. The coproduct on multiple polylogarithms has been used in a variety of contexts [69,26,70,71,72,73]. It serves as a powerful tool to help lift symbols to full functions and to construct functions or identities iteratively in the weight.
Let A denote the Hopf algebra of multiple polylogarithms and A n the weight-n subspace, so that, (3.13) Then, for G n ∈ A n , the coproduct decomposes as, where ∆ p,q ∈ A p ⊗ A q . It is therefore sensible to discuss an individual {p, q} component of the coproduct, ∆ p,q . In fact, we will only need two cases, {p, q} = {n − 1, 1} and {p, q} = {1, n − 1}, though the other components carry additional information that may be useful in other contexts. A simple (albeit roundabout) procedure to extract the coproduct of a generic multiple polylogarithm, G, is reviewed in appendix A. One first rewrites G in the notation of a slightly more general function, usually denoted by I in the mathematical literature. Then one applies the main coproduct formula, eq. (A.15), and finally converts back into the G notation.
Let us discuss how the coproduct can be used to construct identities between multiple polylogarithms iteratively. Suppose we know all relevant identities up to weight n − 1, and we would like to establish the validity of some potential weight-n identity, which can always be written in the form, for some combination of weight-n functions, A n . If this identity holds, then we may further conclude that each component of the coproduct of A n should vanish. In particular, Since this is an equation involving functions of weight less than or equal to n−1, we may check it explicitly. Equation (3.16) does not imply eq. (3.15), because ∆ n−1,1 has a nontrivial kernel. For our purposes, the only relevant elements of the kernel are multiple zeta values, zeta values, iπ, and their products. Through weight six, the elements of the kernel are the transcendental constants, At weight two, for example, we may use this information to write, for some undetermined rational number c, which we can fix numerically or by looking at special limits. Consider the following example for some real positive x ≤ 1, Using eq. (3.19) and simple identities among logarithms, it is easy to check that so we conclude that A 2 = c ζ 2 . Specializing to x = 1, we find c = 2 and therefore A 2 = 2 ζ 2 . Indeed, this confirms the standard inversion relation for dilogarithms. The above procedure may be applied systematically to generate all identities within a given ring of multiple polylogarithms and multiple zeta values. Denote this ring by C and its weight-n subspace by C n . Assume that we have found all identities through weight n − 1. To find the identities at weight n, we simply look for all solutions to the equation, where G i ∈ C n and the c i are rational numbers. Because we know all identities through weight n − 1, we can write each ∆ n−1,1 (G i ) as a combination of linearly-independent functions of weight n−1. The problem is then reduced to one of linear algebra. The nullspace encodes the set of new identities, modulo elements of the kernel K. The latter transcendental constants can be fixed numerically, or perhaps analytically with the aid of an integer-relation algorithm like PSLQ [74]. For the appropriate definition of C, the above procedure can generate a variety of interesting relations. For example, we can choose C = G L I or C = G L II and confirm that there are no remaining identities within these sets.
We may also use this method to express all harmonic polylogarithms with argument u i or 1 − u i in terms of multiple polylogarithms in the set G L I or the set G L II . The only trick is to rewrite the HPLs as multiple polylogarithms. For example, using the uncompressed notation for the HPLs, H a 1 ,...,an (u) = (−1) w 1 G(a 1 , . . . , a n ; u) = (−1) w 1 G a 1 , . . . , a n ; where w 1 is the number of a i equal to one. With this understanding, we can simply take, and then proceed as above to generate all identities within this expanded ring. In all cases, the starting point for the iterative procedure for generating identities is the set of identities at weight one, i.e. the set of identities among logarithms. All identities among logarithms are of course known, but in some cases they become rather cumbersome, and one must take care to properly track various terms that depend on the ordering of the y i . For example, consider the following identity, which is valid for all complex y i , where Arg denotes the principal value of the complex argument. In principle, this identity can be used to seed the iterative procedure for constructing higher-weight identities, which would also be valid for all complex y i . Unfortunately, the bookkeeping quickly becomes unwieldy and it is not feasible to track the proliferation of Arg's for high weight.
To avoid this issue, we will choose to focus on Regions I and II, defined by eqs. (3.10) and (3.11). In both regions, ∆ > 0, so the y variables are real, and the Arg's take on specific values. In Region I, for example, we may write, In the last line, we have rewritten the logarithms in terms of multiple polylogarithms in the set G L I , which, as we argued in the previous subsection, is the appropriate basis for this region. In Region II, the expression for ln u looks a bit different, In this case, we have rewritten the logarithms as multiple polylogarithms belonging to the set G L II .
We now show how to use these relations and the coproduct to deduce relations at weight two. In particular, we will derive an expression for H u 2 = H 2 (1−u) in terms of multiple polylogarithms in the basis G L I in Region I. A similar result holds in Region II. First, we need one more weight-one identity, Next, we take the {1, 1} component of the coproduct, and substitute eqs. (3.25) and (3.27), Finally, we ask which combination of weight-two functions in G L I has the {1, 1} component of its coproduct given by eq. (3.29). There is a unique answer, modulo elements in K, (3.30) We have written a specific value for the coefficient of ζ 2 , though at this stage it is completely arbitrary since ∆ 1,1 (ζ 2 ) = 0. To verify that we have chosen the correct value, we specialize to the surface y u = 1, on which u = 1 and H u 2 = 0. It is straightforward to check that the right-hand side of eq. (3.30) does indeed vanish in this limit.
An alternative way to translate expressions made from HPLs of arguments u i into expressions in terms of the y variables is as follows. Any expression in terms of HPLs of arguments u, v, w may be thought of as the result of applying the integration map to words made from the letters u i and 1 − u i only. For example, where, to verify the final equality straightforwardly, we may choose the contour γ to run from s i = 0 to (s 1 = u, s 2 = v, s 3 = w) sequentially along the s 1 , s 2 , s 3 axes. In the above simple example the second and third parts of the contour are irrelevant since the form to be integrated only depends on s 1 anyway. Then we can change variables from u, v, w to y u , y v , y w by defining and similarly for s 2 , s 3 . Since the result obtained depends only on the end points of the contour, and not the precise path taken, we may instead choose the contour as the one which goes from the origin t i = 0 to the point t 1 = y u , t 2 = y v , t 3 = y w sequentially along the t 1 , t 2 , t 3 axes, as in the discussion around eq. (3.8). Then expression (3.32) yields an expression equivalent to eq. (3.30). Continuing this procedure on to higher weights is straightforward, although the expressions become increasingly complicated. For example, the expression for H w 4,2 has 9439 terms. It is clear that G L I is not an efficient basis, at least for representing harmonic polylogarithms with argument u i . Despite this inefficiency, G L I and G L II have the virtue of spanning the space of hexagon functions, although they still contain many more functions than desired. In the next subsection, we describe how we can iteratively impose constraints in order to construct a basis for just the hexagon functions.

Constructing the hexagon functions
Unitarity requires the branch cuts of physical quantities to appear in physical channels. For dual conformally-invariant functions corresponding to the scattering of massless particles, the only permissible branching points are when a cross ratio vanishes or approaches infinity. The location of branch points in an iterated integral is controlled by the first entry of the symbol; hence the first entry should be one of the cross ratios, as discussed previously. However, it is not necessary to restrict our attention to the symbol: it was argued in ref. [26] that the condition of only having physical branch points can be promoted to the coproduct. Then the monodromy operator M z k =z 0 (which gives the phase in analytically continuing the variable z k around the point z 0 ) acts on the first component of the coproduct ∆ (see appendix A.2), We conclude that if F n is a weight-n function with the proper branch-cut locations, and then F r n−1 must also be a weight-(n − 1) function with the proper branch-cut locations, for every r (which labels the possible letters in the symbol). Working in the other direction, suppose we know the basis of hexagon functions through weight n − 1. We may then use eq. (3.35) and the coproduct bootstrap of section 3.3 to build the basis at weight n.
There are a few subtleties that must be taken into account before applying this method directly. To begin with, the condition that all the F r n−1 belong to the basis of hexagon functions guarantees that they have symbol entries drawn from S u . However, it does not guarantee that F n has this property since the φ r are drawn from the set S y , which is larger than S u . This issue is easily remedied by simply disregarding those functions whose symbols have final entries outside of the set S u .
In pushing to higher weights, it becomes necessary to pursue a more efficient construction. For this purpose, it is useful to decompose the space of hexagon functions, which we denote by H, into its parity-even and parity-odd components, The coproduct can be taken separately on each component, where L + 1 and L − 1 are the parity-even and parity-odd functions of weight one, To construct H ± n , we simply write down the most general ansatz for both the left-hand side and the right-hand side of eq. (3.37) and solve the linear system. The ansatz for H ± n will be constructed from the either G L I or G L II , supplemented by multiple zeta values, while a parametrization of the right-hand side is known by assumption. For high weights, the linear system becomes prohibitively large, which is one reason why it is useful to construct the even and odd sectors separately, since it effectively halves the computational burden. We note that not every element on the right hand side of eq. (3.37) is actually in the image of ∆ n−1,1 . For such cases, we will simply find no solution to the linear equations. Finally, this parametrization of the {n − 1, 1} component of the coproduct guarantees that the symbol of any function in H n will have symbol entries drawn from S u .
Unfortunately, the procedure we just have outlined does not actually guarantee proper branch cuts in all cases. The obstruction is related to the presence of weight-(n − 1) multiple zeta values in the space H + n−1 . Such terms may become problematic when used as in eq. (3.37) to build the weight-n space, because they get multiplied by logarithms, which may contribute improper branch cuts. For example, but the function ζ n−1 ln(1 − u) has a spurious branch point at u = 1. Naively, one might think such terms must be excluded from our ansatz, but this turns out to be incorrect. In some cases, they are needed to cancel off the bad behavior of other, more complicated functions. We can exhibit this bad behavior in a simple one-variable function, It is easy to write down a weight-three function f 3 (u) that satisfies, Indeed, one may easily check that does the job. The problem is that f 3 (u) ∈ H + 3 because it has a logarithmic branch cut starting at u = 1. In fact, the presence of this cut is indicated by a simple pole at u = 1 in its first derivative, The residue of the pole is just f 2 (1) and can be read directly from eq. (3.40) without ever writing down f 3 (u). This suggests that the problem can be remedied by subtracting ζ 2 from f 2 (u). Indeed, forf for which, More generally, any function whose first derivative yields a simple pole has a logarithmic branch cut starting at the location of that pole. Therefore, the only allowed poles in the u iderivative are at u i = 0. In particular, the absence of poles at u i = 1 provides additional constraints on the space H ± n . These constraints were particularly simple to impose in the above single-variable example, because the residue of the pole at u = 1 could be directly read off from a single term in the coproduct, namely the one with ln(1 − u) in the last slot. In the full multiple-variable case, the situation is slightly more complicated. The coproduct of any hexagon function will generically have nine terms, where F is a function of weight n and the nine functions and completely specify the {n − 1, 1} component of the coproduct. The derivative with respect to u can be evaluated using eqs. (2.23) and (3.47) and the chain rule, Clearly, a pole at u = 1 can arise from F 1−u , F yv or F yw , or it can cancel between these terms. The condition that eq. (3.48) has no pole at u = 1 is a strong one, because it must hold for any values of v and w. In fact, this condition mainly provides consistency checks, because a much weaker set of constraints turns out to be sufficient to fix all undetermined constants in our ansatz.
It is useful to consider the constraints in the even and odd subspaces separately. Referring to eq. (2.9), parity sends √ ∆ → − √ ∆, and, therefore, any parity-odd function must vanish when ∆ = 0. Furthermore, recalling eq. (2.11), we see that any odd function must vanish when y i → 1 or when y u y v y w → 1. It turns out that these conditions are sufficient to fix all undetermined constants in the odd sector. One may then verify that there are no spurious poles in the u i -derivatives.
There are no such vanishing conditions in the even sector, and to fix all undetermined constants we need to derive specific constraints from eq. (3.48). We found it convenient to enforce the constraint for the particular values of v and w such that the u → 1 limit coincides with the limit of Euclidean multi-Regge kinematics (EMRK). In this limit, v and w vanish at the same rate that u approaches 1, where x and y are fixed. In the y variables, the EMRK limit takes y u → 1, while y v and y w are held fixed, and can be related to x and y by, This limit can also be called the (Euclidean) soft limit, in which one particle gets soft. The final point, (u, v, w) = (1, 0, 0), also lies at the intersection of two lines representing different collinear limits: In the case at hand, F is an even function and so the coproduct components F y i are odd functions of weight n − 1, and as such have already been constrained to vanish when y i → 1.
(Although the coefficients of F yv and F yw in eq. (3.48) contain factors of 1/ √ ∆, which diverge in the limit y u → 1, the numerator factors 1 − u ∓ (v − w) can be seen from eq. (3.50) to vanish in this limit, cancelling the 1/ √ ∆ divergence.) Therefore, the constraint that eq. (3.48) have no pole at u = 1 simplifies considerably: Of course, two additional constraints can be obtained by taking cyclic images. These narrower constraints turn out to be sufficient to completely fix all free coefficients in our ansatz in the even sector. Finally, we are in a position to construct the functions of the hexagon basis. At weight one, the basis simply consists of the three logarithms, ln u i . Before proceeding to weight two, we must rewrite these functions in terms of multiple polylogarithms. This necessitates a choice between Regions I and II, or between the bases G L I and G L II . We construct the basis for both cases, but for definiteness let us work in Region I.
Our ansatz for ∆ 1,1 (H + 2 ) consists of the 18 tensor products, which we rewrite in terms of multiple polylogarithms in G L I . Explicit linear algebra shows that only a nine-dimensional subspace of these tensor products can be written as ∆ 1,1 (G 2 ) for G 2 ∈ G L I . Six of these weight-two functions can be written as products of logarithms. The other three may be identified with H 2 (1 − u i ) by using the methods of section 3.3. (See e.g. eq. (3.30).) Our ansatz for ∆ 1,1 (H − 2 ) consists of the nine tensor products, which we again rewrite in terms of multiple polylogarithms in G L I . In this case, it turns out that there is no linear combination of these tensor products that can be written as ∆ 1,1 (G 2 ) for G 2 ∈ G L I . This confirms the analysis at symbol level as summarized in table 1, which shows three parity-even irreducible functions of weight two (which are identified as HPLs), and no parity-odd functions.
A similar situation unfolds in the parity-even sector at weight three, namely that the space is spanned by HPLs of a single variable. However, the parity-odd sector reveals a new function. To find it, we write an ansatz for ∆ 2,1 (H − 3 ) consisting of the 39 objects, , and then look for a linear combination that can be written as ∆ 2,1 (G 3 ) for G 3 ∈ G L I . After imposing the constraints that the function vanish when y i → 1 and when y u y v y w → 1, there is a unique solution, The normalization can be fixed by comparing to the differential equation forΦ 6 , eq. (2.24). This solution is totally symmetric under the S 3 permutation group of the three cross ratios {u, v, w}, or equivalently of the three variables {y u , y v , y w }. However, owing to our choice of basis G L I , this symmetry is broken in the representation (3.56).
In principle, this procedure may be continued and used to construct a basis for the space H n for any value of n. In practice, it becomes computationally challenging to proceed beyond moderate weight, say n = 5. The three-loop remainder function is a weight-six function, but, Table 2: Irreducible basis of hexagon functions, graded by the maximum number of y entries in the symbol. The indicated multiplicities specify the number of independent functions obtained by applying the S 3 permutations of the cross ratios.
as we will see shortly, to find its full functional form we do not need to know anything about the other weight-six functions. On the other hand, we do need a complete basis for all functions of weight five or less. We have constructed all such functions using the methods just described.
Referring to table 1, there are 69 functions with weight less than or equal to five. However, any function with no y's in its symbol can be written in terms of ordinary HPLs, so there are only 30 genuinely new functions. The expressions for these functions in terms of multiple polylogarithms are quite lengthy, so we present them in computer-readable format in the attached files. The 30 new functions can be obtained from the permutations of 11 basic functions which we callΦ 6 , F 1 , Ω (2) , G, H 1 , J 1 , K 1 , M 1 , N, O, and Q ep . Two of these functions,Φ 6 and Ω (2) , have appeared in other contexts, as mentioned in section 2. Also, a linear combination of F 1 and its cyclic image can be identified with the odd part of the two-loop ratio function, denoted bỹ V [25]. (The precise relation is given in eq. (B.20).) We believe that the remaining functions are new. In table 2, we organize these functions by their weight and y-grading. We also indicate how many independent functions are generated by permuting the cross ratios. For example, Φ 6 is totally symmetric, so it generates a unique entry, while F 1 and Ω (2) are symmetric under exchange of two variables, so they sweep out a triplet of independent functions under cyclic permutations. The function Q ep has no symmetries, so under S 3 permutations it sweeps out six independent functions. The same would be true of M 1 , except that a totally antisymmetric linear combination of its S 3 images and those of Q ep are related, up to products of lower-weight functions and ordinary HPLs (see eq. (B.50)). Therefore we count only five independent functions arising from the S 3 permutations of M 1 .
We present the {n − 1, 1} components of the coproduct of these 11 basis functions in appendix B. This information, together with the value of the function at the point (1, 1, 1) (which we take to be zero in all but one case), is sufficient to uniquely define the basis of hexagon functions. We will elaborate on these ideas in the next section.

Integral Representations
In the previous section, we described an iterative procedure to construct the basis of hexagon functions in terms of multiple polylogarithms in the y variables. The result is a fully analytic, numerically efficient representation of any given basis function. While convenient for many purposes, this representation is not without some drawbacks. Because S y has one more element than S u , and because the first entry condition is fairly opaque in the y variables, the multiple polylogarithm representation is often quite lengthy, which in turn sometimes obscures interesting properties. Furthermore, the iterative construction and the numerical evaluation of multiple polylogarithms are best performed when the y i are real-valued, limiting the kinematic regions in which these methods are practically useful.
For these reasons, it is useful to develop a parallel representation of the hexagon functions, based directly on the system of first-order differential equations they satisfy. These differential equations can be solved in terms of (iterated) integrals over lower-weight functions. Since most of the low weight functions are HPLs, which are easy to evaluate, one can obtain numerical representations for the hexagon functions, even in the kinematic regions where the y i are complex. The differential equations can also be solved in terms of simpler functions in various limits, which will be the subject of subsequent sections.

General setup
One benefit of the construction of the basis of hexagon functions in terms of multiple polylogarithms is that we can explicitly calculate the coproduct of the basis functions. We tabulate the {n − 1, 1} component of the coproduct for each of these functions in appendix B. This data exposes how the various functions are related to one another, and, moreover, this web of relations can be used to define a system of differential equations that the functions obey. These differential equations, together with the appropriate boundary conditions, provide an alternative definition of the hexagon functions. In fact, as we will soon argue, it is actually possible to derive these differential equations iteratively, without starting from an explicit expression in terms of multiple polylogarithms. It is also possible to express the differential equations compactly in terms of a Knizhnik-Zamolodchikov equation along the lines studied in ref. [19]. Nevertheless, the coproduct on multiple polylogarithms, in particular the {n − 1, 1} component as given in eq. (3.47), is useful to frame the discussion of the differential equations and helps make contact with section 3.
It will be convenient to consider not just derivatives with respect to a cross ratio, as in eq. (3.48), but also derivatives with respect to the y variables. For that purpose, we need the following derivatives, which we perform holding y v and y w constant, We also consider the following linear combination, Using eqs. (2.23) and (4.1), as well as the definition (4.2), we obtain three differential equations (plus their cyclic images) relating a function F to its various coproduct components, Let us assume that we somehow know the coproduct components of F , either from the explicit representations given in appendix B, or from the iterative approach that we will discuss in the next subsection. We then know the right-hand sides of eqs. (4.3)-(4.5), and we can integrate any of these equations along the appropriate contour to obtain an integral representation for the function F . While eq. (4.3) integrates along a very simple contour, namely a line that is constant in v and w, this also means that the boundary condition, or initial data, must be specified over a two-dimensional plane, as a function of v and w for some value of u. In contrast, we will see that the other two differential equations have the convenient property that the initial data can be specified on a single point.
Let us begin with the differential equation (4.4) and its cyclic images. For definiteness, we consider the differential equation in y v . To integrate it, we must find the contour in (u, v, w) that corresponds to varying y v , while holding y u and y w constant. Following ref. [25], we define the three ratios, . Two of these ratios, r and s, are actually independent of y v , while the third, t, varies. Therefore, we can let t parameterize the contour, and denote by (u t , v t , w t ) the values of the cross ratios along this contour at generic values of t. Since r and s are constants, we have two constraints, (4.7) We can solve these equations for v t and w t , giving, Finally, we can change variables so that u t becomes the integration variable. Calculating the Jacobian, we find, There are two natural basepoints for the integration: u t = 0, for which y v = 1 and (u, v, w) = (0, 1, 0); and u t = 1, for which y v = 1/(y u y w ) and (u, v, w) = (1, 1, 1). Both choices have the convenient property that they correspond to a surface in terms of the variables (y u , y v , y w ), but only to a single point in terms of the variables (u, v, w). This latter fact allows for the simple specification of boundary data. For most purposes, we choose to integrate off of the point u t = 1, in which case we find the following solution to the differential equation, (4.10) The last step follows from the observation that The integral representation (4.10) for F may be ill-defined if the integrand diverges at the lower endpoint of integration, u t = 1 or (u, v, w) = (1, 1, 1). On the other hand, for F to be a valid hexagon function, it must be regular near this point, and therefore no such divergence can occur. In fact, this condition is closely related to the constraint of good branch-cut behavior near u = 1 discussed in section 3.4. As we build up integral representations for hexagon functions, we will use this condition to help fix various undetermined constants. Furthermore, if F is a parity-odd function, we may immediately conclude that F (1, 1, 1) = 0, since this point corresponds to the surface y u y v y w = 1. If F is parity even, we are free to define the function by the condition that F (1, 1, 1) = 0. We use this definition for all basis functions, except for Ω (2) (u, v, w), whose value at (1, 1, 1) is specified by its correspondence to a particular Feynman integral.
While eq. (4.10) gives a representation that can be evaluated numerically for most points in the unit cube of cross ratios 0 ≤ u i ≤ 1, it is poorly suited for Region I. The problem is that the integration contour leaves the unit cube, requiring a cumbersome analytic continuation of the integrand. One may avoid this issue by integrating along the same contour, but instead starting at the point u t = 0 or (u, v, w) = (0, 1, 0). The resulting representation is, If F is a parity-odd function, then the boundary value F (0, 1, 0) must vanish, since this point corresponds to the EMRK limit y v → 1. In the parity-even case, there is no such condition, and in many cases this limit is in fact divergent. Therefore, in contrast to eq. (4.10), this expression may require some regularization near the u t = 0 endpoint in the parity-even case. It is also possible to integrate the differential equation (4.5). In this case, we look for a contour where y v and y u y w are held constant, while the ratio y u /y w is allowed to vary. The result is a contour Again, there are two choices for specifying the boundary data: either we set y u /y w = y u y w for which we may take u t = 0 and (u, v, w) = (0, 0, 1); or we set y u /y w = 1/(y u y w ), for which we may take u t = 1 and (u, v, w) = (1, 0, 0). We therefore obtain two different integral representations, and, Here we used the relation, which follows from the observation that √ ∆/v is constant along either integration contour. Finally, we remark that the boundary values F (1, 0, 0) and F (0, 0, 1) must vanish for parity-odd functions, since the points (1, 0, 0) and (0, 0, 1) lie on the ∆ = 0 surface. In the parity-even case, there may be issues of regularization near the endpoints, just as discussed for eq. (4.12).
Altogether, there are six different contours, corresponding to the three cyclic images of the two types of contours just described. They may be labeled by the y-variables or their ratios that are allowed to vary along the contour: {y u , y v , y w , y u /y w , y v /y u , y w /y v }. The base points for these contours together encompass (1, 1, 1), (0, 1, 0), (1, 0, 0) and (0, 0, 1), the four corners of a tetrahedron whose edges lie on the intersection of the surface ∆ = 0 with the unit cube. See fig. 2 for an illustration of the contours passing through the point (u, v, w) = ( 3 4 , 1 4 , 1 2 ).

Constructing the hexagon functions
In this subsection, we describe how to construct differential equations and integral representations for the basis of hexagon functions. We suppose that we do not have any of the function-level data that we obtained from the analysis of section 3; instead, we will develop a completely independent alternative method starting from the symbol. The two approaches are complementary and provide important cross-checks of one another. In section 3.1, we presented the construction of the basis of hexagon functions at symbol level. Here we will promote these symbols to functions in a three-step iterative process: 1. Use the symbol of a given weight-n function to write down an ansatz for the {n − 1, 1} component of its coproduct in terms of a function-level basis at weight n−1 that we assume to be known.
2. Fix the undetermined parameters in this ansatz by imposing various function-level consistency conditions. These conditions are: (a) Symmetry. The symmetries exhibited by the symbol should carry over to the function.
(b) Integrability. The ansatz should be in the image of ∆ n−1,1 . This condition is equivalent to the consistency of mixed partial derivatives.
(c) Branch cuts. The only allowed branch cuts start when a cross ratio vanishes or approaches infinity.
3. Integrate the resulting coproduct using the methods of the previous subsection, specifying the boundary value and thereby obtaining a well-defined function-level member of the hexagon basis.
Let us demonstrate this procedure with some examples. Recalling the discussion in section 3.1, any function whose symbol contains no y variables can be written as products of single-variable HPLs. Therefore, the first nontrivial example occurs at weight three. As previously mentioned, this function corresponds to the one-loop six-dimensional hexagon integral,Φ 6 . Its symbol is given by, It is straightforward to identify the object in brackets as the symbol of a linear combination of weight-two hexagon functions (which are just HPLs), allowing us to write an ansatz for the {2, 1} component of the coproduct, for some undetermined rational number a.
The single constant, a, can be fixed by requiring thatΦ 6 have the same symmetries as its symbol. In particular, we demand thatΦ 6 be odd under parity. As discussed in the previous section, this implies that it must vanish in the limit that one of the y i goes to unity. In this EMRK limit (3.50), the corresponding u i goes to unity while the other two cross ratios go to zero. The right-hand side of eq. (4.18) vanishes in this limit only for the choice a = −2. So we can write, where, confirming the expression given in eq. (2.21). It is also straightforward to verify that eq. (4.18) is integrable and that it does not encode improper branch cuts. We will not say more about these conditions here, but we will elaborate on them shortly, in the context of our next example. Now that we have the coproduct, we can use eqs. (4.4) and (4.5) to immediately write down the differential equations, These derivatives lead, via eqs. (4.10) and (4.14), to the following integral representations: with (u t , v t , w t ) as in eq. (4.8), or with (u t , v t , w t ) as in eq. (4.13). We have set the integration constants to zero becauseΦ 6 is a parity-odd function.
We have now completed the construction of the hexagon basis through weight three. Moving on to weight four, the symbol-level classification reveals one new parity-even function, Ω (2) (u, v, w), and one new parity-odd function, F 1 (u, v, w), as well as their cyclic images. We will discuss the parity-even function Ω (2) (u, v, w) since it exhibits a variety of features that the parity-odd functions lack.
As discussed in section 2, Ω (2) (u, v, w) is an extra-pure function, and as such its symbol has only three distinct final entries, which were given in eq. (2.18), Furthermore, the symbol is symmetric under the exchange of u with v. Taken together, these symmetry properties dictate the form of the {3, 1} component of the coproduct, There are two independent functions in eq. (4.26), Ω (2),u and Ω (2),yu . The symbols of these functions can be read off from the symbol of Ω (2) (u, v, w). Both functions must be valid hexagon functions of weight three. The symbol indicates that Ω (2),u is parity-even and Ω (2),yu is parity-odd. The most general linear combination of parity-even hexagon functions of weight three whose symbol is consistent with that of Ω (2),u is for four arbitrary rational numbers a i . There is only a single parity-odd hexagon function of weight three, so Ω (2),yu is uniquely determined from its symbol, It is not necessarily the case that the right hand side of eq. (4.26) is actually the {3, 1} component of the coproduct of a well-defined function for arbitrary values of the parameters a i . This integrability condition can be formalized by the requirement that the operator annihilate the right hand side of eq. (4.26). To see this, note that (∆ 2,1 ⊗ id) • ∆ 3,1 = ∆ 2,1,1 , and therefore d ∧ d acts on the last two slots, which are just weight-one functions (logarithms). This can be recognized as the familiar symbol-level integrability condition, eq. (2.2), promoted to function-level.
Another way of thinking about the integrability condition is that it guarantees the consistency of mixed partial derivatives. Since there are three variables, there are three pairs of derivatives to check. To illustrate the procedure, we will examine one pair of derivatives by verifying the equation, We have multiplied by an overall factor of √ ∆ for convenience. To simplify the notation, let us define, U ≡ Ω (2),u and V ≡ Ω (2),u | u↔v . (4.32) The algebra leading to the second line may be simplified by using the fact that (1 − w)/ √ ∆ is independent of y w . Similarly, it is straightforward to write down an expression for the right-hand side of eq. (4.30), where we have used the fact that w/ √ ∆ is annihilated by ∂/∂ ln(y v /y u ). As usual, the superscripts indicate the various coproduct components. A special feature of this example is that the functions U and V are built entirely from single-variable HPLs, so it is straightforward to extract these coproduct components using the definitions in appendix A. More generally, the functions may contain non-HPL elements of the hexagon basis. For these cases, the coproduct components are already known from previous steps in the iterative construction of the basis.
The nonzero coproduct components of U are,  Continuing in this way, we can derive similar constraints from the remaining two mixed partial derivative consistency conditions. The result is that a 2 = −1 , and a 3 = 1 . (4.36) Finally, we must impose good branch-cut behavior. As discussed in section 3.4, this constraint can be implemented by imposing eq. (3.52), or, in this case, which implies that a 4 = 0. The one remaining parameter, a 1 , corresponds to an ambiguity that cannot be fixed by considering mathematical consistency conditions. Indeed, it arises from a well-defined weightfour function with all the appropriate symmetries and mathematical properties. In particular, it is the product of ζ 2 with an extra-pure weight-two hexagon function that is symmetric under In general, we would resolve such an ambiguity by making an arbitrary (though perhaps convenient) choice in order to define the new hexagon function. But because Ω (2) (u, v, w) corresponds to a particular Feynman integral, the value of a 1 is not arbitrary, and the only way to fix it is to bring in specific data about that integral. We are not interested in determining the value of a 1 directly from the Feynman integral since this integral has been evaluated previously [25]. Instead, we will be satisfied simply to verify that a consistent value of a 1 exists. From eq. (4.33) we have, Equation (4.39) is consistent with the differential equations of section 4 of ref. [25] only if the function Q φ from that reference (and eq. (2.28)) is related to U and V by, We remark that the antisymmetric combination appearing in eq. (4.40) is related to another function defined in ref. [25] whereZ appears in a derivative of the odd part of the NMHV ratio function (see eq. (B.19)). Following the discussion in section 4.1, the differential equation (4.39) gives rise to the integral representation, While our conventions for generic hexagon functions require the functions to vanish at the boundary value (1, 1, 1), in this specific case we must specify a nonzero value Ω (2) (1, 1, 1) = −6ζ 4 in order to match a prior definition of the function. The differential equation (4.40) gives rise to another integral representation for Ω (2) , where, There is no constant of integration in eq. (4.45) because in this case Ω (2) vanishes at the lower endpoint, Ω (2) (1, 0, 0) = 0 [62,25].
Continuing onward, we construct the remaining functions of the hexagon basis in an iterative fashion, using the above methods. We collect the results through weight five in appendix B. We present the data by the {n − 1, 1} component of the coproduct, plus the constraint that the functions vanish at (u, v, w) = (1, 1, 1) (except for the special case of Ω (2) ). With this information, we can build an ansatz for the three-loop remainder function, as we discuss in the next subsection.

Constructing the three-loop remainder function
In this subsection, we complete the construction of an ansatz for the three-loop remainder function. We use the decomposition (2.29) of the symbol of R 2. We fix as many undetermined parameters in this ansatz as possible by enforcing various mathematical consistency conditions. In particular, (a) We impose extra-purity and symmetry in the exchange of u and v as function-level conditions on the coproduct entries, since these conditions are satisfied at symbol level: In principle, beyond-the-symbol terms do not need to obey the extra-purity relations. At the end of section 5, we will relax this assumption and use the near-collinear limits to show that the potential additional terms vanish.
(b) We demand that the ansatz be integrable. For the multiple polylogarithm approach, this amounts to verifying that there is a weight-six function with our ansatz as the {5, 1} component of its coproduct. For the approach based on integral representations, we check that there are consistent mixed partial derivatives.
(c) We require that the resulting function have the proper branch-cut structure. We impose this constraint by verifying that there are no spurious poles in the first derivatives, just as we did in the construction of the hexagon basis.
After imposing these constraints, there are still nine undetermined beyond-the-symbol parameters. They correspond to well-defined extra-pure hexagon functions of weight six, and cannot be fixed by mathematical consistency conditions.
3. We integrate the resulting coproduct. This result is a weight-six function, R (α 1 ,α 2 ) ep , which depends on the symbol-level constants, α 1 and α 2 , and nine lower-weight functions r 1 , . . . , r 9 , which come multiplied by zeta values. The r i may be expressed in terms of previouslydetermined hexagon functions, while R (α 1 ,α 2 ) ep may be given as an integral representation or explicitly in terms of multiple polylogarithms.
This procedure leaves us with the following ansatz for R 6 : where the c i are undetermined rational numbers, and (4.49) In the following section we will use the collinear limits of this expression to fix α 1 , α 2 and the c i . After fixing these parameters, we can absorb all but the constant terms into a redefinition of R ep . The {5, 1} component of its coproduct is given in appendix C. The final integral representation for R 6 , having fixed also c 10 and c 11 , is given in section 7, eq. (7.1). The final expression in terms of multiple polylogarithms is quite lengthy, but it is provided in a computer-readable format in the attached files.

Collinear limits
In the previous section, we constructed a 13-parameter ansatz for the three-loop remainder function. It has the correct symbol, proper branch structure, and total S 3 symmetry in the cross ratios. In other words, the ansatz obeys all relevant mathematical consistency conditions. So in order to fix the undetermined constants, we need to bring in some specific physical data.
Some of the most useful data available comes from the study of the collinear limit. In the strict collinear limit in which two gluons are exactly collinear, the remainder function must vanish to all loop orders. This condition fixes many, but not all, of the parameters in our ansatz. To constrain the remaining constants, we expand in the near-collinear limit, keeping track of the power-suppressed terms. These terms are predicted by the OPE for flux tube excitations. In fact, the information about the leading discontinuity terms in the OPE [10,11,12] was already incorporated at symbol level and used to constrain the symbol for the three-loop remainder function up to two undetermined parameters [24].
Here we take the same limit at function level, and compare to the recent work of Basso, Sever and Vieira (BSV) [31], which allows us to uniquely constrain all of the beyond-the-symbol ambiguities, as well as the two symbol-level parameters. The two symbol-level parameters were previously fixed by using dual supersymmetry [29], and also by studying the near-collinear limit at symbol level [31], and we agree with both of these determinations.

Expanding in the near-collinear limit
In the (Euclidean) limit that two gluons become collinear, one of the cross ratios goes to zero and the sum of the other two cross ratios goes to one. For example, if we let k 2 and k 3 become parallel, then x 2 24 ≡ (k 2 + k 3 ) 2 → 0, corresponding to v → 0, and u + w → 1. BSV [31] provide a convenient set of variables (τ, σ, φ) with which one can approach this collinear limit. They are related to the (u i , y i ) variables by [33]: where T = e −τ , S = e σ , and F = e iφ .
As T → 0 (τ → ∞) we approach the collinear limit. The parameter S controls the partitioning of the momentum between the two collinear gluons, according to k 2 /k 3 ∼ S 2 , or k 2 /(k 2 + k 3 ) ∼ S 2 /(1 + S 2 ). The parameter F controls the azimuthal dependence as the two gluons are rotated around their common axis with respect to the rest of the scattering process. This dependence is related to the angular momentum of flux-tube excitations in the OPE interpretation.
By expanding an expression in T we can probe its behavior in the near-collinear limit, order by order in T . Each order in T also contains a polynomial in ln T . In general, the expansions of parity-even and parity-odd hexagon functions f even and f odd have the form, BSV predict the full order T 1 behavior of the remainder function [31]. The part of the T 2 behavior that is simplest for them to predict (because it is purely gluonic) contains azimuthal variation proportional to cos 2 φ, i.e. the T 2 F 2 or T 2 F −2 terms; however, they can also extract the T 2 F 0 behavior, which depends upon the scalar and fermionic excitations as well [33]. To compare with this data, we must expand our expression for R (3) 6 to this order. The expansion of an expression is relatively straightforward when its full analytic form is known, for example when the expression is given in terms of multiple polylogarithms. In this case, one merely needs to know how to take a derivative with respect to T and how to evaluate the functions at T = 0. The derivative of a generic multiple polylogarithm can be read off from its coproduct, which is given in appendix A. Evaluating the functions at T = 0 is more involved because it requires taking y u → 1 and y w → 1 simultaneously. However, the limit of all relevant multiple polylogarithms can be built up iteratively using the coproduct bootstrap of section 3.3.
If the expression is instead represented in integral form, or is defined through differential equations, then it becomes necessary to integrate up the differential equations, iteratively in the transcendental weight, and order by order in the T expansion. Recall that for any function in our basis we have a complete set of differential equations whose inhomogeneous terms are lower weight hexagon functions. The change of variables (5.1) and its Jacobian allow us to go from differential equations in the u i or y variables to differential equations in (F, S, T ).
The structure of the T → 0 expansion makes most terms very straightforward to integrate. In eqs. (5.2) and (5.3), T only appears as powers of T , whose coefficients are polynomials of fixed order in ln T . The variable F only appears as a polynomial in cos φ and sin φ, i.e. as powers of F and F −1 . Hence any T or F derivative can be integrated easily, up to a constant of integration, which can depend on S. The S derivatives require a bit of extra work. However, the differential equation in S is only required for the T -and F -independent term arising in the parity-even case, f even 0,0,0 (S). This coefficient is always a pure function of the same transcendental weight as f itself, and it can be constructed from a complete set of HPLs in the argument −S 2 . Thus we can integrate the one required differential equation in S by using a simple ansatz built out of HPLs. There is still one overall constant of integration to determine for each parity-even function, a term that is completely independent of T , F and S. It is a linear combination of zeta values. (The parity-odd functions all vanish as T → 0, so they do not have this problem.) The constant of integration can be determined at the endpoint S = 0 or S = ∞, with the aid of a second limiting line, (u, v, w) = (u, u, 1). On this line, all the hexagon functions are very simple, collapsing to HPLs with argument (1 − u). In the limit u → 0 this line approaches the point (0, 0, 1), which can be identified with the S → 0 "soft-collinear" corner of the T → 0 collinear limit in the parametrization (5.1). Similarly, the S → ∞ corner of the T → 0 limit intersects the line (1, v, v) at v = 0. Both lines (u, u, 1) and (1, v, v) pass through the point (1, 1, 1). At this point, (most of) the hexagon functions are defined to vanish, which fixes the integration constants on the (u, u, 1) and (1, v, v) lines. HPL identities then give the desired values of the functions in the soft-collinear corner, which is enough to fix the integration constant for the near-collinear limit. We will illustrate this method with an example below.
The coefficients of the power-suppressed terms that also depend on T and F , namely f m,n,p (S) in eqs. (5.2) and (5.3) for m > 0, are functions of S that involve HPLs with the same argument −S 2 , but they also can include prefactors to the HPLs that are rational functions of S. The f m,n,p (S) for m > 0 generally have a mixed transcendental weight. Mixed transcendentality is common when series expanding generic HPLs around particular points. For example, expanding Li 2 (1 − x) around x = 0 gives Using an HPL ansatz for the pure S-dependent terms, we use the differential equations to fix any unfixed parameters and cross-check the ansatz. Repeating this process order by order we build up the near-collinear limiting behavior of each element of the basis of hexagon functions as a series expansion.

Examples
In order to illustrate the collinear expansion, it is worthwhile to present a few low-weight examples. We begin with the simplest nontrivial example, the weight-three parity-odd functioñ Φ 6 . SinceΦ 6 is fully symmetric in the u i and vanishes in the collinear limit (like any parity-odd function), its expansion is particularly simple. To conserve space in later formulas, we adopt the notation, The expansion ofΦ 6 is theñ The sign of eq. (5.6), and of the collinear expansions of all of the parity-odd functions, depend on the values of the y variables used. This sign is appropriate to approaching the collinear limit from Region I, with 0 < y i < 1.
Because Ω (2) lacks the symmetries ofΦ 6 , its expansion must be evaluated in multiple channels, and it is substantially lengthier. Through order T 2 , we find, The integral Ω (2) (u, v, w) is symmetric under the exchange of u and v. This implies that the limiting behavior of Ω (2) (v, w, u) can be determined from that of Ω (2) (u, v, w) by exchanging the roles of u and w in the collinear limit. At leading order in T , this symmetry corresponds to letting S ↔ 1/S. This symmetry is broken by the parametrization (5.1) at order T 2 ; nevertheless, the correction at order T 2 is relatively simple, The last independent permutation is Ω (2) (w, u, v). It is symmetric under u ↔ w and vanishes at order T 0 , which together imply that its near-collinear expansion is symmetric under S ↔ 1/S through order T 2 , although that symmetry is not manifest in the HPL representation, We determine these expansions by integrating the differential equations in F , S, and T , as described in the previous subsection. For parity-even functions, it is necessary to fix the constants of integration. Here we present one technique for doing so. Suppose we set S = T in eq. (5.1). Then the limit T → 0 corresponds to the EMRK limit, u = v → 0, w → 1, approached along the line (u, u, 1). As an example, let us consider applying this limit to the expansion of Ω (2) (u, v, w), eq. (5.7). We only need to keep the T 0 terms, and among them we find that the H w terms vanish, L → ln u, and ln T → 1 2 ln u (since u ∼ T 2 ). Therefore, as u → 0 we obtain, The constant of integration, 5 2 ζ 4 , clearly survives in this limit. So, assuming we did not know its value, it could be fixed if we had an independent way of examining this limit.
This independent method comes from the line (u, u, 1), on which all the hexagon functions have simple representations. This can be seen from the form of the integration contour parametrized by v t and w t in eq. (4.8). Setting v = u and w = 1, it collapses to v t = u t , w t = 1 . The integral (4.43) then becomes where ω u (u, u, 1) = [Ω (2),u + (u ↔ v)](u, u, 1) = 2 H u 3 + H u 2,1 + ln u H u 2 + 1 2 ln 3 u .
At the point u = 1, all the H u w = H w (1 − u) vanish, as does ln u = −H u 1 , so we see that eq. (5.14) becomes in agreement with the explicit −6ζ 4 in eq. (5.12). In order to take the limit u → 0, we use HPL identities to reexpress the function in terms of HPLs with argument u instead of (1 − u) : (5. 16) In the limit u → 0, the H w (u) vanish, leaving only the zeta values and powers of ln u, which are in complete agreement with eq. (5.10). In particular, the coefficient of ζ 4 agrees, and this provides a generic method to determine such constants.
In this example, we inspected the (u, u, 1) line, whose u → 0 limit matches the S → 0 limit of the T → 0 expansion. One can also use the (1, v, v) line in exactly the same way; its v → 0 limit matches the S → ∞ limit of the T → 0 expansion.
Continuing on in this fashion, we build up the near-collinear expansions through order T 2 for all of the functions in the hexagon basis, and ultimately for R 6 itself. The expansions are rather lengthy, but we present them in a computer-readable file attached to this document.

Fixing most of the parameters
In section 4.3 we constructed an ansatz (4.48) for R (3) 6 that contained 13 undetermined rational parameters, after imposing mathematical consistency and extra-purity of R ep . Two of the parameters affect the symbol: α 1 and α 2 . (They could have been fixed using a dual supersymmetry anomaly equation [29].) The remaining 11 parameters c i we refer to as "beyond-the-symbol" because they accompany functions (or constants) with Riemann ζ value prefactors. Even before we compare to the OPE expansion, the requirement that R 6 vanish at order T 0 in the collinear limit is already a powerful constraint. It represents 11 separate conditions when it is organized according to powers of ln T , ln S 2 and H w (−S 2 ), as well as the Riemann ζ values. (There is no dependence on F at the leading power-law order.) The 11 conditions lead to two surviving free parameters. They can be chosen as α 2 and c 9 .
Within R ep , the coefficient c 9 multiplies ζ 2 Ω (2) (u, v, w), as seen from eq. (4.49). However, after summing over permutations, imposing vanishing in the collinear limit, and using eq. (2.14), c 9 is found to multiply ζ 2 R (2) 6 . It is clear that c 9 cannot be fixed at this stage (vanishing at order T 0 ) because the two-loop remainder function vanishes in all collinear limits. Furthermore, its leading discontinuity is of the form T m (ln T ), which is subleading with respect to the three-loop leading discontinuity, terms of the form T m (ln T ) 2 . It is rather remarkable that there is only one other ambiguity, α 2 , at this stage.
The fact that α 1 can be fixed at the order T 0 stage was anticipated in ref. [24]. There the symbol multiplying α 1 was extended to a full function, called f 1 . It was observed that the collinear limit of f 1 , while vanishing at symbol level, did not vanish at function level, and the limit contained a divergence proportional to ζ 2 ln T times a particular function of S 2 . It was argued that this divergence should cancel against contributions from completing the α i -independent terms in the symbol into a function. Now that we have performed this step, we can fix the value of α 1 . Indeed when we examine the ζ 2 ln T terms in the collinear limit of the full R

Comparison to flux tube OPE results
In order to fix α 2 and c 9 , as well as obtain many additional consistency checks, we examine the expansion of R 6 to order T and T 2 , and compare with the flux tube OPE results of BSV. BSV formulate scattering amplitudes in planar N = 4 super-Yang-Mills theory, or rather the associated polygonal Wilson loops, in terms of pentagon transitions. The pentagon transitions map flux tube excitations on one edge of a light-like pentagon, to excitations on another, non-adjacent edge. They have found that the consistency conditions obeyed by the pentagon transitions can be solved in terms of factorizable S matrices for two-dimensional scattering of the flux tube excitations. These S matrices can in turn be determined nonperturbatively for any value of the coupling, as well as expanded in perturbation theory in order to compare with perturbative results [31,32]. The lowest twist excitations dominate the near-collinear or OPE limit τ → ∞ or T → 0. The twist n excitations first appear at O(T n ). In particular, the O(T 1 ) term comes only from a gluonic twist-one excitation, whereas at O(T 2 ) there can be contributions of pairs of gluons, gluonic bound states, and pairs of scalar or fermionic excitations. As mentioned above, BSV have determined the full order T 1 behavior [31], and an unpublished analysis gives the T 2 F 2 or T 2 F −2 terms, plus the expansion of the T 2 F 0 terms around S = 0 through S 10 [33].
BSV consider a particular ratio of Wilson loops: the basic hexagon Wilson loop, divided by two pentagons, and then multiplied back by a box (square). The pentagons and box combine to cancel off all of the cusp divergences of the hexagon, leading to a finite, dual conformally invariant ratio. We compute the remainder function, which can be expressed as the hexagon Wilson loop divided by the BDS ansatz [15] for Wilson loops. To relate the two formulations, we need to evaluate the logarithm of the BDS ansatz for the hexagon configuration, subtract the analogous evaluation for the two pentagons, and add back the one for the box. The pentagon and box kinematics are determined from the hexagon by intersecting a light-like line from a hexagon vertex with an edge on the opposite side of the hexagon [31]. For example, if we have lightlike momenta k i , i = 1, 2, . . . , 6 for the hexagon, then one pentagon is found by replacing three of the momenta, say k 4 , k 5 , k 6 , with two light-like momenta, say k ′ 4 and k ′ 5 , having the same sum. Also, one of the new momenta has to be parallel to one of the three replaced momenta: The requirement that k ′ 5 is a null vector implies that The correction term to go between the logarithm of the BSV Wilson loop and the six-point remainder function requires the combination of one-loop normalized amplitudes V n (from the BDS formula [15]), which is finite and dual conformal invariant. There is also a prefactor proportional to the cusp anomalous dimension, whose expansion is known to all orders [75], . Including the proper prefactor, we obtain the following relation between the two observables, Here W hex is BSV's observable (they use the expansion parameter g 2 = λ/(16π 2 ) = a/2) and R 6 is the remainder function. In the near-collinear limit, the correction function X(u, v, w) becomes, Next we apply this relation in the near-collinear limit, first at order T 1 . We find that the T 1 ln 2 T terms from BSV's formula match perfectly the ones we obtain from our expression for R 6 . The T 1 ln T terms also match, given one linear relation between α 2 and the coefficient of ζ 2 R (2) 6 . Finally, the T 1 ln 0 T terms match if we fix α 2 = 7/32, which is the last constant to be fixed. The value of α 2 is in agreement with refs. [29,31]. The agreement with ref. [31] (BSV) is no surprise, because both are based on comparing the near-collinear limit of R Here we give the formula for the leading, order T term in the near-collinear limit of R 6 , after fixing all parameters as just described: The T 2 terms are presented in an attached, computer-readable file. The T 2 terms match perfectly with OPE results provided to us by BSV [33], and at this order there are no free parameters in the comparison. This provides a very nice consistency check on two rather different approaches.
Recall that we imposed an extra-pure condition on the terms in eq. (4.49) that we added to the ansatz for R (3) 6 . We can ask what would happen if we relaxed this assumption. To do so we consider adding to the solution that we found a complete set of beyond-the-symbol terms. Imposing total symmetry, there are 2 weight-6 constants (ζ 6 and (ζ 3 ) 2 ), and 2 weight-5 constants (ζ 5 and ζ 2 ζ 3 ) multiplying ln uvw. Multiplying the zeta values ζ 4 , ζ 3 and ζ 2 there are respectively 3, 7 and 18 symmetric functions, for a total of 32 free parameters. Imposing vanishing of these additional terms at order T 0 fixes all but 5 of the 32 parameters to be zero. We used constraints from the multi-Regge limit (see the next section) to remove 4 of the 5 remaining parameters. Finally, the order T 1 term in the near-collinear limit fixes the last parameter to zero. We conclude that there are no additional ambiguities in R 6 associated with relaxing the extra-purity assumption.

Multi-Regge limits
The multi-Regge or MRK limit of n-gluon scattering is a 2 → (n − 2) scattering process in which the (n − 2) outgoing gluons are strongly ordered in rapidity. It generalizes the Regge limit of 2 → 2 scattering with large center-of-mass energy at fixed momentum transfer, s ≫ t. Here we are interested in the case of 2 → 4 gluon scattering, for which the MRK limit means that two of the outgoing gluons are emitted at high energy, almost parallel to the incoming gluons. The other two gluons are also typically emitted at small angles, but they are well-separated in rapidity from each other and from the leading two gluons, giving them smaller energies.
The strong ordering in rapidity for the 2 → 4 process leads to the following strong ordering of momentum invariants: (6.1) In this limit, the cross ratio u = s 12 s 45 /(s 123 s 345 ) approaches one. The other two cross ratios vanish, In this section, we denote the original cross ratio w byŵ, in order to avoid confusion with another variable which we are about to introduce. The cross ratios v andŵ vanish at the same rate that u → 1, so that the ratios x and y, defined by remain fixed. The variable y in eq. (6.3) should not be confused with the variables y i . In the y variables, the multi-Regge limit consists of taking y u → 1, while y v and y w are left arbitrary.
(Their values in this limit are related to x and y by eq. (3.51).) It is very convenient [39] to change variables from x and y to the complex-conjugate pair (w, w * ) defined by, . (6.4) (Again, this variable w should not be confused with the original cross ratio calledŵ in this section.) This change of variables rationalizes the y variables in the MRK limit, so that As an aside, we remark here that the variables T, S, F in eq. (5.1), used by BSV to describe the near-collinear limit, are closely related to the variables w, w * introduced for the MRK limit. To establish this correspondence, we consider (in this paragraph only) the MRK limit u → 0, v → 0,ŵ → 1, which is related to eq. (6.2) by a cyclic permutation u i → u i−1 , y i → y i−1 . This limit corresponds to the T → 0 limit in eq. (5.1) if we also send S → 0 at the same rate, so that T /S is fixed. Let's rewrite y u from eq. (5.1) as and compare it with the limiting behavior of y v in eq. (6.5). (Comparing y u with y v is required by the cyclic permutation of the u i and y i variables which we need for the two limits to correspond.) then y v in eq. (6.5) correctly matches eq. (6.6). If we start with the variables T, S, F in eq. (5.1), insert the inverse relations to eq. (6.7), and then let S → 0 with w, w * fixed, we can check that all variables approach the values appropriate for the multi-Regge limit u → 0, v → 0,ŵ → 1. The cross-ratioŵ approaches unity as S vanishes, through the relationŵ = (1 + S 2 |1 + w| 2 ) −1 . Finally, we note that the MRK limit interpolates between three different limits: the collinear limit v → 0, corresponding to |w|→ 0; the endpoint of the line (u, u, 1) with u → 0, corresponding to w → −1; and a second collinear limit u → 0, corresponding to |w|→ ∞. Now we return to the u → 1 version of the MRK limit in eq. (6.2). If this limiting behavior of the cross ratios is approached directly from the Euclidean region in which all cross ratios are positive, we call it the EMRK limit (see also eq. (3.50)). In this limit, the remainder function vanishes, as it does in the Euclidean collinear limit discussed in the previous section. However, the physical region for 2 → 4 scattering is obtained by first analytically continuing u → e −2πi u, then taking u → 1, v,ŵ → 0 as above. The analytic continuation generates imaginary terms corresponding to the discontinuity of the function in the u channel, which survive into the MRK limit; in fact they can be multiplied by logarithmic singularities as u → 1.
The general form of the remainder function at L loops in the MRK limit is where the coefficient functions g r , by using a crossing relation from the 3 → 3 channel [40,45].
The coefficient functions in this limit are built out of HPLs with arguments −w and −w * . Only special combinations of such HPLs are allowed, with good branch-cut behavior in the (w, w * ) plane, corresponding to symbols whose first entries are limited to x and y [45]. Such functions may be called single-valued harmonic polylogarithms (SVHPLs), and were constructed by Brown [46].
Using a Fourier-Mellin transformation, Fadin, Lipatov, and Prygarin wrote an all-loop expression for the MRK limit in a factorized form depending on two quantities, the BFKL eigenvalue ω(ν, n) and the impact factor Φ Reg (ν, n) [41]: . (6.10) Here where the cusp anomalous dimension γ K (a) is given in eq. (5.22). By taking the MRK limit of the symbol of the three-loop remainder function, it was possible to determine all of the coefficient functions g (l) r and h (l) r through three loops, up to four undetermined rational numbers, d 1 , d 2 , γ ′ and γ ′′ , representing beyond-the-symbol ambiguities [24]. (Two other parameters, c and γ ′′′ , could be fixed using consistency between the MRK limits in 2 → 4 kinematics and in 3 → 3 kinematics.) One of these four constants was fixed by Fadin and Lipatov [41], using a direct calculation of the NLLA BFKL eigenvalue: γ ′ = −9/2. The remaining three undetermined constants, d 1 , d 2 and γ ′′ , all appear in the NNLLA coefficient g (3) 0 (w, w * ). In ref. [45], the coefficient functions g (3) r (w, w * ) and h (3) r (w, w * ) that appear in the MRK limit (6.9) of R (3) 6 were expressed in terms of the SVHPLs defined in ref. [46]. More specifically, they were rewritten in terms of particular linear combinations of SVHPLs, denoted by L ± w , that have definite eigenvalues under inversion of w and under its complex conjugation. The coefficient function g (3) 0 (w, w * ) then becomes [45]: In the remainder of this section we will describe how to extract the MRK limit of the threeloop remainder function at the full function level. Comparing this limit with eq. (6.13) (as well as the other g , and it will also provide for us the remaining three-loop MRK constants, d 1 , d 2 and γ ′′ .

Method for taking the MRK limit
Let us begin by discussing a method for taking the multi-Regge limit of hexagon functions in general, or of R (3) 6 in particular, starting from an expression in terms of multiple polylogarithms. The first step is to send u → e −2πi u, i.e. to extract the monodromy around u = 0. Owing to the non-linear relationship between the u i and the y i , eq. (2.11), it is not immediately clear what the discontinuity looks like in the y variables. The correct prescription turns out simply to be to take y u around 0. To see this, consider the ∆ 1,n−1 component of the coproduct, which can be written as, There are only three terms, corresponding to the three possible first entries of the symbol.
Using the coproduct formulas in appendix A, it is straightforward to extract the functions u F , v F , and w F for any given hexagon function. These functions capture information about the discontinuities as each of the cross ratios is taken around zero. In particular, since the monodromy operator acts on the first component of the coproduct, we have (c.f. eq. (3.34)), Equation (6.15) is not quite sufficient to deduce M u=0 (F ). The obstruction comes from the fact that all higher powers of (2πi) live in the kernel of ∆ 1,n−1 . On the other hand, these terms can be extracted from the other components of the coproduct: the (2πi) k terms come from the piece of ∆ k,n−k (F ) with ln k u in the first slot. If we write eq. (6.14) in terms of the y i variables, we find, where we have now assumed that we are working in Region I. Equation (6.16) indicates that u F can be extracted uniquely from the terms with G(0; y u ) in the first slot. Similarly, the elements of the full coproduct with ln k u in the first slot are given exactly by the terms with G(0; y u ) k in the first slot. Therefore the discontinuity around u = 0 is the same as the discontinuity around y u = 0. Furthermore, because our basis G L I exposes all logarithms G(0; y u ) (by exploiting the shuffle algebra), the only sources of such discontinuities are powers of G(0; y u ). As a result, we have a simple shortcut to obtain the monodromy around u = 0, M u=0 (F ) = F | G(0;yu) → G(0;yu)−2πi . (6.17) The final step in obtaining the MRK limit is to take y u → 1. This limit is trivially realized on functions in the basis G L I because the only source of singularities is G(1; y u ); all other functions are finite as y u → 1. Writing the divergence in terms of ξ ≡ 1 − u, which approaches 0 in this limit, we take and then set y u = 1 in all other terms. The result of this procedure will be a polynomial in ln ξ whose coefficients are multiple polylogarithms in the variables y v and y w . On the other hand, we know from general considerations that the coefficient functions should be SVHPLs. To translate the multiple polylogarithms into SVHPLs, we use the coproduct bootstrap of section 3.3, seeded by the weight-one identities which follow from eq. (6.5) and from combining eqs. (3.51), (6.4) and (6.5), We obtain, Alternatively, we can extract the MRK limits of the hexagon functions iteratively in the weight, by using their definitions in terms of differential equations. This procedure is similar to that used in section 5 to find the collinear limits of the hexagon functions, in that we expand the differential equations around the limiting region of u → 1.
However, first we have to compute the discontinuities from letting u → e −2πi u in the inhomogeneous (source) terms for the differential equations. For the lowest weight non-HPL function, Φ 6 , the source terms are pure HPLs. For pure HPL functions we use standard HPL identities to exchange the HPL argument (1 − u) for argument u, and again use the Lyndon basis so that the trailing index in the weight vector w in each H w (u) is 1. In this new representation, the only discontinuities come from explicit factors of ln u, which are simply replaced by ln u − 2πi under the analytic continuation. After performing the analytic continuation, we take the MRK limit of the pure HPL functions.
Once these limits are known, we can integrate up the differential equations for the non-HPL functions in much the same fashion that we did for the collinear limits, by using a restricted ansatz built from powers of ln ξ and SVHPLs. The Jacobian factors needed to transform from differential equations in (u, v,ŵ) to differential equations in the MRK variables (ξ, w, w * ), are easily found to be: We compute the derivatives on the right-hand side of these relations using the formula for ∂F/∂u i in terms of the coproduct components, eq. (4.3). We also implement the transformation u → e −2πi u on the coproduct components, as described above for the HPLs, and iteratively in the weight for the non-HPL hexagon functions. When we expand as ξ → 0, we drop all powersuppressed terms in ξ, keeping only polynomials in ln ξ. (In ∂F/∂ξ, we keep the derivatives of such expressions, i.e. terms of the form 1/ξ × ln k ξ.) In our definition of the MRK limit, we include any surviving terms from the EMRK limit. This does not matter for the remainder function, whose EMRK limit vanishes, but the individual parity-even hexagon functions can have nonzero, and even singular, EMRK limits.

Examples
We first consider the simplest non-HPL function,Φ 6 . Starting with the expression forΦ 6 in Region I, eq. (3.56), we take the monodromy around u = 0, utilizing eq. (6.17), (6.23) Next, we take the limit y u → 1. There are no divergent factors, so we are free to set y u = 1 without first applying eq. (6.18). The result is, To transform this expression into the SVHPL notation of ref. [45], we use the coproduct bootstrap to derive an expression for the single independent SVHPL of weight two, the Bloch-Wigner dilogarithm, L − 2 , In the last line we used eqs. (6.20) and (6.21). Lifting eq. (6.25) from coproducts to functions introduces one undetermined rational-number constant, proportional to ζ 2 . It is easily fixed by specializing to the point y v = y w = 1, yielding, which, when compared to eq. (6.24), gives, Let us derive this result in a different way, using the method based on differential equations. Like all parity-odd functions,Φ 6 vanishes in the Euclidean MRK limit; however, it survives in the MRK limit due to discontinuities in the function Ω (1) given in eq. (2.22), which appears on the right-hand side of theΦ 6 differential equation (2.24). The MRK limits of the three cyclic permutations of Ω (1) are given by (6.28) Inserting these values into eq. (2.24) for ∂Φ 6 /∂u and its cyclic permutations, and then inserting those results into eq. (6.22), we find that (6.29) The first differential equation implies that there is no ln ξ term in the MRK limit ofΦ 6 . The second two differential equations imply that the MRK limit is proportional to the Bloch-Wigner dilogarithm,Φ (6.30) Now that we have the MRK limit ofΦ 6 , we can find the limiting behavior of all the coproduct components of Ω (2) appearing in eq. (4.26), and perform the analogous expansion of the differential equations in the MRK limit. For Ω (2) (u, v,ŵ) we obtain, plus the complex conjugate equation for ∂Ω (2) (u, v,ŵ)/∂w * . The solution to these differential equations can be expressed in terms of SVHPLs. One can write an ansatz for the result as a linear combination of SVHPLs, and fix the coefficients using the differential equations. One can also take the limit first at the level of the symbol, matching to the symbols of the SVHPLs; then one only has to fix the smaller set of beyond-the-symbol terms using the differential equations. The result is In this case the constant term, proportional to ζ 3 , can be fixed by requiring vanishing in the collinear-MRK corner where |w| 2 → 0. The last set of terms, multiplying (4π) 2 , come from a double discontinuity. The MRK limit of Ω (2) (ŵ, u, v) is related by symmetry to that of Ω (2) (u, v,ŵ): The final MRK limit of Ω (2) is, where L X = ln ξ + L + 1 . Note that this orientation of Ω (2) has a nonvanishing (indeed, singular) EMRK limit, i.e. even before analytically continuing into the Minkowski region to pick up the imaginary part. On the other hand, there is no surviving double discontinuity for this ordering of the arguments.
As our final (simple) example, we give the MRK limit of the totally symmetric, weight five, parity-odd function G(u, v,ŵ). As was the case forΦ 6 , the limit of G is again proportional to the Bloch-Wigner dilogarithm, but with an extra factor of ζ 2 to account for the higher transcendental weight of G: As usual for parity-odd functions, the EMRK limit vanishes. In this case the double discontinuity also vanishes. In general the MRK limits of the parity-odd functions must be odd under w ↔ w * , which forbids any nontrivial constants of integration. Continuing onward, we build up the MRK limits for all the remaining hexagon functions. The results are attached to this document in a computer-readable format.
6.3 Fixing d 1 , d 2 , and γ ′′ Using the MRK limit of all the hexagon functions appearing in eq. (4.48), we obtain the MRK limit of R 6 . This is a powerful check of the function, although as mentioned above, much of it is guaranteed by the limiting behavior of the symbol. In fact, there are only three rational parameters to fix, d 1 , d 2 and γ ′′ , and they all enter the coefficient of the NNLLA imaginary part, g 0 (w, w * ), given in eq. (6.13). Inspecting the MRK limit of R r (w, w * ) entering the real part. (These can be determined on general grounds using consistency between the 2 → 4 and 3 → 3 MRK limits.) We also agree perfectly with the imaginary part coefficients g 2 at LLA and g 1 at NLLA. Finally, we find for the NNLLA coefficient g 0 (w, w * ) = (6.36) Comparing this result with eq. (6.13) fixes the three previously undetermined rational parameters, d 1 , d 2 , and γ ′′ . We find These three parameters were also the only ambiguities in the expression found in ref. [45] for the two-loop (NNLLA) impact factor Φ (2) Reg (ν, n) defined in ref. [41]. Inserting eq. (6.37) into that expression, we obtain, (6.38) Here Φ (1) Reg is the one-loop (NLLA) impact factor, and E ν,n and E (1) ν,n are the LLA and NLLA BFKL eigenvalues [41,45]. These functions all are combinations of polygamma (ψ) functions and their derivatives, plus accompanying rational terms in ν and n. For example, Additional rational dependence on ν and n enters eq. (6.38) via the combinations We recall that the NNLLA BFKL eigenvalue E (2) ν,n also has been determined [45], up to nine rational parameters, a i , i = 0, 1, 2, . . . , 8. These parameters enter the NNLLA coefficient function g (4) 1 (w, w * ). If the above exercise can be repeated at four loops, then it will be possible to fix all of these parameters in the same way, and obtain an unambiguous result for the NNLLA approximation to the MRK limit.
Finally, we ask whether we could have determined all coefficients from the collinear vanishing of R (3) 6 and the MRK limit alone, i.e. without using the near-collinear information from BSV. The answer is yes, if we assume extra purity and if we also take the value of α 2 from ref. [29]. After imposing collinear vanishing, we have two parameters left: α 2 and the coefficient of ζ 2 R (2) 6 . We can fix the latter coefficient in terms of α 2 using the known NLLA coefficient g 3 . We find that we can fix d 2 and γ ′′ to the values in eq. (6.37), but that α 2 is linked to d 1 by the equation, If we do take α 2 from ref. [29], then the near-collinear limit of our result for R provides an unambiguous test of BSV's approach at three loops, through O(T 2 ).

Final formula for R
(3) 6 and its quantitative behavior Now that we have used the (near) collinear limits to fix all undetermined constants in eq. (4.48) for R 6 , we can write an expression for the full function, either in terms of multiple polylogarithms or integral representations. We absorb the c i r i (u, v) terms in eq. (4.48) into R ep . In total we have, in terms of multiple polylogarithms, valid in Regions I and II, are too lengthy to present here, but they are attached to this document in computer-readable format. To represent R ep as an integral, we make use of its extra purity and similarity to Ω (2) (u, v, w), writing a formula similar to eq. (4.43): with v t and w t as defined in eq. (4.44). Note that the function Q φ in eq. (4.43) is given, via eq. (4.41), as −[Ω (2),u + (u ↔ v)], the analogous combination of coproduct components entering eq. (7.2). The function R u ep is defined in appendix C. We may also define R 6 via the {5, 1} component of its coproduct, which is easily constructed from the corresponding coproducts of R ep in appendix C, and of the product function P 6 . The general form of the {5, 1} component of the coproduct is, ⊗ ln y w . (w, u, v) .

(7.4)
The two independent functions may be written as, Since the {5, 1} component of the coproduct specifies all the first derivatives of R 6 , eqs. (7.5) and (7.6) should be supplemented by the value of R 6 at one point. For example, the value at (u, v, w) = (1, 1, 1) will suffice (see below), or the constraint that it vanishes in all collinear limits.
In the remainder of this section, we use the multiple polylogarithmic and integral representations to obtain numerical values for R  On the line (u, u, 1), the two-and three-loop remainder functions can be expressed solely in terms of HPLs of a single argument, 1 − u. The two-loop function is, while the three-loop function is, We will see that over large swaths of the positive octant, the ratio R 6 does not stray too far from −7. We plot the function R (7.14) Hence the ratio R 6 diverges logarithmically as u → 0 along this line: 6 (u, u, 1) ∼ 1 2 ln u, as u → 0. (7.15) This limit captures a piece of the near-collinear limit T → 0, the case in which S → 0 at the same rate, as discussed in section 5 near eq. (5.10). The fact that R 6 has one more power of ln u than does R (2) 6 is partly from its extra leading power of ln T (the leading singularity behaves like (ln T ) L−1 ), but also from an extra ln S 2 factor in a subleading ln T term. As u → ∞, the leading behavior at two and three loops is, 6 (u, u, 1) = 6097 96 As u → ∞ along the line (u, u, 1), the two-and three-loop remainder functions, and thus their ratio R 6 /R 6 , approach a constant. For the ratio it is: 6 (u, u, 1) ∼ − 50 3 (ζ 3 ) 2 π 4 + 871 972 π 2 = −9.09128803107 . . . , as u → ∞. (7.17) We plot the ratio R 6 on the line (u, u, 1) in fig. 4. The logarithmic scale for u highlights how little the ratio varies over a broad range in u.
The line (u, u, 1) is special in that the remainder function is extra pure on it. That is, applying the operator u(1 − u) d/du returns a pure function for L = 2, 3: The extra-pure property is related to the fact that the asymptotic behavior as u → ∞ is merely a constant, with no ln u terms. Indeed, if one applies u(1 − u) d/du to any positive power of ln u, the result diverges at large u like u times a power of ln u, which is not the limiting behavior of any combination of HPLs in H u .

The line (1, 1, w)
We next consider the line (1, 1, w). As was the case for the line (u, u, 1), we can express the two-and three-loop remainder functions on the line (1, 1, w) solely in terms of HPLs of a single argument. However, in contrast to (u, u, 1), the expressions on the line (1, 1, w) are not extra-pure functions of w.
The two-loop result is, It is not extra pure on this line, because the quantity contains explicit factors of w.

The line (u, u, u)
The symmetrical diagonal line (u, u, u) has the nice feature that the remainder function at strong coupling can be written analytically. Using AdS/CFT to map the problem to a minimal area one, and applying integrability, Alday, Gaiotto and Maldacena obtained the strong-coupling result [52], where φ = 3 cos −1 (1/ √ 4u). The extra constant term −π 2 /12 is needed in order for R (∞) 6 (u, v, w) to vanish properly in the collinear limits [76]. 2 In perturbation theory, the function R (L) 6 (u, v, w) is less simple to represent on the line (u, u, u) than on the lines (u, u, 1) and (1, 1, w). It cannot be written solely in terms of HPLs with argument (1−u). At two loops, using eq. (2.14), the only obstruction is the function Ω (2) (u, u, u), One way to proceed is to convert the first-order partial differential equations for all the hexagon functions of (u, v, w) into ordinary differential equations in u for the same functions evaluated on the line (u, u, u). The differential equation for the three-loop remainder function itself is, with similar differential equations for Ω (2) (u, u, u), H 1 (u, u, u) and J 1 (u, u, u). Interestingly, the parity-even weight-five functions M 1 and Q ep do not enter eq. (7.26).
We can solve the differential equations by using series expansions around three points: u = 0, u = 1, and u = ∞. If we take enough terms in each expansion (of order 30-40 terms suffices), then the ranges of validity of the expansions will overlap. At u = 1, ∆ vanishes, and so do all the parity-odd functions, so we divide them by √ ∆ before series expanding in (u − 1). These expansions, and those of the parity-even functions, are regular, with no logarithmic coefficients, as expected for a point in the interior of the positive octant. (Indeed, we can perform an analogous three-dimensional series expansion of all hexagon functions of (u, v, w) about (1, 1, 1); this is actually a convenient way to fix the beyond-the-symbol terms in the coproducts, by using consistency of the mixed partial derivatives.) At u = 0, the series expansions also contain powers of ln u in their coefficients. At u = ∞, there are two types of terms in the generic series expansion: a series expansion in 1/u with coefficients that are powers of ln u, and a series expansion in odd powers of 1/ √ u with an overall factor of π 3 , and coefficients that can contain powers of ln u. The square-root behavior can be traced back to the appearance of factors of ∆(u, u, u) = (1 − u) √ 1 − 4u in the differential equations, such as eq. (7.26).
The constants of integration are easy to determine at u = 1 (where most of the hexagon function are defined to be zero). They can be determined numerically (and sometimes analytically) at u = 0 and u = ∞, either by evaluating the multiple polylogarithmic expressions, or by matching the series expansions with the one around u = 1.
At small u, the series expansions at two and three loops have the following form: 6 (u, u, u) = − (7.27) while the strong-coupling result is, Note that the leading term at three loops diverges logarithmically, but only as ln 2 u. Alday, Gaiotto and Maldacena [52] observed that this property holds at two loops and at strong coupling, and predicted that it should hold to all orders. At large u, the two-and three-loop remainder functions behave as, 2 ln 3 u + 15 ln 2 u + 6 (6ζ 2 + 11) ln u + 24ζ 3 + 126ζ 2 + 138 6 (u, u, u) = −

48
ζ 6 + ζ 2 3 + 3 π 5 4 u 1/2 while the strong-coupling behavior is, In fig. 6 we plot the two-and three-loop and strong-coupling remainder functions on the line (u, u, u). In order to highlight the remarkably similar shapes of the three functions for small and moderate values of u, we rescale R where R (∞) 6 (1, 1, 1) = π/6 − π 2 /12. A necessary condition for the shapes to be so similar is that the limiting behavior of the ratios as u → 0 is almost the same as the ratios' values at u = 1. From eq. (7.27), the three-loop to two-loop ratio as u → 0 is, R 6 (u, u, u) R (2) 6 (u, u, u) ∼ − 21 5 ζ 2 = −6.908723 . . . , as u → 0, (7.32) which is within 1.5% of the ratio at (1, 1, 1), eq. (7.13). The three-loop to strong-coupling ratio is, R which is again within 1.5% of the corresponding ratio (7.31) at u = 1. The similarity of the perturbative and strong-coupling curves for small and moderate u suggests that if a smooth extrapolation of the remainder function from weak to strong coupling can be achieved, on the line (u, u, u) it will have a form that is almost independent of u, for u < 1.
As mentioned in the introduction, the numerical similarity of two-loop and strong-coupling remainder functions has been explored previously, starting with two-dimensional kinematics at eight points [53,54], and later for the general 2n-point case [57]. The rescaling of the remainder functions used to perform those comparisons is similar to our rescaling in fig. 6. For the six-point case studied in the present paper, ref. [59] has compared the two-loop and strong-coupling rescaled remainder functions along a curve which runs from (u, v, w) = (1/4, 1/4, 1/4) to (1, 0, 0), as well as analytically in the expansion around (1/4, 1/4, 1/4) using conformal perturbation theory, and the results are very similar. The curve runs from the ultraviolet to the infrared region of the renormalization group flow associated with an integrable two-dimensional system. It would be very interesting to perform this comparison with the three-loop remainder function as well, but we will reserve this exercise for future work.
We return now to the line (u, u, u). Whereas all the curves in fig. 6 are very similar for u < 1, they diverge from each other at large u, although they each approach a constant value as u → ∞. The three-to-two-loop ratio at very large u, from eq. (7.29), eventually approaches −1.227 . . ., which is quite different from −7. The three-loop-to-strong-coupling ratio approaches −3.713 . . ., which is very different from −63.4.
On the line (u, u, u), all three curves in fig. 6 cross zero very close to u = 1/3. The respective zero crossing points for L = 2, 3, ∞ are: below. Another way to examine the progression of perturbation theory, and its possible extrapolation to strong coupling, is to use the Wilson loop ratio adopted by BSV, which is related to the remainder function by eq. (5.23). This relation holds for strong coupling as well as weak coupling, since the cusp anomalous dimension is known exactly [75]. In the near-collinear limit, considering the Wilson loop ratio has the advantage that the strong-coupling OPE behaves sensibly. The remainder function differs from this ratio by the one-loop function X(u, v, w), whose near-collinear limit does not resemble a strong-coupling OPE at all. On the other hand, the Wilson loop ratio breaks all of the S 3 permutation symmetries of the remainder function. This is not an issue for the line (u, u, u), since none of the S 3 symmetries survive on this line. However, there is also the issue that X(u, u, u) as determined from eq. (5.24) diverges logarithmically as u → 1.
In fig. 7 we plot the perturbative coefficients of ln[1+W hex (a/2)], as well as the strong-coupling value, restricting ourselves to the range 0 < u < 1 where X(u, u, u) remains real. Now there is also a one-loop term, from multiplying X(u, u, u) by the cusp anomalous dimension in eq. (5.23). We normalize the results in this case by dividing the coefficient at a given loop order by the corresponding coefficient of the cusp anomalous dimension, and similarly at strong coupling. Equivalently, from eq. (5.23), we plot R (L) 6 (u, u, u) γ for L = 1, 2, 3, ∞.
The Wilson loop ratio diverges at both u = 0 and u = 1. The divergence at u = 1 comes only from X and is controlled by the cusp anomalous dimension. This forces the curves to converge in this region. The ln 2 u divergence as u → 0 gets contributions from both X and R 6 . The latter contributions are not proportional to the cusp anomalous dimensions, causing all the curves to split apart at small u. Because X(u, u, u) crosses zero at u = 0.394 . . ., which is a bit different from the almost identical zero crossings in eq. (7.34) and in fig. 6, the addition of X in fig. 7 splits the zero crossings apart a little. However, in the bulk of the range, the perturbative coefficients do alternate in sign from one to three loops, following the sign alternation of the cusp anomaly coefficients, and suggesting that a smooth extrapolation from weak to strong coupling may be possible for this observable as well.

Planes of constant w
Having examined the remainder function on a few one-dimensional lines, we now turn to its behavior on various two-dimensional surfaces. We will now restrict our analysis to the unit cube, 0 ≤ u, v, w ≤ 1. To provide a general picture of how the remainder function behaves throughout this region, we show in fig. 8 the function evaluated on planes with constant w, as a function of u and v. The plane w = 1 is in pink, w = 3 4 in purple, w = 1 2 in dark blue, and w = 1 4 in light blue. The function goes to zero for the collinear-EMRK corner point (u, v, w) = (0, 0, 1) (the right corner of the top sheet). Except for this point, R  be seen clearly from the plot, the w = 1 4 surface actually intersects the w = 1 2 surface near this corner.) 7.5 The plane u + v − w = 1 Next we consider the plane u + v − w = 1. Its intersection with the unit cube is the triangle bounded by the lines (1, v, v) and (w, 1, w), which are equivalent to the line (u, u, 1) discussed in section 7.1, and by the collinear limit line (u, 1 − u, 0), on which the remainder function vanishes.
In fig. 9 we plot the ratio R 6 on this triangle. The back edges can be identified with the u < 1 portion of fig. 4, although here they are plotted on a linear scale rather than a logarithmic scale. The plot is symmetrical under u ↔ v. In the bulk of the triangle, the ratio does not stray far from −7. The only place it deviates is in the approach to the collinear limit, the front edge of  the triangle corresponding to T → 0 in the notation of section 5. Both R 6 contains an extra power of ln T in its vanishing, and so the ratio diverges like ln T . Otherwise, the shapes of the two functions agree remarkably well on this triangle. 7.6 The plane u + v + w = 1 The plane u + v + w = 1 intersects the unit cube along the three collinear lines. In fig. 10 we give a contour plot of R 6 (u, v, w) on the equilateral triangle lying between these lines. The plot has the full S 3 symmetry of the triangle under permutations of (u, v, w). Because R 6 has to vanish on the boundary, one might expect that it should not get too large in the interior. Indeed, its furthest deviation from zero is slightly under −0.07, at the center of the triangle.
From the discussion in section 7.3 and fig. 6, we know that along the line (u, u, u) the twoand three-loop remainder functions almost always have the opposite sign. The only place they have the same sign on this line is for a very short interval u ∈ [0.3325, 0.3343] (see eq. (7.34)). This interval happens to contain the point (1/3, 1/3, 1/3), which is the intersection of the line (u, u, u) with the plane in fig. 10, right at the center of the triangle. In fact, throughout the entire unit cube, the only region where R 6 on the triangle here, but it is easy to verify that it is also uniformly negative in the region  is perfectly continuous across the ∆ = 0 surface. We can also see that the function diverges as w goes to zero, and as u and v go to zero, everywhere except for the two places that this plane intersects the collinear limits, namely the points (u, v, w) = (1/2, 1/2, 0) and (u, v, w) = (0, 0, 1). The line of intersection of the u = v plane and the u + v + w = 1 plane passes through both of these points, and fig. 11 shows that R 6 is very close to zero on this line.
Based on considerations related to the positive Grassmannian [50], it was recently conjectured [77] that the three-loop remainder function should have a uniform sign in the "positive region", or what we call Region I: the portion of the unit cube where ∆ > 0 and u + v + w < 1, which corresponds to positive external kinematics in terms of momentum twistors. On the sur-  face u = v, this region lies in front of the parabola shown in fig. 11. It was already checked [77] that the two-loop remainder function has a uniform (positive) sign in Region I. Fig. 11 illustrates that the uniform sign behavior (with a negative sign) is indeed true at three loops on the plane u = v. We have checked many other points with u = v in Region I, and R also holds in the other regions of the unit cube with ∆ > 0, namely Regions II, III, and IV, which are all equivalent under S 3 transformations of the cross ratios. In these regions, the overall signs are reversed: R only seem to occur where ∆ < 0, and in fact very close to u + v + w = 1.
6 on the plane u + v = 1. This plane provides information complementary to that on the plane u = v, since the two planes intersect at right angles. Like the u = v plane, this plane shows smooth behavior over the ∆ = 0 surface, which intersects the plane u + v = 1 in the parabola w = 4u(1 − u). It also shows that the function vanishes smoothly in the w → 0 collinear limit.

Conclusions
In this paper, we successfully applied a bootstrap, or set of consistency conditions, in order to determine the three-loop remainder function R 6 (u, v, w) directly from a few assumed analytic properties. We bypassed altogether the problem of constructing and integrating multi-loop integrands. This work represents the completion of a program begun in ref. [24], in which the symbol S(R (3) 6 ) was determined via a Wilson loop OPE and certain conditions on the final entries, up to two undetermined rational numbers that were fixed soon thereafter [29].
In order to promote the symbol to a function, we first had to characterize the space of globally well-defined functions of three variables with the correct analytic properties, which we call hexagon functions. Hexagon functions are in one-to-one correspondence with the integrable symbols whose entries are drawn from the nine letters {u i , 1−u i , y i }, with the first entry restricted to {u i }. We specified the hexagon functions at function level, iteratively in the transcendental weight, by using their coproduct structure. In this approach, integrability of the symbol is promoted to the function-level constraint of consistency of mixed partial derivatives. Additional constraints prevent branch-cuts from appearing except at physical locations (u i = 0, ∞). These requirements fix the beyond-the-symbol terms in the {n − 1, 1} coproduct components of the hexagon functions, and hence they fix the hexagon functions themselves (up to the arbitrary addition of lower-weight functions multiplied by zeta values). We found explicit representations of all the hexagon functions through weight five, and of R itself at weight six, in terms of multiple polylogarithms whose arguments involve simple combinations of the y variables. We also used the coproduct structure to obtain systems of coupled first-order partial differential equations, which could be integrated numerically at generic values of the cross ratios, or solved analytically in various limiting kinematic regions.
Using our understanding of the space of hexagon functions, we constructed an ansatz for the function R in the collinear limits fixed all but one of these parameters. The last parameter was fixed using the near-collinear limits, in particular the T 1 ln T terms which we obtained from the OPE and integrability-based predictions of Basso, Sever and Vieira [31]. (The T 1 ln 0 T terms are also needed to fix the last symbol-level parameter [31] independently of ref. [29].) With all parameters fixed, we could unambiguously extract further terms in the near-collinear limit. We find perfect agreement with Basso, Sever and Vieira's results through order T 2 [33]. We have also evaluated the remainder function in the multi-Regge limit. This limit provides additional consistency checks, and allows us to fix three undetermined parameters in an expression [45] for the NNLLA impact parameter Φ (2) Reg (ν, n) in the BFKL-factorized form of the remainder function [41].
Finally, we found simpler analytic representations for R 6 along particular lines in the threedimensional (u, v, w) space; we plotted the function along these and other lines, and on some two-dimensional surfaces within the unit cube 0 ≤ u i ≤ 1. Throughout much of the unit cube, and sometimes much further out from the origin, we found the approximate numerical relation R 6 . The relation has only been observed to break down badly in regions where the functions vanish: the collinear limits, and very near the plane u + v + w = 1. On the diagonal line (u, u, u), we observed that the two-loop, three-loop, and strong-coupling [52] remainder functions are almost indistinguishable in shape for 0 < u < 1.
We have verified numerically a conjecture [77] that the remainder function should have a uniform sign in the "positive" region {u, v, w > 0; ∆ > 0; u + v + w < 1}. It also appears to have an (opposite) uniform sign in the complementary region {u, v, w > 0; ∆ > 0; u + v + w > 1}. The only zero-crossings we have found for either R Our work opens up a number of avenues for further research. A straightforward application is to the NMHV ratio function. Knowledge of the complete space of hexagon functions through weight five allowed us to construct the six-point MHV remainder function at three loops. The components of the three-loop six-point NMHV ratio function are also expected [25] to be weightsix hexagon functions. Therefore they should be constructible just as R was, provided that enough physical information can be supplied to fix all the free parameters.
It is also straightforward in principle to push the remainder function to higher loops. At four loops, the symbol of the remainder function was heavily constrained [45] by the same information used at three loops [24], but of order 100 free parameters were left undetermined. With the knowledge of the near-collinear limits provided by Basso, Sever and Vieira [31,33], those parameters can now all be fixed. Indeed, all the function-level ambiguities can be fixed as well [34]. This progress will allow many of the numerical observations made in this paper at three loops, to be explored at four loops in the near future. Going beyond four loops may also be feasible, depending primarily on computational issues -and provided that no inconsistencies arise related to failure of an underlying assumption.
It is remarkable that scattering amplitudes in planar N = 4 super-Yang-Mills -polygonal (super) Wilson loops -are so heavily constrained by symmetries and other analytic properties, that a full bootstrap at the integrated level is practical, at least in perturbation theory. We have demonstrated this practicality explicitly for the six-point MHV remainder function. The number of cross ratios increases linearly with the number of points. More importantly, the number of letters in the symbol grows quite rapidly, even at two loops [60], increasing the complexity of the problem. However, with enough understanding of the relevant transcendental functions for more external legs [78,79], it may still be possible to implement a similar procedure in these cases as well. In the longer term, the existence of near-collinear boundary conditions, for which there is now a fully nonperturbative bootstrap based on the OPE and integrability [31], should inspire the search for a fully nonperturbative formulation that also penetrates the interior of the kinematical phase space for particle scattering. manuscript. We are also grateful to Nima Arkani-Hamed, Simon Caron-Huot, John Joseph Carrasco, Claude Duhr, Sasha Goncharov, Diego Hofman, Yuji Satoh, Jaroslav Trnka and Cristian Vergu for helpful discussions. LD, MvH and JP thank the Perimeter Institute for hospitality during the completion of this article. This research was supported by the US Department of Energy under contracts DE-AC02-76SF00515 and DE-FG02-92ER40697, and in part by Perimeter Institute for Theoretical Physics. Research at Perimeter Institute is supported by the Government of Canada through Industry Canada and by the Province of Ontario through the Ministry of Economic Development and Innovation.
A Multiple polylogarithms and the coproduct A.1 Multiple polylogarithms Multiple polylogarithms are a general class of multi-variable iterated integrals, of which logarithms, polylogarithms, harmonic polylogarithms, and various other iterated integrals are special cases. They are defined recursively by G(z) = 1, and, G(a 1 , . . . , a n ; z) = z 0 dt t − a 1 G(a 2 , . . . , a n ; t) , where we have introduced the vector notation a n = (a, . . . , a n ).
For special values of the weight vector (a 1 , . . . , a n ), multiple polylogarithms reduce to simpler functions. For example, if a = 0, where S p,q is the Nielsen polylogarithm. More generally, if a i ∈ {−1, 0, 1}, then G(a 1 , . . . , a n ; z) = (−1) w 1 H a 1 ,...,an (z) , where w 1 is the number of a i equal to one. Multiple polylogarithms are not all algebraically independent. One set of relations, known as the shuffle relations, derive from the definition (A.1) in terms of iterated integrals, where w 1 Xw 2 is the set of mergers of the sequences w 1 and w 2 that preserve their relative ordering. Radford's theorem [66] allows one to solve all of the identities (A.4) simultaneously in terms of a restricted subset of multiple polylogarithms {G(l w ; z)}, where l w is a Lyndon word. The Lyndon words are those words w such that for every decomposition into two words w = {u, v}, the left word is smaller (based on some ordering) than the right, i.e. u < v. One may choose whichever ordering is convenient; for our purposes, we choose an ordering so that zero is smallest. In this case, no zeros appear on the right of a weight vector, except in the special case of the logarithm, G(0; z) = ln z. Therefore, we may adopt a Lyndon basis and assume without loss of generality that a n = 0 in G(a 1 , . . . , a n , z). Referring to eq. (A.1), it is then possible to rescale all integration variables by a common factor and obtain the following identity, G(c a 1 , . . . , c a n ; c z) = G(a 1 , . . . , a n ; z) , a n = 0, c = 0 . (A.5) Specializing to the case c = 1/z, we see that the algebra of multiple polylogarithms is spanned by ln z and G(a 1 , . . . , a n ; 1) where a n = 0. This observation allows us to establish a one-to-one correspondence between multiple polylogarithms and particular multiple nested sums, provided those sums converge. In particular, if for |x i |< 1 we define, Li m 1 ,...,m k (x 1 , . . . , x k ) = n 1 <n 2 <···<n k x n 1 1 x n 2 2 · · · x n k k n m 1 1 n m 2 2 · · · n m k Equation (A.7) is easily established by expanding the measure dt/(t − a i ) in eq. (A.1) in a series and integrating. Furthermore, a convergent series expansion for G(a 1 , . . . , a n ; z) exists if |z|≤ |a i | for all i; otherwise, the integral representation gives the proper analytic continuation. The relation to multiple sums endows the space of multiple polylogarithms with some additional structure. In particular, the freedom to change summation variables in the multiple sums allows one to establish stuffle or quasi-shuffle relations, Li m 1 ( x)Li m 2 ( y) = n Li n ( z) . (A.8) The precise formula for n and z in terms of m 1 , m 2 , x, and y is rather cumbersome, but can be written explicitly; see, e.g., ref. [80]. For small depth, however, the stuffle relations are quite simple. For example, Li a (x)Li b (y) = Li a,b (x, y) + Li b,a (y, x) + Li a+b (xy) . (A.9) Beyond the shuffle and stuffle identities, there are additional relations between multiple polylogarithms with transformed arguments and weight vectors. For example, one such class of identities follows from Hölder convolution [80], G(a 1 , . . . , a n ; 1) = n k=0 (−1) k G 1 − a k , . . . , 1 − a 1 ; 1 − 1 p G a k+1 , . . . , a n ; 1 p , (A.10) which is valid for any nonzero p whenever a 1 = 1 and a n = 0. One way to study identities among multiple polylogarithms is via the symbol, which is defined recursively as, S(G(a n−1 , . . . , a 1 ; a n )) = n−1 i=1 S(G(a n−1 , . . . ,â i , . . . , a 1 ; a n )) ⊗ (a i − a i+1 ) −S(G(a n−1 , . . . ,â i , . . . , a 1 ; a n )) ⊗ (a i − a i−1 ) , (A.11) While the symbol has the nice property that all relations result from simple algebraic manipulations, it has the drawback that its kernel contains all transcendental constants. To obtain information about these constants, one needs some more powerful machinery.

A.2 The Hopf algebra of multiple polylogarithms
When equipped with the shuffle product (A.4), the space of multiple polylogarithms forms an algebra, graded by weight. In ref. [67], it was shown how to further equip the space with a coproduct so that it forms a bialgebra, and, moreover, with an antipode so that it forms a Hopf algebra. The weight of the multiple polylogarithms also defines a grading on the Hopf algebra. In the following we will let A denote the Hopf algebra and A n the weight-n subspace, so that, The coproduct is defined most naturally on a slight variant of eq. (A.1), I(a 0 ; a 1 , . . . , a n ; a n+1 ) = a n+1 a 0 dt t − a n I(a 0 ; a 1 , . . . , a n−1 ; t) . (A.13) The two definitions differ only in the ordering of indices and the choice of basepoint. However, as shown in ref. [26], it is possible to reexpress any multiple polylogarithm with a generic basepoint as a sum of terms with basepoint zero. This manipulation is trivial at weight one, where we have, I(a 0 ; a 1 ; a 2 ) = I(0; a 1 ; a 2 ) − I(0; a 1 ; a 0 ) = G(a 1 ; a 2 ) − G(a 1 ; a 0 ) . (A.14) To build up further such relations at higher weights, one must simply apply the lower-weight identity to the integrand in eq. (A.13). In this way, it is easy to convert between the two different notations for multiple polylogarithms. The coproduct on multiple polylogarithms is given by [67], ∆(I(a 0 ; a 1 , . . . , a n ; a n+1 )) = 0<i 1 <···<i k =n I(a 0 ; a i 1 , . . . , a i k ; a n+1 )⊗ k p=0 I(a ip ; a ip+1 , . . . , a i p+1 −1 ; a i p+1 ) .
(A.15) Strictly speaking, this definition is only valid when the a i are nonzero and distinct; otherwise, one must introduce a regulator to avoid divergent integrals. We refer the reader to refs. [67,26] for these technical details.
It is straightforward to check a number of important properties of the coproduct. First, it respects the grading of A in the following sense. If G n ∈ A n , then, ∆(G n ) = p+q=n ∆ p,q (G n ) , (A. 16) where ∆ p,q ∈ A p ⊗ A q . Next, if we extend multiplication to tensor products so that it acts on each component separately, meaning that one may iterate the coproduct in any order and always reach a unique result. This last property allows one to unambiguously define components of the coproduct corresponding to all integer compositions of the weight. Consider G n ∈ A n and a particular integer composition of n, {w 1 , . . . , w k }, such that w i > 0 and k i=1 w i = n. The component of the coproduct corresponding to this composition, ∆ w 1 ,...,w k (G n ), is defined as the unique element of the (k − 1)-fold iterated coproduct in the space A w 1 ⊗ · · · ⊗ A w k . For our purposes it is sufficient to consider k = 2, although other components have been useful in other contexts.
Consider the weight-n function f The monodromy of f (n) around the point z k = z 0 is encoded by the first entry of the symbol, where M z k =z 0 (ln φ i 1 ) is defined in eq. (6.15), and we have ignored higher powers of (2πi) (see section 6). Similarly, derivatives act on the last entry of the symbol, One may trivially extend the definition of the coproduct to include odd ζ values, ∆(ζ 2n+1 ) = 1 ⊗ ζ 2n+1 + ζ 2n+1 ⊗ 1 (A. 24) but including even ζ values and factors of π is more subtle. It was argued in ref. [22,26] that it is consistent to define, ∆(ζ 2n ) = ζ 2n ⊗ 1 and ∆(π) = π ⊗ 1 .

B Complete basis of hexagon functions through weight five
We present the basis of hexagon functions through weight five by providing their {n − 1, 1} coproduct components. For a hexagon function F of weight n, we write, where the nine functions {F u i , F 1−u i , F y i } are of weight n−1 and completely specify the {n−1, 1} component of the coproduct. They also specify all of the first derivatives of F , The other derivatives can be obtained from the cyclic images of eq. (B.2). These derivatives, in turn, define integral representations for the function. Generically, we define the function F by (see eq. (4.10)), We choose F (1, 1, 1) = 0 for all functions except for the special case Ω (2) (1, 1, 1) = −6ζ 4 . Other integral representations of the function also exist, as discussed in section 4.1.
We remark that the hexagon functionsΦ 6 , G, N and O are totally symmetric under exchange of all three arguments; Ω (2) is symmetric under exchange of its first two arguments; F 1 is symmetric under exchange of its last two arguments; and H 1 , J 1 and K 1 are symmetric under exchange of their first and third arguments.

B.4 G
The {4, 1} component of the coproduct of the parity-odd weight five function G can be written as, Furthermore, G is totally symmetric. In particular, G yv (u, v, w) = G yu (v, w, u) , and G yw (u, v, w) = G yu (w, u, v) . Therefore, it suffices to specify the single independent function, N u , (B.42)

B.10 O
The {4, 1} component of the coproduct of the parity-even weight-five function O can be written as,

B.11 Q ep
The {4, 1} component of the coproduct of the parity-even weight-five function Q ep can be written as, where, The three independent functions consist of one parity-odd function, Q yv ep , which is fairly simple, and two parity-even functions, Q v ep and Q w ep , which are complicated by the presence of a large number of HPLs,

C Coproduct of R ep
We may write the {5, 1} component of the coproduct of the parity-even weight-six function R ep as, where, 2) The two independent functions may be written as,