Hexagon Wilson Loop OPE and Harmonic Polylogarithms

A recent, integrability-based conjecture in the framework of the Wilson loop OPE for N=4 SYM theory, predicts the leading OPE contribution for the hexagon MHV remainder function and NMHV ratio function to all loops, in integral form. We prove that these integrals evaluate to a particular basis of harmonic polylogarithms, at any order in the weak coupling expansion. The proof constitutes an algorithm for the direct computation of the integrals, which we employ in order to obtain the full (N)MHV OPE contribution in question up to 6 loops, and certain parts of it up to 12 loops. We attach computer-readable files with our results, as well as an algorithm implementation which may be readily used to generate higher-loop corrections. The feasibility of obtaining the explicit kinematical dependence of the first term in the OPE in principle at arbitrary loop order, offers promise for the suitability of this approach as a non-perturbative description of Wilson loops/scattering amplitudes.


Introduction
The discovery of the duality between Maximally Helicity Violating (MHV) amplitudes and null polygonal Wilson loops in planar N = 4 super-Yang Mills theory [1][2][3] has elucidated remarkable features of its structure, such as dual conformal invariance [4] (see also [5] for a review). The controlled manner in which the latter symmetry is broken, implies that the all-loop behavior of the amplitudes is accurately captured by the BDS ansatz [6] for four and five points, and needs to be corrected by a scalar function of conformally invariant cross-ratios u i , known as the remainder function R n , for n = 6 points and beyond [7].
For the simplest nontrivial case of six points at two loops R 6 (u 1 , u 2 , u 3 ), a long expression involving transcendental functions of many variables, known as multiple (or Goncharov) polylogarithms, was first found on the Wilson loop side in [8]. This was then drastically simplified and reexpressed in terms of classical polylogarithms with the method of symbols in [9]. In extracting higher-loop information, it proves advantageous to consider kinematical limits where simplifications occur, such as the multi-Regge [10,11] and (near-)collinear [12] limit.
The Operator Product Expansion (OPE) approach to null polygonal Wilson loops is precisely an expansion in terms approaching the collinear limit at different paces, each of which receives contributions at any loop order. For the hexagon, which will be the focus of this paper, it predicts that the leading term at weak 't Hooft coupling λ = g 2 Y M N has the form where {τ, σ, φ} is a particular parametrization of {u 1 , u 2 , u 3 }, in which τ → ∞ conveniently describes the limit where two consecutive segments become collinear. As we review in the next section, each term in the expansion corresponds to a different excitation of a colorelectric flux tube, created by the two segments adjacent to the ones becoming collinear, whose energy can be calculated exactly [13] with the help of AdS/CF T integrability (see [14] for a review). The functions f (l) n (σ) above are given in terms of a single Fourier integral, whose precise integrand was found in [12] only for n = l − 1. Still, this information combined with other reasonable assumptions was enough to fix the 3-loop symbol of the hexagon remainder function up to two unknown parameters [15], which were later determined in [16]. Apart from the propagation of the flux tube excitation, what was further necessary for obtaining the integrals for any n, was knowledge of the transition amplitude, or form factor, describing how the excitation is emitted/absorbed at the two sides of the flux tube.
Great progress in this respect was recently made in [17,18], where an all-loop expression for the aforementioned form factor, or 'pentagon transition', was proposed. This formulation, which now holds for any n-gon, again crucially relies on integrability. In particular, it relates the pentagon transition to the S-matrix of excitations on top of the Gubser-Klebanov-Polyakov string [19,20], which is the string dual to the flux tube vacuum. The integrals of the so called 'flattened' part f (l) 0 (σ) were computed up to l = 4 loops by means of an ansatz, and perhaps more importantly, it is also possible to obtain more terms in the OPE expansion in this framework [21]. The data of the leading and subleading OPE contributions were sufficiently strong constraints for determining the full 3-loop hexagon remainder function [22], with the authors of the latter paper stating that this procedure could be extended to four loops and higher.
A practical question that naturally arises in this context, is how one can efficiently evaluate the integrals defining f (l) n (σ) at higher loops. A more conceptual question, is whether there exists a basis of functions large enough to describe the answer at any loop order, especially in light of the conjecture [23], that multiple polylogarithms should form such a basis for the hexagon Wilson loop in general kinematics.
In this work, we address the two aforementioned questions simultaneously. By identifying the building blocks of the leading OPE integrand, reducing the integral to a sum of residues, and employing the technology of Z-sums [24], we prove that at any loop order the σ-dependence of the contribution (1) to R 6 is given by where the c ± are numerical coefficients, and H m 1 ,...,mr (x) are a single-variable subset of multiple polylogarithms known as harmonic polylogarithms (HPLs) [25]. Our proof is algorithmic in nature, which allows us to perform the integrations in principle at any loop order. We implement it in order to determine f (l) n for any n up to l = 6 loops, and for n = l − 1, l − 2 up to l = 12, thereby providing new high-loop predictions for the MHV hexagon remainder function.
Finally, we also analyze a particular component of the NMHV ratio function R [26], for which an OPE framework has also been developed [18,27,28]. This framework is based on efforts to extend the Wilson loop/amplitude duality beyond the MHV case [29][30][31][32][33], and also on the establishment of an interesting triangle of relations of the latter two observables to correlation functions of operators in the stress tensor multiplet of N = 4 theory, in the limit where their consecutive spacetime separations become lightlike [34][35][36][37][38] 2 . We find that at weak coupling where for any l, and the c's are numerical coefficients. We similarly employ our algorithm implementation in order to obtain explicit expressions for all f (l) n up to l = 6, and the n = l, l − 1 terms up to l = 12.
This paper is organized as follows. We begin by reviewing the basic ingredients of the OPE approach for the hexagon Wilson loop in section 2, and present the integral formulas describing the leading term for the MHV remainder function and for the component of the NMHV ratio function. In section 3 we prove that their weak coupling expansion at arbitrary loop order always evaluates to the basis of functions described above. Section 4 focuses on the utility of our proof as a direct evaluation method of these integrals. We first summarize the steps of the algorithm, and then apply it in order to obtain new predictions for the remainder and ratio functions at high loop order. We conclude with an extensive discussion on the implications of our work, and possible directions of further inquiry.
The appendix contains additional information on several functions which were important in our treatment, and most notably harmonic polylogarithms. The results from all our highloop calculations, as well as the Mathematica code used to generate them, are included in seven ancillary files accompanying the version of this paper on the arXiv. 2 In particular, it was shown in [39] that the dimensional regularization of the initial super-Wilson loop proposal [29,30] breaks superconformal symmetry, and hence cannot be in correspondence with scattering superamplitudes. More recent studies of the 1-and 2-loop super-Wilson loop suggest that anomalies appear only in certain components [32], for which the symmetry can be restored either by finite counterterms [31], or by explicitly rewritting them in an invariant form with the help of the superpropagator [33]. In any case the approach [18,27,28] may be thought of as a collinear-limit expansion of the null correlators [34][35][36][37][38].

