Hexagon OPE Resummation and Multi-Regge Kinematics

We analyse the OPE contribution of gluon bound states in the double scaling limit of the hexagonal Wilson loop in planar N=4 super Yang-Mills theory. We provide a systematic procedure for perturbatively resumming the contributions from single-particle bound states of gluons and expressing the result order by order in terms of two-variable polylogarithms. We also analyse certain contributions from two-particle gluon bound states and find that, after analytic continuation to the $2\to 4$ Mandelstam region and passing to multi-Regge kinematics (MRK), only the single-particle gluon bound states contribute. From this double-scaled version of MRK we are able to reconstruct the full hexagon remainder function in MRK up to five loops by invoking single-valuedness of the results.


Introduction
Light-like polygonal Wilson loops have been the subject of much study recently, in particular because of their relation to scattering amplitudes in the planar limit of N = 4 super Yang-Mills theory [1][2][3][4][5][6][7][8][9]. The conformal symmetry of the theory and its associated Ward identity [10] imply that the Wilson loops with four or five sides are essentially trivial, with no interesting dependence on the configuration of the contour. From six points onwards however there is a conformally invariant function of the loop contour which needs to be determined. Here we will focus on the bosonic hexagonal Wilson loop, which corresponds to the six-particle MHV amplitude under the amplitude/Wilson loop duality. The conformally invariant function of the loop is known as the 'remainder function' in this case and is a function of three cross-ratios.
One approach to studying the kinematical dependence of light-like Wilson loops is to make an ansatz based on the analytic behaviour of explicit results obtained for terms in the perturbative expansion such as [11][12][13]. This approach has yielded results up to four loops at six points [14][15][16][17][18] and up to three loops at seven points [19]. The ansätze are given in terms of iterated integrals over words formed from a specified set, or alphabet, of rational functions (known as letters). The alphabets appearing in the explicit two-loop expressions of [13] were observed [20] to be described by the A-coordinates of a class of cluster algebras [21,22] associated to the Grassmannians G(4, n). Taking this observation as an assumption, together with basic analytic information on the locations of possible cuts of the final expressions produces a rather restrictive ansatz for the relevant functions. Indeed, at three loops and seven points, these analytic assumptions were essentially used to replace the dynamical information of the theory, producing a unique result for the symbol of the heptagon remainder function [19].
A second very powerful approach to describing the Wilson loops is based on a type of operator product expansion applied to configurations of null Wilson lines [23][24][25][26][27]. In this approach the Wilson loops are given by an infinite sum over excitations of a light-like flux tube. This sum amounts to an expansion around a collinear limit, and the excitations may also be thought of as insertions of the fields of the theory (gluons, fermions and scalars) on the Wilson line segments that are becoming collinear. The number of excitations or insertions is equal to the order at which these appear in the near-collinear expansion. The spectrum of excitations was calculated exactly in [28] using techniques based on the integrable structure exhibited in the planar N = 4 theory. Apart from the spectrum, a set of overlap functions, called pentagon transitions, are required to calculate the near-collinear expansion. In [29] an all-loop formula was proposed to describe the subset of pentagon transitions consisting of any number of gluons, including multi-particle bound states.
It is interesting therefore to understand how the near-collinear OPE expansion can be resummed into the type of cluster polylogarithmic functions appearing in the perturbative expansions of the null Wilson loops. As a first step towards the resummation of the OPE series, we develop a systematic procedure for summing the contributions to the OPE of single-particle bound states of an arbitrary number of gluons in the weak coupling expansion of the MHV 6-point (or hexagon) remainder function. Since these states have a restricted dependence on the helicities of the gluons, this can really be thought of as a contribution to the 'double scaling' limit of the Wilson loop, where one of the three cross-ratios is taken to zero. In the perturbative regime the double scaling limit is entirely governed by the gluon insertions, i.e. to determine it entirely at each perturbative order one needs all multi-particle bound states of gluons but not the contributions of fermions and scalars [29].
Building on the previous works [30,31], our procedure relies on the technology of nested sums [32], and in particular its nestedsums C++ library realisation [33]. Using these algorithms we are able to show that the OPE expansion for single-particle gluon bound states can always be resummed into two-variable polylogarithmic functions based on a five-letter alphabet. These functions may be described as A 2 polylogarithms in the cluster algebra language [20] or polylogarithms on M 0,5 in the language of [34] or as a subset of the two-dimensional polylogarithms of [35]. In fact we find that only a particular subset of these functions arises, which is consistent with the functions having restricted branch cuts, and indeed with the idea that they are particular limits of hexagon functions [16,17] describing the full six-point remainder function.
Perhaps more importantly, from the knowledge of the single-particle bound state contributions to the double scaling limit we can produce all but the power-suppressed terms of the MHV hexagon in multi-Regge kinematics corresponding to high energy gluon scattering in the 2 → 4 Mandelstam region. The kinematic limit on the cross-ratios is formally the same as the soft limit, in which the remainder function vanishes. To obtain a non-trivial dependence in multi-Regge kinematics, it is necessary to analytically continue the remain-der function to the 2 → 4 Mandelstam region, achieved by going around a singularity where one of the cross-ratios vanishes before taking the limit. In [36] it was observed that this continuation can be carried out order by order in the OPE expansion. Using this property we provide evidence that after analytically continuing and taking the limit, the single-particle bound states we have summed are the only excitations with non-vanishing contributions on a one-dimensional double scaling slice of the full two-dimensional space parameterising multi-Regge kinematics. In other words we find that, after analytic continuation, the contribution from multi-particle bound states of gluons in the double scaling limit is always suppressed when taking the multi-Regge limit.
Under the well-justified assumption that the natural function space in these kinematics is given by the class of single-valued polylogarithms (or SVHPLs) [37], we may then reconstruct the remainder function in the full parameter space, from its knowledge on the line obtained from the double-scaling limit. We have carried out this procedure explicitly up to five loops. We emphasise that in our approach the contribution of the scalars and fermion insertions in the OPE expansion is non-vanishing, even after analytic continuation and taking the Regge limit. Indeed one can already see this in the analysis of [36]. The matter contributions just vanish on the line we are considering, hence we are actually reconstructing their contribution from the gluon bound states by invoking single-valuedness.
Our logic is very much in line with the analysis of [38] who also reconstructed the multi-Regge limit of the hexagon amplitudes in the 3 → 3 Mandelstam region from the OPE contributions of just the single-particle gluon bound states. Indeed, starting with the finite-coupling integral expression for the same OPE contribution that we consider at weak coupling, and performing an additional analytic continuation in Mellin space, the authors of [38] arrived at finite-coupling expressions for the MHV hexagon in multi-Regge kinematics. As this second continuation is fundamentally non-perturbative, it is quite interesting that in our approach we can trade it for single-valuedness directly at the perturbative level.
The paper is structured as follows. We begin in section 2 by describing the OPE expansion of the hexagon Wilson loop in the double scaling limit. Our focus is on the contribution of the single-particle gluon bound state contributions and their resummation order by order in perturbation theory into two-variable polylogarithms of a particular kind, which we perform in section 2.3. We then describe in 2.4 why the class of functions so obtained is consistent with the idea that the full hexagonal Wilson loop remainder function is expressed in terms of hexagon functions.
In section 3.1 we describe the procedure we use to analytically continue the results of the resummation to the double scaling limit of the 2 → 4 Mandelstam region and then take the limit to (double scaled) multi-Regge kinematics. Having obtained a particular limit of the result in multi-Regge kinematics we then describe in section 3.2 how we can complete our expressions to restore the full kinematical dependence in multi-Regge kinematics by demanding that the result is single-valued. The analysis of this section is also supplemented by appendix B, where several two-particle bound states are computed and shown not to contribute to MRK.
Attached to the arXiv submission for this paper are files containing our results for the resummation of the single-particle gluon bound states up to five loops, particular contribu- The collinear limit is indicated by the arrows, and each term in the expansion around it may be mapped to an excitation of an integrable colour-electric flux tube, sourced by the two sides of W adjacent to the ones becoming collinear. tions from two-particle gluon bound states contributing to the double scaling limit and the new results at N 3 LLA and N 4 LLA for the completion of the five-loop remainder function in multi-Regge kinematics.

