Glue-and-Cut at Five Loops

We compute epsilon-expansions around 4 dimensions of a complete set of master integrals for momentum space five-loop massless propagator integrals in dimensional regularization, up to and including the first order with contributions of transcendental weight nine. Our method is the glue-and-cut technique from Baikov and Chetyrkin, which proves extremely effective in that it determines all expansion coefficients to this order in terms of recursively one-loop integrals and only one further integral. We observe that our results are compatible with conjectures that predict $\pi$-dependent contributions.


Introduction
Despite ongoing efforts do develop alternative representations, practical calculations of scattering amplitudes in quantum field theories are still described mostly through Feynman diagrams. As one studies observables at higher orders in perturbation theory in this framework, one needs to evaluate a huge number of corresponding many-loop Feynman integrals. The predominant approaches exploit integration by parts (IBP) identities in momentum space [1,2] to express all required integrals as linear combinations of a much smaller number [3][4][5] of basis integrals, called master integrals (MI). Although this reduction procedure may become challenging in practice, it is thus in principle sufficient to compute and tabulate only one such set of master integrals for a given problem. Due to the presence of divergences in the integrals, calculations are performed in dimensional regularization in d = 4 − 2ε dimensions [6,7], and the master integrals have to be expanded to a sufficiently high order in the ε-expansion.
In this paper we focus on the evaluation of massless propagator-type integrals in momentum space, also known as "p-integrals". The importance of p-integrals is demonstrated by the success of their application to determine arbitrary renormalization group functions up to four loops [1,2], in scalar φ 4 theory up to 6 loops [8,9], and recently to five loops in QCD [10,11]. Alternative effective routes to renormalization group functions include massive vacuum integrals [12][13][14] and position-space techniques [15]. In this comparison, the coefficients of the ε-expansions of p-integrals tend to be rather simple and are often expressible exactly in terms of Riemann zeta values. On the other hand, masslessness introduces infrared divergences, which have to be carefully separated from the sought-after ultraviolet singularities. This is achieved with the R * -operation [16][17][18], which is now in a very mature state [19,20] and provides a technique that expresses L-loop renormalization group functions in terms of poles and finite terms of the ε-expansions of p-integrals with only L − 1 loops.
As a step towards an extension of the currently very good understanding of the master p-integrals [21][22][23] and their IBP reduction [24] at four loops, we compute in this paper a complete set of master p-integrals at five loops, which should be useful in future determinations of renormalization group functions at the six-loop order. Our results are also interesting from a structural point of view and give further evidence for general patterns in the coefficients of ε-expansions of p-integrals, and therefore, about the coefficients of renormalization group functions [25,26]. Concretely, we confirm that finite five-loop p-integrals in d = 4 dimensions are linear combinations of multiple zeta values of transcendental weight at most nine, and that, in a well-defined sense, all π-dependent contributions can be predicted by the π-free terms. We supply our results for the ε-expansions of our master integrals and other results in the ancillary files that supplement this article, the use and content of each file is described in the README file.
All five-loop p-integrals can be obtained by contracting edges of one of 64 maximal sectors 1 , by which we mean diagrams with only cubic vertices and 14 internal propagators, as illustrated in Figure 1. Note that the planar graphs can also be viewed as propagator integrals in position-space, and we completely determined those in previous work [27]. Our calculation presented here confirms these earlier findings for the planar p-integrals, but provides also the expansions for the non-planar integrals.
In order to determine all five-loop p-integrals we follow the algorithm introduced at four loops by Baikov and Chetyrkin [22], which combines the glue-and-cut (GaC) symmetry [35] of massless p-integrals with IBP reductions to bootstrap the determination of master integrals. This elegant algorithm manages to get by without having to actually compute a lot of complicated integrals explicitly. Instead, it produces so many relations between the coefficients of the ε-expansions of the p-integrals, that ultimately all get reduced to coefficients of only a small number of particularly simple integrals. Concretely, the constraints turned out to be so strong that the input of recursively one-loop integrals (see Figure 7), which are easily evaluated to all orders in ε, was sufficient to determine all master integrals at four loops [22].
The bottleneck in this strategy resides mostly in the efficiency of IBP reduction routines, and only recent progress in that field has allowed us to finally apply the GaC strategy at five loops and to fix all p-integrals at that order. We find that it suffices to input the expansions of 21 recursively one-loop master integrals, together with only one further integral. In fact, this single additional input datum can be chosen in form of a very simple product integral, see Figure 8. This demonstrates strikingly that the remarkable power of the GaC method persists at the five loop order, circumventing completely the daunting direct evaluation of a large number of very complicated Feynman integrals.
The paper is organized as follows. In section 2 we review the main ideas of the method and describe how we obtain the relations needed for bootstrapping the p-integrals. In section 3 we describe an algorithm to determine the convergence of vacuum diagrams. In section 4 we describe how we performed the IBP reductions needed while in section 5 we describe how we generated the needed vacuum diagrams and how we constructed and solved the equations needed for the bootstrap procedure. Finally in section 6 we present and discuss our results.