The Wilson Loop OPE
This section serves as a review of the OPE approach for the hexagon Wilson loop, and helps in establishing our notations. In subsection 2.1 we discuss how to take the collinear limit, and outline how at weak coupling the Wilson loop decomposes into terms approaching the limit at different paces, mostly based on [12,13]. Subsections 2.2 and 2.3 focus on the extension and refinement of this approach for the MHV and NMHV hexagon respectively, as presented in [17,18], and also building on [26][27][28].
For the knowledgeable reader, the equations which will form the basis of our subsequent analysis are the definition of the conformally invariant, finite Wilson loop observable (13), and its leading OPE contribution for the MHV case (16) and NMHV case (27), in terms of an all-loop integral.

Kinematics and dynamics in the collinear limit
In order to take the collinear limit of the hexagon, we start by picking two non-intersecting segments, and form a square by connecting them with another two lightlike segments (see figure 1a). We can then fix all of its 16 coordinates, 4 of them from the lightlike constraints, and the rest by conformal transformations. Specifically, the 4-dimensional conformal group SO(2, 4) will have 15 generators, which implies that not only all squares will be conformally equivalent, but also that each given square will be invariant under a subset of 3 transformations.
As we show in figure 1b, a convenient choice will be to arrange the square to lie on the (x 0 , x 1 ) plane, with its points located at the origin, past lightlike infinity (x − = x 0 −x 1 → −∞ with x + = x 0 + x 1 fixed), future lightlike infinity (x + → ∞ with x − fixed), and spacelike infinity (x 1 → ∞ with x 0 fixed). In this case, the three symmetries of the square will be dilatations, boosts in the (x 0 , x 1 ) plane, and rotations in the transverse (x 2 , x 3 ) plane, generated by D, M 01 and M 23 respectively.
For concreteness, we may choose the generators and one-parameter group elements to be whereas D is just the identity matrix, e λD = e λ . It is straightforward to show their action on the square leaves it invariant, also keeping in mind that in two dimensions spacelike infinity is a single point.
Going back to the hexagon, we can parametrize all conformally inequivalent geometries by acting with the symmetries of the square on the cusps below the reference square. In fact,  Figure 1: In (a), the dashed null segments connect two non-intersecting edges of the hexagon, crossing from two of its cusps. They break it up into three squares, with every two adjacent ones forming a pentagon. In (b), using conformal transformations we place O at the origin, and P, F, S at null past, null future, and spacelike infinity respectively. The conformal group element e −τ (D−M 01 ) leaves the middle square invariant, and its action on A and B makes them parallel to x + when τ → ∞.
it will be advantageous to consider the group element exp [−τ (D − M 01 )], as the collinear limit will now simply correspond to τ → ∞. In more detail, This implies that all segments below the middle square will flatten out onto its lower edge, as we illustrate in figure 1b, thus becoming collinear. Depending on the choice of initial hexagon, we will get a slightly different relation between the cross ratios u i and the group coordinates τ, σ, φ parametrizing the symmetries of the square. Following the conventions of [17,18], we will be using 3 3 Note that this parametrization differs from the one used in the initial Wilson loop OPE approach [12].
A great advantage of describing the collinear limit in a matter which takes into account the symmetries of the square, is that it also makes the description of the dynamics more transparent. In particular, we can think of the Wilson loop segments belonging to the middle square as a flux tube sourced by two quarks moving at the speed of light, and decompose the Wilson loop with respect to all possible excitations of this flux tube. These excitations will be eigenstates of the symmetries of the square with eigenvalues E, p, φ, so that schematically we may write where the exponential part describes their propagation, n labels different excitations, and C bot (C top ) denotes the transition amplitude, or overlap, between the initial (final) state of the bottom (top) part of the polygon and the intermediate eigenstate. This picture is reminiscent of expressing the product of two neighboring operators A and B as a sum of local operators, whose scaling dimensions control the dependence of the coefficients on the distance between A and B. Hence the decomposition (9) around the collinear limit has been coined 'Wilson loop OPE'. Since its derivation only relied on the symmetries of the problem, the above formula should hold for any conformal field theory where the flux is conserved. The good news is that in N = 4 super Yang-Mills, the flux tube excitations are in 1-1 correspondence with excitations of an integrable spin chain, with the collinear twist operator [40] D − M 01 as its hamiltonian. The states of the spin chain are single-trace operators, with the vacuum made of a sea of derivatives acting on the complex combination of two scalars of the theory Z, and excited states built by inserting any fundamental field of the theory Φ on the vacuum, e.g. single excitation = tr ZD S 1 where S ∼ S 1 + S 2 0. In fact, minimizing the energy suggests that elementary excitations of the spin chain can only be drawn from the components of Φ which have minimal classical twist eigenvalue ∆ − S = 1. At quantum level the (shifted with respect to the vacuum) twist of the operators (11) receives anomalous contributions due to renormalization, Due to integrability, the 'energy' E can be calculated to all loops [13], and moreover for M elementary excitations we will have E M = M + O(λ). This teaches us that at weak coupling the Wilson loop OPE (9) will be a sum of terms with different integer exponential behaviors e −τ M as τ → ∞, where M is the classical twist of the state, or equivalently the number of elementary excitations it consists of. In this paper, we will focus on the leadingtwist, or single-particle term O(e −τ ), and we will frequently refer to it as the 'leading OPE contribution'.

