PolyLogTools - Polylogs for the masses

We review recent developments in the study of multiple polylogarithms, including the Hopf algebra of the multiple polylogarithms and the symbol map, as well as the construction of single valued multiple polylogarithms and discuss an algorithm for finding fibration bases. We document how these algorithms are implemented in the Mathematica package PolyLogTools and show how it can be used to study the coproduct structure of polylogarithmic expressions and how to compute iterated parametric integrals over polylogarithmic expressions that show up in Feynman integal computations at low loop orders.

An important problem in the practical use of these various classes of functions has been the wealth of functional identities that they satisfy. In many cases these functional identities needed to be derived manually for each case under consideration, hampering an efficient use of the more general classes of polylogarithms. However, in recent years, Goncharov's seminal paper [4] on the Hopf algebra structure of multiple polylogarithms has spurred a flurry of research both on the side of mathematics as well as in physics to better understand the algebraic structures of multiple polylogarithms. This has lead in particular to the realization that the Hopf algebra of multiple polylogarithms and its associated coproduct or coaction allow the derivation of functional identities between MPLs in a purely combinatorial fashion, see, e.g., refs. [7,28]. Based on this, an algorithmic way to compute Feynman integrals was devised [7,30] that makes use of the ability to derive any needed functional equation using the coproduct. Similar algorithms were further developed and automated in ref. [43][44][45][46]. In particular, the package HyperInt has recently been used to evaluate complicated Feynman parameter integrals [47][48][49][50].
In parallel the algebraic properties of multiple polylogarithms have been exploited in more formal developments such as the study of motivic amplitudes [18,19] and the amplitudes / cluster bootstrap in planar N = 4 SYM [21-24, 26, 51-53]. These approaches exploit a particular component of the Hopf algebra of the multiple polylogarithms, the socalled symbol map, to study the algebraic structure of scattering amplitudes. Building on the understanding of the branch-cut structure of the multiple polylogarithms, it was also possible to employ the coaction to build a special class of multiple polylogarithms, called single-valued multiple polylogarithms (SVMPLs) that are free of branch cuts [8,54]. These SVMPLs have played an important role in the study of scattering amplitudes in the so-called multi-Regge Limit [55][56][57][58] and in the calculation and study of certain periods in φ 4 -theory [59]. Inspired by the Hopf-algebra of the multiple polylogarithms and Brown's generalization to a conjectured cosmic Galois group [60], it was shown in refs. [61][62][63] that the coaction of one-loop Feynman integrals can be cast in a remarkably simple and compact form.
All these developments have profited enormously from the mathematical advances in the understanding of the algebraic properties of the multiple polylogarithms. The goal of this paper is to document the Mathematica package PolyLogTools that provides an implementation of many of these algebraic structures. The focus of this implementation is on providing a basis for future exploration and experimentation, as well the ability to compute Feynman integrals that appear in two-loop and three-loop amplitudes [30]. This package in particular is not optimized for the performance that would be required to tackle Feynman integrals at arbitrary high loop orders, but it provides well-tested implementations of all required algorithms and provides the flexibility to adapt it to virtually any application mentioned above. 1 PolyLogTools provides many routines that allow the user to work with multiple polylogarithms in Mathematica; from integrating expressions in terms of polylogarithms through numerical evaluation and calculation of symbols and coproducts to the automated derivation of functional identities. The goal of this paper is to review the relevant mathematics of polylogarithms and to document the routines in PolyLogTools that implement the respective mathematical algorithms.
This paper is organized as follows: In section 2 we describe how to obtain and install the package and its prerequisites. Then, in section 3 we review the basic construction of the multiple polylogarithms and describe how they are realized in PolyLogTools. In section 4 we review the first important algebraic structure that multiple polylogarithms are equipped with, the so-called shuffle algebra, and show how to use it in the package. Next, in section 5, we discuss the degeneration of multiple polylogarithms to multiple zeta values for special values of their arguments and show how they are implemented in PolyLogTools. In section 6 we review the Hopf algebra and coproduct of multiple polylogarithms and its implementation. Then, in section 7 we discuss the relation between the coproduct and the symbol map. In section 8 we discuss how to perform basic calculations in PolyLogTools. In section 9 we review the theory of SVMPLs and discuss the routines dedicated to handling these functions. Afterwards, in section 10 we review a few applications of PolyLogTools in the literature. Finally, in section 11 we discuss the calculation of a Feynman integral from start to finish to illustrate the practical use of PolyLogTools for actual calculations. consistent with the case n = 0 above. In the case where the a i are constants, MPLs are often referred to as hyperlogarithms. MPLs define a very general class of functions that generalize the well-known logarithm and (Nielsen) polylogarithm functions, e.g., for a = 0, G( a n ; z) = 1 n! log n 1 − z a , G( 0 n−1 , a; z) = −Li n z a , G( 0 n−k , a k ; z) = (−1) k S n−k,k z a . (3. 3) The function G(a 1 , . . . , a n ; z) is represented within PolyLogTools by the symbol G[a1,...,an,z]. Here the arguments a1,. . . ,an and z can be any valid Mathematica expression (in practice, we will only deal with cases where the arguments are rational or algebraic expressions). Since G(z) = 1, the function G automatically evaluates to unity if it only has a single argument.
In the case where all the a i are 0 or ±1, MPLs are often referred to as harmonic polylogarithms (HPLs) in the physics literature [42]. HPLs are equal to MPLs, up to a sign, H(a 1 , . . . , a n ; z) = (−1) p G(a 1 , . . . , a n ; z) , a i ∈ {0, ±1} , (3.4) where p denotes the number of elements in a equal to (+1). HPLs are represented within PolyLogTools by the symbols H[a1,...,an,z]. The functions HToG [expr] and GToH[expr] allow one to switch from the H to the G notation inside the Mathematica expression expr. Note that since the Hpl package is loaded automatically with PolyLogTools, also the notation for HPLs from the Hpl package can be used inside PolyLogTools. The functions HToHPL[expr], GToHPL[expr], HPLToH [expr] and HPLToG[expr] allow the user to switch between the notations used by Hpl and PolyLogTools. In particular, the function HPLToG[expr] can be used to convert a HPL in compressed notation to a G in standard notation, for example: There is a second way to define MPLs, using nested sums rather than iterated integrals [2]: Li m 1 ,...,m k (z 1 , . . . , z k ) = 0<n 1 <n 2 <···<n k z n 1 1 z n 2 2 · · · z n k k n m 1 1 n m 2 2 · · · n m k k (3.5) . . .  The multiple polylogarithm G(a 1 , . . . , a n ; z). The indices a i can be numbers, symbolic constants or functions of other variables. H[a1,a2,...,an,x] The harmonic polylogarithm H(a 1 , . . . , a n ; z). The indices a i have to be integers from the set {−1, 0, 1}.
The multiple polylogarithm Li m1,...,mn (z 1 , . . . , z n ) in sum notation. The m i are integers. Conversion functions: Replaces every G in the expression expr with the corresponding H, if the indices of the multiple polylogarithm are integers from the set {−1, 0, 1}.

HToHPL[expr]
Converts every H in expr to the HPL notation.

GToHPL[expr]
Replaces every G in expr with the corresponding HPL, if the indices of the multiple polylogarithm are integers from the set {−1, 0, 1}.

HPLToH[expr]
Replaces every HPL in expr with the corresponding H.
Automatically converts from the a−notation defined by Hpl to the m−notation.

HToG[expr]
Replaces every H in expr with the corresponding G.

HPLToG[expr]
Replaces every HPL in expr with the corresponding G.
Automatically converts from the a−notation defined by Hpl to the m−notation.

GToLi[expr]
Converts every G in expr to the Li notation.

LiToG[expr]
Converts every Li and classical polylogarithm in standard Mathematica notation in expr to the G notation.
The functions LiToG[expr] and GToLi[expr] allow the user to switch between the two representations of MPLs, either as iterated integrals (G) or as nested sums (Li).

The shuffle and stuffle algebras of MPLs
One of the most basic properties of MPLs is that they can be equipped with two algebra structures, one related to the representation of MPLs as iterated integrals -the shuffle algebra -and another one related to the representation as nested series -the stuffle algebra. Most of the functionalities of PolyLogTools rely on the shuffle algebra properties of MPLs. Since we will encounter another shuffle algebra in Section 7, we start this section by reviewing the mathematics of shuffle algebras in general before discussing the implementation of the shuffle algebra of MPLs in PolyLogTools.

