The four-loop remainder function and multi-Regge behavior at NNLLA in planar N=4 super-Yang-Mills theory

We present the four-loop remainder function for six-gluon scattering with maximal helicity violation in planar N=4 super-Yang-Mills theory, as an analytic function of three dual-conformal cross ratios. The function is constructed entirely from its analytic properties, without ever inspecting any multi-loop integrand. We employ the same approach used at three loops, writing an ansatz in terms of hexagon functions, and fixing coefficients in the ansatz using the multi-Regge limit and the operator product expansion in the near-collinear limit. We express the result in terms of multiple polylogarithms, and in terms of the coproduct for the associated Hopf algebra. From the remainder function, we extract the BFKL eigenvalue at next-to-next-to-leading logarithmic accuracy (NNLLA), and the impact factor at NNNLLA. We plot the remainder function along various lines and on one surface, studying ratios of successive loop orders. As seen previously through three loops, these ratios are surprisingly constant over large regions in the space of cross ratios, and they are not far from the value expected at asymptotically large orders of perturbation theory.


Introduction
Recently there has been great interest in studying perturbative scattering amplitudes in N = 4 super-Yang-Mills theory, both for their own sake and as prototypes for the kinds of mathematical functions that will be encountered at the multi-loop level in QCD and other theories. In the planar limit of a large number of colors, N = 4 super-Yang-Mills amplitudes are dual to Wilson loops for closed polygons with light-like edges, and possess a dual conformal symmetry [1,2,3]. This symmetry is anomalous [4], but the anomaly, as well as various infrared divergences, can be removed by factoring out the BDS ansatz [5]. For the case of the maximally-helicity-violating (MHV) configuration of external gluon helicities, the finite remainder function [6] that is left behind is a function only of the dual conformally invariant cross ratios. The first scattering amplitude to have nontrivial cross ratios and a nonvanishing remainder function is the six-point case, corresponding to a hexagonal Wilson loop, for which there are three such cross ratios.
The remainder function is expected to be a pure transcendental function with a transcendental weight 2L at loop order L. Examples of transcendental functions include the logarithm (weight 1), the classical polylogarithms Li k (weight k), and products thereof. A more general class of transcendental functions is provided by iterated integrals [7], or multiple polylogarithms [8,9]. Other types of functions, such as elliptic integrals, can appear in scattering amplitudes. The two-loop equal-mass sunrise integral is elliptic [10], as is an integral entering a particular 10point scattering amplitude in planar N = 4 super-Yang-Mills theory [11]. However, based on a novel form of the planar loop integrand [12], and also a recent twistor-space formulation [13], it is expected that all six-point amplitudes are non-elliptic and can be described in terms of multiple polylogarithms.
The purpose of this paper is to use the six-point amplitude to demonstrate the power of a bootstrap [14,15,16] for scattering amplitudes in planar N = 4 super-Yang-Mills theory. This bootstrap operates at the level of integrated scattering amplitudes, not loop integrands. It imposes physical constraints at this level, in terms of the external kinematics alone, in order to uniquely determine the final answer. The critical assumption is that the amplitude belongs to a certain space of functions that can be identified at low loop order. In the present case it will be a particular class of iterated integrals. Suppose one can enumerate all such functions and characterize their properties in the kinematic limits that are needed to impose the physical constraints. Then one can write an ansatz for the amplitude as a linear combination of the functions with unknown coefficients (which should all be rational numbers). Physical constraints provide simple linear equations relating the coefficients.
If the basic ansatz is correct, then the only other question of principle is whether there is enough "boundary data"; that is, whether one has enough physical constraints to fix all the coefficients. 1 Fortunately, there is a great deal of data indeed. Much of it comes from the operator product expansion (OPE) for Wilson loops, which corresponds to the near-collinear limit of scattering amplitudes. The OPE was first analyzed by Alday, Gaiotto, Maldacena, Sever and Vieira [17,18,19,20]. More recently, even more powerful OPE information has become available via integrability [21,22,23,24]. The application of integrability to the relevant system of flux tube excitations has been pioneered by Basso, Sever and Vieira (BSV) [25,26,27,28]. We will show that when this data is combined with that from the multi-Regge limit [29,30,31,32,33,34,14,35,36], it is exceedingly powerful, uniquely determining the six-point remainder function through at least four loops.
The need for a remainder function beginning at six points and two loops was first identified in the study of the multi-Regge limit [29], and also from direct numerical evaluation of the amplitude and hexagonal Wilson loop at finite values of the cross ratios [6]. (There were also previous indications at strong coupling that a remainder function would be required, at least in the limit of a large number of external legs [37].) The two-loop hexagon Wilson loop integrals were performed analytically in terms of multiple polylogarithms [38], and then simplified dramatically to classical polylogarithms using the notion of the symbol of a transcendental function [39].
Based on the form of the two-loop symbol, it was conjectured [14,15] that for six-point amplitudes to all loop orders the transcendental functions entering the remainder function (and also the next-to-MHV ratio function [15]) should be polylogarithmic functions whose symbols are made from an alphabet of nine letters, corresponding to nine projectively-inequivalent differences z ij of projective variables z i [39]. These letters can also be represented in terms of momentum twistors [40]. For any weight, there are a finite number of such functions. Using the symbol, one can enumerate them all, and then impose physical constraints on a generic linear combination of them. In this way, the symbol for the three-loop six-point remainder function was obtained, up to two undetermined parameters [14] which were fixed [41] using a dual supersymmetry "anomaly" equation [41,42].
However, the symbol does not determine the full function. Lower-weight functions multiplied by constant Riemann ζ values give rise to pure functions but vanish at the level of the symbol. In ref. [16] it was shown how to identify and fix these parameters at the level of the full three-loop remainder function. In this paper, we will follow the same general strategy at four loops.
In fact, two separate strategies were pursued in ref. [16]. One strategy was to pick a particular region in the space of cross ratios, and promote the symbol to an explicit linear combination of multiple polylogarithms. The additional beyond-the-symbol parameters multiply products of Riemann ζ values with multiple polylogarithms of lower weight. Knowledge of the limiting behavior of the multiple polylogarithms on certain boundaries of this region can then be used to impose the physical constraints. A second strategy is to characterize the remainder function by its coproduct. The coproduct is part of the Hopf algebra conjecturally satisfied by multiple polylogarithms [43,44,45]. It has been applied to a number of different physical problems recently [46,47,48,49,50,51]. In particular, the "{k−1, 1}" element of the coproduct of a weight k function specifies all of its first derivatives in terms of weight k − 1 transcendental functions. One can iterate in the weight, and define a candidate remainder function in terms of a set of coupled first-order differential equations. In the limits relevant for the physical constraints, the coupled equations can be solved in terms of a simpler set of transcendental functions, involving harmonic polylogarithms in a single variable [52]. In the present work, we use the multiplepolylogarithm approach to constrain all of the parameters, and both strategies to examine the limiting behavior of the uniquely-determined function.
Besides certain standard symmetry and parity constraints, and the physical constraints to be described shortly, we also impose a constraint on the final entry of the symbol. The final entry should be expressible in terms of only six letters rather than all nine. This constraint comes from a supersymmetric formulation of the polygonal Wilson loop [53] and also from examining the differential equations obeyed by one-loop [54,55,56] and multi-loop integrals [14,15] related to N = 4 super-Yang-Mills scattering amplitudes. The final-entry constraint on the symbol corresponds to a differential constraint we shall impose at function level, which also has a simple description in terms of the coproduct of the function.
The two limiting regions in which we impose physical constraints on the remainder function are the near-collinear limit and the multi-Regge limit. In the near-collinear limit, one of the cross ratios vanishes and the sum of the other two ratios approaches one. Because the remainder function has a total S 3 permutation symmetry under exchange of the three cross ratios, it does not matter which cross ratio we take to vanish. Let's call this variable v for definiteness, and let v = T 2 + O(T 4 ) as v and T → 0. Because there is no remainder function at the five-point level, the six-point remainder function must vanish as v → 0. The precise way in which it vanishes is controlled by the OPE. The first OPE information to be determined [17,18,19,20] concerned the leading-discontinuity terms, which correspond to just the maximum allowed power of ln T (ln L−1 T at L loops). Terms with arbitrary power suppression in T can be determined, as long as they have L − 1 powers of ln T . These terms are dictated by the one-loop anomalous dimensions of the operators corresponding to excitations of the Wilson line, or flux tube. Higherloop corrections to anomalous dimensions and OPE coefficients only generate terms with fewer logarithms of T . At two loops, the leading discontinuity is the only discontinuity, and it suffices to completely determine the remainder function [18]. At three loops [14], and particularly at four loops, more information is required.
Recently, Basso, Sever and Vieira [26,27,28] were able to exploit integrability in order to provide much more OPE information. They partition a generic polygonal Wilson loop into a number of "pentagon transitions" between flux tube excitations. They find that certain bootstrap consistency conditions for the pentagon transitions can be solved in terms of factorizable S matrices for two-dimensional scattering of the flux tube excitations. These S matrices are known for finite coupling, and they can be expanded out in perturbation theory to any desired order. The powers of T in the OPE expansion correspond to the number of flux tube excitations. In their initial papers [26,27], the leading nonvanishing OPE terms, O(T 1 ), were described, corresponding to single excitations. The O(T 1 ) information, combined with the multi-Regge limits and an assumption about the final entry of the symbol, was enough to completely fix the three-loop remainder function [16]. However, it is not enough at four loops. Fortunately, Basso, Sever and Vieira [28] have also been able to determine the contributions to the OPE of two flux-excitation states, and thereby obtain the O(T 2 ) terms. 2 The O(T 2 ) terms from the OPE were found to agree perfectly with those extracted from the three-loop remainder function [16]. Because there were no free parameters in this comparisonall parameters had been fixed at O(T 1 ) -the agreement is a powerful check on the assumptions underlying both approaches. At four loops, we will need to use some of the O(T 2 ) information, supplied to us by BSV, to fix a small number of remaining parameters in the four-loop remainder function -four parameters in the symbol, and then one more at the level of the full function. However, there is considerably more information in the O(T 2 ) OPE expansion, and so the fact that it agrees between our approach and BSV's at four loops is certainly a strong indication that both approaches are correct.
The other physical limit which can be used to constrain the remainder function is the multi-Regge limit. In this limit, two incoming gluons scatter into four gluons that are well separated in rapidity. Whereas the near-collinear OPE limit can be approached in the Euclidean region, this kinematical configuration is in Minkowski space. Coming from the Euclidean region, one first needs to analytically continue to Minkowski space by rotating the phase of one of the cross ratios, let's call it u, by 2π. Then u should be taken to unity at the same rate that the other two cross ratios, v and w, vanish. The analytic continuation in u generates an imaginary part, as well as a real part from a double discontinuity. Both the imaginary and real parts diverge as powers of ln(1 − u) as u → 1. The leading logarithmic approximation (LLA) has a behavior proportional to ln L−1 (1 − u) at L loops, and it is pure imaginary [29,30,31,32].
It has been proposed that factorization in the multi-Regge limit can be extended to subleading logarithmic accuracy [33,35,58]. In ref. [35] the functions that should control the factorization were computed directly through the next-to-leading-logarithmic approximation (NLLA). In ref. [58] a closely-related form of multi-Regge factorization was proposed, based on the hypotheses of rapidity factorization and the completeness of a description in terms of undecorated, null, infinite Wilson lines. In principle, if these hypotheses are true, then the factorization could hold to arbitrary subleading logarithmic order, up to terms that are power-suppressed like O(1−u). In this paper, we will assume that the factorization holds through arbitrary subleading logarithmic accuracy. In practice, our four-loop results are sensitive to at most N 3 LLA. The fact that we find a consistent solution provides evidence in favor of factorization beyond NLLA.
The assumption of factorization makes it possible to bootstrap multi-Regge information from one loop order to the next. That is, the leading-logarithmic behavior of the remainder function is present already at two loops [29,30] and can be used to predict the LLA ln L−1 (1 − u) behavior at three [33] and higher loops [36]. Similarly, the NLLA behavior [33,35] first appears fully at three loops, and can be used to predict the ln L−2 (1 − u) behavior at four and higher loops [36].
The factorization takes place in variables which are related to the original variables by a Fourier-Mellin transform [35]. Two functions control the expansion: the BFKL eigenvalue and the impact factor. Each function has an expansion in the coupling; successive orders in the expansion are needed for higher accuracy in the logarithmic expansion. The N k LLA term in the impact factor makes its first appearance in the remainder function in the ln 0 (1 − u) term at k + 1 loops; whereas the N k LLA term in the BFKL eigenvalue appears one loop order later, at k + 2 loops, accompanied by one power of ln (1 − u).
In ref. [36] it was observed that in the multi-Regge limit the coefficients in the expansion of remainder function in powers of ln(1 − u) are single-valued harmonic polylogarithms (SVHPLs), first introduced by Brown [59]. Based on this observation, techniques for performing the inverse Fourier-Mellin transform were developed, in order to efficiently find the consequences of the N k LLA approximation for the remainder function at a given loop order. Furthermore, part of the program of this paper to determine the four-loop remainder function was carried out there [36]: Several constraints were applied to the relevant space of symbols: S 3 symmetry, parity, the OPE leading discontinuity and the final-entry condition were applied. These constraints left 113 symbol-level parameters undetermined. However, in the multi-Regge limit only one symbol-level parameter, called a 0 , survived. This allowed the NNLLA BFKL eigenvalue and N 3 LLA impact factor to be almost completely constrained at symbol level. At function level, however, there were an additional 26 undetermined rational numbers.
In the continuation of this program in the present paper, we apply additional multi-Regge constraints from NLLA [35] that we did not impose earlier, in order to fix 33 of the 113 remaining symbol-level parameters. Then we match the O(T 1 ) and O(T 2 ) behavior to the OPE [26,27,28], to fix the final 80 symbol-level parameters. We then account for 68 additional beyond-the-symbol parameters, and fix them all using the same OPE information, which we now implement at the level of full functions using the multiple-polylogarithmic representation. With the remainder function uniquely determined, we return to the Minkowski multi-Regge limit and determine the values of the 27 parameters we had previously introduced. This completes the determination of the NNLLA BFKL eigenvalue and N 3 LLA impact factor begun in ref. [36]. We find that the NNLLA BFKL eigenvalue has a very suggestive form that is closely related to the spectrum of anomalous dimensions for flux tube excitations [25].
We then study the quantitative behavior of the four-loop remainder function in various regions, including special lines in the space of cross ratios where it collapses to linear combinations of harmonic polylogarithms of a single variable. We will explore various numerical observations made at three loops in ref. [16] about the sign and constancy of ratios of successive loop orders. We will find that these observations remain true, and are even reinforced at four loops. We will also discuss how close the remainder function at four loops might be, in a certain region, to its expected behavior at large perturbative orders.
The remainder of this paper is organized as follows. In section 2 we describe the construction of the four-loop remainder function. In section 3 we describe its behavior in the multi-Regge limit and extract the NNLLA BFKL eigenvalue and N 3 LLA impact factor. In section 4 we explore the sign of the four-loop remainder function in a certain "positive" region. We plot the ratio of successive loop orders on a two-dimensional surface, and on various lines where its functional form simplifies considerably, as well as discussing expectations for large perturbative orders. Finally, in section 5 we conclude and discuss avenues for future research. We include one appendix on the coproduct representation, and a second one characterizing logarithmic divergences of the remainder function on two particular boundaries of the Euclidean region.
Many of the analytic results in this paper are too lengthy to present in the manuscript. Instead we provide a set of ancillary files in computer readable format.