MHV hexagon
Let us now become more specific, and define the hexagon Wilson loop-related observable which will be most suited for analyzing its OPE. By construction, the middle square always has two of its segments coinciding with the cusps of the hexagon, and hence the bottom and top part of the hexagon which lie outside of it will also be squares. Then, where W is our hexagon, and W , W bot , W top are the Wilson loops defined on the contours of the middle square, and of the pentagons created by joining the middle and lower, and middle and upper squares (see figure 1a). The particular ratio we are considering removes all (cusp-induced) ultraviolet divergences, leaving a finite function of conformal cross ratios. Aside this, it does not cause any loss of information, as the square and pentagon Wilson loops are given by the BDS ansatz [6].
Specializing on the MHV case, from symmetry arguments the single-particle contribution is expected to be bosonic and uncharged under the R-symmetry. This picks out only one out of the twist-1 excitations, the component F +i of the gauge field [12], where the first component is projected on x + and the second component is on the (x 3 , x 4 ) plane, according to the notations discussed in the beginning of section 2.1. As we also mentioned in that section, its all-loop dispersion relation has been found with the help of integrability [13], in parametric form with respect to the Bethe rapidity u. Following the notations of the latter paper, from this point on we will rescale the 't Hooft coupling λ by in terms of which the gauge field dispersion relation reads where the functions γ ∅ ± (2gt) are independent of u, and as we review in appendix A.1, can be obtained iteratively as a Taylor expansion in g 1. In order to obtain information about the hexagon from the analogue of (9) for r, we need however to know the (rescaled) creation/absorption form factors C bot , C top , which also depend on g. Initially, these were determined to the first few orders by comparing with the explicit computation of the hexagon up to two loops [8,9]. Recently however, all-loop expressions for them were proposed, relying again on the integrability of the theory [17,18]. In particular, these form factors, also dubbed as 'pentagon transitions', are related to the S-matrix of excitations on top of the GKP string [19,20], which forms the vacuum of the flux tube. 4 With this ingredient in place, and in the more natural rapidity parametrization, the final formula for the leading (single-particle) OPE contribution for r reads where the measure µ 1 (u) is given by In the last formula, x ± are the Zhukowski variables J i is the i-th Bessel function of the first kind, and the form of the functions Once the observable (16) has been determined, the near collinear limit of the MHV hexagon remainder function R (we drop the index as we will only be dealing with n = 6), parametrizing the part of the Wilson loop W which is not captured by the BDS ansatz W BDS [6], is given by where [41] is defined as in (13), but only including the BDS contribution for each of the polygons. The equality (20) is a consequence of the definition (19), given that the BDS ansatz accurately describes square and pentagon Wilson loops. In more detail, the logarithm of the BDS part of any null polygonal Wilson loop in N = 4 super Yang-Mills, is proportional to the logarithm of the same Wilson loop in the U (1) theory, with the proportionality constant being a fourth of the cusp anomalous dimension Γ cusp . The latter quantity can also be calculated to all loops (see [42] for a review), and is given to the first few orders in the weak coupling by