Shuffle algebras
Consider a (finite) set A, whose elements we will denote as letters, and we call A the alphabet. For concreteness, we will choose here A = {a, b, c, . . .}, though A can be any finite set. We define a vector space V as the vector space formed by all Q-linear combinations of words formed from the letters in A, including the empty word, which we denote simply as 1. The length of a word is defined as the number of letters that the word is made of. Words of length 1 are simply the letters. The concatenation of two words w 1 and w 2 is defined in the obvious way and denoted by w 1 w 2 . V has the structure of a graded vector space, i.e., it admits a direct sum decomposition where V 0 = Q and V n and V >0 denote the subspaces of V spanned by all words of length n and all words of non-zero length respectively. V can be given the structure of a commutative algebra equipped with the shuffle product. The shuffle product assigns to a pair of words (w 1 , w 2 ) the sum of all their shuffles, i.e., the sum of all possible ways of permuting the letters of their union without changing the order of the letters within each word. The shuffle product can also be defined recursively: if α, β ∈ A are letters and w 1 , w 2 ∈ V are words, then the shuffle product of the words αw 1 and βw 2 is defined recursively by 2) and the recursion starts with 1 ¡ w = w ¡ 1 = w. For example, the shuffle product of the words ab and cd is ab ¡ cd = abcd + acbd + cabd + acdb + cadb + cdab . It is easy to check that the shuffle product is associative and commutative. Moreover, it preserves the length, i.e., the shuffle product of two words of lengths n 1 and n 2 is a linear combination of words of length n 1 + n 2 . In this way V becomes a graded algebra, It is often convenient to work with a set of generators for V , i.e., a minimal set of words such that every element of V can be written as a linear combination of products (a polynomial) in these generators. In order to define such a set of generators, we first need to define an ordering < on the set of letters A (for concreteness, we can choose here the lexicographic ordering among the letters, but any other choice of ordering would do).
A theorem by Radford states that the shuffle algebra V is isomorphic to the polynomial algebra (over Q) formed by a subset of words, called Lyndon words [64]. A Lyndon word is a non-empty word that is less (for the chosen ordering <) than any of its proper right factors, i.e., a word w is a Lyndon word if for all non-empty words w 1 and w 2 with w = w 1 w 2 we have w < w 2 . A consequence of Radford's theorem is that for a given ordering every word can be uniquely written as a polynomial in Lyndon words. For example, the word baa is not a Lyndon word for the lexicographic ordering (because baa > aa), so it can be written as a polynomial in the Lyndon words, (4.5)

The shuffle algebra of MPLs
Let us now consider the Q-vector space of all MPLs of the form G( a; z) ending in the same variable z. It can be shown that this vector space forms a graded shuffle algebra, where G( a; z) is identified with the word a = (a 1 , . . . , a n ). The weight of an MPL corresponds to the length of the word. Explicitly, the shuffle product of two MPLs ending in the same variable can be written as where the sum runs over all shuffles of the two weight vectors a and b. It is possible to use PolyLogTools to linearise all shuffle products of MPLs ending in the same argument z by using the function ShuffleG as shown in the following example (cf. eq. (4.3)), Since every word can be represented as a polynomial in Lyndon words for a given ordering, we can also decompose every MPL into a linear combination of products of MPLs whose weight vectors are Lyndon words. This is achieved as follows (cf. eq. (4.5)), The function DecomposeToLyndonWords does not only work on single MPLs, but its argument can be any Mathematica expression. PolyLogTools has tables of Lyndon words hard-coded up to words of length six, and so only MPLs up to weight six are reduced to Lyndon words. The tables of Lyndon words implemented in PolyLogTools have been generated with Sage [65]. The optional argument Alphabet takes as value a list of symbols and defines the alphabet as well as the ordering among the letters (simply the order in which the letters appear inside the list). Its default value is {0,1,-1}. Note that only those MPLs whose weight vectors contain only letters from the Alphabet are decomposed into Lyndon words. Finally, let us make a comment about singularities. The MPL G(a 1 , . . . , a n ; z) diverges whenever z = a 1 , because in that case the integral in eq. (3.1) has an end-point singularity. It may happen that the decomposition into Lyndon words leads to divergent MPLs, even though the input was finite (e.g., take z=a in the above example box). In such a case the decomposition is not performed and a warning is shown. It is often convenient to decompose MPLs not into a Lyndon word basis, but into a set of functions corresponding to words where a given letter does not appear in the last position (except for words that are only made of that letter). For example, in applications it is useful to represent a function in terms of MPLs where the last letter is non-zero, except for powers of G(0; z). Such a representation can always be achieved algorithmically by unshuffling powers of logarithms. For example, we have This operation is implemented via the function ExtractZeroes, which can be used as shown,

The stuffle algebra of MPLs
There is another algebra structure on MPLs coming from their series representation, called stuffle algebra [66]. Since in many applications to Feynman integrals the shuffle algebra structure seems more relevant, we will be brief and only present an example. For a more general discussion, we refer to ref. [66]. Consider a product of two polylogarithms of depth one. We can rearrange the sums in the following way:   Expands all products of G-functions that end in the same argument into a sum of shuffles.

StuffleLi[expr]
Expands all products of Li-functions into a sum of stuffles.

DecomposeToLyndonWords[expr]
Decomposes all G-functions into a sum of products of Lyndon words. The set of letters and their ordering can be passed as a list via the optional argument Alphabet. The default value is {0,1,-1}. MPLs which depend on a letter that is not in the Alphabet are not decomposed.
The same applies to MPLs that would be decomposed into divergent quantities.

ExtractZeroes[expr]
Uses the shuffle algebra to remove all trailing zeroes from all G-functions in expr.

Special values of MPLs
When evaluated at certain special values of the arguments, MPLs often reduce to known transcendental constants. In particular, if HPLs are evaluated at z = ±1, then they reduce to (coloured) multiple zeta values (MZVs), The weight and the depth of (coloured) multiple zeta values are defined in the same way as for MPLs. The series in eq. (5.1) diverges whenever m k = 1. For depth k = 1, MZVs reduce to Riemann's zeta function at positive integer arguments, MZVs and their coloured generalisations (which correspond to some of the m i being negative) are implemented in the Hpl package, where ζ m k ,...,m 1 is represented by the symbol MZV[{mk,...,m1}]. In addition, the Hpl package also knows how to reduce HPLs evaluated at z = ±i to a smaller set of transcendental constants. Since PolyLogTools uses Hpl the reduction rules of HPLs evaluated at z ∈ {±1, ±i} are readily available in PolyLogTools also for the G and H functions. We refer to the documentation of the Hpl package for further details [67,68]. In addition to the values of HPLs evaluated at z = ±i, PolyLogTools can also reduce G( a; z) with a i ∈ {0, ±i} and z ∈ {±1, ±i}, as well as HPLs for some small weight evaluated at z ∈ {±1/2, ±1/3}, to known transcendental constants.
By default, PolyLogTools will automatically express any MPL in terms of transcendental constants whenever possible. It is possible to disable the automatic reduction to transcendental constants by setting The default value of PLT$AutoConvertToKnownConstants is True.
We conclude this section with a comment on the regularisation of MPLs. We have seen that G(a 1 , . . . , a n ; z) diverges whenever z = a 1 . Hence, when evaluating HPLs at z = ±1 the result may be divergent. It is easy to see that all divergencies are logarithmic, and it is useful to introduce a formal quantity lim z→1 G(1; z) = −ζ 1 which acts as a regulator and allows one to keep track of divergences. Indeed, all divergencies show up as powers of ζ 1 , and the dependence on the regulator ζ 1 must cancel for finite quantities. The regulator ζ 1 is represented by HPL [{1},1] in Hpl and PolyLogTools.
There is, however, a subtle point about the regularisation discussed in the previous paragraph. As usual, there is an ambiguity in how to introduce the regulator, related to shifting divergent quantities by finite terms. In practice, one would like to choose a regulator that preserves as many of the algebraic structures as possible. MPLs are equipped with both a shuffle and a stuffle algebra structure. However, it is not possible to preserve both structures at the same time, leading to two different 'schemes' to regulate MPLs and MZVs [69], referred to as shuffle regularisation and stuffle regularisation respectively. For example, using the shuffle product on MPLs, we obtain so that the shuffle-regulated value of G(1, 1; z) at z = 1 is given by 1 2 ζ 2 1 . Instead, using the fact that G(1; z) = −Li 1 (z) as well as the stuffle product, we find so that the stuffle-regulated value of G(1, 1; z) at z = 1 is given by 1 2 ζ 2 1 − 1 2 ζ 2 . We see that the shuffle-and stuffle-regulated values are different. Note that any finite quantity must be independent of the regularisation scheme chosen to compute it. However, care is needed that the same scheme is applied consistently throughout a computation! By default, PolyLogTools uses shuffle regularisation to regulate the G and H functions, because in applications it is often more convenient to preserve the shuffle algebra structure (which is more closely related to the definition of MPLs as iterated integrals). It is possible to switch to stuffle regularisation by changing the value of the variable PLT$ShuffleRegularisation from True to False. We point out that, unlike PolyLog-Tools, the Hpl package uses stuffle regularisation to define the values of MZVs and HPLs at z = ±1. Care is thus needed if the output and functions from this package are used in conjunction with PolyLogTools.

The Hopf algebra of MPLs
An important property of MPLs is that they can be equipped with a coproduct turning them into a Hopf algebra. The Hopf algebra of MPLs has led to the development of novel techniques to deal with Feynman integrals that evaluate to MPLs, and it has seen a multitude of applications over the last few years. In particular, it allows one to control more rigorously identities among MPLs, see ref. [28,70] for a pedagogical introduction of how to derive identities among MPLs using the coproduct. Another important application of the coproduct is the amplitude bootstrap in planar N = 4 Super Yang-Mills [21-24, 26, 51]. Here we will not discuss these applications in detail, but we focus on the implementation of the coproduct on MPLs in PolyLogTools. Since we will discuss different Hopf algebras in subsequent sections, we give a short review of Hopf algebras in general before focusing on the case of MPLs.