The construction 2.1 Hexagon functions
The six-point remainder function R 6 is defined by factoring off the BDS ansatz from the MHV planar amplitude, The BDS ansatz accounts for all of the amplitude's infrared divergences, or ultraviolet divergences in the case of the Wilson loop interpretation. It also absorbs the (related) anomaly in dual conformal transformations. Because R 6 is invariant under such transformations, it can only depend on the dual conformal cross ratios, where x i , 1 ≤ i ≤ 6, denote dual coordinates, related to the external momenta by The six-point remainder function admits the perturbative expansion, starting at two loops, 3 where a = g 2 YM N c /(8π 2 ) is the 't Hooft coupling constant, g YM is the Yang-Mills coupling constant and N c is the number of colors.
The coefficients R (L) 6 (u, v, w) are expected to be pure functions of transcendental weight 2L, i.e., they should be Q-linear combinations of polylogarithmic functions of weight 2L. For this reason, it is convenient to consider the symbol of R (L) 6 (u, v, w). The symbol of a transcendental function f (k) of weight k can most conveniently be defined as follows: if the total differential of f (k) can be written as a finite sum of the form where the φ r are rational functions and the f (k−1) r are transcendental functions of weight k − 1, then the symbol of f (k) can be defined recursively by, (2.5) The six-point remainder function for arbitrary values of the cross ratios is currently known at two [38,39] and three loops [14,16]. One of the main results of this paper is to present the fully analytic answer for the four-loop remainder function R 6 (u, v, w). The construction of the result 3 Beginning at four loops, it is important to specify whether or not R 6 is exponentiated in the definition (2.1), because the two alternative definitions would differ by 1 2 [R (2) 6 ] 2 at this order.
will be performed following closely the ideas of ref. [16], which allow us to bootstrap the four-loop answer without ever inspecting the multi-loop integrand. This bootstrap will be described in the remainder of this section. In ref. [16], a set of polylogarithmic functions called hexagon functions were introduced. Their symbols are built out of the nine letters, where and We sometime also use the labeling u 1 = u, u 2 = v, u 3 = w, y 1 = y u , y 2 = y v , y 3 = y w ). The branch cut locations for hexagon functions are restricted to the points where the cross ratios u i either vanish or approach infinity. In terms of the symbol, this implies [19] that the first entry must be one of the cross ratios u, v, w. In ref. [16], a method based on the coproduct on multiple polylogarithms (or, equivalently, a corresponding set of first-order partial differential equations) was developed that allows for the construction of hexagon functions at arbitrary weight. Using this method, the three-loop remainder function was determined as a particular weight-six hexagon function. In this article, we extend the analysis and construct the four-loop remainder function, which is a hexagon function of weight eight.