NMHV hexagon
Scattering amplitudes involving gluons in other helicity configurations, or other particles of N = 4 super Yang-Mills, are most conveniently described by exploiting its (dual) superconformal symmetry [26], see also [5] for a review. Very briefly, one starts by packaging the particle content of the theory 5 in a single superfield Φ with the help of a Grassmann variable η A , whose index transforms in the fundamental representation of the R-symmetry group SU (4). Namely, all external states of ±1 helicity gluons G ± , ± 1 2 helicity Majorana fermions Γ A ,Γ A , and zero helicity real scalars S AB can be simultaneously described by which in turn allows us to combine all n-point amplitudes in a superamplitude A n (Φ 1 , . . . , Φ n ).
All MHV amplitudes form the part of the superamplitude which has 8 powers of Grassmann variables, starting the the MHV gluon amplitude, On the basis of tree-and 1-loop level amplitude computations, it was argued in [26] that NMHV amplitudes will similarly organize in a homogeneous polynomial of degree 12 in η A i , and more importantly, that they have the same infrared divergence structure as the MHV amplitudes, so that the two superamplitudes are related by where R n is the (dual conformal invariant) NMHV ratio function. Evidently, it will consist of terms involving 4 powers of the Grassmann variables, whose components we can be denoted as R (ijkl) n . Following attempts for generalizing the Wilson loop/scattering amplitude duality beyond the MHV case [29,30], an analogous proposal for the OPE of certain components of the NMHV hexagon Wilson loop was put forth in [27,28]. According to the latter, the dual of the R (i,i+1,j,j+1) 6 component is given by a Wilson loop 6 W (i,i+1,j,j+1) 6 -normalized by bosonic 5 The on-shell fields of N = 4 super Yang-Mills form the CPT self-conjugate, 'doubleton' representation of the superconformal group P SU (2, 2|4) [43]. 6 More precisely, R 6 and W 6 will only differ by O(e −2τ ) terms [18].
squares and pentagons as in (16) -which at tree-level has insertions of a complex scalar field combination at the cusp between segments i, i + 1, and its complex conjugate field at the cusp between segments j, j + 1. The two segments forming the flux tube reference square are chosen to lie between i, i + 1 and j, j + 1 without coinciding with any of them, and |i − j| ≥ 3.
Thus the top and the bottom part of the polygon have a scalar insertion each, and it is also expected that only scalar excitations of the flux tube will propagate to leading order in the OPE expansion. For a particular component of the NMHV ratio function, this leading OPE contribution was also predicted in [18], by a similar analysis of the scalar pentagon transitions. It is the component with the scalar insertions between edges 6-1 and 3-4, and edges 2 and 5 forming the reference frame, for which the OPE reads (we drop the lower index) where γ 0 (u) and p 0 (u) are the anomalous part of the energy and the momentum of a scalar excitation propagating between the bottom and top part of the polygon [13], and the scalar integration measure µ 0 is now given by Here as well J i is the i-th Bessel function of the first kind, γ ∅ ± are the same functions which appeared in the MHV hexagon subsection, and which are reviewed in appendix A.1, and f i (u, v) are functions which are in turn reviewed in appendix A.2. We'll also have the same parametrization (8) of the conformally invariant cross ratios u 1 , u 2 , u 3 in terms of τ, σ, φ.

General Analysis of the Integrals
In section 2, we recalled that the Wilson loop OPE approach yields the leading, O(e −τ ) term in the τ → ∞ collinear limit of the (N)MHV hexagon as a 1-dimensional Fourier integral, to all orders in g. Here we will focus on the weak coupling expansion of the integral, and prove that it can be evaluated at any loop order, in terms of harmonic polylogarithms.

Reduction to sum over residues
Let us start by combining the descriptions of the MHV and NMHV integrals (16), (27) simultaneously 7 , with the help of the definition where η = 1 corresponds to the propagation of a gauge field excitation and describes the MHV case, and η = 0 corresponds to a scalar excitation and describes the NMHV case. Similarly γ, p and the functions f 3 , f 4 will implicitly depend on η, as reviewed in appendix A.2, equations (78), (79) and (81). In obtaining (31), we made use of the identity which may be proven as was done in [13,44] for similar integrals involving Bessel functions 8 , in order to replace e −t → e t in the second line of (17). We will proceed to evaluate the integral (30) order by order at weak coupling g by turning it into a sum over residues. To this aim, we will need to analyze the general dependence of µ, γ and p on the integration parameter u.
Focusing first on the exponential part of µ, and in particular the f i functions contained therein, and defined in (78), it is evident that the dependence on u enters only through the κ,κ vectors (79). As a result of the regular Taylor expansion of the Bessel functions involved (77), at weak coupling these can always be expressed in terms of the following integrals which can be calculated exactly 9 , where ψ (m) (z) is the polygamma function. It is easy to see that in both NMHV (η = 0) and MHV (η = 1) case the possible arguments include z = 1 2 ± iu and z = 1, and for the second argument we further reduce to Riemann zeta functions, For the MHV case, there exists one additional possibility for the argument, z = 3 2 ±iu, and by definition the weak coupling expansion of the f i functions will contain bilinear combinations of these ψ-functions.
A similar analysis can be performed for the remaining exponential part of µ(u) (31), and also for γ(u) and p(u), which shows that we now obtain monomials of the aforementioned polygamma functions with the specific arguments. For γ(u), p(u) this is particularly easy to see due to the alternative expression (81) mentioned in appendix A.2, which relates them to the same building blocks of the f i functions.
Finally, it is straightforward to show that at any loop order in g 1, the x ± -dependent factor of the measure (31) will be a sum of products of inverse powers of u ± i 2 . In more detail, the factor in question is equal to (u 2 + 1 4 ) −1 at g = 0, has a regular Taylor expansion around that point, and due to 10 it has effective expansion parameters g/(u ± i 2 ). Gathering the information we obtained, we deduce that the weak coupling expansion of the integral (30) will be a sum of the general form where we have absorbed all factors that don't depend on u in the coefficients c (which obviously also depend on the different indices l i , m i , r i etc). From the last formula, we arrive at the following important conclusion: The only possible locations where the integrand may have poles are for u = (k + 1 2 )i, k ∈ Z. In particular these poles may come from the denominator (if any), the hyperbolic cosecant, or the polygamma functions at negative integer arguments.
In what follows we will restrict to σ > 0, in which case we can close the contour with a semicircle on the u > 0 plane, whose integral at infinite radius will go to zero due to the decaying exponential. Therefore by Cauchy's residue theorem we will have 10 We have chosen the branch √ u 2 = u, which is equivalent to assuming x(u) starts at O(g 0 ). The opposite branch corresponds to x(u) starting at O(g 2 ), and the two solutions are related by x(u) → g 2 /x(u), as can be seen by the defining relation u = x(u) + g 2 /x(u). From this it follows that the expansions in the two branches of the factor in question only differ by an overall sign, but as pointed out in [18], there already exists a sign ambiguity in µ(u), that needs to be fixed by physical input. and in order to proceed we will need to find an analytic expression for the residues as a function of k. This cannot be achieved by directly Taylor expanding the expressions for the integrands around u = (k + 1 2 )i, as ψ (n) 3 2 + iu and ψ (n) 1 2 + iu develop poles there. Instead, we first employ the reflection formula which allows us isolate the singular terms into elementary functions with known expansions.
In more detail, we first use the recurrence relation with z = 1 2 ± iu, in order to eliminate the ψ (m) ( 3 2 ± iu) factors in (37) 11 . Then, we apply (39) with z = 1 2 + iu, also noting that the cotangent derivative will give a polynomial in cotangents and cosecants, which in terms of u become −i tanh πu and −sechπu respectively 12 In this manner, we have achieved to reexpress (37) as a sum of products of a restricted set of functions, having known expansions for u = (k + 1 2 )i + , small, which we also present below: tanh πu = coth 11 Notice that after expanding, we will not get any new rational factors apart from the already existing u ± i 2 in (37). 12 In more detail, it can be shown that where B 2n are the Bernoulli numbers. A few remarks are in order. First, since the hyperbolic functions become periodic in the imaginary axis, their k-dependence reduces to an overall sign at most. Furthermore, because of the last equation, we will have to determine the residue for k = 0 separately. Aside this value, we may use the above expansions in order to determine the residues in (38) for general positive integer k. Finally, the fact that the residues will only contain ψ-functions with integer arguments, allows us to replace them with generalized harmonic numbers, defined as where γ E = −ψ(1) = 0.577 is the Euler-Mascheroni constant.