Review of the algorithm
All five-loop massless propagator integrals can be represented by families of the form where i denote the loop momenta, and each D i is a quadratic form in these loop momenta and the external momentum p. There are 64 such families, indexed by 1 ≤ k ≤ 64, each corresponding to a different cubic graph with 14 internal lines, and D 1 , . . . , D 14 encode precisely the inverse propagators associated to these lines. The corresponding indices a 1 , . . . , a 14 will always be positive or zero. The remaining six quadratic forms D 15 , . . . , D 20 will only appear with non-positive indices a 15 , . . . , a 20 ≤ 0 and are needed to encode propagator integrals with numerators. See Figure 4 and Table 1 for the example of family k = 46. We work in dimensional regularization with d = 4 − 2ε, so every propagator admits a Laurent expansion whose coefficients c n ∈ R are the numbers we want to determine. Throughout we assume a positive definite (Euclidean) metric and we set p 2 = 1 in our calculations. The dependence on the external momentum p is a power law and thus completely determined by the exponent The algorithm of [22] can be summarized in the following steps: 1. In each family P k , choose a sufficient number of integer seeds a ∈ Z 20 and reduce each of those p-integrals P to (any preferred choice of) master integrals M i : 2. Expand both sides of (4) in ε to express the Laurent coefficients c n (P ) of each seed p-integral P = P k (a; d) in terms of the corresponding coefficients c n (M i ) of the master integrals.
3. Enforce c n (P ) = 0 for all n < −5. 4. Enforce c n (P ) = 0 for all n < 0 in the case of finite p-integrals.
5. Enforce the identities c 0 (P ) = c 0 (P ) for all pairs of p-integrals P and P that are related by the glue-and-cut symmetry.
The glue-and-cut symmetry of massless p-integrals explains why sometimes different finite p-integrals evaluate to the same number in four dimensions. Diagrammatically, this process takes an L-loop p-integral P and glues its two external legs together, to form an (L + 1)-loop vacuum diagram V (see Figure 2 for an example). Conversely, cutting the glued line from V , we recover P . But if we cut another line of V , we may produce a different p-integral P . The glue-and-cut symmetry [22] states that: , and if V is dimensionless, then any other cut P (ε) of V is also finite in d = 4, and P (0) = P (0).