Preliminaries
An operator product expansion (OPE) for light-like Wilson loops was introduced in [23] and refined in many papers [24-27, 29, 39-42]. It describes the near-collinear regime of a particular ratio of light-like Wilson loops, denoted by W. Here we will focus on the hexagonal Wilson loop where the ratio takes the following form, Here W is the hexagonal Wilson loop, while W top and W bot correspond to the pentagonal Wilson loops and W to the square Wilson loop indicated in Fig. 1. The ratio W is finite and conformally invariant and is hence a function of the three available conformal cross-ratios, The OPE describes the ratio W in an expansion around the limit where the two adjacent edges x 56 and x 61 , of the hexagon become collinear so that x 2 15 vanishes. The limit is most conveniently parameterised by the variables {τ, φ, σ} which are related to the conformal cross-ratios as follows u 1 = 1 2 e 2σ+τ sechτ 1 + e 2σ + 2 e σ−τ cos φ + e −2τ , (2. 3) The limit τ → ∞ corresponds to the collinear limit. The quantity W becomes 1 in the τ → ∞ limit, up to exponentially suppressed corrections. The Wilson loop OPE gives a prediction for the form of the exponentially suppressed corrections, namely at l loops Here [x] denotes the integer part of x.
The basic idea behind the OPE is to express the bottom part of the Wilson loop as a coherent sum of excitations of the GKP string. These excitations then propagate to the top part of the Wilson loop where they are absorbed. Both the spectrum of excitations [28], which controls the propagation of states, and overlap functions [25], describing their production and absorption, can be studied at finite 't Hooft coupling using integrability.
The OPE expansion then has the following schematic form It is a sum over intermediate states ψ, weighted by the overlap functions for production and absorption (called 'pentagon transitions') P and a factor due the propagation from bottom to top involving the GKP energy (or 'twist') E ψ , momentum p ψ and helicity m ψ . For the purposes of this paper all the relevant quantities are available in the literature and we will describe them in greater detail in the following subsections. The ratio W is very simply related to the remainder function R 6 via and where the overall coefficient Γ cusp in the latter formula is the cusp anomalous dimension, whose expansion in the coupling is as follows, (2.9) In the next section we will describe the 'double-scaling' limit of W (and hence the remainder function R 6 ) which allows us to consider contributions to the sum over states (2.5) coming from gluons only.