Constraints at symbol level
As in the three-loop case, we begin by constructing the symbol. Referring to the discussion in ref. [36], the symbol may be written as where α i are undetermined rational numbers. The S i are drawn from the complete set of eightfold tensor products (i.e. symbols of weight eight) that satisfy the first-entry condition. They also are required to obey the following properties: 0. All entries in the symbol are drawn from the set 1. The symbol is integrable (i.e. it is the symbol of some function).
2. The symbol is totally symmetric under S 3 permutations of the three cross ratios u i .
3. The symbol is invariant under the parity transformation y i → 1/y i . 4. The symbol vanishes in the collinear limit u 2 → 0, u 1 + u 3 → 1. (The other two collinear limits follow from the S 3 symmetry.) 5. In the near-collinear limit, the symbol agrees with the predictions of the leading discontinuity terms in the OPE [17]. We implement this condition exactly as was done at three loops [14].
6. The final entry of the symbol is drawn from the set Imposing the above constraints on the most general ansatz of all 9 8 possible words will yield eq. (2.9); however, performing the linear algebra on such a large system is challenging. Therefore, it is useful to employ the shortcuts described in refs. [36,16]: the first-and second-entry conditions reduce somewhat the size of the initial ansatz, and applying the integrability condition iteratively softens the exponential growth of the ansatz with the weight. Even still, the computation requires a dedicated method, since out-of-the-box linear algebra packages cannot handle such large systems. We implemented a batched Gaussian elimination algorithm, performing the back substitution with FORM [60], similar to the method described in ref. [61].
As discussed in ref. [36], the factorization formula of Fadin and Lipatov [35] in the multi-Regge limit (see section 3.2) provides additional constraints on the 113 parameters entering eq. (2.9), 7. The symbol agrees with the prediction from BFKL factorization at NLL [35].
We may also apply constraints in the near-collinear limit by matching onto the recent predictions by BSV based on the OPE for flux tube excitations [26], 8. The symbol is in agreement to order T 1 with the OPE prediction of the near-collinear expansion [26,27]. 9. The symbol is in agreement to order T 2 with the OPE prediction of the near-collinear expansion [28,57].
The dimension of the ansatz for the symbol after applying each of these constraints successively is summarized in table 1. In this table, we also provide the corresponding numbers at two and three loops, so that one can appreciate the increased computational complexity of the four-loop problem. It is worth noting that some constraints become even stronger when promoted to function-level properties, not only fixing beyond-the symbol terms, but also implying additional relations on the symbol-level parameters. An example of this was already seen at three loops [14] where, ultimately, only a single free parameter remained to be determined by the O(T ) nearcollinear limit [16].
In ref. [16], the last two constraints were applied at function level to fully determine the three-loop remainder function. In fact, we will soon apply them at function level in the four-loop case as well, but first we will apply them at symbol level in order to determine the constants not fixed by the first seven constraints. For this purpose, it is necessary to expand the symbol S(R  in the near-collinear limit v → 0, u + w → 1. It is convenient to adopt the parametrization of ref. [26] in terms of variables F , S, and T , which are related to the u i and y i variables by, (2.10) The near-collinear limit is the limit T → 0 for fixed F and S.

Expanding the symbol in a limit
We wish to expand symbols and functions in a particular kinematic limit, which in the present case is T → 0. To this end, we formulate the expansion of an arbitrary pure function F (T ) in a manner that can easily be extended to the symbol. The function may contain arbitrary dependence on S and F , which is not shown explicitly. The expansion is not entirely trivial because it will in general contain powers of ln T , as well as powers of T , and some care must be given to keep track of them. Let us explicitly separate the power-law behavior from the logarithmic behavior by writing, then we have and so forth. Now consider a pure function F (T ) for which F (0) = [F (T )] 0 = 0. The function can contain powers of ln T in the expansion around T = 0, as long as they are accompanied by positive powers of T so that the limit as T → 0 vanishes. Because symbols provide information about the derivatives of functions in a convenient way, we write F as the integral of its derivative, to leading order in the expansion around T = 0, (2.14) Owing to the presence of logarithms, it is possible that in evaluating [F (T )] 0 we might generate a pole in T . We let where the first term comes from differentiating explicit ln T factors in F (T ). Then we can write the expansion of the integrand in eq. (2.14) as 1 by again applying eq. (2.14), this time with F → f −1 . Therefore eq. (2.14) defines a recursive procedure for extracting the first term in the expansion around T = 0. The recursion will terminate after a finite number of steps for a pure function. The only data necessary to execute this procedure are the ability to evaluate the function when T = 0, and the ability to take derivatives. Since both of these operations carry over to the symbol, we can apply this method directly to S(R (4) 6 ). To give a flavor of how the recursion works, we expand the symbol in the following way, (2.17) where R i = T is defined to have length one and the A i have length 7 − i. There are terms with up to six consecutive T entries in the final slots. Although we have made explicit the T entries at the back end of the symbol, there may be other T entries hidden inside the A i . Applying eq. (2.14), we obtain, 6 ) that have more factors of T on the back end. Notice that the innermost integrals have no 1/T i in the measure, and as such they will generate terms of mixed transcendentality. The mixed transcendentality is not surprising; indeed it is typical whenever one expands a function of uniform transcendentality to subleading order in a given limit. For example, ln The extension of eq. (2.18) gives the expansion of S(R (4) 6 ) to order T 1 . One can easily generalize this method to extract more terms in the T expansion. To obtain the T n term, we first subtract off the expansion through order T n−1 and divide by T n−1 , yielding a function that vanishes when T = 0. Then we can proceed as above and calculate the T 1 term, which will correspond to the T n term of the original function.
Proceeding in this manner, we obtain the expansion of the symbol of R (4) 6 through order T 2 . To compare this expansion with the data from the OPE, we must first disregard all terms containing factors of π or ζ n , since these constants are not captured by the symbol. We must also convert from the remainder function to the logarithm of the specific Wilson loop ratio considered by BSV. Both expressions are finite and dual conformal invariant, but they differ by a simple additive function: 4 where the cusp anomalous dimension is 20) and the function X(u, v, w) is given by where H u 2 = H 0,1 (1 − u) = Li 2 (1 − u) denotes a harmonic polylogarithm (HPL) [52]. The conventional loop expansion parameter for the Wilson loop, g 2 , is related to our expansion parameter by g 2 = a/2.
Performing the comparison in eq. (2.19) at four loops, we find that the information at order T 1 is sufficient to fix all but four of the remaining parameters. At order T 2 , all four of these constants are determined and many additional cross-checks are satisfied. The final expression for the symbol of R