A short review of Hopf algebras
A coalgebra is a Q-vector space H together with a linear map ∆ : H → H ⊗ H called the coproduct. The coproduct is required to be co-associative, that means A coalgebra must admit another map, called the counit, which does not play any role in the context of the Hopf algebras discussed in this paper. A Hopf algebra is a vector space that is both an algebra and coalgebra, such that the product and the coproduct are compatible, In the previous equation the multiplication of tensors on the right-hand side is defined component-wise, In addition, there is a map S : H → H, called the antipode, which will be defined below. The coalgebras encountered in this paper are graded and connected, i.e., they admit a direct sum decomposition 4) and the coproduct respects the grading, We can iterate the coproduct to obtain tensors with more and more factors. This iteration can be done in different ways, e.g., by iterating on the first or on the second factor of the coproduct. The co-associativity of the coproduct ensures that the different ways of iterating the decomposition into simpler objects give the same result. Since the coproduct respects the grading, it makes sense to define ∆ i 1 ,...,i k as the the part of the iterated coproduct that takes values in H i 1 ⊗ · · · ⊗ H i k . The maps ∆ i 1 ,...,i k : H → H i 1 ⊗ · · · ⊗ H i k satisfy the obvious recursion Let us also discuss the antipode S. In the case of graded and connected Hopf algebras the antipode is uniquely determined by the coproduct to be is the multiplication in H. The antipode is linear, acts trivially on elements of weight 0, S(1) = 1, and preserves the multiplication and the coproduct Let us conclude this section by illustrating these definitions on a concrete example of a Hopf algebra. This example will come back in Section 7 in the context of MPLs. We start from the shuffle algebra V of words from Section 4.1. Every shuffle algebra can be turned into a graded and connected Hopf algebra whose coproduct is given by the deconcatenation of words, where the sum runs over all deconcatenations of the word w, and including the trivial ones where either w 1 or w 2 is the empty word. One can check that this definition satisfies all the properties of a coproduct. The antipode is given by the reversal of words, where |w| denotes the length of w, and the tilde denotes the reversal of words (e.g., if w = abc, thenw = cba).

The Hopf algebra of MPLs
In ref. [4] it was argued that MPLs form a graded and connected Hopf algebra, where the grading comes again from the weight of the functions. The explicit formula for the coproduct is rather involved, and most conveniently expressed not in terms of the MPLs defined in eq. (3.1), but in terms of the iterated integrals I(a 0 ; a 1 , . . . , a n ; a n+1 ) = a n+1 a 0 dt t − a n I(a 0 ; a 1 , . . . , a n−1 ; t) .
It is clear that the functions defined in eqs. (3.1) and (6.11) define the same space of functions. In terms of the functions defined in eq. (6.11), the coproduct on MPLs can be written in the following compact form [4], ∆(I(a 0 ; a 1 , . . . , a n ; a n+1 )) = Here the sum runs over all ordered subsets of a given length k of (a 1 , . . . , a n ). Since the Hopf algebra of MPLs is graded and connected, the antipode is uniquely determined by the coproduct and eq. (6.7), so we do not show it here explicitly. Let us make some comments about the formula for the coproduct in eq. (6.12). First, the individual terms in eq. (6.12) admit a simple combinatorial interpretation in terms of polygons inscribed into a semi-circle [3,4]. It is based on this combinatorial framework that the coproduct on MPLs is implemented in PolyLogTools, and not via the explicit formula in eq. (6.12). We do not review the combinatorial picture based on semi-circles here, but we refer to ref. [70] where worked-out examples can be found. Second, we note that the formula for the coproduct in eq. (6.12) is only valid in the generic case where all the arguments a i are distinct. In the non-generic case individual terms in eq. (6.12) may diverge, and we need to replace the MPLs on the right-hand side of eq. (6.12) with suitably regularised versions. We follow closely refs. [4,28,70], identifying all MPLs on the right-hand side of eq. (6.12) with their shuffle-regularised versions (cf. Section 5), and setting to zero the regulator ζ 1 .
The implementation of the coproduct on MPLs is one of the main features of the PolyLogTools package, because the coproduct is the basis for many applications. The coproduct is called via the function Delta, which takes as argument any valid expression in terms of MPLs of weight up to twelve and returns its coproduct. For example, we have Tensors are represented in PolyLogTools by lists with head CT, i.e., the tensor A 1 ⊗A 2 ⊗ A 3 ⊗ . . . is represented by the symbol CT[A1,A2,A3,...]. 2 Note that Delta can be applied to any valid expression made of polylogarithms, independently of how they are represented (G, H, Li, Log, PolyLog). The only restriction is that no MPLs of weight higher than twelve are allowed in the current implementation. The coproduct also acts non trivially on MZVs (because they are just special values of MPLs), e.g., The symbol CT satisfies the basic properties of a tensor product. In particular, it Moreover, it is linear with respect to rational numbers and the imaginary unit i, e.g., It also automatically evaluates products of tensors (cf. eq. (6.3)), The different components ∆ i 1 ,...,i k of the iterated coproduct can be accessed as in the following example for ∆ 2,1,1 , The different components of the iterated coproduct are evaluated internally using the recursion in eq. (6.6).
At this point we have to make an important comment: MPLs are multivalued functions. It is only MPLs modulo their discontinuities that form a Hopf algebra. The discontinuities are related to taking residues at the simple poles in the integrand in eq. (3.1). Therefore the discontinuities of MPLs are always proportional to iπ. Hence, in order to obtain a Hopf algebra, we need to work modulo iπ. In applications, however, it is important to keep track of powers of π. It is possible to incorporate iπ into the construction by introducing the special rule [7] (see also ref. [28]), 3 ∆(iπ) = iπ ⊗ 1 .

The Lie coalgebra of indecomposables
In applications it can be useful to focus on the most complicated part of an expression. One way to do so is to look at an expression defined modulo products. In particular, one may want to decide if two expressions are equal up to 'simpler' product terms. Since MPLs form a shuffle algebra, it can be hard to decide if an expression is zero up to product terms. For example, one may be tempted to believe that the expression T = G(a, b; z) + G(b, a; z) does not involve any product terms. In this case it is easy to see that the sum in T is precisely a shuffle product, so that T actually vanishes modulo products. In this section we show how we can construct a map whose kernel is precisely generated by all products among MPLs.
If H denotes a graded and connected Hopf algebra, then we define its space of indecomposables Q(H) as the Hopf algebra H modulo all non-trivial products (i.e., products among objects of weight at least one), (6.14) Q(H) is obviously a vector space (because it is the quotient of two vector spaces). It is, however, not an algebra (and thus not a Hopf algebra), because it was defined precisely by removing all products. We now construct a projector P : H → Q(H) which allows one to remove all product terms. We start by recursively defining a linear map R which acts on x ∈ H n , n > 0, by [71] R where m is the multiplication in H and ∆ is the reduced coproduct, defined below eq. (6.7). One can show that the kernel of R is precisely generated by all non-trivial products in The projector P = P 2 is then obtained by correctly normalising R, P (x) = 1 n R(x) for x ∈ H n . For example, we have We see from the previous example that the image of an MPL may contain product terms, despite the fact that the goal was to remove product terms. There is no contradiction: the elements of Q(H) are equivalence classes of MPLs defined modulo product terms. The projector P assigns to a function a canonical representative of its equivalence class modulo products. The product terms ensure that relations among equivalence classes are satisfied. For example, we have In agreement with the fact that G(1, 0; z)+G(0, 1; z) is a shuffle product, and thus vanishes modulo products. The projector P is implemented via the function ProductProjector. For example, we can obtain eq. (6.18) from PolyLogTools, The coproduct on H induces a new structure on the space of indecomposables Q(H), called a Lie coalgebra, equipped with a cobracket δ : where τ (a⊗b) = b⊗a is the map that reverses tensors. The cobracket has recently appeared in physics in the context of scattering amplitudes in planar N = 4 SYM [19,72]. We take the opportunity to make some technical comments: First, the cobracket in eq. (6.19) is really only defined on MPLs modulo their discontinuities, i.e., we have to put to zero of powers of π in both entries of the wedge product. Second, in applications the cobracket in eq. (6.19) is directly defined on the MPLs themselves, i.e., on elements of H, while it is in principle only defined on the equivalence classes living in the space of indecomposables Q(H). This is possible because the map δ in eq. (6.19) has the property δP = δ, i.e., the image of an element of H under δ agrees with the image of its projection on Q(H). For example, the cobracket of the dilogarithm is The cobracket on MPLs is implemented via the function Cobracket, which can be used as shown in the following example (cf. eq. (6.20)), The wedge product is implemented via the symbol CTW. 4 This symbol inherits its properties from the symbol CT defining the tensor product. In particular, it is linear with respect to rational numbers and i, and puts to zero all occuurences of π,  Represents the tensor a 1 ⊗ · · · ⊗ a n .