Z-sums and Harmonic Polylogarithms
Let us now focus on the structure of the residue ofh for different values of the indices s, r, l 1 , l 2 , m 1 , . . . m r and numerical constants c. In more detail, (43) will only contribute numerical factors independent of k. The same will apply for (44), except for an overall (−1) k factor, as the discussion in footnote 12 implies that we will always have odd powers of sechπu. The latter k-dependent factor combines with the exponential coming from (41), which will also contribute the powers of σ. Finally the inverse powers of k, k + 1 come from (45), (46), and the products of harmonic numbers from (42) due to (47).
At this point, we will choose to concentrate on the MHV case (η = 1), and come back to examine what changes for the simpler NMHV case at the end of this subsection. There exist three additional steps we need in order to bring these terms in a form where we will be able to perform the summation over all poles k. First of all, we can always partial fraction in order to get only inverse powers of either k or k + 1, since the degree of the numerator is always smaller than the degree of the denominator, and in particular zero. Secondly, we can perform certain manipulations in order to reduce the number of different arguments in the sums. This procedure, generally known as synchronization [45], is necessary for the efficient evaluation of sums of this type with the help of the computer.
In our case, achieving this goal involves treating the terms with powers of k and k + 1 differently. In particular, we can eliminate the latter in terms of the former by redefining the summation index k = k + 1, and extending the summation range to also include k = 1 by adding and subtracting the term in question. More concretely, and notice that this redefinition of the summation index will turn e −σ → −e σ in (49). In addition, the argument of the harmonic numbers included in these terms will now be k − 1, and we can turn this into the argument of the harmonic numbers multiplying inverse powers of k, by replacing The formula above follows immediately from the definition (48), and is the analogue of (40) for harmonic numbers. Evidently we do not need to partial fraction again, as we only obtain additional powers of the same monomial in k.
Finally, we will make use of the fact that harmonic numbers are the simplest case of a more general set of nested sums, known as Z-sums [24]. These are defined as Z(n; m 1 , . . . , m j ; x 1 , . . . , x j ) = where Z(n) is equal to the unit step function. The sum of m i is known as the weight, or transcendentality, and the number of summations j as the depth. Clearly, harmonic numbers are depth-1 Z-sums, S m (k) = Z(k; m; 1). An important property of these objects, also known as the quasi-shuffle algebra 13 , is that a product of two Z-sums with the same outer summation index can be reexpressed as a linear combination of single Z-sums. This easily follows by splitting the square double summation range into regions with definite index ordering, namely For our purposes, rather than the general recursive procedure for decomposing a product of Z-sums, we will just need the particular case (we set all x i = 1, and for compactness drop both them and the outer summation index n from our notation), In other words we take all permutations that preserve the order of indices of the two Z-rums on the left hand side, and also add the l index to all m i . The above formula can then be used recursively in order to decompose the product of harmonic numbers/depth-1 Z-sums into single harmonic sums in (49). After these three steps, the k-dependent part of all terms in the sum over residues (49) will be itself proportional to a Z-sum, The critical observation, based on [24], is that this particular Z-sum precisely coincides with the series representations of the harmonic polylogarithm H m 1 ,m 2 ,...,m j (−e −2σ ) of Remiddi and Vermaseren [25] 14 ! We have therefore rigorously proven that the integral yielding the single-particle contribution to the Wilson loop OPE of the MHV hexagon (16) will have a weak coupling expansion of the form (36), with its σ dependence always given by where c ± are numerical coefficients. Finally, let us go back and extend the last part of our analysis to the case of the NMHV hexagon component integral, namely (36), (37) with η = 0. We recall that since the latter formula will not contain any ψ (m) ( 3 2 ± iu) or inverse powers of u ± i 2 , the sum over residues (49) will not contain any inverse powers of k, k + 1. Thus we may directly reexpress the harmonic numbers as a linear combination of Z-sums with the help of the the quasi-algebra (55) in (49), and obtain ∞ k=0 ∼ σ s e −σ d dy H 1,m 1 ,...,mr (y)| y=−e −2σ = σ s e −σ 1 + e −2σ H m 1 ,...,mr (−e −2σ ) .
The equality of the last line follows from the definition of HPLs (83), also taking into account the discussion about their "a"-and 'm"-notation before equation (86). We similarly conclude that the weak coupling expansion of the single-particle contribution to the Wilson loop OPE of the NMHV hexagon component (27) will be proportional to (36) with The general, all-loop structure of the hexagon Wilson loop OPE contributions for the MHV (57) and NMHV (59) case agree with the ansätze, based on empirical evidence 15 , used in [18] 16 . Before concluding, we should note that by refining our analysis so as to keep track of the maximum powers of g that can multiply the terms of the integrand (37), it is possible to show that the indices of the terms summing up to h l n (σ) (57), (59) are constrained as follows, where again η = 0, 1 corresponds to the NMHV, MHV case respectively. This of course is in agreement with the fact that all N = 4 amplitudes computed to date have maximal transcendentality 2L at L loops, given also that the Taylor expansion of transcendental functions yields terms with the same or smaller transcendentality.