The gluon contributions and the double-scaling limit
Apart from the collinear limit we described in the previous section, another kinematical limit which will relevant for our discussion is the so-called [24,29] 'double scaling limit': τ, iφ → ∞ , −τ + iφ fixed . (2.10) At the level of the OPE expansion (2.4) it amounts to the subset of contributions with p = 0 and cos mφ → e imφ , which are evidently the only ones surviving the limit. And as far as the cross-ratios are concerned, from (2.3) and (2.10) we deduce that namely the double scaling limit describes a two-dimensional subspace of the kinematics with vanishing u 2 , but general values for the remaining cross ratios.
From the point of view of the OPE, this limit is interesting because a restricted set of relatively simpler flux tube excitations contribute to it. In particular, these are the gluon excitations with positive helicity, given that all fundamental scalar, fermion and gluon excitations have helicity (the charge conjugate to φ in (2.5)) 0,±1/2 and ±1, all of them have twist 1 at weak coupling, and bound states are formed between gluons with the same helicity, with their charges being just the sum of the charges of their constituents. It is important to note that the twists receive corrections at each order in the coupling while the helicities do not. It is therefore important that, in restricting to the gluon contributions only, we are first expanding perturbatively in the 't Hooft coupling g 2 before taking the double scaling limit. That is, in the limit in question 1 where W BDS is obtained from (2.7) after substituting (2.11), and the all-plus gluon contribution is a sum over an arbitrary number of effective particles N , each of which may consist of any number of bound gluons a i , i = 1, . . . , N [29] , In the above formula n k counts the number of bound states made up from a i = k gluons, µ a (u) ≡ µ a (u)e −Ea(u)τ +ipa(u)σ+iaφ , (2.14) and the energy E a , momentum p a , and measure µ a of the a-th gluon bound state, as well as its pentagon transition to the b-th bound state P a|b , are reviewed in appendix A. As a final remark, notice that the weak coupling expansions (A.14)-(A.17) of all the different ingredients of (2.13)-(2.14), imply that the N -particle all-plus gluon state starts contributing at order O(g 2N 2 ). Namely one need only consider the states with N = 1 up to three loops, and N ≤ 2 up to 8 loops.

Resumming all single-particle gluon bound states
In this section, we will present a systematic procedure for perturbatively resumming the contribution of all single-particle gluon bound states to the hexagon Wilson loop/scattering amplitude, namely the N = 1 term in (2.13), which we may rewrite as This procedure is an extension of the method developed in [30,31] for the evaluation of the individual a = 1 and a = 2 terms above. After describing its details, we will apply it in order to obtain explicit for expressions up to l = 5 loops, and also deduce the relevant class of functions describing this contribution to arbitrary loop order. The reader interested in the final result may jump directly to the discussion around eqs. (2.29)-(2.30) and (2.32). We start by expanding the integrand in (2.15) with respect to g, following a well-known procedure that we review towards the end of appendix A. This reveals that up to u-and a-independent factors (e.g. e −τ or powers of τ and σ) which we can factor out, all terms are of the general form From (2.16), we can immediately infer that at any loop order the poles of the integrand are located at u = ±i(a/2 + k), k = 0, 1, 2, . . .. After we restrict to σ > 0, we can close the contour of integration on the u > 0 plane, and trade the integral for a sum over residues by virtue of Cauchy's residue theorem. With the help of the recurrence and reflection relations of the gamma function, as well as analogous equations for polygamma functions arising upon differentiation of (2.17), we arrive to sums over residues which always have the form 2 The next step is to reexpress the polygamma functions in terms of S-or Z-sums [32] via 3 , and replace the products of S-or Z-sums with the same outer summation index with linear combinations thereof, simply by nesting the independent summation ranges of each term in the product 4 . For example, where the first line yields S-sums, and the second line Z-sums.
For reasons that will become apparent very shortly, we will choose to replace ψ (m) (k) by S-and ψ (m) (a + k) by Z-sums respectively. After shifting the summation variable a → j = a + k, partial fractioning in k and expanding, (2.18) splits into terms that look like where we have combined the gamma functions of (2.18) into a binomial coefficient. Quite remarkably, the sum in k can be done for any collection of n i , with the help of algorithm C of [32]. Let us give a simple example where n 1 = 1 and the S-sum is absent, 2 The k = 0 residue is treated separately, as it is the only case where the denominators have poles, and yields simple sums which may be evaluated exactly as in [30]. 3 In particular, S(k; m; 1) = Z(k; m; 1) are the generalised harmonic numbers, ζm the Riemann zeta function, and γE = −ψ(1) = 0.577 the Euler-Mascheroni constant. 4 This is nothing but the quasi-shuffle algebra property of these objects.
so as to convey the basic idea of the algorithm, which is to appropriately manipulate the expression so that the sum in k can be done by means of the binomial theorem, For the example in question, this entails getting rid of the denominator, roughly speaking by rewriting the summand as the integral of its derivative with respect to x, n k=1 n k By changing the integration variable to y = 1 + x, the integrand becomes equal to the sum of the first n terms of a geometric series, so that n k=1 n k where in the last line we used the definition (2.20). By recursively going through the exact same steps, we may similarly obtain n k=1 n k In practice, the simplest way do any sum in k as in eq. (2.23), is to exploit an already existing C++ implementation of the algorithm, as part of the nestedsums library [33] within the GiNaC symbolic computation framework [43] 5 . The relevant command is transcendental_sum_type_C, and we have built an interface that calls it directly from Mathematica, which we use for the remaining manipulations. Combining it also with the Ssum_to_Zsum command, the evaluation returns Z-sums with outer summation index j − 1.
Then, after employing quasi-shuffle algebra relations generalising (2.22) in order to eliminate the resulting products of Z-sums in favour of their linear combinations, we may immediately evaluate the remaining sum in j in (2.23), in terms of multiple polylogarithms, see appendix C for more details on their definition and properties. The final result is most conveniently expressed in terms of the variables 28) and for example at one loop we have (we denote H m 1 ,...,m j (x) ≡ Li m 1 ,...,m j (x, 1, . . . , 1) the subclass of harmonic polylogarithms [46]), (2.29) and at two loops, As we mentioned at the beginning of this section, we have computed W (l) 1 up to l = 5 loops, and the remaining results may be found in the accompanying file allboundstates.m so as to avoid clutter.
We stress that the procedure we have described can be applied in principle at any loop order, and in fact we can make a precise statement about the particular class of multiple polylogarithms in which W (l) 1 lies: From (2.23) it is evident that the leftmost argument of the MPLs (2.27) will always be x. Furthermore, by inspecting the steps of algorithm C in [32], we observe that it can only generate Z-sums with arguments y, 1/y and 1, such that y and 1/y always appear in alternating order with y leftmost, if one removes all arguments equal to 1. Combining these two requirements, we infer that only MPLs of the form can appear in our results. The basis of MPLs (2.31) we have obtained may be expressed more transparently in terms of G-functions, whose definitions we review in appendix C. More specifically, the equivalent Li-and G-function representations of MPLs are related by (C.6), and it also proves advantageous to rescale all resulting G-functions so that the rightmost argument becomes x as in (2.28), by virtue of the identity (C.7).
In this manner, we finally arrive at the following important conclusion: Apart from explicit factors of σ = − 1 2 log(x(1 − y)) and τ ∼ − 1 2 log u 2 (in the double scaling limit), the OPE contribution of all single-particle bound states W (l) 1 is expressed in terms of G(a 1 , . . . , a n ; x) , a i ∈ {0, 1, 1/y} with a n = 0 , and a j = 1/y if a 1 = . . . = a j−1 = 0 , (2.32) at any loop order l. This basis of MPLs has dimension 3 n−1 at weight n 6 , and in fact it turns out to be part of the two-dimensional harmonic polylogarithms (2dHPL) of Gehrmann and Remiddi, which were first introduced in the computation of four-point functions of three on-shell and one off-shell leg at two loops [35].
In the next section, we will explore whether these results are consistent with the fundamental assumption of the hexagon amplitude bootstrap [14][15][16][17][18] for the class of functions describing six-particle scattering in general kinematics.
Calculating the contribution of multi-particle bound states to the double scaling limit is considerably harder than computing the single-particle bound state contributions. We will leave the resummation over such states to future work. However we can say something about the OPE expansion of certain contributions and this is presented in appendix B.

