A Novel Algorithm for Nested Summation and Hypergeometric Expansions

We consider a class of sums over products of Z-sums whose arguments differ by a symbolic integer. Such sums appear, for instance, in the expansion of Gauss hypergeometric functions around integer indices that depend on a symbolic parameter. We present a telescopic algorithm for efficiently converting these sums into generalized polylogarithms, Z-sums, and cyclotomic harmonic sums for generic values of this parameter. This algorithm is illustrated by computing the double pentaladder integrals through ten loops, and a family of massive self-energy diagrams through $O(\epsilon^6)$ in dimensional regularization. We also outline the general telescopic strategy of this algorithm, which we anticipate can be applied to other classes of sums.


Introduction
Hypergeometric functions appear ubiquitously in physics-including as they do special cases such as Legendre polynomials and Bessel functions-and in particular are known to appear in dimensionally-regularized perturbative quantum field theory.This ubiquity is due in part to the extensive class of second-order differential equations that hypergeometric functions solve, which has been the subject of much dedicated research.While many aspects of these functions are well understood as a result of this research, they can still prove unwieldy in practical calculations.
Generalized (or multiple) polylogarithms also appear in many places in quantum field theory.While they are less general than hypergeometric functions, they are correspondingly under better control; in particular, a great deal of progress has been made leveraging motivic aspects of these functions, considered as iterated integrals on the moduli space of the Riemann sphere with marked points [1][2][3][4].For instance, all functional relations between generalized polylogarithms can in principle be exploited with the use of the coaction [5,6], as applied in [7][8][9][10].There also exist public codes for the efficient numerical evaluation of polylogarithms, for instance in GiNaC [11,12].
In some cases, the hypergeometric functions that appear in physics can be expressed as infinite sums of polylogarithms, allowing us to leverage this polylogarithmic technology.This occurs, for instance, in the case of one-loop dimensionally regulated Feynman integrals [13,14].For finite values of the dimensional regularization parameter , these integrals can be expressed in terms of hypergeometric functions, while their series expansion around small values of can be written in terms of polylogarithms.While the sum representation of hypergeometric functions is a useful starting point for performing this expansion, explicitly converting the expansion coefficients into polylogarithms can prove to be a nontrivial task.Great progress in this respect has been made in [15], where several summation algorithms, covering a large number of cases, were presented.Even so, these algorithms are not always capable of carrying out the series expansion around symbolic integer indices of hypergeometric functions. 1In the context of dimensionally-regularized Feynman integrals these symbolic integers essentially correspond to propagators raised to generic powers.
In this paper, we extend the work of [15] by presenting an algorithm for explicitly evaluating series expansions of Gauss hypergeometric functions taking the form in terms of well-known classes of functions.Importantly, since we place only two conditions on the integer parts of the hypergeometric indices, (1.1) is a function of an arbitrary symbolic integer, which we will generally denote α, on top of the complex argument x.Hence, our expansion may prove especially useful when hypergeometric functions of the form (1.1) appear as part of a larger expression in which the symbolic integer α is summed over, or simply because this expansion can be evaluated a single time and then used for different values of α. 2  More generally, we show how sums of the form N n=1 x n n p (n + α) q Z m 1 ,...,m d (n−1|y 1 , . . ., y d )Z r 1 ,...,r h (n+α−1|z 1 , . . ., z h ) , (1.2) 1 Similarly, although there are several algorithms for treating sums that depend on symbolic integer parameters [16][17][18], in general these are not applicable when the summand depends on continuous parameters, such as the arguments of hypergeometric functions. 2It is also well known that all 2F1 functions whose arguments differ by integer shifts may be expressed in terms of a basis of two functions in this family.In this sense, our algorithm could be thought of as the explicit expansion not only of the basis, but also of the coefficients of any other function in the family.
can be efficiently evaluated for symbolic values of {x, y 1 , . . ., y d , z 1 , . . ., z h } ∈ C and {N, α} ∈ Z ≥0 , where (as we will review in section 2) Z-sums are given by Z m 1 ,...,m d (N |x 1 , . . ., x d ) = (1.3) In particular, the sums in (1.2) evaluate to generalized polylogarithms and Z-sums when N is infinite, and cyclotomic harmonic sums [19] for generic N .The sums appearing in the expansion (1.1) correspond to the special case in which x i = y i = m i = 1 and N → ∞.
We again highlight that the novelty of our algorithm in comparison with the algorithms of [15] is that we allow α to be an arbitrary symbolic integer.To achieve this generalization, we rely on telescoping identities.In their simplest guise, these take the form and can be used to compute Φ(α) if the quantity ∆Φ(α) is simpler than the original sum.We will make use of a generalized identity of the form (1.4), in which the analog of ∆Φ(α) has lower depth-that is, involves fewer nested sums-than the original sum.This will allow us to leverage a recursion in the depth of sums taking the form (1.2), and in this manner express them in terms of nested sums that depend on α and the other symbolic parameters.
In principle, the algorithm of [20] is also capable of converting the sums in (1.2) into nested sums for generic values of α.Nevertheless, we find that this algorithm proves overly computationally expensive in the cases we consider, due to the fact that it generates spurious term-wise divergences at intermediate steps in the calculation.These divergences require regularization, and cancel out in the final result.For sums of the form (1.2), we therefore view our approach as simpler and more efficient.
In order to illustrate the utility of our algorithm, we apply it to two examples.First, we consider the class of double pentaladder integrals introduced in [21].A compact generating function for these integrals was derived in [22], which gives rise to a sum representation involving products of 2 F 1 functions that fall into the class (1.1).Using our algorithm, we explicitly evaluate these integrals in terms of generalized polylogarithms through ten loops.Second, we consider dimensionally-regularized one-loop self-energy diagrams for generic masses and propagator powers.In the limit of zero external momentum, these integrals become expressible in terms of 2 F 1 functions [14], and our algorithm can be used to simultaneously expand families of integrals that have different propagator powers around 4 − 2 dimensions.We carry out the expansion and resummation of one such family of self-energy diagrams through O( 6).
This paper is organized as follows.In section 2, we begin by reviewing the aspects of Z-sums and generalized polylogarithms that will be relevant for our analysis.We then initiate the series expansion of Gauss hypergeometric functions taking the form (1.1), and deduce the types of nested sums they give rise to.In section 3, we present the general strategy of our algorithm, and state the recursion it gives rise to for sums of the form (1.2).Finally, we derive this recursion in a simplified setting, for compactness of presentation.Sections 4 and 5 deal with the applications of our algorithm to the double pentaladder integrals and the massive self-energy integrals, respectively.In section 6 we conclude.We also include two appendices, in which we derive the most general form of our algorithm and present an example involving cyclotomic harmonic sums.
2 Z-sums, Polylogarithms, and 2 F 1 Functions We begin by reviewing the types of sums that arise when 2 F 1 functions are expanded around integer values of their indices.Although the coefficients of these expansions are known to be expressible in terms of generalized polylogarithms around fixed integer values [37], they can evaluate to the (more general) class of Z-sums [15] around generic symbolic integers.We review this class of sums, and then describe how they relate to generalized polylogarithms and 2 F 1 functions.