Algorithm
In the previous section, we presented a general proof for the exact basis of harmonic polylogarithms, including the dependence of their coefficients and arguments on the kinematical data, which is suitable for describing the leading OPE contribution of the hexagon Wilson loop at any order in the weak coupling expansion g 1. Furthermore, it is perhaps evident that our proof forms an algorithmic process which allows us to directly compute the Fourier integrals (36) and obtain explicit expressions for the relevant part of the hexagon in terms of the kinematical variables (8). We stress again that this method of computation can be applied at arbitrary loop order l, subject to restrictions in computational power.
Let us now summarize the steps of the algorithm, which facilitate its implementation in any computer algebra system. We start with the weak coupling expansion of the integrand (30), which as we reviewed in section 2, can be performed with the help of the definitions of µ, γ, p, and the integral (33). Then, we 1. Replace all ψ-functions appearing with ψ (n) 1 2 − i p 2 with the help of (39) and (40). 2. Replace u = (k + 1 2 )i + and find the residue for general k ≥ 1, with the help of (41)- (46). 3. Replace ψ (m−1) (k + 1) → H m (k) = Z(k, m, 1) by means of (47).
4. Partial fraction the inverse powers of k, k + 1, if any. 5. For the k ≥ 1 sum over residues, synchronize the arguments by applying (50) to the terms with k + 1 powers, and (51) to terms with k powers.
6. Reduce products of Z-sums to linear combinations of single Z-sums by recursively applying (55).

New predictions: MHV hexagon
With the steps of the algorithm set in place, we can now proceed with its implementation, in order to determine the leading OPE contribution of the hexagon remainder function to high loop order. We remind the reader that this contribution will be given in terms of the h by virtue of equations (16), (20), (22), (23) and (36).
We have included the results, together with a Mathematica code that generates them, in ancillary text files accompanying the version of this paper on the arXiv. The code is fully general, and can be used in principle for any given value of l. In particular, we have used it to obtain h (l) n (σ) for all allowed values of l − 1 ≥ n ≥ 0, up to l = 6, thereby providing information about the 5-and 6-loop hexagon for the first time.
Generating results at higher loops from this code is only a matter of computational power and optimization, and to illustrate this we also calculated the n = l − 1, l − 2 components up to l = 12. As far as the efficiency of the code is concerned, we should stress that even without particular attention to optimization, the computation of all 4-loop leading-twist terms takes about 10 seconds on a portable computer, and all 5-loop terms about a minute! Let us start by briefly mentioning what is already known about the hexagon in the nearcollinear limit, and in the first few orders at weak coupling. In [18], the n = 0, also dubbed "flattened" part of the hexagon was computed up to 4 loops under the assumption of an ansatz for the general structure of the expression. More specifically, the free parameters of the ansatz were determined by comparing its Taylor expansion with a finite number of terms in the sum over residues, which the Fourier integral (36) reduces to 17 .
In the previous section, we proved that the structure of this ansatz is correct at any loop order, which thus places the aforementioned calculations on a firmer setting. As a further consistency check, we first aimed to reproduce the results reported in [18] for the h (l) 0 (σ) functions 18 , n = 1, . . . , 4 with the help of our direct computation method.
More specifically, we compared with the expressions contained in the Mathematica file Functionshf.nb accompanying the latter paper, and found indeed agreement. Due to the many functional identities between harmonic polylogarithms, the two expressions are not identical, but they can be brought in the same form with the help of the HPL package [49,50]. To this end, we use property (90) in order to change the argument of the HPLs to e −2σ , and also replace all powers of σ → − 1 2 H 0 (e −2σ ) multiplying them. Then, we employ the command HPLProductExpand in both expressions, which eliminates any products of HPLs, in favor of their linear combinations.
Next, we moved on to determine all non-flattened leading-twist contribitions h (l) n (σ), n = 0 for l ≤ 4. As all HPLs in the expressions have argument −e −2σ , we will omit it for compactness. At two loops we have at three loops 19 + 2 20 + π 2 − 24σ + 4σ 2 H 1,1 − 4(2σ + 1)H 2,1 + 48(1 − σ)H 1,1,1 + 48H 1,1,1,1 . 17 See also [48] for an application of the same method in the multi-Regge limit. 18 In the notations of [18], h (l) 0 (σ) → f l (σ). 19 Recently [22] appeared, which calculates the 3-loop hexagon in general kinematics, and also specializes to the near-collinear kinematics considered here. We have checked that the results of the two calculations agree if we take into account that the expansion parameter in the latter reference is 2g 2 . See also [51] for a calculation of the analogue of h n (σ), as a function of τ, σ. Colors of the visible spectrum denote different values of h (6) , increasing from blue to red. The function is always positive, and monotonically increasing and decreasing in τ and σ respectively.
Evidently the length of the expressions grows quite fast with the loop order, due to the increase not only in the number of terms in the integrand, but also in the quasi-shuffle algebra decomposition (55) for higher powers of harmonic numbers.
Finally, we computed all 5-and 6-loop leading-twist contributions to the OPE of the MHV hexagon, as well as the h