Hexagon functions in the double-scaling limit
An important finding of the previous section, is that a particular OPE contribution (2.15) to the hexagon remainder function R 6 (u 1 , u 2 , u 3 ) (2.6) is always expressed in the basis (2.32) at any order in the weak coupling expansion. More specifically, this contribution is the only one that survives in the double scaling limit (2.11) up to 3 loops, and is supplemented by the N = 2 particle gluonic contributions up to 8 loops, as we noted at the end of section 2.1. Given that R 6 is believed to lie in the space of hexagon functions for general kinematics, here we would like to explain why our results are indeed consistent with the double scaling limit u 2 → 0 of hexagon functions.
First of all we note that hexagon functions are just a special subset of all iterated integrals based on the nine-letter alphabet {u i , 1 − u i , y i } for i = 1, 2, 3 described in [14]. In other words they are a subset of polylogarithms on M 0,6 , the moduli space of six points on a Riemann sphere. It is simple to see that in the double scaling limit u 2 → 0 these reduce to iterated integrals on the five-letter alphabet 7 {u 1 , u 3 , 1 − u 1 , 1 − u 3 , 1 − u 1 − u 3 }, or in terms of the variables x and y defined above, on the alphabet {x, y, 1 − x, 1 − y, 1 − xy}. Thus the double scaling limit will give polylogarithms on M 0,5 . A basis for such polylogarithms is given in terms of products of iterated integrals of the form G(m y ; y)G(m x ; x) where the weight vector m y is made of elements drawn from the set {0, 1} and the weight vector m x is drawn from the set {0, 1, 1 y }. Here we have chosen a specific contour of integration over integrable words in the five letters such that one performs the integration along the y-axis first and then in the x-direction. We could of course have chosen to do it the other way round, yielding a basis with x and y exchanged.
Hexagon functions, as defined in [16], are a special subset of polylogarithms obeying restrictions on their branch cuts, namely that the locations of the branch points in the Euclidean region occur only at the boundary u i = 0 or u i = ∞. At the level of the symbols of hexagon functions this is reflected in the fact that the initial entries are only the u i and not any of the other nine letters. When taking the double scaling limit, the allowed branch point locations take the form In other words there should be no branch cuts around y = 0 or x = 1 at all and the branch cuts around y = 1 must match those around x = 0.
To make the branch cuts around x = 0 explicit we may use the shuffle relations to extract any trailing zeros from the weight vectors of the G-functions with argument x. In other words we reexpress the G-functions with argument x as linear combinations of Gfunctions which have no trailing zeros in their weight vectors and explicit factors of log x. We may then replace the explicit logarithms of x with logarithms of the product x(1 − y) at the cost of redefining the G-functions with argument y for each element of the basis. However once we have done so we are no longer allowed any further branch cuts at y = 0 or at y = 1. Thus we must have no further factors of the form G(m y ; y) at all.
Thus we find that a basis compatible with the constraints on the branch points at y = 0, 1 is given by all products of the form where there are no trailing zeros in the weight vector m x . Of course finding a function in the above basis does not necessarily satisfy the remaining constraint, namely that there should be no branch point at x = 1. Only special linear combinations of functions in the above basis will also satisfy this additional constraint. We have checked that the combinations appearing in the results derived above do indeed obey the property of having no branch points at x = 1.
Finally let us mention that, as described in (2.32), we never have a weight vector in the G-functions which begins with a string of zeros (of any length) followed by a 1 y . We have verified, by analysing the known four-loop remainder function [17], that this property holds also for the full R (4) 6 in the double scaling limit. In other words it holds not only for the single-particle gluon bound state contributions but also for the two-particle ones. This property appears to be an additional constraint on the form of the remainder in the double scaling limit which, at least in the case of the single particle states, we can observe as coming from the structure of the OPE.