Antipode[expr]
Computes the antipode of the symbol expr.

ProductProjector[expr]
Applies the projector P to expr.

Symbols in PolyLogTools
While the coproduct on MPLs is very useful in applications, it often contains too much information. It is often useful to focus on a piece of the coproduct which is easier to work with, albeit at the price of losing some information. Such a quantity is the symbol [4,6,18,28,73,74], which can be defined as the maximal iteration of the coproduct (modulo iπ), The symbol contains the same information as the maximal iteration of the coproduct. It is very easy to work with, because its entries are MPLs of weight one, i.e., ordinary logarithms. For this reason it is customary to drop the log-signs when talking about the symbol, e.g., Besides this notational difference, there is no difference between the symbol S and the maximal iteration of the coproduct ∆ 1,...,1 .
PolyLogTools contains a function that allows the user to turn a maximal iteration of a coproduct into a symbol (which effectively only amounts to removing the log-signs), e.g., There is also a function which readily combines the maximal iteration of the coproduct with the function ToSymbol: Hence, whenever an entry in a list with head CiTi is ±1, it automatically evaluates to zero, We note here that the alphabet returned by GetSymbolAlphabet does not necessarily consist of irreducible polynomial letters, but it simply collects all the entries in the symbol. In order to obtain an alphabet of irreducible letters, the GetSymbolAlphabet can be composed with the SymbolExpand function, There are various equivalent definitions of symbols in the literature. An important feature of the implementation of the symbol map in PolyLogTools is that it does not put to zero symbol tensors that contain a constant letter, e.g., Keeping constant letters in the symbol sometimes provides valuable information on the underlying function, cf., e.g., refs. [28,74].
PolyLogTools contains another implementation of the symbol map S based on dissections of decorated rooted polygons attached to MPLs [74]. It can be called through the function ComputeSymbol. The result is fully equivalent to applying the function ToSymbol to the maximal iteration of the coproduct, e.g., While the two implementations are fully equivalent, we mention here that (based on experience) the implementation of ComputeSymbol is usually faster for weights less than four, whereas SymbolMap works more efficiently for higher weights. 6 The symbol map assigns to an MPL expression a symbol tensor. The question then naturally arises if any symbol tensor constructed from a certain alphabet can be the symbol of a function. It turns out that this is not the case in general, however, but the symbol tensor S = is the symbol of a function, i.e., there is a function F such that S(F ) = S, if and only if it satisfies the integrability condition for every consecutive pair of indices (i p , i p+1 ), 1 ≤ p < n and ∧ denotes the usual wedge product of differential forms. For every value of p, the integrability condition translates into a system of linear constraints on the coefficients c I . PolyLogTools can generate these constraints for a given symbol S via the command IntegrablityCondition[ S, p ]. We emphasise that PolyLogTools does not attempt to solve the integrability constraints, because these linear constraints usually give rise to very large linear systems whose solution may require dedicated algorithms and/or specialised software.

The shuffle Hopf algebra of symbols
The symbol map is not only linear, but it also preserves the multiplication of MPLs and maps a product of MPLs to the shuffle product of their symbol tensors, S(x · y) = S(x) ¡ S(y) .

(7.7)
In Section 6.1 we have seen that every shuffle algebra is naturally a Hopf algebra, with the coproduct and antipode given by the deconcatenation and reversal of words, cf. eqs. (6.9) and (6.10). It follows that the shuffle algebra B of all symbol tensorsintegrable or not -forms a Hopf algebra with the deconcatenation coproduct and antipode. It is then interesting to ask what is the relation between the coproduct ∆ on MPLs and the deconcatenation coproduct ∆ dec on symbol tensors. The shuffle algebra B, however, is too large, because the image of the symbol map is the subspace of B consisting of integrable symbol tensors, which we denote by H 0 B. It is easy to check that H 0 B is a Hopf subalgebra of B. Moreover, the symbol map preserves the coproducts (and also the antipodes), i.e., S maps the coproduct of an MPL to the deconcatenation coproduct of its symbol, and similarly for the antipode, S dec (S(x)) = S(S(x)) . (7.8) The deconcatenation coproduct and antipode are implemented in the code via the functions DeltaDeconcatenation and AntipodeDeconcatenation, which act on symbol tensors as follows, Since H 0 B is a graded and connected Hopf algebra, we can apply the results of Section 6.3 and study its Lie coalgebra of indecomposables. First, following eq. (6.15) we define a linear map ρ whose kernel is precisely generated by all shuffles, ρ(a 1 ⊗ . . . ⊗ a n ) = n a 1 ⊗ . . . ⊗ a n − ¡(id ⊗ ρ)∆ dec (a 1 ⊗ . . . ⊗ a n ) . The previous definition can be cast in the equivalent form [75,76] ρ(a 1 ⊗ . . . ⊗ a n ) = ρ(a 1 ⊗ . . . ⊗ a n−1 ) ⊗ a n − ρ(a 2 ⊗ . . . ⊗ a n ) ⊗ a 1 . (7.10) From this map we define a projector Π 2 = Π by [74] Π(a 1 ⊗ . . . ⊗ a n ) = 1 n ρ(a 1 ⊗ . . . ⊗ a n ) , (7.11) and the cobracket on symbols is given by (cf. eq. (6.19)) The projector Π and the deconcatenation cobracket δ dec are implemented via the functions ProductProjector and CobracketDeconcatenation. They have the same properties as their analogues acting on MPLs described in Section 6.3, so we will not describe them here in more detail.

The cobracket conjecture, classical and Nielson polylogarithms
In ref. [18] a conjectural criterion on the symbol of a function was presented that allows one to determine if a function of weight four can be expressed in terms of classical polylogarithms only, and it was applied to the analytic result for the two-loop six-point remainder function in planar N = 4 Super Yang-Mills of ref. [77,78]. The criterion can be concisely stated in terms of the cobracket δ dec : if T is a polylogarithmic function of weight four such that its symbol satisfies δ dec,2,2 (S(T )) = 0 , (7.13) then T can be expressed in terms of classical polylogarithms only. The cobracket on symbol tensors of length n = 4 has a particularly simple form, 14) It is of course interesting to speculate how this criterion extends beyond weight four. Naive speculation would lead to the guess that if T is a polylogarithmic expression weight n such that δ dec,p,n−p (S(T )) = 0 for all 1 < p < n − 1, then T can be expressed (possibly up to terms in the kernel of S) in terms of classical polylogarithms only. This extension, however, seems to be too naive, because already at weight five it is not known how to express the Nielsen polylogarithm S 2,3 (x) in terms of classical polylogarithms only, even though δ dec,2,3 (S(S 2,3 (x))) = 0. Here S n,p (x) denotes the Nielsen polylogarithm [79], Instead, a folklore conjecture (which can be proven in some cases [80]) states that 7 If T is a polylogarithmic expression of weight n such that δ dec,p,n−p (S(T )) = 0 for all 1 < p < n − 1, then T can be expressed (possibly up to terms in the kernel of S) in terms of Nielsen polylogarithms only.
It would be interesting to have a sharper extension of the criterion of ref. [18], which allows one to determine if a function can be expressed in terms of classical polylogarithms only even at higher weight. In the following we present such a conjecture. As a starting point, we note that the cobrackets in eqs. (6.19) and (7.12) are related by In the remainder of this section we give evidence for the following conjecture: If T is a polylogarithmic expression of weight n such that δ p,n−p (T ) = 0 for all 1 < p < n − 1, then T can be expressed in terms of classical polylogarithms only.
In the remainder of this section we present some evidence for this conjecture. In particular, we show through weight six relations that relate different Nielsen polyogarithms up to classical polylogarithms. First we note that the conjecture is always true up to weight three, in agreement with the fact that up to weight three all MPLs can be expressed in terms of classical polylogarithms only [1,81,82]. More generally, it is known that S 1,n−1 (x) and S n−1,1 (x) can always be expressed in terms of classical polylogarithms, and indeed we have δ p,n−p (S 1,n−1 (x)) = δ p,n−p (S n−1,1 (x)) = 0 , for all 1 < p < n − 1 . (7.17) Next let us discuss the Nielsen polylogarithm S 2,2 (x). It is easy to check that δ 2,2 (S 2,2 (x)) = 0. We can indeed express S 2,2 (x) in terms of classical polylogarithms only, At weight five, we find for the first time Nielsen polylogarithms that have a non-trivial cobracket, Since classical polylogarithms of weight five lie in the kernel of δ 2,3 , we conclude from the previous equations that S 2,3 (x) and S 3,2 (x) cannot be expressed in terms of classical polylogarithms alone. Note that, since S(ζ 3 ) = 0, we have, δ dec,2,3 (S(S 2,3 (x))) = δ dec,2,3 (S(S 3,2 (x))) = 0 , (7.20) and so the cobracket computed from the symbol alone does not allow one to reach this conclusion. However, we learn something more from eq. (7.19): we see that δ 2,3 (S 2,3 (x)) = δ 2,3 (S 3,2 (x)), and so the conjecture implies that the two functions are equal up to classical polylogarithms. Indeed, we find the following relation, Finally, let us analyse weight six. We find δ 2,4 (S 2,4 (x)) = δ 2,4 (S 4,2 (x)) = δ 2,4 (S 3,3 (x)) = 0 , and ) .