New predictions: NMHV hexagon
Similarly, we employed our algorithm implementation in order to compute all O(e −τ ) contributions in the OPE of the NMHV hexagon component (27) up to 6 loops, and all of its leading and subleading powers of τ up to 12 loops. Also in this case, our results agreed with the previous indirect computation of [18], contained in the file Fnl.nb attached to the latter paper, up to 3 loops.
Of our new results, for reasons of space we will only mention here three out of the five single-particle contributions at 4 loops, saving the rest for the attached ancillary files. After we rescale the overall σ-dependent factor in (59), and also shift the loop index so as to measure the total powers of g 2 in (27), the aforementioned component of the NMHV ratio function will be given by

Conclusions and Outlook
In this paper, we explored the implications of the recent conjecture [17,18], which adds significant new ingredients to the OPE approach to null polygonal Wilson loops, for the case of the MHV and NMHV hexagon. We found that the integrability underlying this proposal, and manifesting itself in the presence of polygamma functions in the energy and higher conserved charges of the flux tube 20 , heavily constrains the conformally invariant functions 20 In [52], the hamiltonian density of the Heisenberg XXX, or SL(2), spin chain of spin s was found to be a digamma function of the operator measuring the total spin of the adjacent sites. A similar expression for the hamiltonian density of the P SU (2, 2|4) spin chain, giving rise to the complete 1-loop dilatation operator of N = 4 SYM, was found in [53] by appropriate identification of an SL(2) subsector. The generalization to parametrizing the hexagon: It implies that the leading OPE contribution for these functions is expressed in terms of harmonic polylogarithms with specific kinematical dependence, at any loop order in the weak coupling expansion! This general result is consistent with the evidence based on a particular "d log" form for the all-loop integrand, and presented in [23], that a basis of multiple (or Goncharov) polylogarithms is sufficient for describing the MHV and NMHV hexagon. In particular, harmonic polylogarithms are a single-variable subset of multiple polylogarithms. The basis of harmonic polylogarithms we found also agrees with the ansätze assumed in [18] in order to compute (parts of) the leading OPE contribution up to 3 loops for the NMHV case, and up to 4 for the MHV, thus further supporting their validity.
Moreover, our analysis has important consequences in the practical front of extracting data for the hexagon at higher orders in the weak coupling expansion. Starting from the Fourier integral form of the aforementioned OPE contribution as predicted by [18], our proof of general structure forms an algorithm for the direct computation of these integrals, in principle at any loop order. By implementing this algorithm in Mathematica, we were able to obtain new results for the full O(e −τ ) term in the near-collinear limit of the MHV remainder function and NMHV ratio function up to 6 loops, and the part of this term with the leading and subleading powers of τ up to 12 loops. We include these results in the attached ancillary files, together with an implementation of the algorithm, which can be readily used to evaluate higher-order contributions as well.
Input from the Wilson loop OPE has been crucial, together with other reasonable assumptions, for the recent determination of the full 3-loop MHV remainder function in general kinematics [22]. The authors of the latter paper report that the same methodology may be used for obtaining the remainder function at 4 loops and higher, and we hope that the data presented in this paper will again prove useful in that respect. The same applies for the NMHV ratio function, which is currently known to 2 loops [54]. Furthermore, it has been argued in [22], that the hexagon MHV remainder and NMHV ratio functions are described by hexagon functions, a more restricted class of combinations of multiple polylogarithms with proper branch cuts. It would be an interesting consistency check to investigate whether their near-collinear limit reduces to the basis (57) at arbitrary loop order.
There is a number of exciting open questions which require further inquiry. Reconstruction of the full hexagon Wilson loop from its OPE would require obtaining all terms with higher powers of e −τ in the near-collinear limit, which are contributions of multiparticle states propagating in the flux tube. Indeed, it is possible to extend the current framework in order to incorporate 2-particle excitations [21] and higher, and one would ideally like to find a general basis/method to calculate the resulting integrals for these terms as well. Note that even though the dimensionality of the integrals will be equal to the excitation number, data from the 3-loop hexagon remainder function show that the e −2τ terms are described by harmonic polylogarithms as well [22]. Naïvely similar multidimensional integrals, which are however not expected to only yield HPLs, appear in the near-collinear limit of the heptagon. the all-loop dilatation operator relies on a deformation of the spectral parameter, which loosely speaking is responsible for the appearance of polygamma functions.
For the latter, very little is known beyond the 2-loop total differential [55,56] and motivic avatar structure [57] for the MHV case, and symbol for the NMHV case [16].
Aside the near-collinear limit, a similar all-loop integral formula for the hexagon remainder function also exists in the multi-Regge limit [10,11]. Based on evidence that the relevant class of functions for describing this limit are certain combinations of HPLs that are single-valued on the complex plane, the authors of [48] were able to extract the (nextto-)leading-logarithmic part of the aforementioned formula to high loop order. Given the resemblance of the integral formulas in the near-collinear and multi-Regge limit, it should be possible to extend our method in order to rigorously prove the appropriateness of the single-valued HPLs, and possibly the conjecture [58] generalizing the results of [48] for the leading-lorarithmic part to all loops. More broadly, the similarities between the two limits seem to be quite extensive, and suggest that a physical picture based on integrability should exist for the multi-Regge limit as well. For example, the main quantities entering the master formula of [10,11] only depend on polygamma functions up to the currently known order.
Finally, our work may prove useful for the computation of other observables of N = 4 super Yang-Mills theory. In [59,60], a very interesting class of infrared-finite weighted crosssections was introduced and extensively studied. These measure the total flow of charge registered by detectors positioned at spatial infinity, for an initial state generated by the scalar half-BPS operator acting on the vacuum. For the choice of charges considered in the latter papers, and for two detectors measuring them at different directions simultaneously, these were shown to be related to 4-point correlation functions of components of the stresstensor multiplet. In particular, they are given by the convolution of the Mellin transform of the Euclidean correlators with a coupling-independent 'detector kernel'. The resulting inverse Mellin integral was computed from its known ingredients at one loop, and a generalization of our method may be applicable for performing the computation at higher loops, given the relation between Fourier and Mellin transforms, We expect that again harmonic polylogarithms will play a role in describing the dependence of these double flow correlations on the single variable parametrizing them, the angle between the two detectors.