The soft limit, analytic continuation and multi-Regge kinematics
Having obtained expressions for the single-particle gluon bound state contribution to the double scaling kinematics we may now ask about relevant physical limits of our expressions. Since our aim will be to obtain the six-particle scattering amplitude in multi-Regge kinematics we will briefly review the relevant kinematic limits. We recall that the soft limit and the multi-Regge limit are formally identical in that we send one cross-ratio (u 3 say) to one and the other two (u 1 and u 2 ) to zero so that the ratios are fixed. In the formulas above we have introduced the variables w and w * which parametrise the possible remaining dependence in the limit. In the soft limit we expect the remainder function to simply vanish. However, after we analytically continue the amplitude to the 2 → 4 Mandelstam region (obtained by continuing around the singularity at u 3 = 0 via u 3 → u 3 e −2πi ), the above limit is non-trivial and we expect a series of divergent logarithms whose coefficients are functions of the remaining variables w and w * , (3. 2) The limit (3.2) has been studied in many papers (see for example [48][49][50][51][52][53][54][55][56][57][58]). In our preceding discussion from section 2 we were considering the double-scaling limit where one of the cross-ratios, u 2 was taken to zero. This means that starting from the double scaling limit we actually only have access to multi-Regge kinematics in the regime where w * → 0 with w fixed (or, switching the helicities of the gluon bound states, to the regime w → 0 with w * fixed). We refer to this regime as double-scaled multi-Regge kinematics.
We may relate the variables w and w * to the variables σ, τ and φ via, Note also that the variable x introduced to parametrise the double scaling limit is related to w via x = −w and that therefore w * = −xe −2iφ and hence terms power suppressed in the double scaling limit correspond to terms power suppressed in w * . Now let us turn to the properties of the expressions we have obtained in the various physical limits discussed above. An obvious property, manifest from the form of the OPE expansion is that all W (l) vanish in the strict collinear limit τ → ∞. Likewise, the expressions we have obtained vanish in the soft limit x → 0 with y fixed. All of the G-functions of argument x are power suppressed as x → 0 and there is one such G in each term of our results.
The other soft limit y → −∞ with x fixed is related by the symmetry u 1 ↔ u 3 to the first. Our expressions also vanish in this limit. To see this however requires taking the limit y → −∞ on each of the G-functions appearing in our expressions. It is instructive to do this explicitly as we will need to perform exactly the same operation when taking the Regge limit.
In order to carry out the limit y → −∞ on our expressions it is first useful to express them in such a way that no G-functions are divergent in the limit. To arrange this we use the shuffle relations to unshuffle any zeros or 1 y appearing at the end of the weight vector so that each G-function either has a weight vector ending in a one or there are no ones on the weight vector at all. In the latter case we may rescale the argument by y so that these G-functions have weight vectors with entries 0 and 1 and argument xy (i.e. they are HPLs of argument xy). For example we write  The limit y → −∞ may now be taken straightforwardly using e.g. the HPL mathematica package [59] to give the limits of the G-functions with arguments xy as y → −∞. In this way we find that all our resummed single-particle bound state contributions to the double scaling limit vanish in the soft limit y → −∞.
The very same limit is required to analyse the multi-Regge limit of the double scaling limit of W. However to obtain a non-vanishing contribution we must first analytically continue to the Mandelstam region by passing round the branch but starting at u 3 = 0. This is easily achieved since all of the cuts at u 3 = 0 are manifest in the form of our result since they appear as explicit logarithms of x(1 − y) = u 3 /u 1 . We simply need to continue these so that log x(1 − y) → log x(1 − y) − 2πi. In so doing we obtain the result for W in the Mandelstam region in the double scaling limit.
Having analytically continued our results we may then go to multi-Regge kinematics by taking exactly the same soft limit as above, namely y → −∞ with x fixed. This time of course we obtain a non-vanishing result, dependent on x = −w. Finally, we express our results in a way that will be convenient for us to complete our expression in double scaled multi-Regge kinematics to the full multi-Regge kinematics. In order to do this we note that in the double scaling limit we have We may use these relations to rewrite all contributions in terms of the variables w, w * and (1 − u 3 ) which parametrise multi-Regge kinematics. In fact we also rewrite log(ww * ) = log(−w) + log(−w * ) and neglect the divergent logarithm log(−w * ). This is because these divergent logarithms are tied to the finite parts by the completion to single-valued polylogarithms that we describe in the next section.
The above discussion has focussed on the contribution of the single-particle gluon bound states to double-scaled multi-Regge kinematics. We have also analysed certain contributions from two-particle bound states, as described in appendix B. We find that after analytic continuation and taking the limit these contributions are all power suppressed. Hence, based on this preliminary analysis it seems that the single-particle states give the only non-vanishing contribution to double-scaled multi-Regge kinematics.