Z-sums
We give just a brief review of Z-sums, introducing notation and recalling the properties that will prove useful in later sections.The reader is referred to [15] for more details.Starting from any integer N and the initial definition Z-sums are defined recursively in terms of pairs of variables m i ∈ Z + and 3) The depth of each Z-sum is defined to be its number of summation indices d, and its weight is defined to be w = d i=1 m i .In general, we will adopt an abbreviated notation in which bold characters indicate multi-indices, for instance which leaves the depth of each sum implicit.We will also adopt the notation that the first entry of primed multi-indices has been dropped, for instance m = m 2 , . . ., m d , and denote the number of indices in m by |m|.
The Z-sums obey a stuffle algebra as a consequence of the ability to split up unordered sums into nested ones.More precisely, a product of Z-sums that have the same upper summation limit N (but which don't necessarily have the same depth) may be expressed in terms of Z-sums of higher depth by recursively applying until the depth of the original Z-sums has been reduced to zero (after which new sums are built up using (2.2)).For example, we can reexpress It is also interesting to note that this stuffle algebra has an associated coalgebra, and that together these structures form a Hopf algebra; however, the associated coproduct is not the coproduct usually encountered in Feynman integral calculations, which is associated with their mixed Tate Hodge structure [4] (see for instance [38]).
Identities also exist between sums with different summation bounds.For instance, relations between Z-sums with different upper bounds follow directly from (2.2), namely for M, N ∈ Z + and |m| > 0. Equation (2.7) is particularly useful for synchronizing a product of two Z-sums with different upper summation bounds, as it can be applied (iteratively) to replace one of the Z-sums with a linear combination of (products of) Z-sums with shifted summation bounds.
Additionally, the sums one encounters often have different lower summation bounds than allowed in (2.2); in particular, we will see below that expansions of gamma functions that appear in the denominator of hypergeometric functions more naturally give rise to S-sums, which are closely related to Z-sums [15].Specifically, these sums satisfy which can be compared to (2.3).The S-sums have similar algebraic properties to the Z-sums, but we prefer to work in terms of Z-sums as they are more directly related to generalized polylogarithms.S-sums can be converted to Z-sums iteratively using which separates out the contribution to these sums in which pairs of summation indices are equal; this leaves summation bounds that fit the definition of Z-sums, at the cost of introducing a sum of lower depth (making it clear that this recursion terminates).Notice that we have used double primes to indicate dropping the first two indices in a multi-indix (for example x = x 3 , . . ., x d ).
An important sub-class of the Z-sums are the Euler-Zagier sums [39], which occur when all x i = 1; similarly, S-sums with all arguments evaluated at unity reduce to the harmonic sums [36].We abbreviate both cases using S m (N ) ≡ S m (N |1, . . ., 1) . (2.11) These sums appear naturally in the expansion of the gamma function and its reciprocal, and therefore also in the expansion of hypergeometric functions.

Generalized Polylogarithms
Many of the quantities encountered in quantum field theory can be expressed entirely in terms of generalized polylogarithms.This corresponds to expressing these quantities in terms of Z-sums in which N → ∞, as lim reproduces the normal definition of generalized polylogarithms (note the reversal of indices and arguments) [15,35,[40][41][42][43].The depth and (transcendental) weight of generalized polylogarithms coincide with the definitions they were given above as Z-sums.Generalized polylogarithms form a closed subalgebra within the Z-sums, and thereby inherit the algebraic properties of the larger space.Moreover, they can be given an integral definition where z, x i ∈ C, and where the second definition accounts for the cases in which the first n indices are zero (as the general definition diverges in these cases).This definition is related to the sum definition by , . . ., 0, . . ., 0 and allows polylogarithms to be analytically continued outside of the region of convergence |x i | < 1.This representation also gives rise to a new set of identities analogous to the stuffle identities, corresponding to the ability to triangulate unordered integration ranges (coming from products of polylogarithms) into sums over iterated integrals; these go by the name of shuffle identities.We refer the interested reader to [44] for further details.
In some of our examples, we will also make use of the harmonic polylogarithms (HPLs) [35].These are functions of a single argument z, and correspond to restricting x i ∈ {0, 1, −1} in (2.13).The standard notation is given by where p is the number of indices x i that are 1.It can be useful to express functions in terms of HPLs when possible, due to the existence of several dedicated packages for their analytic and numerical evaluation [45][46][47][48].

Gauss Hypergeometric Function Expansions
Having introduced some of the necessary machinery in the previous subsections, let us now proceed with the expansion of the hypergeometric function (1.1) mentioned in the introduction.Our starting point will be the series definition of the Gauss hypergeometric function, This sum converges for |x| < 1, as can be seen by taking the n → ∞ limit of consecutive terms in the series.
We will be interested in expanding the indices of Gauss hypergeometric functions around for integer k, l, and m.To this end, we make use of the identity When we need to expand gamma functions in the denominator, we also employ so as to place all nested sums in the numerator.With the help of these replacements, it is easy to see that the factor outside of the sum over n in (2.16) readily evaluates to Z-and S-sums in the limit (2.17) (note also that the factors of Γ(1 + ) will cancel out of the overall expression).Therefore, the only nontrivial terms requiring evaluation in (2.16) take the form in the expansion.Ideally, we would like to be able to evaluate sums of the type (2.20) for general, symbolic values of all three integers k, l, and m.However, the presence of gamma functions in the above formula places significant obstacles on the telescoping methods mentioned in the introduction.For this reason, in what follows we specialize to a one-dimensional subspace of this 3D lattice, such that the aforementioned gamma functions reduce to polynomials in n.That is, we impose the constraint that k − m and l are fixed integers , which in turn allows us to replace (where if k > m we simply exchange k and m) and similarly for k → l, m → 1. 3 If this replacement results in a polynomial in n in the numerator, we may use the differential operator introduced in [15], in order to reexpress each monomial in n as With this replacement, we can first evaluate the sum on the right-hand side, and then differentiate the result to evaluate the original sum.The differentiation of nested sums has already been implemented in software such as HarmonicSums [19,[23][24][25][26][27][28][29][30][31][32][33][34][35][36].
After transforming the S-sums in (2.20) to Z-sums with the help of (2.9), applying (2.7) to synchronize them, making the replacements (2.22), and finally partial fractioning the denominator with respect to n, we find that only a single nontrivial class of sums needs to be evaluated: where α is the symbolic integer that remains after imposing the constraint (2.21).One may readily check that these sums are a special case of eq.(1.3) for x i = y i = m i = 1, and N → ∞.
In the next section, we present an algorithm for evaluating the more general class of sums, and therefore also for evaluating the sums (2.25) that appear in the 2 F 1 expansions we are considering.

A Telescopic Nested Summation Algorithm
We begin this section by outlining the main idea behind our nested summation algorithm.The reader interested in its statement and the outline of its proof may jump directly to subsections 3.2 and 3.3, respectively.

Strategy of the Algorithm
Let Φ(α|x) be a sum of depth d, which depends on a continuous variable x ∈ C and a symbolic, non-negative integer α.Our objective is to find a systematic cancellation between the terms in this sum that either decreases the summation range or the sum depth.
Inspired by the telescoping identity (1.4), we begin by writing the following ansatz for this sum: where P is known while Q(µ) and R are to be determined.Expanding this ansatz and inspecting the coefficients of Φ(µ|x) for different values of µ, we see that we must have for (3.1) to be consistent.Solving these constraints and plugging the solution back into (3.1),we arrive at where we have defined The form of (3.5) is analogous to equation ( 17) in [20], albeit much less general.However, what we lose in generality we gain in simplicity and computational efficiency.The strategy is then to determine whether ∆Φ(µ|x) can be reduced to simpler sums than the original Φ(α|x).This requires analyzing the specific form of the sums under consideration; in many cases, no reduction will occur (although in some of these cases, a different form of the generalized telescoping identity may work).In the class of sums we focus on in this paper, we will see that the sum ∆Φ(µ|x) has lower depth than Φ(µ|x), allowing us to telescopically recurse.We now turn to that analysis.

Statement of the Algorithm
Having motivated sums of the form (2.25) by considering the expansion of Gauss hypergeometric functions, we now broaden the scope of our analysis to the more general class of sums4 S p,q m;r (α, N |x; y; z) where {x, x 1 , . . ., x d , y 1 , . . ., y h } ∈ C and {N, α} ∈ Z ≥0 are all allowed to be symbolic.We here present an algorithm for converting any sum S p,q m;r (α, N |x; y; z) with fixed values for {p, q} ∈ Z ≥0 and {m 1 , . . ., m d , r 1 , . . ., r h } ∈ Z + into nested sums that depend on the remaining symbolic parameters.
Our algorithm takes the form of a recursion, in which the application of the relation decreases the overall depth of the Z-sums appearing on the right-hand side of equation (3.7) in each iteration.The boundary term appears at each step, but can be converted into Z-sums using the techniques of [15].Equation (3.8) can be applied until the depth of both Z-sums are zero (using the prescription that Z-sums of negative depth vanish), whereby we are left with the sum We then proceed to convert the sums over µ that appeared at each step in the recursion into nested sums.When N is taken to be infinite, these sums can be evaluated as Z-sums, again using the methods of [15].For generic finite N , these sums evaluate to the more general class of cyclotomic harmonic sums [19].Once these sums have been converted into nested sums, we are left with sums that can be carried out explicitly for fixed values of p and q.
Altogether, this recursion converts the sum over n in (3.7) into a linear combination of Z-sums or cyclotomic harmonic sums with coefficients that depend on α, N , x, y, and z.
In the rest of this section, as well as in sections 4 and 5, we will focus on sums for which N is infinite; we provide more details for sums with generic N in the appendix.For the cases we focus on here, we abbreviate S p,q m;r (α|x; y; z) ≡ lim V p,q m;r (α|x; y; z) ≡ lim A number of simplifications occur in this limit.For instance, the last line in (3.9) vanishes, and many of the Z-sums that appear can be converted into generalized polylogarithms using (2.12).However, there always remain Z-sums that cannot be expressed as polylogarithms, as their upper summation bound only depends on α.This can be seen, for instance, in the the terminal sum (3.10), which becomes in this limit.For any specific choice of α, these remaining Z-sums evaluate to rational functions of their arguments.
(3.14) and then At this point, the recursion terminates in an expression of the form (3.13), and we need to convert the sums over ν and µ into Z-sums.In general, this involves some combination of partial fractioning, reindexing sums, and applying (2.7) to shift the upper bound of existing Z-sums.Applying these techniques, it is not hard to show that and that the boundary contribution is given by Using these results, the sum over µ in (3.14) can be carried out using the same techniques, whereby we find S 0,1 1;1 (α|x; To evaluate the boundary contribution V 0,1 1;1 (α|x; y 1 ; z 1 ) we need to be slightly more careful.Evaluating the sums over i in (3.9) for these indices and arguments, we find where we have additionally shifted the summation index n → n−1 to put the denominator in a form that fits the definition of Z-sums.Since the summand on the right hand side is zero when n = 1, we can change the lower summation bound back to 1.To increase the upper summation bound of Z m (n−2), we would like to use equation (2.7); however, we need to emend this relation so that it remains valid when n = 1.Doing so, we have where θ(k) is the Heaviside function, which is equal to 1 when k ≥ 0 and 0 otherwise.The Heaviside function makes clear that this term must vanish when n = 1 (since all other terms in (3.20) vanish), despite the ambiguity of both its numerator and denominator evaluating to zero.Substituting this relation into (3.19) and converting the expression into Z-sums like above, one finds Putting this all together, we find As emphasized above, this expression depends on not only polylogarithms but also Z-sums with upper summation bound α.However, for any specific value of α, this expression reduces to a linear combination of generalized polylogarithms with rational coefficients that depend on x, y 1 , and z 1 .
More complicated sums of this class can be evaluated following the same strategy, but require an increasing amount of algebra and quickly become tedious.In the ancillary files included with this paper, we provide a Mathematica package that can be used to apply this algorithm to any sum in the class (3.7), up to computer memory and time limitations.

Proof for Euler-Zagier Sums
In this section we illustrate the derivation of (3.8) by proving it in the N → ∞ limit for the case of Euler-Zagier sums, which corresponds to the restriction y = 1, . . ., 1 and z = 1, . . ., 1 in (3.7).This simplifies the notational clutter without altering the telescopic strategy, which can also be applied in the general case.We prove the more general result in appendix A.
We abbreviate the sums we focus on in this section by S p,q m;r (α|x) ≡ S p,q m;r (α|x; 1, . . ., 1; 1, . . ., 1) x n n p (n + α) q Z m (n−1)Z r (n+α−1) . (3.24) A generic sum of this form can be split into a linear combination of cases in which either p or q is zero by partial fractioning: Our strategy will be to find telescopic recursions for S p,0 m;r (α|x) and S 0,q m;r (α|x) separately, after which the case of general p and q can be treated by plugging these recursive formulas into (3.25).
Telescoping S p,0 m;r (α|x) To find a telescopic recursion on S p,0 m;r (α|x), we consider the shifted sum x n n p (n + α) r 1 Z m (n−1)Z r (n+α−1) = S p,0 m;r (α|x) + S p,r 1 m;r (α|x) , where in the second line we have used identity (2.7) to shift the upper summation index of Z r (n + α).In principle, we should treat the |r| = 0 case separately, since (2.7) can't be applied to depth-zero sums; however, it is easy to see that S p,0 m;∅ (α|x) is independent of α, and that (3.28) correspondingly gives the correct answer as long as we adopt the prescription that Z-sums with negative depth evaluate to zero.
Comparing to (3.6), we see that the parameter P in our ansatz is in this case zero, and we have ∆S p,0 m;r (α|x) ≡ S p,0 m;r (α + 1|x) − S p,0 m;r (α|x) (3.29) = S p,r 1 m;r (α|x) .Importantly, the sums on the right hand side of this equation are all strictly simpler than the original; the terms in the sum over µ all involve a Z-sum of one lower depth, and the last term does not depend on α.Note, however, that the terms in the sum over µ no longer satisfy q = 0; thus, while we can again split these sums up into cases in which either p or q is zero as in (3.25), we will need to be able to reduce the depth of sums with nonzero q in order to achieve a genuine recursion in depth.
Telescoping S 0,q m;r (α|x) Following the same strategy as above, we consider where we have shifted the summation index by n → n−1 relative to the definition (3.24).Since the summand on the right hand side of (3.32) is zero when n = 1, we can change the lower summation bound back to 1.To increase the upper summation bound of Z m (n−2), we again use equation (3.20).This gives us x S 0,q m;r (α|x) − S m 1 ,q m ;r (α + 1|x) − δ |m|,0 (α + 1) q Z r (α) . (3.34) where we have shifted n → n+1 in the sum involving the Heaviside function to get the expression into this form.
Comparing to (3.6), we see that P = −1 and that ∆S 0,q m;r (α|x) ≡ S 0,q m;r (α + 1|x) − Therefore, our ansatz (3.5) gives the relation where we have already converted one of the sums over µ into a Z-sum.As in (3.31), all the terms on the right hand side are simpler than the original sum; the sums in the first line either have one lower depth or don't depend on α, and the terms in the second line no longer involve a sum over n.
The Closed Recursion for S p,q m;r (α|x) Substituting equations (3.31) and (3.37) into (3.25),we obtain the recursion relation which is the simplified version of (3.8) that applies to Euler-Zagier sums in cases where N → ∞.We have collected all the boundary terms into Like the more general recursion (3.8), the relation (3.38) can be applied iteratively until the depth of both Z-sums are zero, terminating in the sum (3.10).
4 Application I: The Double Pentaladder Integrals Ω (L) Let us now explore the power of our algorithm by considering some applications.Our first example is a family of integrals first introduced in [21] in the context of N = 4 supersymmetric Yang-Mills theory-the double pentaladder, or Ω (L) , integrals.These integrals consist of (L−2)-loop box ladders capped on either end by a pentagon loop that comes with a numerator factor.These diagrams are depicted in figure 1.Their numerator renders them infrared finite and parity even, and they depend on the kinematic variables where  In [49,50], it was shown that Ω (L) is related to Ω (L−1) by a second-order differential equation.This differential equation was solved for generic values of the coupling g 2 in [22] by where is a hypergeometric function that has been normalized by the factor such that F ν (1) = 1.Note that F ν (x) depends on the coupling g 2 , although we have left this dependence implicit.By virtue of Cauchy's residue theorem, the generating function (4.7) can put in an equivalent sum representation Expanding this expression around small coupling with the use of (2.16), one can derive a sum representation for each of the perturbative double pentaladder integrals Ω (L) .As noted in [22], it is advantageous to carry out this expansion in two steps: first with respect to the small parameter and then by expanding with respect to g.In [22], the resulting sum representation for Ω (L)  was evaluated in terms of generalized polylogarithms in the u → 1 and w → 0 limits through eight loops.However, the techniques used there were not sufficiently powerful to evaluate the sum in general kinematics.Using (4.11) to replace g with , we see that the hypergeometric function in takes the form (1.1), and thus falls into the class of hypergeometric expansions our algorithm can handle.Leveraging this fact, we now proceed to evaluate the sum representation of Ω (L) in general kinematics.

Evaluating the hypergeometric function
We begin by converting the building blocks F −iα (x) into Z-sums.Using the techniques discussed in subsection 2.3, it is possible to express these functions as where we have denoted the relevant specialization of the sums (3.24) by in order to avoid notational clutter.We highlight that in (4.13) implicitly depends on α and g 2 , as per the definition (4.11).The α → 0 case of (4.13) may be obtained by taking the smooth limit → ig, as can be seen from (4.11).In this case, the sum over n may be readily evaluated with the techniques of [15], where after using stuffle relations to combine the product of Z-sums it can be evaluated in terms of HPLs, consistent with equation (2.15). 5In this manner we observe that where H ∅ (x) = 1 and the constants C l are proportional to Bernoulli numbers B l , For example, the first few orders of this quantity are given by We have checked formula (4.17) up to O(g 24 ).
Moving on to the α = 0 case, it is easy to show that the general recursion formula (3.8), restricted to N → ∞ as in (3.38), simplifies to The boundary term similarly reduces to and can be expressed in terms of HPLs.Finally, the recursion terminates when both i and j are zero, at which point (4.15) evaluates to Due to the fact that depends on g 2 , we need to evaluate S i,j (α|x) for all i + j ≤ L − 1 to compute Ω (L) .We have explicitly carried out these sums for all i + j ≤ 11.These results are easy to verify numerically for various values of α by truncating the sum over n in equation (4.14).

4.2
Resumming Ω (L) through L = 10 loops Having evaluated the expansion of F −iα (x) in terms of Z-sums, we now carry out the outermost sum in equation (4.10).Denoting the perturbative expansion of this function as we may exploit the x ↔ y and z ↔ 1/z symmetry of Ω (L) by adopting the notation as well as by introducing the building blocks where a denotes the integer part of a.The sum representation of Ω (L) in (4.10) may then be written as where and finally as Note that the variables r and R used here differ from similarly-named variables used in [22] by a minus sign.To evaluate the double pentaladder integral in terms of generalized polylogarithms at a given loop order, we thus proceed as follows: Expressions through six loops are provided as ancillary computer files with the arXiv submission of this article.Readers interested in the results for seven, eight, nine, or ten loops may contact the authors.Finally, we note that similar sum representations for three additional integrals were derived in [22].Our evaluation method works equally well in these cases.However, these integrals are related to the derivatives of Ω (L) , so explicit polylogarithmic representations for them can be derived more efficiently by starting from the already-resummed representation of Ω (L) .

Comparison to Existing Results in the Literature
Before closing this section, let us compare our results for Ω (L) with the existing results in the literature.As mentioned above, the authors of [22] were able to resum Ω (L) in the u → 1 limit and the w → 0 limit.We have compared to these results by taking the appropriate limits of our expressions through four loops, and find complete agreement.
The double pentaladder integral Ω (L) was also computed in general kinematics via direct integration through four loops in [51].The comparison of our results with the expressions presented there is a bit subtle, as these representations are manifestly real in different regions.The sum representations of Ω (L) in (4.29) converges when x, y, |r|, |R| ≤ 1, and is thus manifestly real when z ∼ 1 and x, y 1.This corresponds to the neighborhood of the point u = v = w = 1. 6Conversely, the expressions given in [51] are manifestly convergent in the so-called positive region, which corresponds to positive values of the variables f 1 , f 2 , and f 3 used in that paper.In particular, the point u = v = w = 1 corresponds in those variables to the limit f 1 , f 3 → −1, f 2 → ∞, well outside of the positive region.
To compare these representations, we analytically continue the expressions for Ω (L) in terms of f 1 , f 2 , and f 3 out of the positive region to the neighborhood of the point u = v = w = 1.To do so, we note that when f 1 , f 2 , and f 3 are all small and positive, u v, and w are also all positive.Thus, the required analytic continuation path, which changes the signs of f 1 and f 3 , entirely exists within the Euclidean region.This implies that the signs of the cross-ratios u, v, and w should not change, which turns out to be true only if f 1 and f 3 are analytically continued in the opposite direction.After sending f 1 → f 1 e ±iπ and f 3 → f 3 e ∓iπ , we indeed find that the resulting expression is manifestly real near u = v = w = 1.This then allows us to confirm that the two representations match in this region.
Self-energy diagram with masses M 1 and M 2 , and propagator powers ν 1 and ν 2 .

Application II: Self-energy Diagram
Another natural use of the resummation algorithm presented in section 3.2 is the expansion of one-loop integrals in dimensional regularization.In this section we illustrate how this can be done for families of massive self-energy diagrams with different propagator powers.
We depict a generic massive self-energy diagram with propagator powers ν 1 and ν 2 in figure 2. This diagram was evaluated in [14], where it was expressed in terms of Appell's F 4 function.It was also studied there in various kinematic limits, including the limit of zero external momentum, Q 2 → 0. In that limit, the diagram can be written in terms of Gauss hypergeometric functions.For M 1 > M 2 , it is given by where x = and D = 4 − 2 .For M 1 < M 2 the expression is the same, but with the two masses exchanged.
The function 2 ) depends on two integers, the powers of the self-energy diagram propagators ν 1 and ν 2 .Thus, to make use of the methods of section 3 we must constrain ν 1 and ν 2 to depend on only a single symbolic integer.For concreteness, we impose the relation although we could have chosen any integer on the right hand side.After imposing the restriction (5.2), we expand the 2 ) around small values of using (2.16).This gives us Using the methods described in sections 2 and 3, f (ν) can be converted into Z-sums for generic integer values of ν.We have explicitly carried out this calculation through O( 6) (which took a reasonable time on a laptop).The result through O( 3) are given by where ) ) ) (5.9) (5.12) Unlike the sums contributing to the double pentaladders, f (ν|x) is not of uniform transcendental weight, and includes contributions ranging from weight 0 to weight n at O( n ).

