<figure><img src="image_1.png" /></figure>

HAL
open science

Glue-and-cut at five loops

Alessandro Georgoudis, Vasco Goncalves, Erik Panzer, Raul Pereira,

Alexander V. Smirnov, Vladimir A. Smirnov

► To cite this version:

Alessandro Georgoudis, Vasco Goncalves, Erik Panzer, Raul Pereira, Alexander V. Smirnov, et al.. Glue-and-cut at five loops. Journal of High Energy Physics, 2021, 2021 (09), pp.098. 10.1007/JHEP09(2021)098. hal-03217628

HAL Id: hal-03217628

https://hal.science/hal-03217628v1

Submitted on 22 Aug 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Glue-and-cut at five loops

Alessandro Georgoudis,$^{a,b,d}$ Vasco Goncalves,$^{c,d}$ Erik Panzer,$^e$ Raul Pereira,$^f$ Alexander V. Smirnov$^{g,h}$ and Vladimir A. Smirnov$^{h,i}$

$^a$Laboratoire de physique de l'Ecole normale supérieure, ENS, Université PSL, CNRS, Sorbonne Université, Université Paris-Diderot, Sorbonne Paris Cité, 24 rue Lhomond, 75005 Paris, France

$^b$Department of Physics and Astronomy, Uppsala University, Box 516, SE-751 20 Uppsala, Sweden

$^c$Centro de Física do Porto e Departamento de Física e Astronomia, Faculdade de Ciências da Universidade do Porto, Rua do Campo Alegre 687, 4169-007 Porto, Portugal

$^d$ICTP South American Institute for Fundamental Research, Instituto de Física Teórica, Universidade Estadual Paulista, Rua Dr. Bento T. Ferraz 271, 01140-70, Sao Paulo, Brazil

$^e$Mathematical Institute, University of Oxford, OX2 6GG, Oxford, U.K.

$^f$School of Mathematics and Hamilton Mathematics Institute, Trinity College Dublin, Dublin, Ireland

$^g$Research Computing Center, Moscow State University, 119992 Moscow, Russia

$^h$Moscow Center for Fundamental and Applied Mathematics, 119991 Moscow, Russia

$^i$Skobeltsyn Institute of Nuclear Physics of Moscow State University, 119992 Moscow, Russia
E-mail: alessandro.georgoudis@phys.ens.fr, vasco.dfg@gmail.com, erik.panzer@maths.ox.ac.uk, raul@maths.tcd.ie, asmirnov80@gmail.com, smirnov@theory.sinp.msu.ru

**ABSTRACT:** We compute $\varepsilon$-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.

**KEYWORDS:** Renormalization Regularization and Renormalons, Scattering Amplitudes, Perturbative QCD

ARXIV ePRINT: 2104.08272
# Contents

<table><tr><td>1 Introduction</td><td>1</td></tr><tr><td>2 Review of the algorithm</td><td>3</td></tr><tr><td>3 Convergence</td><td>6</td></tr><tr><td>4 Reductions</td><td>8</td></tr><tr><td>5 Details of our calculation</td><td>9</td></tr><tr><td>6 Results</td><td>13</td></tr></table>

## 1 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–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\epsilon$ dimensions [6, 7], and the master integrals have to be expanded to a sufficiently high order in the $\epsilon$-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 $\phi^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–14] and position-space techniques [15]. In this comparison, the coefficients of the $\epsilon$-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–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 $\epsilon$-expansions of p-integrals with only $L-1$ loops.
<figure><img src="image_3.png" /><figcaption>**Figure 1.** These are three examples of maximal sectors for five-loop p-integrals. The diagrams have only cubic vertices and 14 internal propagators. The rightmost sector is an example of a non-planar integral.</figcaption></figure>