Completion to full multi-Regge kinematics from single-valuedness
Here we would like to explain why deriving the result in multi-Regge kinematics from the double scaling limit is sufficient to reconstruct the full result in multi-Regge kinematics. We know that in the multi-Regge limit of the remainder function in the 2 → 4 Mandelstam region we will obtain an expansion in powers of the divergent logarithm which we make take to be log(1 − u 3 ) 8 , The coefficients of the divergent logarithms are separated into an imaginary part g (l) r and a real part h (l) r . For each l and r both g (l) r and h (l) r are single-valued polylogarithms (or SVHPLs) [37].
The expressions we have obtained from the analytic continuation of the double-scaling limit allow us to obtain the limit of each g (l) r and h (l) r as w * → 0 with w held fixed 9 . We also discard any divergent logarithms of the form log w * as we take this limit, keeping only the finite term. From this data we may reconstruct the full dependence of all contributions to the real and imaginary parts in multi-Regge kinematics by invoking single-valuedness of the result.
In order to explain this point we briefly recall the construction of single-valued polylogarithms given in [60]. We begin with Knizhnik-Zamolodchikov equation in a single complex variable z, d dz Here e 0 and e 1 are two free non-commuting generators. We take the solution L 0 of (3.8) normalised so that as z → 0 we have L 0 (z) ∼ L an 0 (z)z e 0 where L an 0 (z) is analytic around z = 0 and L an 0 (0) = 1. The solution can be represented as a sum over all words in e 0 and e 1 with coefficients which are harmonic polylogarithms (i.e. regularised iterated integrals over d log z and d log(1 − z) away from z = 0), We may now consider the analytic continuation of the solution around the singular points at z = 0 and z = 1. Under an analytic continuation around z = 0 (via z → ze 2πi ), the monodromy of the solution L 0 is explicit, given the asymptotics as z → 0, (3.14) The monodromy of L 0 around z = 1 may be obtained by transporting to z = 1, where the monodromy of L 1 is explicit, and then back again, Now, to construct single-valued polylogarithms in a complex variable z, we consider a second 'primed' alphabet e 0 and e 1 and form the following series in all four letters e 0 , e 1 , e 0 e 1 , The second factor is built on the primed alphabet while the symbol '∼' means that the words in e 0 and e 1 are reversed with respect to the ones in e 0 and e 1 appearing in L 0 (z). Now the series L(z,z) will be single-valued if it has no monodromy around z = 0 or z = 1. The monodromy around z = 0 is given by (3.17) The series will be unchanged if e 0 = e 0 .

(3.18)
Similarly the monodromy around z = 1 is given by where Φ is built on the letters e 0 and e 1 . The series is unchanged if In other words, at a given weight, there are exactly as many single-valued poylogarithms in z andz as there are harmonic polylogarithms in a single variable z built on d log z and d log(1 − z). Given the asymptotics of L 0 defined above, the limit of each L m (z,z) asz → 0 (now treating z andz as independent complex variables), is a series in divergent logarithms logz with the coefficient of the finite term being simply H m (z). In other words there is a unique completion of a given harmonic polylogarithm H m (z) to a single-valued polylogarithm L m (z,z) such that the finite term in the limitz → 0 is H m (z). This fact was already noticed and found to be very useful in [61] in reconstructing the full form of the threeloop 'Easy' integral appearing in the correlation function of four stress-tensor multiplets in N = 4 super Yang-Mills theory.
Thus from the above discussion it is clear that to complete our expressions given in terms of HPLs with argument (−w) into single-valued expressions with the correct limit, we simply replace each instance of H with L. In this way we obtain the full expressions   r dependent on w and w * . We have verified that our results up to five loops reproduce the known expressions derived in [17,37]. In addition we have obtained the remaining five-loop terms g and g (5) 0 . We also obtain h (5) 0 but the real parts may be simply related to the imaginary ones so this does not constitute independent data. We give plots of the new data g (5) 1 and g (5) 0 along the diagonal line w = w * in Fig. 2 and Fig. 3.