Dimensionless here means that the additional edge in
In fact, as long as this condition is fulfilled, the symmetry holds for arbitrary values of the indices a and for arbitrary ε: P (ε) = P (ε) [30]. Indeed, any p-integral P obtained from cutting V , is equal to the residue of the vacuum integral V at ω ε (V ) = 0. 2 However, note that a 0 will depend on ε, so in specializing to ε = 0 we gain that a 0 becomes an integer, and hence P (0) will again be a p-integral with integer indices.
This symmetry generates an array of identities between c 0 -coefficients of different pintegrals. Once they are reduced to master integrals, we get a linear system of equations which in general allows us to relate the expansions of some master integrals in terms of the others. The hope is that the remaining undetermined coefficients can all be computed without difficulty. The implementation of the algorithm for L-loop p-integrals starts then with the enumeration of (L + 1)-loop vacuum diagrams V . Those are written as with D i a linear combination of scalar products i · j , which denotes either an inverse propagator (a i > 0) or a numerator (a i < 0). The superficial degree of divergence in four dimensions is then According to Theorem 1, we consider only diagrams with ω 0 (V ) = 0, which means that they must contain at least 2(L + 1) propagators. The six-loop maximal sectors (diagrams with cubic vertices) have 15 denominators D 0 , . . . , D 14 (corresponding to the 15 edges in the graph), which means that in that case the candidates for cutting are obtained by adding a 15 + . . . + a 20 = 3 numerators. For graphs with less denominators (a i = 0 for one or more i ≤ 14), we pick accordingly less numerators. Apart from thus assuring ω 0 (V ) = 0, we need also to restrict to diagrams which have no subdivergences, so that the p-integrals resulting from cutting are all finite. Our algorithm for checking the convergence properties of vacuum diagrams will be explained in section 3.
The following step is the cutting of the valid vacuum diagrams in all possible ways. The obvious identities are extracted by cutting along the existing denominators. However, one can also obtain constraints by splitting any of the vertices with valence k > 3, see Figure 3. Starting for instance with a valence-k vertex, we can produce two vertices of valences (k −j +1) and (j + 1), as long as j, (k − j) ≥ 2. These vertices are then connected by a new propagator (with index a i = 0) which we can cut, thus producing an additional identity.
Finally, the last step is the IBP reduction (4) to master integrals. At five loops this is quite an arduous task, further details on the strategy followed are given in section 4. Once the reduction is performed and all master integrals have been identified, one just has to perform the symbolic ε expansion (2), impose finiteness of all p-integrals obtained from cutting, and equate their finite orders. Note that we also impose c n (P ) = 0 for all n < −5, since even for divergent p-integrals, at five-loops they can at worst have poles of order five.

Convergence
Consider a scalar p-integral, that is, some P k (a; d) as per (1) where a 15 = · · · = a 20 = 0 (no numerators). Such a p-integral is convergent if, and only if, we have ω ε > 0 and the conditions  ω ε if g connects the external legs and 0 otherwise, are fulfilled for all non-empty proper subsets g {1, . . . , 14} of the edges of the defining graph [36,37]. Here, we denote by L g the loop number (first Betti number) of the subgraph g. Alternatively, consider the vacuum graph V obtained by gluing the external legs of P into an additional edge, with a 0 chosen such that ω ε (V ) = 0. Then the conditions above can equivalently be stated more symmetrically as for all non-empty subgraphs g ⊆ {0, . . . , 14} of V . It is not necessary to check all of these conditions separately, however, as many of them turn out to be redundant. It suffices to verify the constraint only for those subgraphs g that have the property that both g, and its quotient V /g, are biconnected [38,39].
This provides an efficient way to select the finite integrals from a list of scalar p-integrals: For each of the 14 cubic graphs V with six loops and 15 edges corresponding to a top-level vacuum integral, we compute once and for all the list of all subgraphs g with the property that g and V /g are biconnected. This yields a short list (in the worst case there are 107) of inequalities e∈g a e > d/2 · L g , which can be used to quickly determine if a scalar p-integral obtained by cutting some edge of V is finite or not. For the remaining six-loop vacuum integrals (those with less than 15 edges), we proceed analogously.
With this control on scalar integrals, we can also determine the finiteness of integrals with numerators. Indeed, integrals with numerators can be viewed as linear combinations of scalar integrals in higher dimensions [40]. Concretely, let us introduce Schwinger parameters x i for each denominator D i in (1). Then, the quadratic form 5 ) determines a 5 × 5 matrix A, a 5-vector B and a scalar C (all of these depend on x). Define the polynomials U := det A and F : The parametric representation [37,41] of the integral family (1) then takes the form It is apparent that, after carrying out the derivatives with respect to x j for the numerators (a j < 0), and bringing the integrand on a common denominator, we can write P k (a; d) = n b n (ε)P k (a (n) ; d ) as a finite linear combination of p-integrals in some raised dimension propagators ISPs numerators graphs finite integrals   15  6  56  14  347  14  7  28  22  248  13  8  8  32  206  12  9  1 17 15 Table 2: Counts of non-isomorphic vacuum graphs (see Figure 6) by number of edges (propagators). The third column gives the number (10) of logarithmically divergent numerators.
The convergence of the integral (7) is equivalent to the convergence of each constituent P k (a (n) ; d ). We therefore apply the convergence criteria for scalar integrals as discussed above, in order to decide whether or not a given p-integral P k (a; d) is finite, i.e. convergent in d = 4 dimensions. In fact, we carry out the entire analysis via (6) on the level of the vacuum integrals. See Table 2 for statistics of the results. A detailed description will be given in section 5.
The convergence of such an integral implies that c j (P k (a; d)) = 0 for all singular coefficients (j < 0) in the ε-expansion (2). Furthermore, we know also that the c 0 of such finite p-integrals agree for all cuts of the same vacuum integral. These constraints form the starting point of the GaC algorithm.