As a step towards an extension of the currently very good understanding of the master p-integrals [21–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 $\varepsilon$-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 $\pi$-dependent contributions can be predicted by the $\pi$-free terms. We supply our results for the $\varepsilon$-expansions of our master integrals and other results in the supplementary material attached to 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,<fnref n="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.

A range of methods to compute single-scale massless Feynman integrals has been developed over time [28]. Relatively recent advances and applications include in particular dimensional recurrence relations [21, 29], integration over Schwinger parameters [30–33] and graphical functions [34].

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 [2] 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 $\varepsilon$-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

<footnote n="1">¹A *sector* encompasses all integrals in a family $P_k(a; d)$ where the set $\{j: a_j = 0\}$ is fixed. Since setting $a_j = 0$ for an inverse propagator ($j < 15$) corresponds to contraction of the associated edge in the underlying graph of the family, different families can share common subsectors (different graphs may become isomorphic after several contractions).</footnote>
(see figure 7), which are easily evaluated to all orders in $\epsilon$, 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.

## 2 Review of the algorithm

All five-loop massless propagator integrals can be represented by families of the form

$$P_k(a_1, \dots, a_{20}; d) := \int \frac{d^d \ell_1}{\pi^{d/2}} \dots \int \frac{d^d \ell_5}{\pi^{d/2}} \frac{1}{D_1^{a_1} \dots D_{20}^{a_{20}}} \quad (2.1)$$

where $\ell_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 \le k \le 64$, each corresponding to a different cubic graph with 14 internal lines, and $D_1, \dots, D_{14}$ encode precisely the inverse propagators associated to these lines. The corresponding indices $a_1, \dots, a_{14}$ will always be positive or zero. The remaining six quadratic forms $D_{15}, \dots, D_{20}$ will only appear with non-positive indices $a_{15}, \dots, a_{20} \le 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\epsilon$, so every propagator admits a Laurent expansion

$$P(d) = (p^2)^{-\omega_{\epsilon}(P)} \sum_{n \in \mathbb{Z}} c_n(P) \epsilon^n \quad (2.2)$$

whose coefficients $c_n \in \mathbb{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

$$\omega_{\epsilon} = \sum_{j=1}^{20} a_j - 5(2 - \epsilon). \quad (2.3)$$
<table><caption>**Table 1.** Inverse propagators $D_1, \dots, D_{14}$ and numerators $D_{15}, \dots, D_{20}$ for integral family $P_{46}$.</caption><thead><tr><th>D<sub>1,...,4</sub></th><th>D<sub>5,...,8</sub></th><th>D<sub>9,...,12</sub></th><th>D<sub>13,...,16</sub></th><th>D<sub>17,...,20</sub></th></tr></thead><tbody><tr><td>ℓ<sub>1</sub><sup>2</sup></td><td>ℓ<sub>3</sub><sup>2</sup></td><td>(ℓ<sub>1</sub> + ℓ<sub>3</sub> + ℓ<sub>4</sub> + ℓ<sub>5</sub> - p)<sup>2</sup></td><td>(ℓ<sub>1</sub> + ℓ<sub>4</sub>)<sup>2</sup></td><td>ℓ<sub>2</sub> · ℓ<sub>4</sub></td></tr><tr><td>(ℓ<sub>1</sub> - p)<sup>2</sup></td><td>ℓ<sub>4</sub><sup>2</sup></td><td>(ℓ<sub>1</sub> + ℓ<sub>3</sub> + ℓ<sub>4</sub> - p)<sup>2</sup></td><td>(ℓ<sub>2</sub> - ℓ<sub>3</sub> + p)<sup>2</sup></td><td>ℓ<sub>2</sub> · ℓ<sub>5</sub></td></tr><tr><td>ℓ<sub>2</sub><sup>2</sup></td><td>(ℓ<sub>3</sub> + ℓ<sub>4</sub>)<sup>2</sup></td><td>(ℓ<sub>2</sub> - ℓ<sub>3</sub> - ℓ<sub>5</sub> + p)<sup>2</sup></td><td>ℓ<sub>1</sub> · ℓ<sub>3</sub></td><td>ℓ<sub>3</sub> · ℓ<sub>5</sub></td></tr><tr><td>(ℓ<sub>2</sub> + p)<sup>2</sup></td><td>ℓ<sub>5</sub><sup>2</sup></td><td>(ℓ<sub>1</sub> - ℓ<sub>2</sub> + ℓ<sub>3</sub> + ℓ<sub>4</sub> + ℓ<sub>5</sub> - p)<sup>2</sup></td><td>ℓ<sub>2</sub> · ℓ<sub>3</sub></td><td>ℓ<sub>4</sub> · ℓ<sub>5</sub></td></tr></tbody></table>

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 \in \mathbb{Z}^{20}$ and reduce each of those p-integrals $P$ to (any preferred choice of) master integrals $M_i$:

$$P(d) = r_1(P, d)M_1(d) + \dots + r_N(P, d)M_N(d) \quad (2.4)$$

2. Expand both sides of (2.4) in $\epsilon$ 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:

**Theorem 1.** If the p-integral $P(\epsilon)$ is finite in $d = 4$ dimensions ($\epsilon = 0$), and if $V$ is dimensionless, then any other cut $P'(\epsilon)$ of $V$ is also finite in $d = 4$, and $P(0) = P'(0)$.

Dimensionless here means that the additional edge in $V$ is assigned the index $a_0 = d/2 - \omega_{\epsilon}(P)$, such that we have $\omega_{\epsilon}(V) = 0$ where

$$\omega_{\epsilon}(V) := \sum_{j=0}^{20} a_j - 6(2 - \epsilon).$$

In fact, as long as this condition is fulfilled, the symmetry holds for arbitrary values of the indices $a$ and for arbitrary $\epsilon$: $P(\epsilon) = P'(\epsilon)$ [35]. Indeed, any p-integral $P$ obtained from
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Starting with the five-loop p-integral on the left and glueing the external lines (dashed), we obtain the six-loop vacuum diagram in the center. This diagram has 12 propagators and no numerator, leading to a vanishing superficial degree of divergence. Since there are no subdivergences, we can cut a different propagator (dotted) to produce a distinct five-loop p-integral with the same value in four dimensions.</figcaption></figure>

cutting $V$, is equal to the residue of the vacuum integral $V$ at $\omega_{\epsilon}(V) = 0$.<fnref n="2" /> However, note that $a_0$ will depend on $\epsilon$, so in specializing to $\epsilon = 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 p-integrals. 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

$$ \int \frac{d^d \ell_1 \cdots d^d \ell_{L+1}}{D_0^{a_0} \cdots D_n^{a_n}}, $$

with $D_i$ a linear combination of scalar products $\ell_i \cdot \ell_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

$$ \omega_0 = \sum_{k=0}^{n} a_k - 2(L+1). \qquad (2.5) $$

According to theorem 1, we consider only diagrams with $\omega_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, \dots, 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} + \dots + a_{20} = 3$ numerators. For graphs with less denominators ($a_i = 0$ for one or more $i \le 14$), we pick accordingly less numerators. Apart from thus assuring $\omega_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,

<footnote n="2">²To make sense of a vacuum Feynman integral, it is in principle necessary to introduce some kinematic dependence, e.g. by adding some external legs, or assigning masses to some lines. However, the residue of $V$ at $\omega = 0$ does not depend on any of these choices.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** We split the quartic vertex into a pair of cubic vertices connected by a new propagator (dashed). Cutting this propagator yields yet another finite p-integral.</figcaption></figure>

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) \ge 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 (2.4) to master integrals. At five loops this is quite an arduous task, further details are given in section 4. Once the reduction is performed and all master integrals have been identified, one just has to perform the symbolic $\epsilon$ expansion (2.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 [16].

### 3 Convergence

Consider a scalar p-integral, that is, some $P_k(a; d)$ as per (2.1) where $a_{15} = \dots = a_{20} = 0$ (no numerators). Such a p-integral is convergent if, and only if, we have $\omega_{\epsilon} > 0$ and the conditions

$$ \sum_{e \in g} a_e - \frac{d}{2} \cdot L_g > \begin{cases} \omega_{\epsilon} & \text{if } g \text{ connects the external legs and} \\ 0 & \text{otherwise,} \end{cases} $$

are fulfilled for all non-empty proper subsets $g \subsetneq \{1, \dots, 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 $\omega_{\epsilon}(V) = 0$. Then the conditions above can equivalently be stated more symmetrically as

$$ \sum_{e \in g} a_e > \frac{d}{2} \cdot L_g \qquad (3.1) $$

for all non-empty subgraphs $g \subsetneq \{0, \dots, 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
<table><caption>**Table 2.** Counts of non-isomorphic vacuum graphs (see figure 6) by number of edges (propagators). The third column gives the number (5.1) of logarithmically divergent numerators.</caption><thead><tr><th>propagators</th><th>ISPs</th><th>numerators</th><th>graphs</th><th>finite integrals</th></tr></thead><tbody><tr><td>15</td><td>6</td><td>56</td><td>14</td><td>347</td></tr><tr><td>14</td><td>7</td><td>28</td><td>22</td><td>248</td></tr><tr><td>13</td><td>8</td><td>8</td><td>32</td><td>206</td></tr><tr><td>12</td><td>9</td><td>1</td><td>17</td><td>15</td></tr></tbody></table>

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 $\sum_{e \in g} a_e > d/2 \cdot 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 (2.1). Then, the quadratic form

$$x_1 D_1 + \cdots + x_{20} D_{20} = \ell^\top A \ell + 2 B^\top \ell + C$$

in the loop momenta $\ell = (\ell_1, \dots, \ell_5)$ determines a $5 \times 5$ matrix $A$, a $5$-vector $B$ and a scalar $C$ (all of these depend on $x$). Define the polynomials $\mathcal{U} := \det A$ and $\mathcal{F} := \mathcal{U}(B^\top A^{-1} B - C)$. The parametric representation [37, 41] of the integral family (2.1) then takes the form

$$P_k(a; d) = \Gamma(\omega_{\varepsilon}) \left( \prod_{j: a_j > 0} \int_0^{\infty} \frac{x_j^{a_j-1} dx_j}{\Gamma(a_j)} \right) \left( \prod_{j: a_j \le 0} \left[ -\frac{\partial}{\partial x_j} \right]_{x_j=0}^{-a_j} \right) \frac{\delta(1-x_1)}{\mathcal{U}^{d/2-\omega_{\varepsilon}} \mathcal{F}^{\omega_{\varepsilon}}}. \quad (3.2)$$

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) = \sum_n b_n(\varepsilon) P_k(a^{(n)}; d')$ as a finite linear combination of p-integrals in some raised dimension $d' - d \in 2\mathbb{Z}_{\ge 0}$ that are scalar ($a_j^{(n)} = 0$ whenever $a_j < 0$), with raised indices $a_j^{(n)} \in a_j + \mathbb{Z}_{\ge 0}$ on the denominators ($a_j > 0$).

The convergence of the integral (3.2) 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 (3.1) 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 $\varepsilon$-expansion (2.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.
<figure><img src="image_6.png" /><figcaption>Figure 4. The integral families for which bad denominators survive the change of basis. The diagram on the left shows family $P_{46}$.</figcaption></figure>

## 4 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,

$$0 = \int \frac{d^d \ell_1}{\pi^{d/2}} \cdots \frac{d^d \ell_L}{\pi^{d/2}} \sum_{j=1}^{L} \frac{\partial}{\partial \ell_j^\mu} \frac{k_j^\mu}{D_1^{a_1} \cdots D_n^{a_n}}, \quad (4.1)$$

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–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.<fnref n="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,<fnref n="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 (2.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

<footnote n="3">³See also alternative competitive code presented in [48, 49].</footnote>

<footnote n="4">⁴There could be other simpler sectors for which this happens but this simplification was not necessary to obtain the reduction.</footnote>
<figure><img src="image_7.png" /><figcaption>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.</figcaption></figure>

results we then find the new relation, which has the form

$$ P_{46}(1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 0, 0, 0, 0, 0; d) = \frac{1}{3d-11} \left[ \begin{aligned} & (8d-28)P_{46}(0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 0, 0, 0, 0, 0; d) \\ & -(5d-25)P_{46}(1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 0, 0, 0; d) \\ & +(d-5)P_{46}(0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 0, 0, 0; d) \end{aligned} \right] + \dots \quad (4.2) $$

where dots stand for 37 master integrals with less than eleven positive indices. The complete relation can be found in the supplementary 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.<fnref n="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.

## 5 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 $\omega_0$, defined in (2.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

<footnote n="5">⁵Within FIRE, one can use the function FindRules for obtaining the mapping.</footnote>
<figure><img src="image_8.png" /><figcaption>Figure 6. Examples of vacuum diagrams at six loops with 12, 13, 14 and 15 denominators (reading from left to right and top to bottom). The larger the number of propagators, the larger the number of numerators we must include to ensure vanishing $\omega_0$.</figcaption></figure>

nevertheless, since each diagram with $(12+x)$ propagators produces

$$ (9-x)_x/x! \qquad (5.1) $$

distinct numerators, with $(a)_x$ the Pochhammer symbol. However, since the IBP reductions are very demanding at this loop order, we neglected all vacuum diagrams with 15 denominators and degree 3 numerators, and checked only the convergence properties of the remaining $889 = 28 \cdot 22 + 8 \cdot 32 + 1 \cdot 17$ integrals in table 2, finding that $469 = 248 + 206 + 15$ of them are free of divergences.

The attentive reader will wonder how it is possible to fix the expansions of master p-integrals 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 $(\ell_i - \ell_j)^2$, and a dot basis formed by $\ell_i \cdot \ell_j$ (for example see $D_{15}, \dots, D_{20}$ in table 1), with $\ell_1, \dots, \ell_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 (2.1) contain master integrals in the top sector (that means $a_i > 0$ for all propagators $1 \le i \le 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

$$ \frac{1}{\varepsilon G(1, 1)} \int \frac{d^d \ell}{\pi^{d/2}} \frac{1}{\ell^{2a} (p - \ell)^{2b}} = (p^2)^{d/2-a-b} \frac{G(a, b)}{\varepsilon G(1, 1)}, \qquad (5.2) $$

with G defined as

$$ G(a, b) = \frac{\Gamma(a + b - d/2)\Gamma(d/2 - a)\Gamma(d/2 - b)}{\Gamma(a)\Gamma(b)\Gamma(d - a - b)}. \qquad (5.3) $$

Note that we used the freedom in the definition of dimensional regularization to normalize integrations according to the G-scheme [53]. Using the functional equation $\Gamma(a + 1) = a\Gamma(a)$, the $\varepsilon$-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

$$ \frac{\Gamma(2 - \varepsilon - a)}{\Gamma(a)} / \frac{\Gamma(1 - \varepsilon)}{\Gamma(1)} = \exp \left[ \sum_{k=2}^{\infty} \frac{\zeta(k)}{k} \left( (\varepsilon + a - 1)^n - (1 - a)^n - \varepsilon^n \right) \right]. \qquad (5.4) $$
<figure><img src="image_9.png" /><figcaption>Figure 7. Three examples of recursively one-loop master integrals. Each loop variable can be integrated one by one, and in the end we obtain an exact expression for each of these integrals as a function of the spacetime dimension d.</figcaption></figure>

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 (5.2) 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

$$ \zeta(n_1, \dots, n_r) = \sum_{1 \le k_1 < \dots < k_r} \frac{1}{k_1^{n_1} \dots k_r^{n_r}} \in \mathbb{R}_{>0}, \qquad (5.5) $$

where the indices $n_1, \dots, n_{r-1} \ge 1$ and $n_r \ge 2$ denote integers, and their sum $n_1 + \dots + 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 $\zeta(3, 5)$, that conjecturally cannot be reduced in this way, and can hence never be obtained from the expansion of (5.3), since the latter only contains Riemann zeta values from (5.4). As $\zeta(3, 5)$ has long been known to appear at five-loop order in massless propagators [54, 55], 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 (5.2). 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 $\epsilon^{-5}$, one finds that some of the coefficients $r_i$ in the reductions (2.4) also contain poles in $\epsilon$. Requiring the cancellation of those (spurious) poles of order $\epsilon^{-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 one-loop 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 unde-
<figure><img src="image_10.png" /><figcaption>**Figure 8.** This five-loop master integral is a product of two-and three-loop master integrals. The knowledge of its expansion, together with the recursively one-loop integrals, is sufficient to determine the expansions of all remaining master integrals.</figcaption></figure>

termined 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,<fnref n="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 (5.2), the missing information is thus contained in known higher order corrections to the two-loop master integral [55, 56].

## 6 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 $\epsilon$. As predicted in [22] this type of integrals are not enough to fix completely the solution. In fact the system of equations obtained then fixes the expansions of the remaining 260 master-integrals, up to transcendental weight 9, in terms of a single master integral. This master integral can then be chosen as the product of planar lower-loop diagrams, one of which can be written as the difference of two $_3F_2$ hypergeometric functions [55, 57–59] and fixes the multiple zeta value contribution.

As we have mentioned in section 5, with our choice of master integrals,<fnref n="7" /> some coefficients $r_i(P, 4-2\epsilon)$ of the IBP reductions (2.4) occasionally have poles at $\epsilon = 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. A general master integral will have the following expansion:

$$M_i(\epsilon) = \sum_{n=-5}^{q_i} c_n(M_i)\epsilon^n + \sum_{n>q_i} c_n(M_i)\epsilon^n,$$

<footnote n="6">⁶At four loops, this was done in [21].</footnote>

<footnote n="7">⁷Complete details are provided in the supplementary material.</footnote>
where the first terms with $n \le q_i$ are fully fixed via GaC equations in terms of zeta values $\{\zeta(k), \zeta(k, j)\}$ of transcendental weight $\le 9$. The higher orders $n > q_i$ are not fully determined, but nevertheless the GaC constraints provide several relations among these higher order coefficients between different master integrals. We have incorporated this information in our results by solving these equations, so that some of the higher order coefficients $c_{>q_i}(M_i)$ are expressed as linear combinations of coefficients $c_{>q_j}(M_j)$ of other master integrals. We find that these GaC constraints on the higher orders are sufficient to determine $c_{\le 0}(P)$ for any P-integral: after IBP reduction of P to our basis of master integrals, all remaining undetermined degrees of freedom arising from higher order coefficients $c_{>q_j}(M_j)$ cancel out from the coefficients $c_{\le 0}(P)$.

Given that the undetermined coefficients cancel out in the IBP reductions of p-integrals and that we are free to choose any basis of master integrals, we followed the procedure indicated in [60] and constructed a new basis such that the IBP reduction coefficients never exhibit poles in $\epsilon$ 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.<fnref n="8" /> The basis thus constructed, and its relation to our original basis, is provided in the supplementary material.

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 $\phi^4$ computation [9], were all equally confirmed here. For the finite p-integrals obtained as cuts of $\phi^4$ four-point graphs, our results agree with the periods determined in [61].

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 $\epsilon$-expansions of all master integrals up to transcendental weight 9, and moreover all coefficients $c_n(P)$ with $n \le 0$ of all p-integrals we considered, can be written as polynomials (with rational coefficients) in the following five quantities:

$$ 
\begin{aligned}
\hat{\zeta}(3) &= \zeta(3) + \frac{3\epsilon}{2}\zeta(4) - \frac{5\epsilon^3}{2}\zeta(6) + \frac{21\epsilon^5}{2}\zeta(8), & \hat{\zeta}(5) &= \zeta(5) + \frac{5\epsilon}{2}\zeta(6) - \frac{35\epsilon^3}{4}\zeta(8), \\
\hat{\zeta}(7) &= \zeta(7) + \frac{7\epsilon}{2}\zeta(8), & \hat{\zeta}(3,5) &= \zeta(3,5) - \frac{29}{12}\zeta(8) - \frac{15\epsilon}{2}\zeta(4)\zeta(5), \\
\text{and } \zeta(9). & & &
\end{aligned} 
\quad (6.1) 
$$

The same observation and $\epsilon$-dependent transformation was found previously in the context of position space integrals [27], up to a redefinition of $\hat{\zeta}(3,5)$ following [25] that is

<footnote n="8">⁸It is interesting to notice that the number of replaced master integrals is considerably less than the number of integrals with undetermined higher-order coefficients, which points to a common structure of the spurious divergent poles in a specific dimension.</footnote>
purely a matter of convenience.<fnref n="9" /> This implies that, up to the $\epsilon$-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 $\epsilon$-orders) of polynomials in

$$ \zeta(3), \ \zeta(5), \ \zeta(7), \ \zeta(3,5) - \frac{29}{12}\zeta(8), \ \text{and} \ \zeta(9). \qquad (6.2) $$

The first such connection, tying $\zeta(4) = \frac{\pi^4}{90}$ to $\zeta(3)$ at three loops, was explored in [62]. At the four-loop level, the correlation of $\zeta(6) = \frac{\pi^6}{945}$ to $\zeta(3)$ and $\zeta(5)$ was settled in [22], and extensions of such relations and $\epsilon$-dependent transformations like (6.1) 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 $\hat{\zeta}(3,5)$. These findings suggest that the rich number-theoretic structure [63] of convergent ($\epsilon=0$) p-integrals persists also in some form on the level of $\epsilon$-expansions, which is a domain whose mathematical foundations are still under construction [64]. A systematic understanding of these $\epsilon$-dependent transformations of transcendentals has so far been obtained only in the case of Riemann zeta values [26].

These constraints on the $\epsilon$-expansions of p-integrals have striking implications for renormalization group functions, in form of the no-$\pi$ theorem [65], which predicts the $\pi$-dependent terms (involving even zeta values) in terms of the $\pi$-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 $\beta$-functions or anomalous dimensions. For example, they can also appear as boundary conditions in the context of differential equations when computing higher-point integrals [66–68].

## Acknowledgments

The authors would like to thank Kostja Chetyrkin for useful discussions regarding the glue-and-cut method. AG also thanks Tiziano Peraro for providing the reduction of one of the families of p-integrals, using the method of [69]. Part of the computations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at Uppmax. V.G. is funded by FAPESP grant 2015/14796-7, CERN/FIS-NUC/0045/2015 and Simons Foundation grant 488637 (Simons Collaboration on the Non-perturbative bootstrap). The work of A.G. is supported by the Knut and Alice Wallenberg Foundation under grant 2015.0083 and by the French National Agency for Research grant ANR-17-CE31-0001-02. A.G. would like to thank FAPESP grant 2016/01343-7 for funding part of his

<footnote n="9">⁹In terms of $\hat{\zeta}(3,5)' = \zeta(2,6) - \zeta(5,3) - 3\epsilon\zeta(4)\zeta(5) + \frac{5\epsilon}{2}\zeta(3)\zeta(6)$ as defined in [27], the precise relation is $\hat{\zeta}(3,5)' = \frac{3}{5}\hat{\zeta}(3,5) + \hat{\zeta}(3)\hat{\zeta}(5)$.</footnote>
visit to ICTP-SAIFR in March 2019 where part of this work was done. R.P. was funded under SFI grant 15/CDA/3472. The work of A.S. and V.S. was supported by Russian Ministry of Science and Higher Education, agreement No. 075-15-2019-1621. E.P. is funded as a Royal Society University Research Fellow through grant URF\R1\201473. A.G. and V.S. are grateful to the organizers of the Paris Winter Workshop: The Infrared in QFT (March 2020) where we understood how the present project could be completed.

**Open Access.** This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.

## References

[1] F.V. Tkachov, *A Theorem on Analytical Calculability of Four Loop Renormalization Group Functions*, Phys. Lett. B 100 (1981) 65 [INSPIRE].

[2] K.G. Chetyrkin and F.V. Tkachov, *Integration by Parts: The Algorithm to Calculate $\beta$-functions in 4 Loops*, Nucl. Phys. B 192 (1981) 159 [INSPIRE].

[3] A.V. Smirnov and A.V. Petukhov, *The Number of Master Integrals is Finite*, Lett. Math. Phys. 97 (2011) 37 [arXiv:1004.4199] [INSPIRE].

[4] R.N. Lee and A.A. Pomeransky, *Critical points and number of master integrals*, JHEP 11 (2013) 165 [arXiv:1308.6676] [INSPIRE].

[5] T. Bitoun, C. Bogner, R.P. Klausen and E. Panzer, *Feynman integral relations from parametric annihilators*, Lett. Math. Phys. 109 (2019) 497 [arXiv:1712.09215] [INSPIRE].

[6] C.G. Bollini and J.J. Giambiagi, *Dimensional Renormalization: The Number of Dimensions as a Regularizing Parameter*, Nuovo Cim. B 12 (1972) 20 [INSPIRE].

[7] G. 't Hooft and M.J.G. Veltman, *Regularization and Renormalization of Gauge Fields*, Nucl. Phys. B 44 (1972) 189 [INSPIRE].

[8] H. Kleinert, J. Neu, V. Schulte-Frohlinde, K.G. Chetyrkin and S.A. Larin, *Five loop renormalization group functions of $O(n)$-symmetric $\phi^4$-theory and $\epsilon$-expansions of critical exponents up to $\epsilon^5$, Phys. Lett. B 272 (1991) 39* [Erratum ibid. 319 (1993) 545] [hep-th/9503230] [INSPIRE].

[9] M.V. Kompaniets and E. Panzer, *Minimally subtracted six loop renormalization of $O(n)$-symmetric $\phi^4$ theory and critical exponents*, Phys. Rev. D 96 (2017) 036016 [arXiv:1705.06483] [INSPIRE].

[10] P.A. Baikov, K.G. Chetyrkin and J.H. Kühn, *Five-Loop Running of the QCD coupling constant*, Phys. Rev. Lett. 118 (2017) 082002 [arXiv:1606.08659] [INSPIRE].

[11] F. Herzog, B. Ruijl, T. Ueda, J.A.M. Vermaseren and A. Vogt, *The five-loop $\beta$-function of Yang-Mills theory with fermions*, JHEP 02 (2017) 090 [arXiv:1701.01404] [INSPIRE].

[12] T. van Ritbergen, J.A.M. Vermaseren and S.A. Larin, *The Four loop $\beta$-function in quantum chromodynamics*, Phys. Lett. B 400 (1997) 379 [hep-ph/9701390] [INSPIRE].

[13] M. Czakon, *The Four-loop QCD $\beta$-function and anomalous dimensions*, Nucl. Phys. B 710 (2005) 485 [hep-ph/0411261] [INSPIRE].
[14] T. Luthe, A. Maier, P. Marquard and Y. Schröder, Complete renormalization of QCD at five loops, JHEP 03 (2017) 020 [arXiv:1701.07068] [INSPIRE].

[15] O. Schnetz, Numbers and Functions in Quantum Field Theory, Phys. Rev. D 97 (2018) 085018 [arXiv:1606.08598] [INSPIRE].

[16] K.G. Chetyrkin and V.A. Smirnov, $R^*$-operation corrected, Phys. Lett. B 144 (1984) 419 [INSPIRE].

[17] K.G. Chetyrkin and F.V. Tkachov, Infrared $R$-operation and ultraviolet counterterms in the MS-scheme, Phys. Lett. B 114 (1982) 340 [INSPIRE].

[18] K.G. Chetyrkin, Combinatorics of $R$, $R^{-1}$, and $R^*$-operations and asymptotic expansions of Feynman integrals in the limit of large momenta and masses, Tech. Rep. MPI-Ph/PTh 13/91, Max-Planck-Institut für Physik, Werner-Heisenberg-Institut, Munich, Germany (1991) [arXiv:1701.08627] [INSPIRE].

[19] F. Herzog and B. Ruijl, The $R^*$-operation for Feynman graphs with generic numerators, JHEP 05 (2017) 037 [arXiv:1703.03776] [INSPIRE].

[20] R. Beekveldt, M. Borinsky and F. Herzog, The Hopf algebra structure of the $R^*$-operation, JHEP 07 (2020) 061 [arXiv:2003.04301] [INSPIRE].

[21] R.N. Lee, A.V. Smirnov and V.A. Smirnov, Master Integrals for Four-Loop Massless Propagators up to Transcendentality Weight Twelve, Nucl. Phys. B 856 (2012) 95 [arXiv:1108.0732] [INSPIRE].

[22] P.A. Baikov and K.G. Chetyrkin, Four Loop Massless Propagators: An Algebraic Evaluation of All Master Integrals, Nucl. Phys. B 837 (2010) 186 [arXiv:1004.1153] [INSPIRE].

[23] A.V. Smirnov and M. Tentyukov, Four Loop Massless Propagators: a Numerical Evaluation of All Master Integrals, Nucl. Phys. B 837 (2010) 40 [arXiv:1004.1149] [INSPIRE].

[24] T. Ueda, B. Ruijl and J.A.M. Vermaseren, Forcer: a FORM program for 4-loop massless propagators, PoS LL2016 (2016) 070 [arXiv:1607.07318] [INSPIRE].

[25] P.A. Baikov and K.G. Chetyrkin, Transcendental structure of multiloop massless correlators and anomalous dimensions, JHEP 10 (2019) 190 [arXiv:1908.03012] [INSPIRE].

[26] A.V. Kotikov and S. Teber, Landau-Khalatnikov-Fradkin transformation and the mystery of even $\zeta$-values in Euclidean massless correlators, Phys. Rev. D 100 (2019) 105017 [arXiv:1906.10930] [INSPIRE].

[27] A. Georgoudis, V. Goncalves, E. Panzer and R. Pereira, Five-loop massless propagator integrals, arXiv:1802.00803 [INSPIRE].

[28] A.V. Kotikov and S. Teber, Multi-loop techniques for massless Feynman diagram calculations, Phys. Part. Nucl. 50 (2019) 1 [arXiv:1805.05109] [INSPIRE].

[29] R.N. Lee, Space-time dimensionality $D$ as complex variable: Calculating loop integrals using dimensional recurrence relation and analytical properties with respect to $D$, Nucl. Phys. B 830 (2010) 474 [arXiv:0911.0252] [INSPIRE].

[30] E. Panzer, On the analytic computation of massless propagators in dimensional regularization, Nucl. Phys. B 874 (2013) 567 [arXiv:1305.2161] [INSPIRE].

[31] F. Brown, The Massless higher-loop two-point function, Commun. Math. Phys. 287 (2009) 925 [arXiv:0804.1660] [INSPIRE].
[32] E. Panzer, *Algorithms for the symbolic integration of hyperlogarithms with applications to Feynman integrals*, *Comput. Phys. Commun.* **188** (2015) 148 [arXiv:1403.3385] [INSPIRE].

[33] A. von Manteuffel, E. Panzer and R.M. Schabinger, *A quasi-finite basis for multi-loop Feynman integrals*, *JHEP* **02** (2015) 120 [arXiv:1411.7392] [INSPIRE].

[34] O. Schnetz, *Graphical functions and single-valued multiple polylogarithms*, *Commun. Num. Theor. Phys.* **08** (2014) 589 [arXiv:1302.6445] [INSPIRE].

[35] S.G. Gorishnii and A.P. Isaev, *On an Approach to the Calculation of Multiloop Massless Feynman Integrals*, *Theor. Math. Phys.* **62** (1985) 232 [INSPIRE].

[36] S. Weinberg, *High-energy behavior in quantum field theory*, *Phys. Rev.* **118** (1960) 838 [INSPIRE].

[37] V.A. Smirnov, *Analytic tools for Feynman integrals*, Springer Tracts Mod. Phys. **250** (2012) 1 [INSPIRE].

[38] E.R. Speer, *Ultraviolet and infrared singularity structure of generic Feynman amplitudes*, *Ann. Inst. H. Poincaré Phys. Theor.* **23** (1975) 1 [INSPIRE].

[39] E. Panzer, *Hepp's bound for Feynman graphs and matroids*, arXiv:1908.09820 [INSPIRE].

[40] O.V. Tarasov, *Connection between Feynman integrals having different values of the space-time dimension*, *Phys. Rev. D* **54** (1996) 6479 [hep-th/9606018] [INSPIRE].

[41] N. Nakanishi, *Graph theory and Feynman integrals*, vol. 11 of Mathematics and its applications, Gordon and Breach, New York (1971).

[42] S. Laporta, *Calculation of master integrals by difference equations*, *Phys. Lett. B* **504** (2001) 188 [hep-ph/0102032] [INSPIRE].

[43] A.V. Smirnov and F.S. Chuharev, *FIRE6: Feynman Integral REDuction with Modular Arithmetic*, *Comput. Phys. Commun.* **247** (2020) 106877 [arXiv:1901.07808] [INSPIRE].

[44] R.N. Lee, *Presenting LiteRed: a tool for the Loop InTEgrals REDuction*, arXiv:1212.2685 [INSPIRE].

[45] A. von Manteuffel and C. Studerus, *Reduze 2 — Distributed Feynman Integral Reduction*, arXiv:1201.4330 [INSPIRE].

[46] J. Klappert, F. Lange, P. Maierhöfer and J. Usovitsch, *Integral reduction with Kira 2.0 and finite field methods*, *Comput. Phys. Commun.* **266** (2021) 108024 [arXiv:2008.06494] [INSPIRE].

[47] A.V. Smirnov and V.A. Smirnov, *How to choose master integrals*, Nucl. Phys. B **960** (2020) 115213 [arXiv:2002.08042] [INSPIRE].

[48] J. Usovitsch, *Factorization of denominators in integration-by-parts reductions*, arXiv:2002.08173 [INSPIRE].

[49] J. Boehm, M. Wittmann, Z. Wu, Y. Xu and Y. Zhang, *IBP reduction coefficients made simple*, *JHEP* **12** (2020) 054 [arXiv:2008.13194] [INSPIRE].

[50] A.V. Smirnov and V.A. Smirnov, *FIRE4, LiteRed and accompanying tools to solve integration by parts relations*, *Comput. Phys. Commun.* **184** (2013) 2820 [arXiv:1302.5885] [INSPIRE].

[51] The Sage Developers, *SageMath, the Sage Mathematics Software System (Version 8.3)*, (2018), http://www.sagemath.org.
[52] A. Georgoudis, K.J. Larsen and Y. Zhang, *Azurite: An algebraic geometry based package for finding bases of loop integrals*, Comput. Phys. Commun. **221** (2017) 203 [arXiv:1612.04252] [INSPIRE].

[53] K.G. Chetyrkin, A.L. Kataev and F.V. Tkachov, *New Approach to Evaluation of Multiloop Feynman Integrals: The Gegenbauer Polynomial x Space Technique*, Nucl. Phys. B **174** (1980) 345 [INSPIRE].

[54] D.J. Broadhurst, *Massless scalar Feynman diagrams: five loops and beyond*, Tech. Rep. OUT-4102-18, Open University, Milton Keynes (1985) [INSPIRE].

[55] D.J. Broadhurst, *Exploiting the 1.440 Fold Symmetry of the Master Two Loop Diagram*, Z. Phys. C **32** (1986) 249 [INSPIRE].

[56] I. Bierenbaum and S. Weinzierl, *The Massless two loop two point function*, Eur. Phys. J. C **32** (2003) 67 [hep-ph/0308311] [INSPIRE].

[57] D.I. Kazakov, *Multiloop Calculations: Method of Uniqueness and Functional Equations*, Teor. Mat. Fiz. **62** (1984) 127 [INSPIRE].

[58] D.T. Barfoot and D.J. Broadhurst, $Z_2 \times S_6$ Symmetry of the Two Loop Diagram, Z. Phys. C **41** (1988) 81 [INSPIRE].

[59] D.J. Broadhurst, J.A. Gracey and D. Kreimer, *Beyond the triangle and uniqueness relations: Nonzeta counterterms at large N from positive knots*, Z. Phys. C **75** (1997) 559 [hep-th/9607174] [INSPIRE].

[60] K.G. Chetyrkin, M. Faisst, C. Sturm and M. Tentyukov, *epsilon-finite basis of master integrals for the integration-by-parts method*, Nucl. Phys. B **742** (2006) 208 [hep-ph/0601165] [INSPIRE].

[61] O. Schnetz, *Quantum periods: A Census of $\phi^4$-transcendentals*, Commun. Num. Theor. Phys. **4** (2010) 1 [arXiv:0801.2856] [INSPIRE].

[62] D.J. Broadhurst, *Dimensionally continued multiloop gauge theory*, hep-th/9909185 [INSPIRE].

[63] E. Panzer and O. Schnetz, *The Galois coaction on $\phi^4$ periods*, Commun. Num. Theor. Phys. **11** (2017) 657 [arXiv:1603.04289] [INSPIRE].

[64] F. Brown and C. Dupont, *Lauricella hypergeometric functions, unipotent fundamental groups of the punctured Riemann sphere, and their motivic coactions*, arXiv:1907.06603 [INSPIRE].

[65] P.A. Baikov and K.G. Chetyrkin, *The structure of generic anomalous dimensions and no-$\pi$ theorem for massless propagators*, JHEP **06** (2018) 141 [arXiv:1804.10088] [INSPIRE].

[66] J.M. Henn, A.V. Smirnov, V.A. Smirnov and M. Steinhauser, *A planar four-loop form factor and cusp anomalous dimension in QCD*, JHEP **05** (2016) 066 [arXiv:1604.03126] [INSPIRE].

[67] B. Eden and V.A. Smirnov, *Evaluating four-loop conformal Feynman integrals by D-dimensional differential equations*, JHEP **10** (2016) 115 [arXiv:1607.06427] [INSPIRE].

[68] R. Brüser, C. Dlapa, J.M. Henn and K. Yan, *Full Angle Dependence of the Four-Loop Cusp Anomalous Dimension in QED*, Phys. Rev. Lett. **126** (2021) 021601 [arXiv:2007.04851] [INSPIRE].

[69] T. Peraro, *FiniteFlow: multivariate functional reconstruction using finite fields and dataflow graphs*, JHEP **07** (2019) 031 [arXiv:1905.08019] [INSPIRE].