(7.23)
In agreement with our conjecture, we find that we can express S 4,2 (x) and S 3,3 (x) up to classical polylogarithms in terms of S 2,4 (x) and S 2,4 (1 − x). The explicit expressions are and (7.25) Computes the symbol of expr as the maximal iteration of the coproduct.

ComputeSymbol[expr]
Computes the symbol of expr from dissections decorated polygons.

ToSymbol[expr]
If expr is a maximal iteration of a coproduct, then ToSymbol turns it into a symbol expression, by removing the log-signs. CiTi[a1,...,an] Represents the symbol tensor a 1 ⊗ . . . ⊗ a n .

SymbolExpand[sym]
Maximally factors all the polynomial symbol letter and uses the additivity of the symbol sym to expand them out.

SymbolFactor[sym]
Attempts to combine terms into products and ratios in the symbol. Works best if the coefficients of the symbol are integers.

GetSymbolAlphabet[sym]
Returns the symbol alphabet of the symbol sym.

IntegrabilityCondition[sym, p]
Returns the constraint from the integrability condition applied to the factors p and p + 1 in the symbol sym.

AntipodeDeconcatanation[sym]
Computes the deconcatenation antipode of the symbol sym.

ProductProjector[sym]
Applies the projector Π to the symbol sym and projects it to the space of indecomposables.

CobracketDeconcatenation[sym]
Computes the cobracket attached to the deconcatenation coproduct of the symbol sym.

Working with PolyLogTools
So far we have only discussed the implementation of the basic objects -MPLs, their coproduct and symbols -into PolyLogTools, and we have discussed their very basic usage (e.g., how to compute the coproduct or the symbol of an MPL). In this section we describe a set of functions that can be used at runtime to manipulate expressions involving MPLs in Mathematica. Note that since PolyLogTools automatically loads Hpl, all functions from this package are also available. For a description of these functions we refer to refs. [67,68].

Manipulating expressions
We start by describing a collection of functions that are useful to manipulate MPL expressions inside Mathematica. A valid MPL expression is a polynomial built out of the functions G, H (or HPL) and Li (and the built in Mathematica functions PolyLog and Log), as well as the transcendental constants MZV, Zeta and Pi (as well as certain special symbols like HPLs6 defined by Hpl [67,68]). All the functions described in this section can be applied to any valid MPL expression. When working with big expressions, it is often useful to focus on certain subexpressions, or to read out the elementary building blocks (e.g. MPLs) the expression is composed of. The function GetGs returns a list of all Gs in the expression. Analogously, GetGArguments returns the set of arguments of the G, H and HPL objects in the given expression, while GetGIndices returns the entries of the weight vectors of the G, H and HPL objects. For example, we have The previous functions allow the user to focus on a specific subset of an expression. There are other functions which allow one to manipulate an expression without changing its value. In Section 4 we have already encountered the functions ShuffleG, StuffleLi, DecomposeToLyndonWords and ExtractZeroes. They all fall into this category. Since they have already been discussed in Section 4, we will not discuss them here.
We can use eq. (8.1) to reduce all MPLs with a n = 0 to the form G( a; 1). If a n = 0, we can extract the trailing zeroes using the shuffle algebra properties before rescaling the last argument to unity. In this way we can always express any expression in terms of MPLs of the form G( a; 1) and G( 0; z), a n , z = 0. This operation is implemented in PolyLogTools via the function NormalizeG, e.g., In applications the arguments of MPLs are often complicated rational or algebraic functions. It is possible to instruct PolyLogTools to simplify the arguments of the G functions without touching the rest of the expression. The utility function GArgumentSimplify simplifies the arguments of all G, Li, PolyLog and Log objects in an expression leaving the rest of the expression unchanged.
Also the coefficients multiplying the MPLs are often functions rather than constants. It can be useful to collect the coefficients of all transcendental quantities (i.e., MPLs and MZVs). While this can in principle be done using the built-in Mathematica function Collect, experience shows that this function is rather slow when acting on large expressions. The PolyLogTools function GatherTranscendentals, which is built around the Mathematica function GatherBy, is usually much faster. It is used as shown in the following example where the functions f k (z) are analytic in a neighbourhood of z = z 0 and can be expanded into a Laurent series. In the case of MPLs the value of n is bounded by the weight. In Returns a set of all Gs in expr.

GetGArguments[expr]
Returns the set of all arguments of the Gs in expr.

GetGIndices[expr]
Returns the set of all entries in the weight vectors of the Gs in expr.

GetWeightTerms[expr,n]
Returns the weight n terms in expr.

GArgumentSimplify[expr]
Simplifies the arguments of all G, Li, PolyLog and Log functions in expr.

GatherTranscendentals[expr]
Groups terms in expr by transcendental object.

GCoefficientSimplify[expr]
Same as GatherTranscendentals[expr], but applies the Simplify function to the coefficients. Table 5.
the following we discuss how to compute the expansion of MPLs of the form G( a; z) where the weight vector is independent of z around the point z = 0. In that case the expansion can be obtained in an algorithmic way which has been implemented into PolyLogTools.
Expansions around other points z 0 = 0 can be obtained by letting z = z 0 − z and working out the functional equations that map MPLs of the form G( a; z 0 − z ) to those of the form G( a ; z ). In this way the problem is reduced to finding the expansion around z = 0. Finding the required functional equations can often be done in an algorithmic way as well, and we refer for example to refs. [28,42,67,70,83] and to Section 8.5. Let us now focus on a single MPL of the form G(a 1 , . . . , a n ; z), where we assume that the a i are independent of z. This function has a logarithmic singularity at z = 0 if and only if a n = 0. In Section 4 we have seen that we can use the shuffle algebra properties to write any MPL as a linear combination of products of G(0; z) = log z and MPLs where the last element in the weight vector is non-zero. The result of this operation (implemented in PolyLogTools via the function ExtractZeroes) is precisely the decomposition of G(a 1 , . . . , a n ; z) in eq. (8.2) into powers of logarithms multiplied by functions that are analytic at the origin. We have in this way reduced the problem to finding the series expansion of MPLs with a n = 0. We can then use eqs. (3.5) and (3.6) to obtain the expansion around z = 0 (a i = 0), where the coefficients on the right-hand side are written in terms of Z-sums [84], which can be thought of as a variant of harmonic numbers [85] depending on additional variables, The argument of ExpandPolyLogs can be any expression made out of objects for which the Series function can obtain a series representation, as well as MPLs G( a; z) where the weight vector a is independent of z. ExpandPolyLogs automatically uses ExtractZeroes internally, to extract the logarithmic terms as seen in the above example.

Differentiation and integration
Two of the basic operations on MPLs are differentiation and integration. In this section we discuss how these operations are implemented in PolyLogTools. Let us start by discussing derivatives. From the definition in eq. (3.1) it is easy to see that MPLs satisfy the differential equation ∂ x G(a 1 , . . . , a n ; x) = 1 x − a 1 G(a 2 , . . . , a n ; x) .
The previous equation is only true assuming that the weight vector a is independent of x. In applications one often encounters situations where one needs to compute partial derivatives with respect to the arguments of the weight vector. While it is possible to obtain closed formulae for these partial derivatives [3], in PolyLogTools partial derivatives of MPLs are evaluated with the help of the coproduct. Indeed, the coproduct of a derivative ∆ and the differential operator ∂ x are related through the identity [28,86] Using the fact that the Hopf algebra of MPLs is graded and connected, one can obtain the following formula for the derivative of an MPL expression F of uniform weight n, On the right-hand side of this equation the derivative only acts on MPLs of weight one, i.e., ordinary logarithms, for which the derivatives are easily computed.
Partial derivatives are implemented in PolyLogTools through the function DG. This function takes two arguments. The first argument is the expression whose derivative one wishes to compute, and the second argument is the variable with respect to which one differentiates. Moreover, this function satisfies all the basic properties of a derivative (linearity and Leibniz rule) in its first argument, and it evaluates to the built-in Mathematica implementation of the derivative whenever the first argument does not depend on G, H, HPL or Li functions. When acting on an MPL, it evaluates its derivative with respect to the second argument by means of eq. (8.7). The following example illustrates the usage of the function DG, An important property of MPLs is that they behave nicely under integration. Consider the vector space generated by all functions of the form R(z) G( a; z), where the weight vector is independent of z and R(z) is a rational function. This vector space is closed under taking primitives, i.e., for every function f (z) in that vector space there is a function F (z) in the same space such that ∂ z F = f . The function F can be found in an algorithmic way using integration by parts (see, e.g., ref. [87]). This algorithm is implemented in PolyLogTools through the function GIntegrate. Its first argument is the integrand f and the second argument the variable z. Whenever the weight vector depends on the integration variable z, or if the integrand involves non-rational algebraic functions of z, the algorithm does not converge and returns an unevaluated expression. We emphasise that this is not a limitation of our implementation, but in these cases the space of MPLs may not be enough to perform all the integrals and more general classes of functions, e.g. of elliptic type, may be required. In applications one is usually not directly interested in the primitive of a function, but in its integral over a certain range. One can easily evaluate definite integrals of MPL expressions by first computing a primitive and then evaluating the primitive at the endpoints of the integration range. Doing so often results in spurious singularities that cancel in the final result. For example, consider the integral We can formally do the integral by computing a primitive and evaluating it at z = b and z = 0, This last expression is ill-defined because G(b, . . . ; b) is divergent. In Section 4 we have seen how we can use the shuffle algebra properties to replace divergent MPLs by their shuffle regularised versions, In cases where the integral has not only spurious logarithmic end-point singularities but also poles, replacing MPLs by their shuffle regularised version is not sufficient, and one needs to carefully expand the MPLs around the pole. This can for example be done using the ExpandPolyLogs function described in the Section 8.2.