A Review of Useful Functions
A. 1 The γ ∅ ± functions For completeness, we review the functions γ ∅ ± [13] which enter in the expressions (15) and (28), for the anomalous dimensions and momenta of the flux tube excitations propagating between the two parts of the polygon.
These are defined as where the coefficients γ ∅ n depend on g and obey The last two equations can be solved perturbatively in g 1, with the help of the following Taylor explansion of the i-th Bessel function of the first kind J i (z) around z = 0, In any case, from the latter formula it is evident that for small t, γ ∅ ± have a regular Taylor expansion.

A.2 The f i functions
In this appendix, we review the computation of the f i functions [18] appearing in the integration measure (17), (29) of the leading (N)MHV hexagon OPE contribution.
The functions in question are defined as where the κ,κ are vectors with elements with η = 0 for the NMHV case and η = 1 for MHV case, and Q, M are matrices independent of u. In particular, Q has matrix elements, Q ij = δ ij (−1) i+1 i, and M is related to another matrix K, The Taylor expansion of the Bessel functions J i (77) implies that K ij starts at order O(g i+j ).
Hence even though the vectors and matrices entering (78) are infinite-dimensional (equivalently they are an infinite triple sum), if we wish to obtain f i up to O(g n ), we can truncate them to their first i, j = 1, . . . , n − 1 components. The building blocks of the f i functions can also be used in order to calculate all remaining ingredients of the leading OPE contribution for the (N)MHV hexagon. Namely the anomalous dimension and momentum of the single-particle flux tube excitations can be written as γ(u) = 4g (Q · M · κ(u)) 1 , p(u) = 2u − 4g (Q · M ·κ(u)) 1 , and the cusp anomalous dimension as The expressions (81), which also implicitly depend on η by virtue of (79), are equivalent to the ones used in the main text (15), (28), and can be derived from them.

A.3 Harmonic Polylogarithms
In this appendix, we review the definition and basic properties of harmonic polylogarithms (HPL) [25], also based on [50,61]. For x ∈ (0, 1) and a i = {−1, 0, 1}, harmonic polylogarithms are defined as H(x) = 1 , H(a 1 , . . . , a n ; x) = 1 n! log n x if a 1 = . . . a n = 0 , x 0 dyf a 1 (t)H(a 2 , . . . , a n ; t) otherwise, where the auxiliary functions f a entering the case of the recursive definition are given by Precisely because of the shuffle algebra, not all HPLs of a given weight will be algebraically independent, as they can be expressed in terms of other HPLs of the same weight, and products of lower weights. In any case it is always possible to construct a basis of algebraically independent HPLs, whose indices form Lyndon words [62], although we will refrain from exploiting this property here.
Finally, an argument transformation property for HPLs with nonzero indices in the "m"notation, which we will be making use of, is that H m 1 ,...,m k (−x) = (−1) k H −m 1 ,...,−m k (x) , for m k = 0 .
This follows simply from the definition (83), after we change the sign of the integration variables.