Comparison with BFKL and the approach of Basso, Caron-Huot and Sever
From the reconstructed expressions derived in the previous section one can compare to the BFKL formula for the amplitude in multi-Regge kinematics, The formula expresses the amplitude in multi-Regge kinematics as a Fourier-Mellin transform of a factorised expression involving the BFKL eigenvalue ω(ν, n) and impact factor Φ reg (ν, n) which encode all the kinematical dependence of the amplitude. The other quantities in the above formula are given by (3.23) In principle, from our reconstructed expressions, we can then find Φ reg (ν, n) and ω(ν, n) so that they are consistent with the perturbative expansion of the amplitude in multi-Regge kinematics that we reconstructed from the double-scaling limit after analytic continuation. However, since an all-order form for these quantities was obtained in [38] following a different (though similar) logic, we can simply compare our expressions for the amplitude in multi-Regge kinematics with those obtained from their formula. Doing so we find perfect agreement up to five loops 10 .
Notice that the reconstruction argument of this section essentially amounts to the claim that the full BFKL expression can be reconstructed just from the knowing the residue at the first pole on the positive ν axis, assuming that the perturbative expression is given in terms of single-valued polylogarithms in w and w * . Certainly, under this assumption, the simplest way to evaluate any given expression perturbatively is simply to calculate the integral above around the first pole and resum the result as a single-variable function of w. Then one may complete this boundary information to a single-valued polylogarithm as described above.
It is very interesting that both here and in [38] the crucial ingredients were the singleparticle bound states of gluons appearing in the OPE expansion. We stress that the other states do contribute in multi-Regge kinematics. It is just that we are able to ignore the contributions from the fermions and scalars from the beginning by going to the double scaling limit and then, by extension of our observations, assume that the multi-particle gluon bound states drop out when going to double-scaled multi-Regge kinematics. The contributions of the missing states can then all be reconstructed by appealing to singlevaluedness. Note that the issue of identifying the n = 0 term in the sum (3.22) does not arise in our approach; it follows along with all the others when completing the double-scaled MRK expressions to the full MRK ones. If we were able to fully justify the fact that the multi-particle gluon bound states drop out from double-scaled MRK we would have proven that the full expression can be reconstructed from the single-particle gluon bound states alone. As the authors of [38] stress, very similar assumptions have to be made also in their approach. It would be very interesting indeed if some combination of the two arguments could be used to fully prove the BFKL formula for the hexagon in MRK.

A Dispersion relation, measure and pentagon transitions for gluons
Here we review the building blocks of the Wilson loop OPE, coming from all gluonic excitations of the corresponding flux tube [29], in slightly adjusted notation. Although these formulas hold at finite coupling, they are particularly suited for the weak coupling expansion, as we will comment below, and display to leading order at the end of this section.
The energy and the momentum of the a-th gluon bound state, a ∈ Z * with a > 0 corresponding to positive and a < 0 to negative helicity gluons, are given by where Q is a matrix with elements Q ij = δ ij (−1) i+1 i, M is related to another matrix K, J i is the i-th Bessel function of the first kind, and κ,κ are vectors with elements The subscripts in (A.1) denote that we are only taking the first component of the vector inside the brackets. Moving on, the measure of the a-th gluon bound state, and its pentagon transition to a b-th bound state are given by and in addition where for the last three formulas we have resticted a, b > 0, and under the same restriction the remaining cases may be obtained by Finally, Although the vectors (A.3) and matrices (A.3) are infinite-dimensional, and hence the matrix products (A.1) and (A.6) should really be thought of as infinite-dimensional sums at finite coupling, it is possible to show that at weak coupling K ij starts at order O(g i+j ).
Thus if we wish to obtain the weak coupling expansion up to order O(g 2l ), we may simply truncate all summations up to their first i, j = 1, 2 . . . , 2l − 1 terms.
In order to see this, and also perform the expansion in practice, one starts with the Taylor series of the Bessel functions, and then computes integrals in t by means of where ψ (m) (z) is the polygamma function, For example, the leading order expressions for all necessary ingredients of the Wilson loop OPE that we obtain in this manner, are for a, b > 0 More details about the weak coupling expansion may be found in appendix E of [26], as well as appendix A.2 of [29]).
B Remarks on the two-particle gluon bound state contributions Although we leave the resummation of all two-particle gluon bound states as an exciting open question for future work, in this appendix we shall employ the technology we have developed so far in order to compute the contribution of several such states in the collinear limit (part B.1). The main motivation here is to examine whether these states survive when we transition from collinear to multi-Regge kinematics, and complement the analysis of section 3 (part B.2) . We gather evidence up to 5 loops that this is not the case, thus allowing us to reconstruct the hexagon remainder function in the latter kinematics from the single-particle gluon bound states, resummed in subsection 2.3.