Conclusions
In this paper, we have presented a new algorithm for converting sums of the form (1.2) into nested sums for symbolic values of α, N , x, y i , and z i .For generic values of N , these sums evaluate to cyclotomic harmonic sums, while in the N → ∞ limit they reduce to Z-sums.This algorithm allows us, in particular, to resum the expansion coefficients of Gauss hypergeometric functions in cases where this expansion is around integer indices that depend on a symbolic integer parameter α.While these hypergeometric expansion coefficients cannot be expressed solely in terms of generalized polylogarithms for generic values of α, the additional Z-sums that appear reduce to rational functions of their arguments for any specific value of α.
We have illustrated how this new resummation technology can be applied to the calculation of Feynman integrals in two examples.In the first, we considered the kinematic expansion of the double pentaladder integrals derived in [22], and explicitly evaluated this sum in terms of generalized polylogarithms through ten loops.This represents a significant advance over the previous state of the art, which leveraged the integral representation of these diagrams and became computationally infeasible beyond four loops [51].In the second example, we showed how this technology can be used to simultaneously expand massive one-loop self-energy diagrams with different propagator powers in dimensional regularization.
Arguably the biggest shortcoming of our algorithm when applied to the expansion of Gauss hypergeometric functions is the requirement that the symbolic integer appear with the same coefficient in the first and third indices, and be absent from the second index.This ensures that no binomial coefficients occur in the expansion of the hypergeometric function.It is our hope that the general telescoping strategy outlined in section 3.1 can also be leveraged to convert sums that depend on symbolic binomial coefficients into Z-sums and polylogarithms, but we leave this to future work.In the meantime, we highlight that the telescopic recursion (3.8) is significantly more general than its instantiation in the Gauss hypergeometric case; we anticipate it will find further applications in perturbative quantum field theory computations.
One direction in which new nested resummation technology would prove fruitful is in the Pentagon Operator Product Expansion (POPE) representation of amplitudes in planar N = 4 supersymmetric Yang-Mills [52][53][54][55][56][57][58].The POPE represents amplitudes-or rather, their dual null polygonal Wilson loops-at finite coupling in terms of an expansion in flux tube states propagating across the Wilson loop.This representation can be expanded around small coupling, resulting in infinite sum representations for perturbative amplitudes.Building on earlier work [59,60], the resummation of these expressions was initiated in [61], and further continued in [62][63][64][65].Due to the appearance of binomial coefficients, however, this resummation cannot be carried out systematically at higher orders.We are optimistic that a telescoping strategy will be a viable way forward in these cases.
A great deal is known about the analytic properties of Feynman integrals, and it would be interesting to understand the implications of these properties for sum representations of these integrals.For instance, Feynman integrals are expected to obey the Steinmann relations [66][67][68][69], and the double pentaladder integrals are even known to obey an extended set of Steinmann relations to all orders [22,70].However, it is not clear how this property is encoded in the representation of the double pentaladder integrals in (4.10).It would also be interesting to find a sum representation of loop-level amplitudes in planar N = 4 supersymmetric Yang-Mills that make the observed positivity properties of these amplitudes manifest [71,72], or that make connections with the cluster-algebraic properties of these amplitudes [73][74][75][76][77][78][79][80][81][82].
Plugging these values into (3.10),we can express the terminal sum as This sum over k cannot be carried out in terms of Z-sums.Rather, we are required to make use of cyclotomic harmonic sums, defined by where x = x 1 , . . ., x d is a multi-index of depth d, and S(N ) = 1.As with the Z-sums, cyclotomic harmonic sums satisfy a large number of identities, such as stuffle and synchronization identities; they can also be given an iterated integral representation.We refer the reader to [19] for more details.After carrying out the conversion Z 1 (α+k−1|z 1 ) = Z 1 (α−1|z 1 ) + z α−1 S {1,α−1,1} (z 1 ; k), the sum over k in (B.3) can be evaluated to give Putting this all together, we find It is easy to check that this expression numerically reproduces (3.16) (with y 1 set to 1) for large values of N .

Figure 1 .
Figure 1.The double pentaladder integral Ω (L) , built out of an (L−2)-loop box ladder capped on each end by a pentagon loop.The dashed lines represent specific numerator factors.