Constraints at function level
We now turn to the problem of promoting the symbol to a function. In principle, the procedure is identical to that described in ref. [16]; indeed, with enough computational power we could construct the full basis of hexagon functions at weight seven (or even eight), and replicate the analysis of ref. [16]. In practice, it is difficult to build the full basis of hexagon functions beyond weight five or six, and so we briefly describe a more efficient procedure that requires only a subset of the full basis.
To begin, we construct a function-level ansatz for the {5, 1, 1, 1} components of the coproduct of R . The ansatz is a four-fold tensor product whose first slot is a weight-five function and whose last three slots are logarithms. The symbols of the weight-five functions can be read off of the symbol of R (4) 6 , by clipping off the last three entries. They can then be identified with functions in the weight-five hexagon basis. Therefore we can immediately write down, 6 ] s i ,s j ,s k are the most general linear combinations of weight-five hexagon functions with the correct symbol and correct parity. There will be many arbitrary parameters, all of which are associated with ζ values multiplying lower-weight functions.
Many of these parameters can be fixed by demanding that s i ∈Su [R 6 ] s i ,s j ,s k be the {5, 1} component of the coproduct for some weight-six function for every choice of j and k. This is simply the integrability constraint, discussed extensively in ref. [16], applied to the first two slots of the four-fold tensor product in eq. (2.22). We also require that each weight-six function has the proper branch cut structure; again, this constraint may be applied using the techniques discussed in ref. [16]. Finally, we must guarantee that the weight-six functions have all of the symmetries exhibited by their symbols. For example, if a particular coproduct entry vanishes at symbol level, we require that it vanish at function level as well. We also demand that the function have definite parity since the symbol-level expressions have this property.
After imposing these mathematical consistency conditions, we will have constructed the {5, 1} component of the coproduct for each of the weight-six functions entering ∆ 6,1,1 (R (4) 6 ), as well as all the integration constants necessary to define corresponding integral representations (see section 4 of ref. [16]). There are many undetermined parameters, but they all correspond to ζ values multiplying lower-weight hexagon functions, so they cannot be fixed at this stage.
It is also also straightforward to represent ∆ 6,1,1 (R (4) 6 ) directly in terms of multiple polylogarithms in a particular subspace of the Euclidean region, called Region I in ref. [16]: The fact that the y i are all real and between 0 and 1 facilitates a representation in terms of multiple polylogarithms, as discussed in ref. [16]. This region is also of interest because it corresponds to positive external kinematical data in (2, 2) signature.
To this end, we now describe how to integrate directly the {n − 1, 1} component of the coproduct of a weight-n function in terms of multiple polylogarithms. The method [7,8,63] is very similar to the integral given in eq. (3.8) of ref. [16], which maps symbols directly into multiple polylogarithms. Instead of starting from the symbol, we start from the {n − 1, 1} coproduct component, and therefore we only have to perform one integration, corresponding to the final iteration of the n-fold iterated integration in eq. (3.8) of ref. [16]. As discussed in ref. [16], we are free to integrate along a contour that goes from the origin t i = 0 to the point t i = y i sequentially along the directions t u , t v and t w . The integration is over ω = d ln φ with φ ∈ S y , where S y is the set of 10 letters in the y i variables [16]. The integrand is a combination of weight-(n − 1) multiple polylogarithms in Region I. Together, these two facts imply that the integral may always be evaluated trivially by invoking the recursive definition of multiple polylogarithms, 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 p ; z) = ln p z p! .
(2.24) (Many of the properties of multiple polylogarithms are reviewed in appendix A of ref. [16].) Applying this method to the case at hand, we obtain an expression for ∆ 6,1,1 (R 6 ) in terms of multiple polylogarithms in Region I. Again, we enforce mathematical consistency by requiring integrability in the first two slots, proper branch cut locations, and well-defined parity. We then integrate the expression using the same method, yielding an expression for ∆ 7,1 (R (4) 6 ). Finally, we iterate the procedure once more and obtain a representation for R (4) 6 itself. At each stage we keep track of all the undetermined parameters. Any parameter that survives all the way to the weight-eight ansatz for R (4) 6 must be associated with a ζ value multiplying a lower-weight hexagon function with the proper symmetries, branch-cut locations, and the function-level analog of the final-entry condition. There are 68 such functions. The counting of parameters is presented in table 2.
It is straightforward to expand our 68-parameter ansatz for R 6 in the near-collinear limit. Indeed, the methods discussed in ref. [16] can be applied directly to this case. We carried out this expansion through order T 3 , though even at order T 1 the result is too lengthy to present here. The expansion is available in a computer-readable format from [62].
Demanding that our ansatz vanish in the strict collinear limit fixes all but ten of the beyondthe-symbol constants. Consistency with the OPE at order T 1 , corresponding to contributions of single (gluonic) flux-tube excitation, fixes nine of the ten remaining constants. The final constant k MZVs of weight k Functions of weight 8 − k Total parameters 2 6 after imposing all mathematical consistency conditions. is fixed at order T 2 , corresponding to double flux-tube excitations, as well as twist-two boundstate contributions [28]. The rest of the data at order T 2 provides many nontrivial consistency checks of the result.
In slightly more detail, we can characterize the contributions at a given order in the T expansion by their dependence on F , or equivalently on an azimuthal angle φ, introduced by letting F = e iφ . As discussed in ref. [28], the F dependence is correlated with the helicity of the excitations. The order T 1 term in the near-collinear expansion of the L-loop remainder function always has the form, where c k (S) is a linear combination of HPLs [52] of the form H m (−S 2 ), m i ∈ {0, 1}, multiplied by simple rational functions of S. The weight of the HPLs is at most 2L − k, but can be lower, in accordance with the mixed transcendentality of the T expansion mentioned earlier.
The single powers of F and F −1 correspond to the helicity ±1 gluonic excitations, which have equal contributions due to parity. The expansion of the three-loop remainder function at this order was given explicitly in ref. [16]. At order T 2 , the expansion has the form, where d k (S) and f k (S), like c k (S), are linear combinations of HPLs multiplied by rational functions of S (more complicated ones than appear in c k (S). The terms in eq. (2.26) that have the F ±2 prefactors come entirely from gluonic excitations -either pairs of single excitations, or the contribution of a twist-two gluonic bound-state, either of which can have helicity ±2; whereas the T 2 F 0 terms can come from excitations of pairs of gluons, fermions or scalars [28]. All of the constraints at order T 2 that were needed to fix the five parameters at that stage (four parameters at symbol level and one beyond-the-symbol) came from matching the T 2 F 2 contributions. Hence the comparison of the T 2 F 0 terms, which tests the scalar and fermion contributions as well as gluonic ones, was completely rigid, with no free parameters.
In practice the comparison to the OPE predictions was done after expanding the functions of S in an expansion around S = 0. For the T 2 F 2 comparison we matched the terms through S 20 ; for the T 2 F 0 comparison, through S 10 . Certainly higher orders could be matched if desired; on the OPE side this just amounts to evaluating more residues in the complex rapidity plane [26,27,28]. In some cases one can also perform the residue sums to all orders, see e.g. ref. [64].
The final expression for R (4) 6 in terms of multiple polylogarithms in Region I is available from [62] in a computer-readable format. We also provide a coproduct-based description of it; see appendix A.

.1 Fixing constants at four loops
In the limit of multi-Regge kinematics (MRK), the cross ratios u 1 , u 2 and u 3 approach the values held fixed 5 . While the remainder function vanishes in Euclidean MRK, this is no longer the case once it is analytically continued to a different Riemann sheet, according to u 1 → e −2πi |u 1 | [29].
On this Riemann sheet we can write, The LLA series of coefficients has n = L − 1. The coefficients h L−1 (w, w * ) are known to all orders in perturbation theory [36,65]. At NLLA (n = L − 2), results for the coefficients g L−2 (w, w * ) have been given up to nine loops [33,35,14,36].
At NNLLA (n = L − 3), only the three-loop coefficients are known [33,14,16]. In ref. [36], the four-loop coefficients at NNLLA and N 3 LLA, g 1 (w, w * ) and g (4) 0 (w, w * ), respectively, were heavily constrained and their functional form was completely determined, up to 27 rational numbers a i , b j , i ∈ {0, . . . , 8}, j ∈ {1, . . . , 18}. As mentioned in the introduction, a 0 is a parameter that enters at the level of the symbol. The remaining 26 parameters are beyond-thesymbol; they appear with ζ values multiplying them. Since we have now a complete and unique analytic expression for the four-loop remainder function in general kinematics, the coefficients g (4) 1 and g (4) 0 can be extracted by using the techniques described in ref. [16]. Appendix A gives a brief description of how the coproduct representation of R (4) 6 may be used for this purpose. In this way, we find expressions for the two previously-undetermined MRK coefficients at four loops, and g (4) 0 (w, w * ) = − (3.5) The functions L ± m appearing in these expressions are single-valued harmonic polylogarithms (SVHPLs) [59]. They appear in a basis defined in ref. [36], which diagonalizes the Z 2 × Z 2 action of inversion and conjugation of the variables (w, w * ).
The expressions above match with those of eqs. (7.14) and (7.15) of ref. [36], provided that the constants defined in that reference take the values, 2 (w, w * ) , K are the L-loop coefficients of the cusp anomalous dimension defined in eq. (2.20), and the lower loop g coefficients are given in ref. [36].

The NNLL BFKL eigenvalue and N 3 LL impact factor
The functions g 1 (w, w * ) and g (4) 0 (w, w * ), in turn, determine the NNLLA BFKL eigenvalue and N 3 LLA impact factor, through a master equation [35], where recalling that L − 0 = ln|w| 2 and L + 1 = 1 2 ln(|w| 2 /|1 + w| 4 ). The BFKL eigenvalue ω(ν, n) and the impact factor Φ Reg (ν, n) can be expanded perturbatively, Reg (ν, n) + O(a 4 ) . (3.14) We remark that an alternate version of the master equation has recently been found in ref. [58]. In contrast to eq. (3.11), the denominator ν 2 + n 2 /4 contains an additional term proportional to the square of the cusp anomalous dimension. It also lacks the explicit Regge pole contribution (the cos πω ab term), although this contribution can be recovered by evaluating the n = 0 term and ν = 0 residue in the integral at finite coupling. Then the two factorization forms become equivalent, up to a different definition of the impact factor. In this paper, we will continue to use the form (3.11).
The first two nontrivial orders in the expansion of the BFKL eigenvalue and the impact factor were known previously [35,30,32,16,36], where ψ(z) = d dz ln Γ(z) is the digamma function, ψ(1) = −γ E is the Euler-Mascheroni constant, and V and N are given by, with D ν ≡ −i∂ ν ≡ −i ∂/∂ν. After expanding the master equation (3.11) to the relevant order in a and ln(1 − u), one has to match the resulting combinations of SVHPLs in (w, w * ) against the inverse Fourier-Mellin transforms of suitable functions of ν and n. This was carried out in ref. [36], in terms of the then-undetermined a i and b i constants. Inserting the values (3.6) and (3.7) into the respective expressions, we obtain, ν,n (when the latter becomes available), will permit an evaluation at N 3 LLA -assuming that the factorization continues to hold at this order.
In ref. [36] it was observed that E ν,n in eq. (3.20) has a nonvanishing limit ν → 0 (after setting n = 0), even though E ν,n and E ν,n vanish in this limit [35]. This limit of E ν,n held independently of all the constants in eqs. (3.6) and (3.7), which were unknown at that time. The reason it was independent of the constants was that the four-loop remainder function was required to vanish in the collinear corner of the MRK limit, |w| 2 → 0. This limit in the (w, w * ) plane in turn controls the n = 0, ν → 0 limit of the BFKL eigenvalue ω(ν, n). In ref. [58], the general constraints imposed by collinear triviality of the remainder function were derived at finite coupling, and eq. (3.22) was obtained as a byproduct.

NNLL coefficient functions at five loops
The MRK factorization implicit in the master equation lets us bootstrap higher-loop coefficients in the MRK limit. We simply insert the results for the BFKL eigenvalue and the impact factor through NNLLA into the master equation (3.11). We then use the techniques of ref. [36] to perform the inverse Fourier-Mellin transform from (ν, n) space back to (w, w * ) space. This transform is facilitated by having a complete basis of SVHPLs at the appropriate transcendental weight. The inverse Fourier-Mellin transform leads to double sums, which can either be summed explicitly, or truncated and then matched to a Taylor expansion of the SVHPL basis. In this way we can obtain explicit expressions for R 6 in MRK at NNLLA, just as was done at LLA and NLLA in ref. [36]. These data will be important in order to help constrain the functional form of the remainder function at higher loop orders.
As an example, we present here the result for the five-loop six-point remainder function at NNLLA. For the imaginary part, we find, g (3.23) The corresponding results at LLA and NLLA were given in ref. [36].

(3.24)
Finally, we give the N 3 LL real-part coefficient, which is related to eq. (3.23) and to imaginary parts at lower logarithmic or lower loop orders, h 1 (w, w * ) 1 (w, w * ) .

(3.25)
Although the real parts are related by analyticity to the imaginary parts, they still provide useful additional constraints on ansätze for the remainder function.

Connection between BFKL and the flux tube spectrum?
We conclude this section by noting that the result for the BFKL eigenvalue at NNLLA suggests an intriguing connection between the BFKL eigenvalues E ν,n , E ν,n , and E ν,n and the weak-coupling expansion of the energy E(u) of a gluonic excitation of the GKP string as a function of its rapidity u, given in ref. [25]. First we rewrite the expressions for E ν,n , E (1) ν,n , and E (2) ν,n explicitly in terms of ψ functions and their derivatives, where ξ ± ≡ 1 ± iν + |n| 2 . Next, we keep only the pure ψ (and ζ) terms, dropping anything with a V or an N , ν,n ψ only = 1 8 Now we compare this formula to equation (4.21) of ref. [25] for the energy E(u) of a gauge field ( = 1) and its bound state ( > 1), where g 2 = a/2 is the loop expansion parameter, s = 1 + /2, and Neglecting the constant offset at order a 0 (the classical operator scaling dimension), eq. (3.29) matches perfectly with eq. (3.28) at order a 1 and a 2 , provided that we identify, The correspondence continues to order a 3 if we also drop the term 24 ζ 3 ψ (+) 1 (s−1, u). It would be very interesting to understand the origin of this correspondence, and whether there is a physical meaning to the operation of dropping all terms with a N or a V . We leave this question to future work and return our attention to the quantitative behavior of the four-loop remainder function.

Quantitative behavior
In this section we investigate the quantitative behavior of the four-loop remainder function in the Euclidean region where all three cross ratios are positive. It will prove particularly instructive to plot the ratios of successive loop orders, R 6 . It was observed in ref. [16] that the former ratio was quite stable along large portions of a line and a two-dimensional surface where it was examined. We will find that the stability of such ratios extends to four loops, i.e. to the latter ratio, and to a number of different lines and one two-dimensional surface, as long as the cross ratios are not too large or too small. We will also examine certain limiting behavior analytically, where it can sometimes shed light on the remarkable stability of the ratios. Finally, we will discuss how perturbation theory is doing with respect to the approach to large orders.

Region I
While the full function R (4) 6 is too lengthy to be shown here, its representation in terms of multiple polylogarithms can easily be evaluated numerically in Region I, defined in eq. (2.23), using GiNaC [66,67]. In table 3, we show the value of the four-loop remainder function for five reference points. In addition, in fig. 1 we plot the ratio R  Table 3: Numerical evaluation of the four-loop remainder function at a selection of points in Region I.
For u = v, the boundary of Region I in the interior of the Euclidean region is defined by ∆(u, u, w) = 0, where ∆ is given in eq. (2.8); this parabola w = (1 − 2u) 2 is shown as the red line in the plot. We restrict the plot to stay slightly away from the boundaries of the Euclidean region, taking u, w > 0.06. At these boundaries, R 6 (u, u, w) diverges logarithmically, order by order in perturbation theory, whenever one of the cross ratios becomes very small and the other one is kept finite. At a given loop order, the degree of the logarithmic divergence is one power lower when w → 0 with u fixed, than it is for the opposite case when u → 0 with w fixed:  are quite lengthy, so we do not show them here. We list the results for the coefficients of just the leading logarithmic divergence up to four loops in appendix B. Because the leading logarithm increases by one with each additional loop, the ratios plotted in fig. 1 diverge like a single logarithm as either boundary is approached. However, the leading logarithms in the numerator and denominator of the ratio are far from dominant at the boundaries of the plot where u or w = 0.06. If one keeps all subleading logarithms, and neglects the power-suppressed terms, one gets quite close to the exact numerical value of the ratio at either boundary of the plot.
It was recently conjectured [68] that the remainder function should have a uniform sign in Region I, which corresponds to the kinematic regime of positive external momentum twistor kinematics. Recent formulations of the planar scattering amplitude loop integrand [69,12,70] lead to manifestly positive integrands in this region. On the other hand, an infinite subtraction is required to pass to the remainder function. Nevertheless, it was observed that this conjecture indeed holds at two loops [68] and also at three loops [16]. Given that R (3) 6 is negative in Region I [16], it is obvious from fig. 1 that R (4) 6 has a uniform (positive) sign in Region I, at least on the surface u = v. In total we checked more than 1000 points in Region I, both on and off the u = v surface; for all points checked, the value of R (4) 6 is positive, in agreement with the conjecture. In the rest of this section we focus on the remainder function restricted to certain onedimensional subspaces where the functional form simplifies drastically. These lines may prove useful in trying to find a form for the remainder function that is valid to all loop orders, i.e. at finite coupling, beyond what is presently known in the OPE limit [26,27,28]. The first line we discuss has one endpoint which intersects the OPE limit. Perhaps this proximity could allow the knowledge of the OPE limit to anchor such a finite-coupling construction. The other two lines never approach the OPE limit, although they have other interesting properties.

The line (u, u, 1)
As noted in ref. [16], the two-and three-loop remainder functions can be expressed solely in terms of harmonic polylogarithms (HPLs) of a single argument, 1 − u, on the line (u, u, 1), and we use the notation H u m ≡ H m (1 − u). The same is true at four loops, although the resulting expression is rather lengthy. It can be obtained by taking the limit of the general coproduct representation described in appendix A onto the line (u, u, 1). The quantity ∆ defined in eq. (2.8) vanishes on this line. As a consequence, all parity-odd functions vanish on the line too. The derivatives of weight n parity-even functions can be expressed using eq. (A.2) in terms of parity-even and parity-odd coproduct components of weight n − 1. The vanishing of the parity-odd functions as one approaches the line is fast enough that they can be neglected in computing the derivative along the line. Then one obtains from eq. (A.2), which is easily integrated in terms of the functions H u m , given that the coproduct components F u , F v , F 1−u and F 1−v are also expressible in this form.
Because the four-loop expression is still rather lengthy, in order to save space we first expand all products of HPLs using the shuffle algebra. The resulting "linearized" representation will have HPL weight vectors m consisting entirely of 0's and 1's, which we can interpret as binary numbers. Finally, we can write these binary numbers in decimal, making sure to keep track of the length of the original weight vector, which we write as a superscript. For example, 11 . (4.3) In this notation, R 6 (u, u, 1) and R 6 (u, u, 1) read, 6 (u, u, 1) = −3h [6] 1 + 5h [ − 5h [8] 159 + 6h [8] 161 − 11h [8] 163 − 3h [8] 165 + 3h [8] 167 − 4h [8] 171 − 4h [8]  + 9h [8] 193 − 25h [8] 195 − 9h [8] 197 + 27h [8] 199 − 2h [8] 201 + 9h [8] 203 + 2h [8] 205 − 23h [8] 207 + 2h [8] 209 − h [8] 213 − 8h [8] 215 + 2h [8] 217 − 3h The remainder function R 6 (u, v, w), as a function of three variables, satisfies a differential constraint, corresponding to the final-entry condition imposed on the symbol. As discussed in appendix A, this means that the {7, 1} components of the coproduct obey R This property of the partial derivatives does not necessarily extend to the ordinary derivatives along a generic line. However, from eq. (4.2) it is easy to see that it must hold along the line (u, u, 1), where it implies that dR (L) It is easy to check that the property (4.7) holds for the expressions for R   In fact, the ratios are also similar away from this point, as can be seen in fig. 2. The logarithmic scale for u highlights how little the ratios vary over a broad range in u, as well as how the u-dependence differs minimally between the successive ratios.
We also give the leading term in the expansion of R  We note the intriguing observation that the maximum-transcendentality piece of the u 1 ln 0 u term is proportional to the four-loop cusp anomalous dimension, 219 K . In fact, the corresponding pieces of the two-and three-loop results, given in ref. [16], can be checked to similarly correspond to − 1 4 γ K and − 1 4 γ K . In the limit u → 0, the line (u, u, 1) touches the end of the collinear line v = 0, u + w = 1. So one could ask where the cusp anomalous dimension seen in eq. (4.11) originates in the nearcollinear limit from the OPE perspective. Actually, it is not there at all in the limiting behavior S → 0, T → 0 of the Wilson loop ratio W hex employed in refs. [26,27,28]. To see this, first recall from eq. (2.10) that to leading order in T , u = S 2 /(1 + S 2 ), v = T 2 , and w = 1/(1 + S 2 ). Hence the line (u, u, 1) for u → 0 matches the S → 0, T → 0 limit, after making the identification u = S 2 = T 2 , to leading order. Now let's inspect the additive term 1 8 γ K (a) X(u, v, w) in eq. (2.19) relating R 6 to ln W hex . The function X(u, v, w) defined in eq. (2.21) is suppressed by a power of u in this limit, X(u, u, 1) = 2 u + O(u 2 ), (4.12) as u → 0. This limiting behavior has the precise form and value to cancel the − 1 4 γ K (a) · u in R 6 (u, u, 1) in passing to ln W hex (a/2) via eq. (2.19).
Suppose, however, that we look at the other end of the collinear line v = 0, u + w = 1; namely the line (1, u, u) as u → 0. This line matches the S → ∞, T → 0 limit, with the identification u = 1/S 2 = T 2 to leading order. The S 3 permutation symmetry of the remainder function implies that R 6 (1, u, u) = R 6 (u, u, 1). However, the function X has a different behavior in this limit, (4.13) The logarithmic term implies that in the S → ∞, T → 0 limit the cusp anomalous dimension is visible in the OPE. The difference between the two limits (or more generally, the lack of symmetry of the Wilson loop ratio) is related to changing the "framing" of the hexagonal Wilson loop, by making the other possible choice of pentagons and box to remove the ultraviolet divergences. This change of frame always involves the cusp anomalous dimension [28]. It may be useful to study the limiting behavior of the (u, u, 1) and (1, u, u) lines in more detail, as an avenue along which the OPE might potentially be resummable at finite coupling. Comparing eq. (4.11) with the corresponding results for R 6 and R 6 [16], we see that the ratios R (L) both diverge logarithmically as u → 0 along this line: 6 (u, u, 1) R 6 (u, u, 1) (4.14) The slight difference in these coefficients is reflected in the slight difference in slopes in the region of small u in fig. 2.
As u → ∞, the leading behavior at four loops is, (4.16) Just like at two and three loops, R 6 (u, u, 1) approaches a constant as u → ∞. Comparing with eq. (7.17) of ref. [16], we find

The line (u, 1, 1)
Next we consider the line (u, 1, 1), which, due to the total S 3 symmetry of R 6 (u, v, w), is equivalent to the line (1, 1, w) discussed in ref. [16]. As was the case at two and three loops, we can express R Using eq. (4.8), it is easy to check that none of these functions satisfies a property like eq. (4.7), where the derivative is expressed in terms of a single pure function multiplied by a rational prefactor. The reason is related to the nonvanishing contributions of the parity-odd functions in the coproduct representation. At both large and small u, these functions all diverge logarithmically. At two and three loops, this was observed in ref. [16]. At four loops, we find at small u, R (4) 6 (u, 1, 1) = 1 24 6 (u, 1, 1) R 6 (u, 1, 1) 6 (u, 1, 1) R 6 (u, 1, 1) 6 (u, 1, 1) R 6 (u, 1, 1) 6 (u, 1, 1) R

The line (u, u, u)
At strong coupling, using the AdS/CFT correspondence, gluon scattering amplitudes can be computed in the semi-classical approximation by minimizing the area of a string world-sheet propagating in AdS 5 × S 5 [2]. The world-sheet boundary conditions depend on the scattering kinematics. The amplitude has the generic form, 6 where λ = g 2 YM N c = 8π 2 a. As discussed in refs. [72,16], on the symmetrical diagonal line (u, u, u), the remainder function at strong coupling can be written analytically, where φ = 3 cos −1 (1/ √ 4u). The simplicity of this formula motivates us to evaluate the four-loop remainder function on the line (u, u, u), as we did earlier at two and three loops [16].
In perturbation theory, the function R (L) 6 (u, u, u) cannot be written solely in terms of HPLs with argument (1−u). However, it is possible to use the coproduct structure to derive differential equations which may be solved by using series expansions around the three points u = 0, u = 1, and u = ∞. This method was applied in ref. [16] at two and three loops, and here we extend it to the four-loop case.
The expansion around u = 0 takes the form, The leading term at four loops diverges logarithmically, but, just like at two and three loops, the divergence appears only as ln 2 u. This is another piece of evidence in support of the claim by Alday, Gaiotto and Maldacena [72] that this property should hold to all orders in perturbation theory. Because of this fact, the ratios R 6 (u, u, u) R The ratios R 6 (u, u, u)/R 6 (u, u, u) and R 6 (u, u, u)/R 6 (u, u, u) approach constants in the limit u → ∞,  In contrast to the expansions around u = 0 and u = ∞, the expansion around u = 1 is regular, (4.31) We take 100 terms in each expansion, around 0, 1 and ∞ and piece them together to obtain a numerical representation for the function R (4) 6 (u, u, u) that is valid along the entire line. In the regions of overlap, we find agreement to at least 15 digits. In fig. 4, we plot the ratios R (L) 6 (u, u, u)/R (L−1) 6 (u, u, u) for a large range of u. The spike in the plot is not a numerical instability; it occurs because the denominators in the respective ratios go through zero at a slightly different point from the numerators, around u = 1/3.
As noted in ref. [16], the two-and three-loop remainder functions vanish along the line (u, u, u), very close to the point u = 1/3. More precisely, it was found that the vanishing relation for two and three loops, respectively. The point (u, v, w) = (1/3, 1/3, 1/3) is special because it is where the line (u, u, u) pierces the plane u + v + w = 1. This plane passes through all three of the lines marking the collinear limits (v = 0, u + w = 1; and cyclic permutations thereof). Because R 6 (u, v, w) vanishes on all three lines, one might expect it to vanish close to the equilateral triangle that is bounded by them, which lies in the plane u + v + w = 1. Indeed, that is what is seen at three loops [16]. In this paper, we will not evaluate the four-loop remainder function on this triangle, but we can verify that the zero-crossing point remains close to u = 1/3. The precise zero-crossing value at four loops is u With respect to the three-loop value in eq. (4.32), the zero-crossing point has shifted slightly further away from u = 1/3. As can be seen from fig. 4, R 6 (u, u, u) actually crosses zero in a second place, at a very large value of u,ũ This phenomenon does not happen at two or three loops: R 6 (u, u, u) and R 6 (u, u, u) have unique zero crossings, at the values given in eq. (4.32). Aside from the zero-crossing neighborhood,  4 shows excellent agreement between the two successive ratios for relatively small u, say u < 1000. For large u, the ratios approach constant values that differ by a factor of about −17.6 (see eq. (4.30)).
In fig. 5, we plot the two-, three-, and four-loop and strong-coupling remainder functions on the line (u, u, u). In order to compare their relative shapes, we rescale each function by its value at (1, 1, 1). The remarkable similarity in shape that was noticed at two loops [73] 7 and at three loops [16] clearly persists at four loops, particularly for the region 0 < u < 1.
As discussed in ref. [16], 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. Comparing eq. (4.28) to eq. (4.10), we find, R 6 (u, u, u) R  These ratios are indeed quite close to 1, despite their complicated representations in terms of ζ values. The agreement is slightly better for the double ratio between four and three loops, than it is for the one between three and two loops. We can also compute similar double ratios involving the perturbative and strong coupling coefficients, R (∞) 6 (u, u, u) R   The ratio between the two-loop and strong-coupling points is exactly 1, while the corresponding ratios for three and four loops deviate slightly from one. The deviations increase as L increases, suggesting that the shapes of the weak-coupling curves on the line (u, u, u) are getting slightly further from the shape of the strong coupling curve, at least for small L. This observation is also evident in fig. 5 at large u.
Let us conclude this section by making a comment on hexagon functions on the line (u, u, u). It is easy to check that on this line we have 38) and the symbol of R 6 (u, u, u) has all its entries drawn form the set {y, Φ 2 (y), Φ 3 (y)}, where Φ 2 (y) = 1 + y and Φ 3 (y) = 1 + y + y 2 (4.39) denote the second and third cyclotomic polynomials. It follows then that R 6 (u, u, u) can be entirely expressed through iterated integrals over d ln forms with cyclotomic polynomials as arguments. This class of iterated integrals is a generalization of HPLs, called cyclotomic HPLs, and was studied in detail in ref. [77]. Note that this observation only follows from the entries in the symbol, and is by no means restricted to four loops. As a consequence, we conclude that on the line (u, u, u) hexagon functions, and thus the six-point remainder function, can always be expressed in terms of cyclotomic HPLs.

Approach to large orders
In the previous subsections, we have found that the portion of the line (u, u, u) with 0 < u < 1 leads to quite constant ratios of successive loop orders L. We can also ask what this ratio should become as L → ∞. Most quantum field theories have a zero radius of convergence for their perturbative expansions; that is, the series are asymptotic. There are two generic reasons for this: renormalons and instantons, each of which leads to factorial growth of perturbative coefficients. However, planar N = 4 super-Yang-Mills theory is free from both of these phenomena. Because it is conformally invariant, the beta function vanishes and there are no renormalons. Because the number of colors N c is very large, at fixed 't Hooft coupling λ instantons are exponentially suppressed as N c → ∞ by a factor of exp(−8π 2 /g 2 YM ) = exp(−8π 2 N c /λ). Hence we should expect the perturbative expansion to have a finite radius of convergence r. The ratio r corresponds to a growth rate of successive perturbative coefficients c (L) , which approaches a constant as L becomes large, In eq. (4.40) we have assumed an alternating series, which holds for R (L) 6 for L = 2, 3, 4 throughout Region I and on the lines (u, u, 1) and (u, 1, 1), and for L = 2, 3 throughout almost all of the unit cube 8 .
There is another quantity, closely related to the scattering amplitude, which we can use as a simple benchmark for assessing large order behavior. That quantity is the cusp anomalous dimension. Its perturbative expansion can be computed to all orders using the exact formula of Beisert, Eden and Staudacher [23]. Using this formula, we give the ratio of successive loop orders in table 4. At very large loop orders, the ratio approaches −8, corresponding to a radius of convergence of 1/8 when using the loop expansion parameter a. (In terms of g 2 = a/2, the radius of convergence is 1/16.) However, the approach to this asymptotic value is quite slow. Table 4 also shows the two nontrivial ratios currently available for the remainder function at (u, v, w) = (1, 1, 1), as representative of the fairly constant region (u, v, w) = (u, u, u) with 0 < u 1. We also give values for the three available ratios for the Wilson loop ratio evaluated at two interior points, u = 1 4 and u = 3 4 . (The Wilson loop ratio diverges at (u, v, w) = (1, 1, 1).) There is an extra ratio available for the Wilson loop because its one-loop value is nonzero, due to the function X(u, v, w) appearing in eq. (2.19).
An optimist would say that the remainder-function ratios exhibit a precocious approach to the expected asymptotic value of −8: the cusp anomalous dimension ratio does not reach −6.5 until eight loops. A pessimist would say that the trend is the wrong way: the ratio for L = 4 is further from −8 than is the ratio for L = 3. On the other hand, the Wilson loop ratios are actually approaching −8 monotonically. For both u = 1 4 and u = 3 4 , they appear to be converging more quickly to −8 than is the cusp anomalous dimension.
It is worth remarking that in Region I for u = v, the region shown in fig. 1, the ratio R lies between −6.6 and −7 over the entire region shown. More generally, sampling 1352 points 8 As noted in ref. [16], there is a small region surrounding the plane u + v + w = 1 in which R 6 and R     in Region I, including ones with u = v, the ratio is always between −6.60 and −8.67. Clearly a computation of the remainder-function ratio at the next loop order, R 6 , would be very illuminating in this regard.

Conclusions
In this article, we presented the four-loop remainder function, which is a dual-conformally invariant function that describes six-point MHV scattering amplitudes in planar N = 4 super Yang-Mills theory. The result was bootstrapped from a limited set of assumptions about the analytic properties of the relevant function space. Following the strategy of ref. [14], we constructed an ansatz for the symbol and constrained this ansatz using various physical and mathematical consistency conditions. A unique expression for the symbol was obtained by applying information from the near-collinear expansion, as generated by the OPE for flux tube excitations [26]. The symbol, in turn, was lifted to a full function, using the methods described in ref. [16]. In particular, a mathematically-consistent ansatz for the function was obtained by applying the coproduct bootstrap described in ref. [16]. All of the function-level parameters of this ansatz were fixed by again applying information from the near-collinear expansion.
The final expression for the four-loop remainder function is quite lengthy, but its functional form simplifies dramatically on various one-dimensional lines in the three-dimensional space of cross ratios. While the analytic form for the function on these lines is rather different at two, three, and four loops, a numerical evaluation shows that they are in fact quite similar for large portions of the parameter space, at least up to an overall rescaling. On the line where all three cross ratios are equal, an analytical result at strong coupling is available. The perturbative coefficients are very similar in shape to the strong-coupling one, particularly in the region where the common cross ratio is less than one. This agreement suggests that an interpolation from weak to strong coupling may depend rather weakly on the kinematic variables, at least on this one-dimensional line.
Given the full functional form of the four-loop remainder function, it is straightforward to extract its limit in multi-Regge kinematics. This information allowed us to fix all of the previously undetermined constants in the NNLLA BFKL eigenvalue and the N 3 LLA impact factor. Although we used some multi-Regge factorization information as input, the fact that we found a solution consistent with all the OPE data suggests that factorization does hold beyond NLLA. We also observed an intriguing correspondence between the BFKL eigenvalue and the energy of a gluonic excitation of the GKP string. It would be very interesting to better understand this correspondence.
There are many avenues for future research. For example, it would be interesting to try to understand the correspondence between the integrated results found here (and at three loops) and the types of multi-loop integrals that appear in recent formulations of the planar multi-loop integrand [69,12,70].
In implementing the kind of bootstrap used here beyond the six-point case, it is important to have a good understanding of the relevant space of functions from results at low loop order. Progress is being made on this front [53,78], most recently through the introduction of cluster coordinates [79] and cluster polylogarithms [80].
In principle, the methods used in this work could be extended to five loops and beyond. The primary limitation is computational power and the availability of boundary data, such as the near-collinear limit, to fix the proliferation of constants. It is remarkable that a fully nonperturbative formulation of the near-collinear limit now exists. Ultimately, the hope is that the full analytic structure of perturbative scattering amplitudes, as exposed here through four loops for the six-point case, might in some way pave the way for a nonperturbative formulation for generic kinematics.

A Sample coproducts
As mentioned in section 2.4, the construction of a complete set of hexagon functions at weight five [16] facilitated the construction of R (4) 6 at function level in the present paper. We could identify the symbols of all the coefficients [R (4) 6 ] s i ,s j ,s k of the {5, 1, 1, 1} coproduct, ∆ 5,1,1,1 (R (4) 6 ) in eq. (2.22), with linear combinations of the functions constituting the weight-five basis, modulo the ζ value ambiguities listed in table 2. Besides facilitating the construction, writing the {5, 1, 1, 1} coproduct elements in terms of weight-five hexagon functions also provides a compact way to define the final answer for R (4) 6 . Essentially we are specifying the function via its derivatives. In this appendix, we will list a few of the coproduct elements of R (4) 6 to give a flavor for this description, although they are still too lengthy to list all of them here. We will provide the complete set as a computer-readable file.
First, though, we briefly review the connection between the coproduct and derivatives of hexagon functions [16]. A hexagon function F of weight n has a {n − 1, 1} coproduct component of the form, where the nine functions {F u i , F 1−u i , F y i } are of weight n−1. The first derivatives of F , in either the u i variables or the y i variables, are simple linear combinations of these coproduct elements: Derivatives with respect to v, w, y v and y w can be obtained from the cyclic images of eq. (A.2). As discussed extensively in ref. [16], the derivatives can be used to define various integral representations for F , which can be evaluated numerically. It is also possible to integrate the differential equations analytically in various kinematical limits. For example, in the MRK limit, the appropriate variables are (ξ, w, w * ), where ξ ≡ 1 − u is vanishing and (w, w * ) are defined via eq. (3.2). The differential equations in the MRK variables are [16], Using eq. (A.2) and its cyclic images, we find that the w derivative can be rewritten directly in terms of the coproduct elements as, This differential equation can be integrated up systematically in terms of SVHPLs. The MRK limiting behavior of all the weight-five hexagon functions was given in ref. [16]. These results give directly the MRK limits of all the independent elements ∆ 5,1,1,1 (R (4) 6 ). Then we can integrate up eq. (A.4) in order to get the MRK behavior of all the ∆ 6,1,1 (R (4) 6 ) elements, integrate once more to get the limiting behavior of the ∆ 7,1 (R (4) 6 ) elements, and integrate a final time to get the desired MRK behavior of R (4) 6 itself. How many coproduct components have to be specified? Thanks to the S 3 permutation symmetry of R (4) 6 (u, v, w) and the differential constraint corresponding to the final-entry condition, the number is manageable. First of all, there are only two independent {7, 1} coproduct elements, R u and R yu , (A.5) where we have suppressed the subscript 6 and superscript (4) to avoid clutter in subsequent equations. The final-entry constraint becomes for the coproduct. The S 3 symmetry implies that the other elements can be obtained by permuting the two elements given in eq. (A.5), R v (u, v, w) = R u (v, w, u), R w (u, v, w) = R u (w, u, v), R yv (u, v, w) = R yu (v, w, u), R yw (u, v, w) = R yu (w, u, v). (A.7) There are 11 independent {6, 1, 1} coproduct elements: R u,u , R 1−u,u , R yu,u = R u,yu , R 1−u,yu , R yu,yu , R v,u , R 1−v,u , R yv,u , R v,yu , R 1−v,yu , R yv,yu . (A.8) The counting is as follows: Using the cyclic symmetry, the last entry can be rotated to be u, 1 − u or y u . However, the final-entry condition at function level eq. (A.6) says that a last entry of 1 − u can be exchanged for a last entry of u, at the price of a minus sign. There is still a residual flip symmetry, exchanging v ↔ w, which allows the next-to-last entry to be forbidden from being w, 1 − w or y w . That counting leaves 12 possibilities; however, we also find that R yu,u = R u,yu , which presumably follows from integrability. Here we will give the {5, 1, 1, 1} coproduct elements that allow the construction of R u,u . In fact, the {5, 1, 1, 1} coproduct entries allow us to construct the total derivative of R u,u , so we need to supplement them with a constant of integration, which we specify at the point (u, v, w) = (1, 1, 1): Now, using the residual v ↔ w flip symmetry for R u,u , the six independent elements required to specify R u,u are: R u,u,u , R 1−u,u,u , R v,u,u , R 1−v,u,u , R yu,u,u and R yv,u,u . The parity-odd elements R yu,u,u and R yv,u,u are much simpler to represent, because the basis of weight-five parity-odd functions is much smaller than the parity-even basis. They are given by, R yu,u,u = 1 128 −3 H 1 (u, v, w) + H 1 (v, w, u) + H 1 (w, u, v) + 1 4 11 [J 1 (u, v, w) + J 1 (v, w, u)] The four parity-even elements are given by,