Reductions
In order to apply the bootstrap procedure explained in section 2, we have to reduce the integrals involved to MI using IBP relations. These relations are generated by setting to zero a total derivative, where k j is constructed as a linear combination of the external and internal momenta. With enough of these relations one can utilize the Laporta algorithm [42] to row reduce the associated linear system. There exist several available computer codes to perform such reductions [43][44][45][46]. In this work we applied the current version of FIRE [43] (coupled with LiteRed [44]) with some features which are private at the moment but will soon be made public, as usually. After revealing primary master integrals, we used the recently developed code [47], which can be used to improve a given basis of master integrals by getting rid of so-called bad denominators from the coefficients of the IBP reduction 3 . By definition, a factor in the denominator is 'good' if it is either a linear function of d and independent of kinematic variables, or solely a polynomial of the kinematic variables and thus independent of d. Any other configuration is considered 'bad'. Applying this code to the most complicated sectors we obtained a new basis which is free of bad denominators in all but three of the families 4 , which are depicted in Figure 4. Since the resolution of this issue is similar for all three families, we shall focus here only on the leftmost diagram of Figure 4, which corresponds to family P 46 defined by (1) for the inverse propagators and numerators shown in Table 1.
The fact that some bad denominators survive might be explained by the existence of a hidden relation between master integrals. Unfortunately, the recipes based on symmetries presented in [50] did not help us find them. In order to reveal this hidden relation we analyzed IBP reductions for a set of integrals with all but one of the exponents a i set to 0 or 1, while the remaining exponent is equal to 2. By running FIRE on integrals with thirteen positive indices in two different ways, with the option no_presolve and without it, we then obtained two equivalent yet distinct reductions. By equating the corresponding results we then find the new relation, which has the form P 46 (1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 0 where dots stand for 37 master integrals with less than eleven positive indices. The complete relation can be found in the ancillary file ExtraRelation.m.
After taking into account this additional relation and using the option rules we observe that all the bad denominators disappear and the IBP reduction becomes faster and requires less RAM. The corresponding relations in the other two families of Figure 4 can be obtained from this one, since the relevant subsector is shared by all three families 5 .
By construction, this additional relation follows from IBP relations and symmetries which come into the game with LiteRed [44] so that there is nothing mysterious in it. However, it looks reasonable to try to describe this and similar relations in a more natural language. Figure 5: The diagrams on the left and center are rejected since they would produce p-integrals with tadpoles after the cutting procedure. The condition for removal is the existence of a zero-momentum edge (manifest in the center diagram but revealed on the left only after blow-up of the higher-valence vertex). Meanwhile the rightmost diagram is rejected because it produces integrals with a double propagator.

Details of our calculation
Using the computer program SAGE [51], we found 102 six-loop vacuum diagrams with at least 12 propagators, 99 of which are connected. Removing graphs with zero-momentum edges or double propagators (see Figure 5 for more details), we are left with 85 diagrams which are viable for cutting, some of which are depicted in Figure 6.
In order to have zero superficial degree of divergence ω 0 , defined in (5), the diagrams with (12 + x) propagators must be accompanied by x numerators. The complete basis of scalar products has 21 elements at this loop order, therefore the numerator can be any polynomial of degree x built from the (9 − x) irreducible scalar products (ISPs) available. In practice we avoided the complex search for all generic convergent numerators, and instead we only tested convergence for each monomial. We have many available integrals nevertheless, since each 4 There could be other simpler sectors for which this happens but this simplification was not necessary to obtain the reduction. 5 Within FIRE, one can use the function FindRules for obtaining the mapping.
diagram with (12 + x) propagators produces The attentive reader will wonder how it is possible to fix the expansions of master pintegrals with 14 denominators if we cut only vacuum diagrams with 14 denominators or less. The trick resides in the splitting of higher-valence vertices, which effectively increases the number of propagators by one. Note also that we checked convergence for numerators which are written as products of ISPs, hence the number of finite integrals in Table 2 depend on our choice of ISPs. We considered two options, a difference basis given by ( i − j ) 2 , and a dot basis formed by i · j (for example see D 15 , . . . , D 20 in Table 1), with 1 , . . . , 5 the loop momenta. The analysis at lower loops showed that the latter produced more convergent integrals, and therefore more equations, and so we considered the dot basis for the bootstrap at five loops.
We then systematically cut each of the 469 finite vacuum diagrams in all possible ways, either through an existing propagator, or by splitting a higher-valence vertex, as explained in section 2. In the end we obtained 7647 equations relating different cuttings of the vacuum diagrams. Each cut produces a single finite p-integral with a given numerator which is inherited from that of the vacuum diagram, but one needs to embed the integral into one of the 64 five-loop maximal sectors (diagrams with cubic vertices). Therefore, one must decompose the numerator into the basis of the chosen maximal sector, which effectively produces a linear combination of many divergent p-integrals that need to be reduced.
Once we run the IBP reductions for each of the 64 families, we obtain a set of master integrals for each of them. However, some of the integrals are common to several families. For example, the watermelon diagram inevitably appears in the reductions of all 64 maximal sectors. While it is easy to find relations between subsectors of different families, one must be careful when a given sector has more than one master integral. The reduction routine outputs in that case master integrals with a double denominator, but the precise prescription for its positioning depends on the family upon consideration. In order to avoid redundancies one must therefore set a convention for the sector that provides the preferred integral, and then map to it all variants arising from other families with a different double-denominator position.
After all such considerations, we find that at five loops there are 281 master integrals from 245 distinct sectors. There are 30 sectors with 2 master integrals, and also 3 sectors containing 3 master integrals each. Only 16 of the integral families P k from (1) contain master integrals in the top sector (that means a i > 0 for all propagators 1 ≤ i ≤ 14), but 12 of these families do contain more than one master integral in the top sector.
In principle it would have been possible to miss some master integrals, since the maximal sectors of a few integral families are not probed by any of the integrals produced through the cutting procedure. However, we have used Azurite [52] to confirm that there are no extra master integrals in any of those sectors. Furthermore, and since we applied the glue-and-cut algorithm to integrals with at most 2 numerators, we have also used Azurite to reduce all maximal sectors with degree 3 and verified the absence of new master integrals.
It is clear that the linear system of equations will impose a myriad of relations between the expansions of the different master integrals. However, the key point is whether we can fix all expansions in terms of coefficients that we can easily evaluate. At four loops this turned out to be possible, with all master integrals given in terms of the expansions of recursively one-loop integrals. This class of diagrams can be evaluated exactly by a recursive application of the one-loop bubble integration with G defined as Note that we used the freedom in the definition of dimensional regularization to normalize integrations according to the G-scheme [35]. Using the functional equation Γ(a + 1) = aΓ(a), the ε-expansion of G(a, b)/G(1, 1) consists, at any order, of rational linear combinations of Riemann zeta values: For example, near a = 1, the expansion can be obtained from Meanwhile, at five loops there are 21 recursively one-loop master integrals appearing in the IBP reductions of p-integrals (see Figure 7 for some examples). Equation (11) can then be used to obtain their expansions explicitly, so that we provide some inputs to the linear system of equations. In our case, however, these integrals cannot provide sufficient information to fix all expansions up to transcendental weight 9. Namely, at five loops there is for the first time a contribution from a multiple zeta value (MZV) of weight 8. MZVs are defined as ζ(n 1 , . . . , n r ) = where the indices n 1 , . . . , n r−1 ≥ 1 and n r ≥ 2 denote integers, and their sum n 1 + · · · + n r is referred to as the transcendental weight. Up to weight 7, all such MZV can be written as rational linear combinations of products of Riemann zeta values (r=1). But at weight 8, such relations still leave a 1-dimensional quotient space, generated e.g. by ζ (3,5), that conjecturally cannot be reduced in this way, and can hence never be obtained from the expansion of (12), since the latter only contains Riemann zeta values from (13). As ζ (3,5) has long been known to appear at five-loop order in massless propagators [53,54], it is thus clear that we will need further inputs, beyond recursively one-loop integrals.
At this point the linear system of equations obtained from the cutting procedure has been written in terms of the master integrals. We input the expansions of all recursively  one-loop integrals, which we can evaluate explicitly according to equation (11). However, before we consider the relations introduced by different cuttings, we can already extract some constraints from an analysis of the IBP identities themselves. While five-loop p-integrals diverge at most as ε −5 , one finds that some of the coefficients r i in the reductions (4) also contain poles in ε. Requiring the cancellation of those (spurious) poles of order ε −6 and worse, fixes already 563 coefficients in the expansions of the master integrals. Then we impose that all 8116 p-integrals obtained from distinct cuttings of the finite vacuum diagrams are finite, which further fixes 1300 coefficients. Finally, looking at the 7647 equations which relate the finite orders of those p-integrals, we are able to fix another 259 coefficients.
In total we were able to fix 2122 coefficients using the expansions of the 21 recursively oneloop integrals, which themselves contain only 264 coefficients. In the end, the expansions of all master integrals up to transcendental weight 9 depend only on 2 undetermined coefficients. As explained above, they encode the dependence on the expected MZV of weight 8. Fortunately, these undetermined coefficients are present in the expansions of many integrals, including 21 which are products of lower-loop masters. All lower-loop master integrals were computed up to transcendental weight twelve 6 , and so we can easily determine the expansions of all five-loop master integrals which are a product of lower-loop integrals.
In fact, the last 2 undetermined coefficients can be extracted from the expansion of the p-integral depicted in Figure 8. Using (11), the missing information is thus contained in known higher order corrections to the two-loop master integral [54,55].

Results
In this work we have demonstrated that the glue-and-cut algorithm still provides an extremely powerful way of determining master p-integrals at the five-loop order. Following the strategy developed at lower loops, we input the expansions of 21 recursively one-loop integrals, which can be evaluated easily to arbitrary orders in ε. The system of equations obtained fixes the expansions of the remaining 260 master-integrals, up to transcendental weight 9, in terms of a single master integral, which can be chosen as the product of planar lower-loop diagrams.
As we have mentioned in section 5, with our choice of master integrals, 7 some coefficients r i (P, 4 − 2ε) of the IBP reductions (4) occasionally have poles at ε = 0. These poles are called spurious poles in [22]. Such poles do not necessarily cause a problem, since the GaC constraints determine each master integral M i all the way up to the order q i where terms of transcendental weight 9 first appear. However, and unlike the four-loop case, we observe instances where the maximal power of a spurious pole is greater than the order q i to which the master integral is fixed. While these higher orders are not fully determined, the GaC constraints are still sufficient because we find that all remaining degrees of freedom arising from higher order coefficients c >q i (M i ) effectively cancel out from the coefficients c ≤0 (P ) of any p-integral P .
Furthermore, and since we are free to choose any basis of master integrals, we followed the procedure indicated in [56] and constructed a new basis such that the IBP reduction coefficients never exhibit poles in ε with degree higher than q i . This new basis is obtained by selecting a master integral M i for which the power of the spurious poles surpasses q i , and replacing it with one of the p-integrals P whose coefficient r i (P, d) realizes the maximal pole. This procedure is iterated, until all such spurious poles are resolved. 8 The basis thus constructed, and its relation to our original basis, is provided in the supplementary files.
We checked that the expansions of the planar integrals precisely match the results obtained in the context of position-space integrals [27]. Similarly, the expansions of all non-planar product-type integrals agree with those obtained from lower-loop master integrals. Finally, a few orders in the expansions of some genuine five-loop non-planar integrals which were known from the Fourier transform of position-space results [27] and a φ 4 computation [9], were all equally confirmed here. For the finite p-integrals obtained as cuts of φ 4 four-point graphs, our results agree with the periods determined in [57].
We wish also to emphasize that the equations obtained with this method do not discriminate between planar and non-planar graphs, and they often involve p-integrals from both categories. In that way, it seems highly unlikely that errors could arise only for the non-planar graphs for which we provide novel results.
Interestingly, we find that the ε-expansions of all master integrals up to transcendental weight 9, and moreover all coefficients c n (P ) with n ≤ 0 of all p-integrals we considered, can be written as polynomials (with rational coefficients) in the following five quantities: The same observation and ε-dependent transformation was found previously in the context of position space integrals [27], up to a redefinition ofζ (3,5) following [25] that is purely a matter of convenience. 9 This implies that, up to the ε-order where the master integrals reach transcendental weight 9, the contributions involving even zeta values are not independent, but in fact determined completely by the coefficients (in lower ε-orders) of polynomials in The first such connection, tying ζ(4) = π 4 90 to ζ(3) at three loops, was explored in [58]. At the four-loop level, the correlation of ζ(6) = π 6 945 to ζ(3) and ζ(5) was settled in [22], and extensions of such relations and ε-dependent transformations like (15) up to eight loops (transcendental weight 13) were proposed in [25]. Our work can be viewed, alongside [27], as a further indication towards the validity of this picture at five loops (transcendental weight 9), where multiple zeta values enter for the first time, asζ (3,5). These findings suggest that the rich number-theoretic structure [59] of convergent (ε = 0) p-integrals persists also in some form on the level of ε-expansions, which is a domain whose mathematical foundations are still under construction [60]. A systematic understanding of these ε-dependent transformations of transcendentals has so far been obtained only in the case of Riemann zeta values [26].
These constraints on the ε-expansions of p-integrals have striking implications for renormalization group functions, in form of the no-π theorem [61], which predicts the π-dependent terms (involving even zeta values) in terms of the π-free terms.
Given its success so far, it seems reasonable to expect that the glue-and-cut method should in principle be applicable also at six loops. However, it requires the ability to reduce p-integrals to a basis of master integrals, and at six loops that task does not seem feasible with any of the currently available tools for IBP reduction.
Finally, the p-integrals presented in this paper have different applications, beside the aforementioned computation of β-functions or anomalous dimensions. For example, they can also appear as boundary conditions in the context of differential equations when computing higher-point integrals [62][63][64].