Numerical evaluation of MPLs
It is important to be able to evaluate MPL expressions, and for this reason a considerable effort has been put by the community in developing fast and reliable computer libraries for the numerical evaluation of (some classes of) MPLs, see for example refs. [67,68,83,[88][89][90][91][92][93][94].
Given the vast amount of publicly available codes for the evaluation of MPLs, Poly-LogTools does not have its own numerical routines for the evaluation of MPLs, but it relies on the Hpl and GiNaC [95] packages. Since the Hpl package is loaded together with PolyLogTools, HPL functions can be evaluated numerically as described in the manual of that package [67,68], and we will not discuss it here any further.
GiNaC is a C++ library for numerical and symbolic computations in high-energy physics. It contains an implementation of the algorithms of ref. [83] for the numerical evaluation of MPLs. The command-line utility ginsh allows the user to start an interactive shell-version of GiNaC. We refer to the GiNaC manual for more details [96].
PolyLogTools contains an interface that allows the user to evaluate valid MPL expressions using the interactive ginsh environment from within Mathematica. The input expression is transformed into a string that is a valid GiNaC expression. This string is written to a temporary file that is then piped through the ginsh shell command via the Run function in Mathematica. The result returned by ginsh is written to a temporary file. The content of this file is imported back into PolyLogTools and returned by the Ginsh function. The temporary files are then deleted from the disc. It is possible to instruct PolyLogTools not to delete the temporary files by setting the option Debug to True (the default is False) as shown in the following example: In Let us make an important comment at this point. MPLs are multi-valued functions, and so care is needed when evaluating MPLs to obtain the correct numerical value. When calling GiNaC, PolyLogTools does not make any assumption on the branch cuts of the functions and it solely relies on the choices made by GiNaC. PolyLogTools merely translates an expression into a format that can be processed by GiNaC, and it is up to the user to make sure that he or she understands the branch cut structures of the functions that are evaluated before calling GiNaC. Closely related to numerical evaluation is fitting numerical constants using the PSLQ algorithm [97]. PolyLogTools provides the user with two possible ways to fit transcendental constants. The function RunPSLQ is based on the implementation of the PSLQ algorithm in Mathematica by P. Bertok. The source code of this implementation is publicly available from the Wolfram Mathematica library [98] and is shipped together with PolyLogTools. In more recent versions, Mathematica provides a built-in routine to find integer relations via the FindIntegerNullVector function. The PolyLogTools function PSLQFit relies on this implementation.
Let us illustrate the use of these functions with an example. The transcendental constant G(0, 1; 1/2) = −Li 2 (1/2) can be expressed as a linear combination of log 2 2 and π 2 . The coefficients of this linear combination can be found using the PSLQ algorithm from the numerical value of G(0, 1; 1/2), The last argument of the RunPSLQ and PSLQFit functions is the number of digits that should be taken into account in the fit. It is usually advisable to evaluate the fit with several different numbers of digits, in order to ensure that the results are stable and that the identified rational coefficients are not accidental.

Fibration bases
An important part of applying MPLs to the computation of Feynman integrals is determining a combination of MPLs whose symbol matches a given (integrable) symbol tensor. While in general this is a highly complicated task for which no algorithmic solution is known, under certain conditions on the symbol alphabet it is possible to determine an MPL expression with the given symbol in an algorithmic way. In this section we describe such a criterion, and the corresponding algorithm implemented into PolyLogTools. Consider a symbol alphabet A. We assume that all letters are non-constant rational functions in a set of variables x 1 , . . . , x m . Without loss of generality we may assume that all letters are irreducible polynomials p i over Z.
Next, we assume that all letters are linear in one of the variables, which we choose to be x 1 . Then all the letters in A can be written in the form Following ref. [99], we define the set A x 1 as the set consisting of all irreducible non-constant polynomial factors in a i , b i and a i b j − a j b i . If all the polynomials in A x 1 are linear in one of the variables, say x 2 , we can iterate the procedure and construct the set A x 1 ,x 2 . We say that a symbol is linearly reducible [99] if we can find an ordering of the variables (for example the natural ordering (x 1 , . . . , x m )), which allows us to iterate this procedure until the end, i.e., we can find at every step a variable x k+1 such that all polynomials in A x 1 ,...,x k are linear in x k+1 . Given an integrable symbol tensor S that is linearly reducible with respect to the ordering (x 1 , . . . , x m ), one can find a function of the form such that S(F ) = S and the weight vector a k,i k only involves rational functions in the variables (x k+1 , . . . , x m ) [99]. We call the MPLs that appear on the right-hand side of eq. (8.13) a fibration basis for the alphabet A and the ordering (x 1 , . . . x m ). The function F can be constructed in an algorithmic way, cf. [30,43,45,46,99,100]. The algorithm implemented in PolyLogTools via the function FiberSymbol follows closely the one described in ref. [30]. The function FiberSymbol takes two arguments. The first one is an integrable and linearly reducible symbol tensor, and the second one is the list of variables with the chosen ordering. This is illustrated in the following example: In If the input symbol tensor is not integrable or not linearly reducible with respect to the ordering in the second argument, then the algorithm fails and a warning will be shown. So far we have only described how to use fibration bases to find a function whose symbol matches a given symbol tensor. In applications it is often useful to write a given MPL expression in terms of a fibration basis. More precisely, consider an MPL expression F with the property that all MPLs have a linearly reducible symbol alphabet with respect to a same ordering of the variables. Then the user can use PolyLogTools to write F in terms of the fibration basis with respect to that ordering. We illustrate this on the following example: Unlike the output of FiberSymbol, the output of ToFibrationBasis involves also terms proportional to MZVs that are not detected by the symbol. These terms are reconstructed using the coproduct combined with the numerical fitting technique described in ref. [30], where the MZVs are reconstructed by evaluating constant combinations of MPLs numerically at a single point. So far this fitting technique has been implemented through weight six, and therefore the use of ToFibrationBasis is limited to expression of weight up to five. Since MPLs are multi-valued functions, the result of the fit can depend crucially on the numerical value used in the fit. By default, each variable is assigned a random value in the range It is strongly recommended that the user always chooses the numerical value in a way which reflects the applications which he or she has in mind (e.g., in applications the variables are often related to kinematic variables, which take values in a certain range). This is particularly important when the symbol contains letters such as x − y that could lead to polylogarithms that can develop imaginary parts depending on the relative values used for fitting the constants. For large expressions, the reduction to a fibration basis can take a considerable amount of time. It it possible to save intermediate results to the disc by setting the option Save -> "file.m". Intermediate results can later be read in via the option Input -> "file.m".
The function ToFibrationBasis is one of the most important features of PolyLog-Tools. In particular, fibration bases play an important role in the computation of Feynman parameter integrals via direct integration, using ideas similar to those described in refs. [30,43,45,46,99,100]. Indeed, if the set of polynomials in the denominator of an integrand is linearly reducible with respect to a given ordering, then we can perform all the integrations one-by-one in that order. At each step one can use the function ToFibrationBasis to write all MPLs in the integrand in a form where the integration variable appears in the last argument. Once that is achieved, the integral can be performed easily in an algorithmic way using the GIntegrate function. Another important application of fibration bases is the derivation of transformation formulae for MPLs. As an example, consider the function G(−1, −y; x), and imagine we want to obtain its series expansion around y = 1/2 with 0 < x < 1. In Section 8.2 we have seen how to use the function ExpandPolyLogs to expand MPLs of the form G( a; z) around z = 0. Hence, if we let y = 1/2 − z and we pass to a fibration basis with respect to the ordering {z, x}, we can easily obtain the desired expansion. The corresponding piece of code is shown below: Let us conclude by mentioning that the criterion of linear reducibility described above is only a sufficient, but not necessary, condition to find a function that matches a given symbol tensor. The sets A x 1 ,...,x k provide an upper bound on the singularities that can appear after the variables x 1 , . . . , x k have been integrated out. There are more refined criteria that allow one to obtain a more refined bound on the singularities, and therefore to find fibration bases for larger classes of functions [5,99]. These more refined criteria, however, have not yet been implemented into PolyLogTools.

Single-valued MPLs
Just like the logarithm function, MPLs are multi-valued functions. It is possible to define a variant of MPLs that are real-analytic and single-valued, while preserving most of their algebraic properties. The price to pay is that these functions are no longer holomorphic, i.e., they depend explicitly on the complex conjugated variables. More precisely, single-valued MPLs are combinations of MPLs and their complex conjugates such that all discontinuities cancel. The simplest example of such a function is the single-valued version of the logarithm, which is simply given by the logarithm of the absolute value of its argument, log |z| 2 = log z + logz . (9.1) It is possible to determine the combinations that lead to single-valued MPLs in an algorithmic way. This was done for the first time in ref. [101], where the single-valued versions of HPLs whose weight vectors only contain 0's and 1's have been defined by constructing single-valued solutions to a unipotent differential equation. This construction was extended in ref. [102] to general classes of hyperlogarithms. In refs. [8,86] it was shown how to define general single-valued MPLs and also MZVs. Single-valued MPLs are not just of interest in pure mathematics, but they also appear in loop computations. In particular, they show up in the computation of Feynman integrals with massless propagators and three off-shell external legs [29], as well as in the computation of conformal four-point functions in four dimensions [59,[103][104][105][106]. In addition, single-valued versions of MPLs are known to describe the multi-Regge limit of scattering amplitudes in planar N = 4 Super Yang Mills [55,107] and the high-energy limit of the dijet cross section in QCD [56,108]. They also appear in the analytic result for the three-loop corrections to the soft anomalous dimension [109,110].
In the remainder of this section we present the implementation of single-valued MPLs in PolyLogTools. If expr involves only MPLs of the form G( a; x), with a independent of x, then expr is expanded into a series around x = 0 up to order n. Ginsh [expr, list] Uses GiNaC to evaluate expr numerically. list is a replacement list that specifies the numerical values that should be assigned to all the variables in expr. Ginsh has an option Debug. If set to True (which is the default), the temporary files containing the GiNaC code are not deleted from the disc. RunPSLQ [num, list, prec] Uses the PSLQ implementation by P. Bertok with precision prec to express num as a rational linear combination of the quantities in list. The argument num must be real. PSLQFit [num, list, prec] Uses the FindIntegerNullVector function to express num as a rational linear combination of the quantities in list. The argument num can be either real or complex. The integer prec specifies that all floats should be interpreted as having a precision of prec digits.

DG[expr, x]
Computes the partial derivative of expr with respect to x. GIntegrate [expr, x] Computes the primitive of expr with respect to x. Only expressions that contain rational functions of x and MPLs of the form G( a; x), with a independent of x are allowed inside expr.

ShuffleRegulate[expr]
Replaces all divergent MPLs in expr by their shuffleregularised version.

Single-valued MPLs in PolyLogTools
In refs. [8,86,101,102] (see also ref. [55]) a map s was constructed which associates to G( a; z) its single-valued version G( a; z). This map can be given explicitly in terms of the coproduct and the antipode on MPLs (see Section 6), whereS is related to the the antipode of the complex conjugate of x, up to a sign, and |x| is the weight of x. The map s is obviously linear, and it also preserves the multiplication, The image of any MPL expression under s is single-valued. Note that s does not only act on functions, but also on numbers. In particular, it sends to zero all powers of π, and acts Returns a combination of MPLs in a fibration basis with respect to the variables and the ordering specified in list whose symbol matches sym. The symbol sym is assumed integrable and linearly reducible.

ToFibrationBasis[expr, list]
Returns expr in a fibration basis with respect to the variables and the ordering specified in list. This function has several options decribed below. FitValue Option of ToFibrationBasis. A list of replacements which specify the numerical values of the all the variables used in the numerical fit. The default is Automatic, which assigns random values between 0 and 1 to each variable. Save Option of ToFibrationBasis. The value is a string (default: the empty string). If the string is non-empty, then intermediate results of ToFibrationBasis are written to the file whose name is the specified string.

Input
Option of ToFibrationBasis. The value is a string (default: the empty string). If the string is non-empty, then the file specified by the string is read in. ProgressIndicator Option of ToFibrationBasis. If set to True, a dynamically updated text indicates how many MPLs still need to be converted to the fibration basis. on odd MZVs of depth one in a simple way, s(π) = 0 and s(ζ 2n+1 ) = 2ζ 2n+1 . (9.5) The single-valued MPLs G(a 1 , . . . , a n ; z) are represented in PolyLogTools by the symbols cG[a1,...,an,z]. The map s is implemented as the function SV, which can act on any MPL expression and returns its single-valued version. We illustrate this with the following example, It is possible to work with the cG functions in very much the same way in PolyLog-Tools as with ordinary MPLs. All functions to manipulate MPL expressions described in Sections 4 and 8.1 can be used in the same way with G and cG functions. In particular, since the map s respects the multiplication (cf. eq. (9.4)) single-valued MPLs form a shuffle algebra, and we can decompose them into a Lyndon word basis, extract trailing zeroes, etc., just like for the G functions. It is also possible to differentiate and integrate single-valued MPLs. The map s commutes with holomorphic differentiation [86,101,102], ∂ z s = s∂ z , so that we can compute derivatives in the same way as for ordinary MPLs. In particular, the function DG can be used just like for G functions to compute the holomorphic derivatives of single-valued MPLs. It is also possible to evaluate single-valued MPLs by means of the Ginsh function, for example, While the map s commutes with holomorphic differentiation, the situation is slightly more subtle when it comes to integration. Indeed, let us compute the following (holomorphic) primitive, and use the fact that single-valued MPLs can be written as ordinary MPLs for which we know how to compute primitives. We find We see that the primitive of a single-valued function, if computed in a naive way, is not single-valued, and the difference is precisely the antiholomorphic function G(0, 1;z). We should keep in mind, however, that the primitive of a function is not unique. In particular, we can add to a holomorphic primitive any antiholomorphic function, and so we can find a single-valued primitive, albeit not the one that we would have naively constructed. This is a general fact: similarly to ordinary MPLs, the space of single-valued MPLs (multiplied by rational functions) is closed under taking primitives. The single-valued primitive of a singlevalued MPL expression can be computed in PolyLogTools via the cGIntegrate function. This function works in the same way as its non-single-valued analogue GIntegrate described in Section 8.3, and so we do not describe it here any further. In Section 6.2 we have seen that ordinary MPLs modulo their discontinuities form a Hopf algebra. Since we have a map that assigns to an MPL its single-valued version, it is natural to ask if the single-valued MPLs themselves form a Hopf algebra. This is indeed the case, and single-valued MPLs form a graded and connected Hopf algebra with a coproduct ∆ sv and antipode S sv given by the same formulas as the coproduct ∆ and the antipode S on MPLs of Section 6.2, except that all MPLs need to be replaced by their single-valued versions. The coproduct ∆ sv and the antipode S sv are implemented in PolyLogTools via the functions DeltaSV and AntipodeSV. These functions are completely analogous to the functions Delta and Antipode of Section 6.2.
Since both coproducts ∆ and ∆ sv are computed via very similar formulas, we find it important to quickly compare the two coproducts. We do this on the example of the function G(0, 1; z) = G(0, 1; z) + G(1;z) G(0; z) + G(1, 0;z) .
The coproduct ∆ sv acts on single-valued MPLs through the same formula as ∆ on ordinary MPLs. We the find ∆ sv (G(0, 1; z)) = G(0, 1; z) ⊗ 1 + 1 ⊗ G(0, 1; z) + G(1; z) ⊗ G(0; z) . (9.8) In particular, by construction, ∆ sv only involves single-valued MPLs. The coproduct ∆, instead, acts on the MPLs on the right-hand side of eq. (9.7). We find We see that for ∆ only the first entries in the coproduct are single-valued. The two coproducts are thus genuinely different.
Let us now discuss single-valued fibration bases. Since the map s preserves the algebra structure, single-valued MPLs satisfy the same relations as their non-single-valued analogues (with all factors of π removed). This implies in particular that single-valued MPLs can be expressed in terms of a fibration basis under the same conditions and in the same way as ordinary MPLs, simply by acting with s on the corresponding relation among ordinary MPLs. For this reason, the function ToFibrationBasis can be applied in the same way to single-valued MPLs as to ordinary ones as described in Section 8.5.
We conclude by mentioning that more general classes of single-valued MPLs show up in Feynman integral computations, where the symbol letters are neither holomorphic nor anti-holomorphic [29,56,59,104,106]. While it is also possible to use PolyLogTools to construct these functions using the algorithm of ref. [106] (cf., e.g., ref. [29,56,104]), this algorithm is not implemented in PolyLogTools in an automated way. We mention, however, that these functions can be obtained in an automated way from the Maple package Hyperlog Procedures [105].

The Lie coalgebra of clean single-valued functions
Since single-valued MPLs form a graded and connected Hopf algebra, we can use the results of Section 6.3 and construct a projector P sv to the Lie coalgebra of indecomposables of The single-valued MPL G(a 1 , . . . , a n ; z). cC[a1,...,an,z] The clean single-valued MPL C(a 1 , . . . , a n ; z).

SV[expr]
Applies the map s to expr. cGToG [expr] Replaces all cG functions in expr by ordinary MPLs and their complex conjugates.

cCTocG[expr]
Replaces all cC functions in expr by single-valued MPLs.

cCToG[expr]
Equivalent to cGToG[cCTocG [expr]]. cGIntegrate [expr, x] Computes the single-valued primitive of expr with respect to x. Only expressions that contain rational functions x and MPLs of the form G( a; x), with a independent of x are allowed inside expr.

AnitpodeSV[expr]
Applies the antipode S sv to expr.

ProductProjectorSV[expr]
Applies the projector P sv to expr.
this Hopf algebra. The projector is constructed recursively following eq. (6.15) using the coproduct. We first define the map where x has weight n and we let P sv (x) = 1 n R sv . In ref. [71] the following clean single-valued MPLs have been defined, C(a 1 , . . . , a n ; z) = P sv (G(a 1 , . . . , a n ; z)) . (9.11) The images of single-valued MPLs under the projector P sv have very interesting properties. In particular, they satisfy clean functional relations, i.e. the same relations as the G functions, but with all product terms removed [71]. The clean single-valued functions are represented by the symbols cC[a1,...,an,z], and the projector P sv is called ProductProjectorSV. The clean single-valued functions can be expressed in terms of singlevalued MPLs. This can be achieved via the function cCTocG. Moreover, all the functions from Sections 4, 8.1 and 8.5 can also be applied to the clean versions. Finally, there is a cobracket δ sv = (P sv ⊗ P sv )(1 − τ )∆ sv acting on these functions. The cobracket is implemented via the function CobracketSV.

Validation
PolyLogTools has already been applied to many computations involving MPLs that have led to publications in peer-reviewed journals, both in physics and in mathematics. In this section we review these applications. The reason for doing this is twofold. First, these applications show the versatility of the code, and they can serve as examples to the variety of problems to which PolyLogTools can be applied. Second, these examples serve at the same time as validation of the code.
PolyLogTools was applied for the first time in ref. [28], where the coproduct on MPLs was used to simplify the two-loop amplitudes for the production of a Higgs boson in association with three partons [111].
The most prominent application of PolyLogTools is the computation of the N 3 LO corrections to Higgs production in gluon-fusion [33,35], where the package was used extensively to compute boundary conditions for differential equations [30,112] and to evaluate phase space integrals for single-emission contributions [31,113]. It was also used in ref. [114] to evaluate convolution integrals that appear in the mass-factorisation formulas beyond LO. Further phenomenological applications to Higgs physics include the computation of the two-loop amplitudes for Higgs production in the Standard Model Effective Field Theory [115] and the top-Yukawa contributions to bbH production [116].
PolyLogTools was successfully used to compute integrated counterterms for the CoLoRFulNNLO subtraction method [117][118][119]. It was also used in the context of Lattice QCD in the computation of certain one-loop integrals that appear in the massive momentum-subtraction scheme [120]. More recently, the package was used to manipulate MPLs that appear in the computation of four and five-gluon scattering using the numerical unitarity approach [121][122][123][124].
The package was not only applied to higher-order computations relevant to collider physics, but it has also seen various more formal applications. In particular, it has played a crucial role in the analytic computation of the three-loop soft anomalous dimension matrix [109]. It has also been used to compute certain two-loop triangles [29] and conformal four-point functions [104] in terms of single-valued MPLs. Single-valued MPLs are known to show up in multi-Regge kinematics, and PolyLogTools was also successfully applied in that context, both in planar N = 4 SYM [55,57,58,107] and in QCD [56,108]. Recently, it was used to manipulate analytic expressions for two-loop form factors in N = 4 SYM [125,126]. The package was also used to compute and analyse the cuts of one and two-loop integrals [61,127,128] and to study their relationship to the coaction of Feynman integrals [62,63]. Finally, (a private extension of) PolyLogTools was used to evaluate various Feynman integrals that evaluate to elliptic generalisations of MPLs [129][130][131][132].
PolyLogTools was not only used to perform research in theoretical and mathematical physics, but it has led to new results in pure mathematics. It was applied in refs. [11,12] to study reduction identities for MPLs to lower depth, and in refs. [9,10] to study some properties of MPLs of weight four, in particular to derive a novel functional relation for MPLs of weight four involving more than 100 terms. Finally, it was used in ref. [71] to define and study some of the clean single-valued MPLs reviewed in Section 9.2.

An example calculation
In the following we present a sample computation using PolyLogTools that illustrates the most important features of the package. We compute the one-loop four-mass scalar box integral shown in fig. 11. This integral is finite and furthermore one of the simplest and most prominent examples of dual-conformal integrals appearing for example in N = 4 SYM. We stress that these properties are by no means prerequisites for the use of PolyLogTools, in particular the package has been used to compute many infrared-divergent and non-dualconformal integrals as pointed out in the previous section. However, this integral allows us to focus on the features of PolyLogTools and elaborate on some potential obstructions in the calculation and how to avoid them using PolyLogTools. This example is also described in the notebook integration manual.nb distributed with the source code of the package. The four-mass box integral can be written as ( , a)( , b)( , c)( , d) . (11.1) Here we have used the notation (a, b) = (x b − x a ) 2 that is inspired by the embedding space formalism. The four dual points a, b, c, d are non-light-like separated, so that we can define the two finite dual-conformal cross ratios u = (a, b)(c, d) (a, c)(b, d) and v = (a, d)(b, c) (a, c)(b, d) , (11.2) as well as the Gram determinant While it is possible to linearize the Gram determinant by introducing new variables, we will refrain from doing so here, in order to illustrate how to handle the appearance of square roots in the calculation. The starting point for our calculation is the Feynman parametrization of the loop integral, that we write as [50], Since none of the resulting polylogarithms have a 0 in their last index, we can immediately set α 1 to zero to see that the integral vanishes on the lower boundary. To derive the result at the upper bound, we need to map +∞ to zero, which we can do by letting α 1 → 1/t. Now we can instruct PolyLogTools to derive the functional identities that are required to express the primitive as a function of t using (11.5) We see that there are no logarithmically divergent polylogarithms with argument t and so we could just set t = 0 to obtain the value of the integral at the upper bound. It is usually advisable to instead expand the polylogarithms to zeroth order around the desired point. This is more robust as it automatically enforces the cancellation of spurious logarithmic divergences or poles. This can be achieved using Here we see now a potential obstruction that is not automatically resolved by PolyLog-Tools: The denominator is quadratic in the next integration variable. In general this can be a fundamental obstruction, in the sense that that it may not be possible to evaluate the integral in terms of MPLs alone. In this case, however, there is only one integration left to do, so the result should be expressible as polylogarithms with algebraic arguments.
PolyLogTools cannot automatically resolve this integral though, as it expects to be able to partial fraction the input into terms with only linear denominators. We can however manually solve the quadratic equation in the denominator and use it to reexpress the denominator in factorized form. At this point we are in principle done and could perform the next integral. We have to be careful to hide the square root from Mathematica as the routines for partial fractioning will otherwise restore the original quadratic form. Here we can just identify the root with the gram determinant ∆ 6 that we introduced before and thus define the new denominator as In [7]:= De = 2*(1/2*(1-u+v-d6)+a [2]) *(1/2*(1-u+v+d6)+a [2]); Now the next integration can be performed by invoking The resulting primitive is a pure function of weight two that we refrain from spelling out here. Again we need to evaluate the primitive at the boundaries of integration. The lower boundary is immediately found using: In [9]:= at0 = ExpandPolyLogs[P,a [2],0,0] Out [9]:= 0 In order to obtain the value at infinity, we once again need to invert the argument of the polylogarithms and use ToFibrationBasis to derive the required functional equations. Rather than calling ToFibrationBasis on the entire expression, we can also call it on a list of polylogarithms to obtain the individual functional equations. This can be particularly useful when we have a basis of functions for which we want to derive functional equations for that can then be stored. It is also useful if we want to inspect the individual functional equations to make sure that our choice of FitValue has not introduced any spurious iπ terms. In the present case there is the potential for such terms, as the polylogarithms contain spurious branchpoints for u − v = 0 and as such the ordering of the variables that is employed when numerically fixing the branches of the logarithms becomes relevant. When we derive the functional equations in the following way we can easily verify that the obtain expressions maintain manifest reality:

Conclusion
In this article we have reviewed the mathematical foundations of multiple polylogarithms and have documented their implementation in the Mathematica package PolyLog-Tools. Multiple polylogarithms have become ubiquitous in many areas of high-energy physics and the algorithms presented in this paper and implemented in PolyLogTools provide a powerful and flexible way to handle expressions involving multiple polylogarithms. PolyLogTools has already served as the backbone of many recent calculations in high energy physics. Its public release accompanying this paper will enable even more studies of the mathematical structure of scattering amplitudes and calculations of multiloop amplitudes and cross sections. It also serves as a well-tested reference implementation for more specialized and optimized implementations that might be developed to handle particular situations.