B.1 Evaluation of the integrals
The two-particle gluon bound state OPE contribution corresponds to keeping only the N = 2 term in eq. (2.13), which we may rewrite as with a, b the number of gluons contained in each of the two bound states (also equal to their twist and helicity), δ ab a Kronecker delta function, and all ingredients of the OPE integral on the right-hand side contained in appendix A. With the help of (A.14)-(A.17), at 4 loops the integral becomes where we have replaced four of the gamma functions coming from (A.17) with the help of the identity which may be derived from (2.17). For specific values of a, b, and depending on whether a/2, b/2 and (a + b)/2 are integer or half-integer, we may similarly replace the remaining gamma functions by virtue of , for n integer. Finally, we may factorise all u− and v−dependence coming from the trigonometric/hyperbolic functions, e.g.
From these considerations, it is clear that the only inseparable dependence on u − v will come from the products in (B.4) appearing in the denominator of (B.1). Given that at higher loops the integrand (B.1) is dressed by sums of products of polygamma functions with arguments 1 + a 2 ± u, a 2 ± u, 1 + b 2 ± v and b 2 ± v, this statement holds true independently of the loop order.
As was first observed in [31], for the case a = b = 1 these (u − v)-dependent terms in the denominator cancel out, and therefore W (l) 2 [1,1] always reduces to a sum of two factorised 1-fold integrals, which may be evaluated along the lines of [30]. Here, we find that exactly the same phenomenon occurs also for a = 2 and b = 1.
Although for a + b ≥ 4 we start dealing with genuine 2-fold integrals, these can still be done by summing over residues, at least at 4 loops. The main subtlety is that we need to treat separately not only the poles at u = ia/2 and v = ib/2, but also the ones where These special poles will lead to simple sums, whereas the remaining ones give rise to double sums, which may be brought to a form similar to the one we encountered in 2.3. Sparing the reader the rest of the details, all in all we arrive at expressions for W (l) 2 [1,1] and W (l) 2 [2,1] up to l = 5 loops, as well W 2 [3,1] and W 2 [2,2] , which are included in the accompanying file boundstatesN2.m.
Before moving on to examine the contribution of these states in MRK, let us end this section by noting a peculiar property of the two twist-4 contributions: Each one separately gives rise to polylogarithms with symbol letters S, 1+S 2 but also 1−S 2 , namely HPLs with both positive and negative weights. At first this may seem to contradict the expectation discussed in section 2.4, that perturbatively R 6 , and hence also W, are described by a particular 9-letter alphabet, since the latter only reduces to the letters S and 1 + S 2 in the collinear limit. Quite remarkably however, in their sum W 2 [2,2] /2, also taking into account the symmetry factor in B.1, all terms containing the letter 1 − S 2 cancel out, and consistency is restored. This perhaps suggests that it is only the combination of states with the same charges and particle number N that has a physical significance in the OPE.

B.2 Contribution to MRK
In section 3.1, we reviewed the analytic continuation which is necessary for obtaining a nontrivial result in MRK, and specialised its form in a regime where it overlaps with the double scaling limit, where our resummed single-particle gluon contribution lives. For the individual two-particle gluon bound states that we computed in the previous part of this appendix, we will need a similar specialisation to a region overlapping with the collinear limit. This was first considered in [62], based on the initial Wilson loop OPE approach [23]. Following the refinement of the latter approach, this procedure has been revisited and extended in [36], which we now briefly review.
In [36] it was also observed that at two loops, the analytic continuation for S > 1 and the collinear limit expansion commute. In other words, first analytically continuing R (l) 6 and then expanding the near collinear limit yields the same result as first expanding, and then analytically continuing term by term. This behaviour was then conjectured to hold to all loops, and in fact it can be proven for all polylogarithms on M 0,6 . In other words if one assumes that R (l) 6 is described by hexagon functions, then the commuting of the two procedures follows.
Taking the aforementioned property as granted, it is then easy to obtain the collinear limit expansion of the analytically continued R (l) 6 . For all (sums of) positive helicity gluon states with given particle number and twist a that we have considered (as discussed in the last paragraph of section B.1), we find that their kinematical dependence is always a linear combination of terms of the form For the rational part in e iφ and S, clearly the value at the endpoint of the analytic continuation does not depend on the path, so that we can simply set e iφ → −e iφ , S → −S, and in fact all extra minus signs always cancel out. The HPLs with positive weights also remain unchanged, since they only have branch cuts for their argument between 1 and infinity, and the path (B.7) never crosses them for S > 0. Therefore analytically continuing our expressions amounts to the almost trivial replacement log S → log S + iπ , or σ → σ + iπ . (B.9) The final step is to pass from the analytically continued collinear to multi-Regge kinematics.
As was first noticed in [16], and can be readily verified from (2.3) and (3.1), in the variables we're currently using MRK corresponds to T → 0, S → 0 with r = T /S fixed. The ratio r has already appeared in (3.3), and it is in fact via this procedure that the relation with the w, w * variables of the MRK appearing on the left-hand side of the latter formula can be established. We thus arrive at the following connection between the two kinematics regions: Starting with the analytically continued collinear limit expansion (around T = 0), expanding additionally around S = 0 and replacing T = rS lands us on the multi-Regge limit expansion (around S = 0), where we are additionally expanding around r = 0. It is specifically the non-power suppressed O(S 0 ) term that has principally been the focus of BFKL analysis, and also of our paper here.
For that reason, we will take the strict S → 0 limit of our expressions (up to large logarithms), under which we see that, very interestingly, all but the k = 0 part of all terms (B.8) already drops out. Further examining the k = 0 term for W 2 [2,2] /2, we see that it vanishes in all cases. We deem this as very compelling evidence that 2-particle gluon bound states do not contribute in MRK in general, since this behaviour remain unchanged across loop orders, and also for qualitatively quite different types of integrals.

C Multiple Polylogarithms
In this appendix we review the definitions and several properties of multiple polylogarithms (MPLs), which will be useful throughout the main text. More detailed expositions may be found for example in [63][64][65].
Multiple polylogarithms may be defined as nested sums 11 Li m 1 ,...,m k (x 1 , . . . , x k ) ≡ The nested sum definition (C.1) is more suited for the evaluation of Mellin-Bernes type integrals by residues, as we did in section 2. However there also exists an alternative, global definition of MPLs in terms of (regularised) iterated integrals, as