Cutting massless four-loop propagators

Among the unitarity cuts of 4-loop massless propagators two kinds are currently fully known: the 2-particle cuts with 3 loops corresponding to form-factors, and the 5-particle phase-space integrals. In this paper we calculate master integrals for the remaining ones: 3-particle cuts with 2 loops, and 4-particle cuts with 1 loop. The 4-particle cuts are calculated by solving dimensional recurrence relations. The 3-particle cuts are integrated directly using 1->3 amplitudes with 2 loops, which we also re-derive here up to transcendentality weight 7. The results are verified both numerically, and by showing consistency with previously known integrals using Cutkosky rules. We provide the analytic results in the space-time dimension 4-2{\epsilon} as series in {\epsilon} with coefficients being multiple zeta values up to weight 12. In the ancillary files we also provide dimensional recurrence matrices and SummerTime files suitable for numerical evaluation of the series in arbitrary dimensions with any precision.


Introduction
Inclusive physical observables like total scattering cross sections and related quantities are naturally defined within the perturbation theory in terms of cut Feynman integrals. Particularly, particle decay cross sections at the level of N 3 LO in massless QCD require the knowledge of cuts of four-loop massless propagator-type integrals (two-point functions).
While the optical theorem allows one to calculate massless total cross sections [1,2] with only the knowledge of the discontinuities of the propagators-without calculating each cut separately-the knowledge of the cuts is necessary when they are used as a building block in calculations of massive processes, exclusive processes, or subtraction terms. Examples of such cases include [3], where a subset of four-loop massless propagator cuts were used in the large-mass expansion procedure needed for the boundary condition of differential equations of B decay master integrals. Also [4], where cuts of three-loop propagators were used to develop infrared subtractions scheme for exclusive 2-jet production at NNLO. And finally [5,6], where cuts of three-loop massless propagators were used in boundary conditions for differential cross section master integrals.
The particular use case for cuts of four-loop propagators we have in mind is the extraction of NNLO time-like splitting functions from a semi-inclusive decay cross section at N 3 LO. The time-like splitting functions are currently known from the space-like case via an analytic continuation procedure, which leaves one of the terms in quark-gluon and gluon-quark NNLO terms undetermined [7][8][9]. A direct calculation is needed to fix those terms. As discussed in [6,10] this direct calculation requires the knowledge of decay cross sections differential in the energy fraction of one of the outgoing partons-a calculation for which the master integrals are not yet known, but can be determined via the differential equations method [11,12], as long as the boundary conditions are known to fix the integration constants. The cuts of four-loop propagators can be used for these boundary conditions by noting that a differential cross section integrated over its parameters must give precisely the fully inclusive one, so integrating over a semi-inclusive master integral will relate it to the cuts of four-loop propagators, providing enough information to fix the integration constants.
These considerations motivate us to calculate all the cuts of four-loop massless propagators. Their 2-particle cuts (three-loop form factors) are already known from [13][14][15], 5-particle cuts are known from [16], and in this article we shall complete this knowledge by calculating the master integrals for the remaining 3-and 4-particle cuts. To this end we shall use dimensional recurrence relations [17,18] as well as direct phase-space integration. At the end we shall gather the values of all the cut master integrals-old and new-and present all of them as ε-series in the space-time dimension of 4 − 2ε, with coefficients being multiple zeta values (MZVs) [19] up to transcendentality weight 12. As an intermediate step of our method we shall also present a recalculation of 1→3 two-loop amplitudes up to weight 7 (building upon the known weight-4 results from [20,21]).

What are cut integrals?
To calculate a total cross section of an (off-shell or massive) particle decay, one needs to integrate the probability density of the final state over its phase space in all possible configurations, σ ∼ n dPS n p 1 , . . . , p n |S|q 2 , where S is the scattering matrix, and the phase-space element is defined as Once the scattering amplitude p 1 , . . . , p n |S|q is expanded in perturbation theory as a sum of Feynman diagrams, p 1 , . . . , p n |S|q = A (1) 1→n + A (2) 1→n + · · · = + + . . . , (1.3) expanding the modulus squared in eq. (1.1) gives rise to phase-space integrals of the form 1→n A (2) 1→n * + · · · = dPS n *

Notation for the integrals
Throughout the paper we shall define our cut integrals in d space-time dimensions as where dPS n is the same as eq. (1.2); l and l are the loop momenta; k and k are the propagator momenta, being linear combinations of l, l , and the cut momenta p i ; and the signs of the i0 prescription depend on whether the propagator is to the left or to the right of the cut. Note that the i0 prescription is only relevant for the propagators involved in the loops, it does not matter for the rest of them.
In practical terms, it is often convenient to factor out a common prefactor from this definition, and only work with the normalized integrals without it. For an integral in d dimensions with n cut lines, m L loops to the left of the cut, m R loops to the right, and p propagators we shall factor out the full n-particle phase space PS n , m L + m R one-loop bubbles, and the full q 2 dependence as follows: with the bubble B given by , (1.8) and the full n-particle phase space by PS n = . . . (1.9) Note that both normalization factors-B and PS n -are dimensionless, so the power of q 2 directly corresponds to the dimensionality of the integral.
This normalization removes all q 2 dependence along with any imaginary numbers from J, making the normalized integrals real functions of d, and removing the distinction between the "left" and the "right" amplitudes (this distinction is fully captured by the prefactors). It also cancels the surface UV divergencies of the integrals; this, for example, makes one-loop integrals finite when d is high enough to suppress the IR divergences. Finally, it normalizes all of the following integrals to unity, simplifying dimensional recurrence relations: VVVV 1 from Table 2, VVVR 1 from Table 13, VVRV 3 from Table 14, VRRV 1 from Table 5, VVRR 4 from Table 4, VRRR 16 from Table 3, and RRRR 1 from Table 15. Topology H. Propagators in the indicated order:  Topology N. Propagators: Topology L. Propagators:  To identify the master integrals for the cuts one can try to solve the integration-by-parts (IBP) relations [26], but a simpler procedure turns out to be sufficient: it is enough to construct all possible 2-, 3-, 4-, and 5-particle cuts of these 28 VVVV integrals, and remove the symmetric duplicates from the obtained set. Further application of IBP reduction reveals no additional linear relations between the obtained cut integrals-we have checked this by using a combination of Fire [27] and LiteRed [28] during the computation of dimensional recurrence relations. Moreover, we believe that the set of integrals obtained this way constitutes a complete basis, because by the optical theorem any full cross section can be expressed via the discontinuity of the VVVV integrals, and any discontinuity of those can be expressed as a linear combination of these cut integrals (we shall construct these relations explicitly in Section 7).   Table 3: VRRR, 4-particle cut master integrals. The index " * " denotes cut propagators.

Direct integration over the phase space
The four-particle phase-space is quite complicated [29][30][31]. The parametrization of dPS 4 from eq. (1.2) in terms of the scalar products has 5 degrees of freedom, and the following form: where ∆ 4 is the Gram determinant, Because the shape of the integration region is given by Θ(∆ 4 ) Θ(s ij ), with ∆ 4 being a polynomial of the third degree, integrating over this region in the general case is a challenge. Parametrizations such as the "tripole parametrization" [29,30] exist that remap this shape onto a hypercube, but they do so at the expense of introducing square roots in the mapping from s ij to the new parameters, which only allows one to complete the integration analytically if the said roots drop out from the integrand.
Still, for simpler integrals such an approach is viable. Integrating VRRR 1 using the tripole parametrization, we find Similarly for VRRR 4 we get (3.5) In [3] a number of other VRRR integrals were calculated in a similar fashion. In the general case though, this approach does not apply, and thus we shall turn to computing these integrals differently: by solving the dimensional recurrence relations.

An overview of dimensional recurrence relations
Any Feynman integral I can be cast into a parametric representation (e.g. Feynman) in which the space-time dimension d becomes a free parameter. This is true for cut integrals as well, where a Baikov-like representation can be used (we shall use it in Section 3.3). With an insight gained from [17], it is then possible to shift d by ±2, and express the resulting integral I(d ± 2) as a linear combination of integrals in space-time dimension d. Thus, a set of master integrals with shifted dimension J i (d ± 2) can be expressed as a linear combination of the same master integrals in space-time dimension d via the IBP reduction, giving us dimensional recurrence relations (DRR). Using the terminology of [28], the "lowering" DRR take the form We have constructed such DRR systems for all the cut configurations and VVVV using LiteRed with Fire. For all of them the matrices M turn out to be triangular, and eq. (3.6) can be easily split into the homogeneous and inhomogeneous parts, The triangular form here is guaranteed by the fact that only a single master integral is present in each sector, and thus no coupled blocks form. The general solution of this recurrence relation system can be represented as is a particular solution, determined by integrals from the lower sectors J j<i ; • ω i (d) is an arbitrary periodic function, such that ω i (d + 2) = ω i (d); this function cannot be determined from the DRR relations alone, and needs to be fixed separately.
The triangular form greatly simplifies the construction of the homogeneous solution compared to the case of coupled blocks (explored for example in [32,33]). For the diagonal entries M ii of the form the homogeneous solution can be immediately constructed as Both forms of the factors are acceptable, and generally one should choose one or the other depending on where it would be most convenient to have the poles and zeros of H(d) located. The function HomogeneousSolution from the Dream package [34] automates this construction.
The particular solution R(d) can be constructed as an infinite sum, (3.11) where the direction of the summation is chosen depending on which one converges. With the help of Dream this sum can be evaluated numerically as a series in ε with arbitrary precision-as long as the integrals from lower sectors J j<i are known, of course. Normally this is a quickly converging geometric sum, and a precision of thousands of digits can be easily achieved.
With H and R being well understood, the most difficult part of solving eq. (3.7) is finding the periodic function ω(d), which plays the same role as integration constants play in the solution of differential equations.

Solving DRR for VRRR integrals
To fix ω(d) we shall loosely follow the "dimensional recurrence and analyticity" method (DRA) from [18]. As we shall soon see, only a simple incarnation of it is needed for the VRRR integrals, but we shall return to the full method in Section 8.
The overarching idea is to find restrictions on the possible forms that ω(d) is allowed to take by analyzing its analytic properties in the whole complex plane. Two essential sources of information are used for this: the location of the poles in d of ω(d), and its asymptotic behavior in the limit Im d → ±∞.
To perform such an analysis, let us first rewrite eq. (3.8) in the form Because ω(d) is periodic with a period of 2, we can restrict the analysis to a stripe in the complex plane where Re d ∈ (d 0 , d 0 + 2] or Re d ∈ [d 0 , d 0 + 2). It is useful to choose this stripe such that J, H −1 , and R have as few poles on it as possible. The poles of J(d) are particularly to be avoided.
Conveniently, all VRRR integrals normalized according to eq. (1.7) are finite in the stripe (6,8]. All H i can be chosen according to eq. (3.10) to not have any zeros on this stripe too. To analyze the behavior of R i constructed via eq. (3.11), the knowledge of the previous integrals J j<i is required, having which we can evaluate R i (d) numerically via Dream. Proceeding with the solution steps for each J i , and plotting the numerical values of R i (d) for d ∈ (6,8], every time we find that R i (d) is smooth and finite-which is expected, because in eq. (3.11) we have an infinite sum with finite coefficients that converges geometrically.
Because none of the integrals J(d), H −1 (d), or R(d) have poles when Re d ∈ (6,8], it follows that ω(d) is free from poles in d in the whole complex plane.
Next we shall turn to the investigation of the asymptotic behavior of ω(d) in the limit of Im d → ±∞. For this we shall investigate each component of eq. (3.12) separately.

Asymptotic behavior of VRRR integrals
To investigate the behaviour of VRRR integrals at Im(d) → ±∞, a parametric representation is needed. To construct it, let us first rewrite eq. (1.6) in the form where l is the loop momentum, A L is the loop part of the amplitude (i.e. the propagators that contain l), k i are some linear combinations of the cut momenta p j , A T is the tree part (i.e. the propagators that do not contain l), and J is the normalized integral.
The loop amplitude A L entering this form can be parametrized with Feynman parameters, giving us For the phase-space element dPS 4 the direct parametrization from eq. (3.2) can be used. Specifically, to analyze the normalized integral J, this form with a factorized PS 4 will be useful: To analyze the asymptotic behaviour of |J(d)| in the limit Im d → ∞, it is enough to note that eqs. (3.14) and (3.15) contain only two kinds of structures involving d: the Gamma functions Γ(α + βd) and the powers x αd . The asymptotic behavior of Γ(z) follows from the Stirling formula, As long as the base of a power is positive, its modulus does not depend on Im d at all, This suits our case, because the base of the power in eqs. (3.15) and (3.14) are all in fact always positive. Combining all together, it follows that |J(d)| is bounded in the limit by

Asymptotic behavior of the homogeneous solutions
All the diagonal elements of the DRR matrix for VRRR integrals have the form of eq. (3.9), with all a k ∈ [0, 7/3]. For this reason it is possible to construct the homogeneous solution via eq. (3.10), and to avoid both poles and zeros at d > 14/3 by choosing the first form of the factors. Then, applying the asymptotic of the Gamma functions to this construction, for each VRRR master integral we find that where α depends on the integral, but is always a small rational number, α ∈ [0, 5/2]. It is important to note the absence of factors like e α|Im d| in this asymptotic that could appear from eq. (3.16): it turns out that they cancel each other.

Asymptotic behavior of the particular solutions
Next, let us look at the particular solution R(d) in the form of the infinite series from eq. (3.11). The bounds established so far tell us that each term in that series is asymptotically bounded by C(Re d) |Im d| β , for some values of β. To claim that the series as a whole is bounded similarly too, it is enough to show that a series composed of these bounds converges. This is in fact the case, because the dependence of the series terms on Re d approaches an exponential irrespective of the value of Im d, and the whole series converges geometrically.

Asymptotic behavior of the periodic function
Because J, H, and R all have the same form of asymptotic behavior, from the definition of ω(d) in eq. (3.12) it follows that ω(d) is bounded similarly, for some γ. Furthermore, because ω(d) is periodic, we can view it as a function of z, In terms of z, the limit Im d → +∞ corresponds to z → 0, and Im d → −∞ corresponds to z → ∞. Because ω(d) has no poles in d, ω(z) viewed as a function on the Riemann sphere can only have poles at 0 and ∞. Moreover, its growth at both of these limits is bounded by |Im d| γ ≈ |ln z| γ , which grows slower than any non-zero power of z. Then, representing ω(z) as its Taylor series, these constraints mean that only the z 0 term is allowed. In other words, this means that ω(z) can only be a constant.

Fixing the constants
Now that we have determined that all ω i for VRRR master integrals are constants, what is left is fixing them. This is just one constant per integral, so it is sufficient to for example calculate the leading term of ε-expansion of each integral. To make this easy, one can use the following observation: for a VRRR master integral with n loop propagators, the superficial degree of divergence of the loop part becomes zero when d = 2n (meaning that the integral begins to diverge logarithmically in the UV region), and importantly no IR divergences are present in this d as well. Changing d to 2n − 2ε regulates the UV divergence via a single ε pole, and being an UV pole, it does not depend on any masses in the diagram. Therefore, one can just as well insert some mass m into the loop without affecting the pole. Then, applying the large mass expansion [35] to the massive diagram factorizes it into a massive one-loop vacuum bubble (equal to the massive loop with external legs amputated) and a 4-particle phase-space integral (equal to the original integral with the loop shrinked into a dot), while still not changing the pole.
For VRRR 31 this process can be illustrated as follows: all with d = 10 − 2ε. The vacuum bubble here was evaluated via 24) and the 4-particle phase-space integral was reduced to the masters from [29] using IBP and dimensional recurrence relations.
From here it is possible to determine ω 31 by inserting eq. (3.23) into the definition of ω from eq. (3.12) along with a high-precision numerical series for R 31 (d) calculated by Dream.
Finally, the same procedure needs to be performed for all other VRRR master integrals. Interestingly, we have found that all ω i are identically zero, except for ω 1 , ω 4 , and ω 16 , which all correspond to simple integrals.
This concludes the calculation of ω i for VRRR master integrals. Together with H i from eq. (3.10) and R i from eq. (3.11) this gives us the full solution to the DRR eq. (3.6) in terms of nested infinite sums. With the help of Dream or SummerTime these sums can be evaluated as a series in ε around arbitrary d to any desired precision. We have done so for d = 4 − 2ε with 4000 digits of precision and have restored the analytical answers via the PSLQ [36] algorithm in the basis of MZVs [19] up to weight 12 (see Appendix C for the precise definition of this basis). The results of this restoration up to weight 6 can be found in Section A.1. The full results up to weight 12 as well as the explicit expressions for the sums in SummerTime format can be found in the ancillary files, as described in Appendix D.

Cross-checks
One way to cross-check the obtained results is to recalculate them numerically. To this end the form of eq. (3.13) can be used along with the Feynman parametrization of the loop amplitude from eq. (3.14) and the tripole parameterization of the phase-space [29,30]. Because VRRR integrals contain only one loop, their ultraviolet divergences manifest themselves only in the prefactor of the Feynman parameterization. The infrared divergences also disappear if we look at the series around d ≥ 6. For this reason this parameterization can be integrated directly via the standard Monte-Carlo numerical integration methods (we use the Vegas algorithm [37] implemented in Cuba [38]) in d = 6 − 2ε, and then lowered back to d = 4 − 2ε via DRR-avoiding the need for methods like sector decomposition.
We have performed this numerical integration and can report that for each VRRR integral the results for the first 3 orders in ε match with the analytic answers within 1% accuracy.
Finally, a cross-check based on Cutkosky rules will be described in Section 7.

2-loop 3-particle cut integrals (VVRR, VRRV)
There are 22 VVRR and 9 VRRV master integrals; 27 in total, if one omits the duplicates between these two sets. Both of these sets are depicted in Table 4 and Table 5.
Ideally we would like to calculate VVRR master integrals by solving the dimensional recurrence relations just as we did in Section 3. The difficulty here lies in the fact that while for VRRR only one constant per integral was needed to uniquely fix the periodic function ω(d), for VVRR integrals dozens might be needed. Therefore we shall postpone solving DRR-until Section 8. Fortunately the 3-particle phase space is considerably simpler than the 4-particle one, and we can return to the idea of direct integration over it.
In principle, the 3-particle phase-space volume element from eq. (1.2) can be parameterized in terms of the kinematic invariants in the following way: and the integration volume given by the δ-function and the condition s ij ≥ 0 is simple enough that it is often possible to integrate eq. (1.6) directly, as long as both amplitudes entering the integral are known. Specifically, one way to make the integration volume explicit is In practice many of the amplitudes are only known as series in ε, and a series expansion operation does not necessarily commute with integration. To illustrate this issue, let us   Table 5: VRRV, 3-particle cut master integrals with 1 loop on each side of the cut.
consider the following integral: The amplitude in the integrand here is a single-scale integral (the precise value is given by eq. (B.7)), and just by dimensional analysis must be proportional to q 2 s 12 If we were to first integrate eq. (4.4) directly and then expand the result in ε, we would get (4.6) In contrast, expanding both the amplitude and the phase-space element in ε first using multiplying it by the expansion of the phase-space element from eq. (4.2), and then integrating order-by-order in ε would result in divergences corresponding to the limit s 12 → 0 in each order of the series: This should not be surprising, seeing that the true ε-expansion in eq. (4.6) starts with ε −1 , while the integrand in eq. (4.10) starts with ε 0 , so integrating it order-by-order can never give the expected series. In other words, the integral in eq. (4.4) is not infrared finite, and the pole in ε regulating this infinity cannot be obtained by integration in 4 dimensions.
Our solution to this problem stems from the fact that every Feynman integral is free from infrared divergences if the dimension of space-time is high enough. This can be seen already from eq. (4.2): higher d results in higher powers of the (s 12 s 13 s 23 ) factor, which can eventually compensate any singularity of the integrand in the infrared region. In particular, at d = 6 all our 3-particle master integrals are infrared-finite. Thus, we can overcome the divergence in eq. (4.10) by integrating the ε-series not around d = 4 − 2ε, but rather around d = 6 − 2ε. Once this is done, we can use dimensional recurrence relations to restore the series in d = 4 − 2ε.
This procedure requires the knowledge of 1→3 amplitudes at 1 and 2 loops as series in ε, expanded around d = 6 − 2ε, and DRR for the phase-space master integrals themselves.
Note that there is nothing magical about the basis at d = 6 − 2ε specifically: any IR-finite basis would be sufficient for our purposes.

1-loop 1→3 amplitudes
Up to relabeling of p i , all the 1→3 amplitudes at 1 loop fit into the box topology: (4.11) There are only four master integrals in this topology, with only two being meaningfully distinct: the bubble and the box with one off-shell leg.  In addition to the master integrals, a triangle with two off-shell legs also appears in the VRRV integrals. It can be found via IBP as (4.12)

2-loop 1→3 amplitudes
The master integrals for these amplitudes were first calculated in [20,21]. The results there are provided as series in ε with coefficients in terms of "2dHPLs", a subclass of multiple polylogarithms, up to transcendentality weight 4. This turns out to be insufficient for our needs, because the ε-finite parts of the VVVV integrals are known to contain ζ 7 , and thus we need the amplitudes to at least that weight as well. Thus, a re-derivation of these master integrals is required.
The overall idea of the method is to write down the differential equation system for the master integrals in external kinematic invariants, solve it, and determine the integration constants by enforcing regularity conditions on the solution. Let us do this step by step.

Topologies and master integrals
We start by determining IBP topologies that cover the master integrals. Up to relabeling of p i , three topologies are sufficient: one planar ("PA"), and two non-planar ("NA" and "NB"). These are depicted in Table 9.
Solving the IBP relations in these topologies using LiteRed [28] we arrive at 18 master integrals in PA, 22 in NA, and 29 in NB. Note that there is an overlap between these three sets of master integrals, but this will not bother us.  Table 9: Generic topologies for the 1→3 master integrals at two loops.  There are only three independent external kinematic invariants for the 1→3 amplitudes: p 2 12 , p 2 13 , and p 2 23 (p ij ≡ p i + p j ), with the incoming energy q 2 being the sum of all three. This allows us to extract one of them as a dimensionful prefactor, with a function of two dimensionless ratios remaining: where k is determined from dimensionality (+ d 2 for each loop, −1 for each propagator), (4.14) Solution via differential equations in the ε-form  To solve these systems as series in ε, it is convenient to transform them into the combined ε-form, where the dependence on d of both M matrices is simultaneously factorized by a basis change: where J i is the set of master integrals in the ε-basis, related to the original basis I i via a transformation matrix T , and S ij are matrices of the form  In particular our S matrices only contain these six denominators (the "alphabet"): which correspond to singularities at limits s ij = 0, and s ij = 1.
Once the ε-form is constructed, writing down the solution as a series in ε becomes easy.
First, start at some arbitrary point, for example (y, z) = (0, 0), and fix the value of J ≡ {J i } at the point as a series of constants, where k 0 is some arbitrary starting order of the series, chosen high enough to cover the highest ε-pole of the integrals.
Then, integrate along the z axis using the right side of eq. (4.16) to obtain where W is the fundamental solution matrix of the corresponding differential equation system constructed as a path-ordered exponential along the specified path,

Solution as multiple polylogarithms
Because the S matrices have the structure of eq. (4.18), the iterated integrals in the general solution will result in terms of the form This is a class of functions known as "multiple polylogarithms" [39], "Goncharov polylogarithms" (GPLs), or "hyperlogarithms". If the set of parameters is restricted to {0, 1, −1}, these correspond to a well known class of "harmonic polylogarithms" [40]. The set {0, 1, 1 − x, −x} corresponds to "2d harmonic polylogarithms" [20].
Note that for a consistent definition of GPLs, one needs to introduce a special case for Because we have chosen to first integrate along the (0, 0) → (0, z) direction, and our alphabet is limited to eq. (4.19), eq. (4.23) will result in each J i having the form where C are some constants, and the parameters of G are taken from sets w z ∈ {0, 1, 1 − z, −z}, and w c ∈ {0, 1}. (4.28) Finding the combined ε-form The possibility and usefulness of transforming the differential equations into the form of eq. (4.16) were first pointed out in [41]. Later in [42] an algorithm was described that constructs such transformations (in the form of the matrices T ij ) for differential equation systems in one variable. The algorithm was later updated in [43] and [44, section 8]. Two implementations of that algorithm followed, in the form of Fuchsia [45] and Epsilon [46] tools. Other approaches include Canonica [47] with its own algorithm based on a rational ansatz for T ij explicitly designed for differential equations in multiple variables, and an approach based on finding integrals with constant leading singularities [48].
While not being described in the original paper [42], the same algorithm treating singlevariable case can be reused for the case of multiple variables too. Here is how: 1. Look at the first system from eq. (4.15), the differential equation system in y, and construct a transformation I i = T (I,y) ij J (y) j that reduces it to an ε-form as described in [42]. The variable z stands as a free parameter during this reduction.
2. Write down the differential equation system in the second variable, z, for the new basis J  This is an outline of the approach that Fuchsia uses to find ε-form transformations for multi-variate systems, 1 and how it was possible for us to compute the three T (I) ij from eq. (4.17).

Fixing the integration constants
To determine all the integration constants, it is sufficient to use the following three conditions. i0 prescription is specified here explicitly to fix the analytic continuation of the ln(−1) terms that will appear if one is to expand this prefactor in ε. Practically speaking, because all the integrals are proportional to this factor, there is no need to expand it: it is better to separate it out into a common prefactor, and only work with the remaining part of the integral, which is free of imaginary numbers.
2. All of our integrals are massless, and therefore must only have discontinuities at limits s ij → 0 and nowhere else. On the other hand, the differential equations we are solving also have poles at s ij → 1, as the list of denominators in eq. (4.19) demonstrates. Thus, requiring that the apparent discontinuities of the general solution at s ij → 1 vanish will generate nontrivial identities between the integration constants. This requirement can be written down by separating the terms proportional to ln(1 − s ij ), and enforcing that the coefficient in front of them vanishes in the limit.
3. The planar integrals only have discontinuities at limits where adjacent momenta go to zero. For the PA topology this means that it should be regular at s 12 → 0 (i.e. y + z → 1), as long as q 2 = 0. Similarly for the planar integrals from other topologies. Here again we are looking at the logarithmic terms like ln(s 12 ), enforcing the cancellation of the coefficients in front of them in the limit.
To apply the regularity conditions above one needs to separate the terms proportional to ln k (1 − s ij ), and require that the coefficient of each is exactly zero in the limit s ij → 1.
For the limit y → 1, to separate the divergent logarithms, it is enough to employ the GPL shuffle relations to rewrite every G(1 . . . , − → w y ; y) in eq. (4.27) into a product of the divergent factor G(1 . . . ; y) and a part finite at y → 1. For the limit z → 1 the same cannot be done directly on eq. (4.27), because z appears in the parameter list of G(. . . ; y). Instead, we can rewrite eq. (4.27) into the reverse form, where w y ∈ {0, 1, 1 − y, −y}, (4.30) and factor our logarithmic terms from that.
Such a rewrite can be achieved by the quite general technique of recursively differentiating the GPL and then integrating it back: Practically speaking, a much more general incarnation of this technique is implemented in the fibrationBasis routine from HyperInt [49]. In our experience it is powerful enough to perform these transformations up to weight 9. At weight 10 the recursion results in too many terms so that memory and performance optimizations become necessary. For this reason we had to implement this transformation manually in C++ with the help of GiNaC [50].
Note that applying the regularity conditions above to the weight-10 expansion of eq. (4.23) only fixes the constants C up to weight 8. These results are quite large (megabytes of text), and we are providing them in machine-readable format in the ancillary files, as described in Appendix D.
Interestingly, weight 8 amplitudes are insufficient to compute all of the VVRR integrals to weight 8 as well: because of an apparent cancellation during DRR, some of the weight 8 information is lost, and VRRR integrals 20-23 can only be obtained to weight 7. Still, this is sufficient for practical needs, because weight 7 corresponds to terms of the order ε 1 or higher. Moreover, after cross-checking these results we shall improve upon them in Section 8.

Cross-checks
One way to cross-check the obtained 3-particle cut integrals is to calculate them numerically. This can be conveniently done by using sector decomposition as implemented in Fiesta4 [51], all that is required is a parametrization of the integrals suitable for the SDEvaluateDirect function (Fiesta does not construct such parametrizations automationcaly for cut integrals).
To that end, the Feynman parametrization can be used for the loop parts of the integrals, and the naive parameterization from eq. (4.2) can be used for the phase-space parts.
Proceeding this way we were able to calculate the first 3 terms of the ε-expansion for each VVRR and VRRV integral numerically in d = 6 − 2ε and d = 8 − 2ε. Note that we were unable to do the same for the integrals in d = 4 − 2ε: due to fairly high pole orders, Fiesta has identified thousands of sectors in some cases, and the whole calculation effectively froze. Still, a numerical match within 1% of accuracy for both 6 − 2ε and 8 − 2ε gives us confidence in both the 4 − 2ε results and the DRR matrices.
Additionally, we have compared our results to the weight-6 series for VVRR 16 , VVRR 8 , and VVRR 9 reported in [3], and find them to match fully.

3-loop 2-particle cut integrals (VVVR, VVRV)
There are 22 distinct 2-particle cuts in total, all depicted in Table 13 and Table 14. These integrals correspond to 3-loop form factors, and have been calculated in [13][14][15]. We list these integrals here for completeness, and as a reference to the order in which their values are presented in the ancillary files (as described in Appendix D).  Table 13: VVVR, 2-particle cut master integrals with 3 loops on one side of the cut.

5-particle phase-space integrals (RRRR)
There are 31 master integrals for the 5-particle cuts. All of these have been calculated in [16], and we list them here for completeness. These integrals will be important for the Cutkosky relations in Section 7.

Relations from Cutkosky rules
By this point we already know all the cuts of VVVV integrals to at least weight 7. A good test for the consistency of all these results is the optical theorem (Cutkosky rules, to be more precise), which relates the discontinuity of the VVVV integrals to its cuts.
The general optical theorem comes from the requirement of unitarity of the scattering matrix S, S † S = 1.
Introducing the transition matrix T as For a decay of a single particle with momentum q, rewriting this relation in terms of the transition amplitudes produces (7.4) This is the optical theorem. Note that once the transition amplitude is expanded in terms of individual Feynman diagrams, the right-hand side will consist precisely of our cut integrals as in eq. (1.6).
Cutkosky rules [52,53] provide a stronger form of this relation that holds not only for the whole transition amplitude q| iT |q , but for each individual Feynman diagram F that comprise it too, where the sum goes over all possible cuts of the diagram, each cut being a partition into two sides, with the right-hand side complex-conjugated, and the propagators between sides set on shell, exactly like in eq. (1.6).
To write down these relations for the VVVV master integrals in this convenient form, we'll augment our integrals with Feynman rules stemming from a simple scalar field theory with the Lagrangian of the form The momenta-space Feynman rules corresponding to this Lagrangian are An additional prescription for cut integrals is this: every vertex and propagator to the right side of the cut needs to be conjugated.
Note that the values of λ n are not important for us here, because a cut of a diagram will have the same overall λ factor as the initial diagram, which will thus factor out from eq. (7.5).
With this in mind, writing down eq. (7.5) for each Feynman diagram corresponding to a VVVV master integral, and mapping the cuts onto our master integrals, we obtain the following relations: Inserting the values of our cut integrals into these relations, we find that they all hold precisely. This concludes our cross-check.

Dimensional recurrence relations for 3-particle cut integrals
In Section 4 we have postponed solving the dimensional recurrence relations for 3-particle cut master integrals because we did not have enough information to fix all the ω i . Since that time we have obtained multiple terms of the ε expansion, as well as the Cutkosky relations. This information combined will be enough to fix ω i and solve the DRR. Such a solution will upgrade our knowledge of 3-particle cut integrals from weight-7 series to weight-12, and more generally will provide expressions that can be evaluated in arbitrary d to arbitrary ε order numerically with any desired precision.

DRA method by example: VRRV
We have already gone through the process of solving DRR in the simpler case of VRRR integrals in Section 3. Let us now consider the more complicated cases. To summarize, following the "dimensional recurrence and analyticity" (DRA) method described in [18], to solve DRR in the form of eq. (3.7) for a given integral, one needs to: 1. Solve all integrals in subsectors.

Choose a semi-open stripe of width 2 in the complex plane, such that
, and restrict the analysis to this stripe. 7. Fix the constants in the ansatz by expanding eq. (3.12) near the poles, and/or other considerations.
Let us illustrate these steps using a simple, but non trivial case of VRRV 8 from Table 5.
The VRRV family of integrals is considerably simpler than the VVRR, because the loop part of these integrals factorizes into a product of two one-loop integrals, for which arbitrary-d expressions are known (see Section 4.1), and which simplifies the pole analysis and makes a numerical evaluation of the integrals in arbitrary d easy via Monte-Carlo integration.
We shall restrict our analysis to the strip of Re d ∈ (6,8], because all normalized VRRV integrals are finite on it. The reason is that IR divergences are suppressed at d ≥ 6, and surface UV divergences are cancelled by the normalization of eq. (1.7). We can verify this finiteness numerically via Monte-Carlo integration, which can be performed by parametrizing the three-particle phase-space element entering eq. (1.6) via eq. (4.2), and using the analytical expressions for the remaining two one-loop amplitudes found in Section 4.1. The result for VRRV 8 is presented on Figure 1a; it is finite as expected.

The homogeneous solution H(d)
Next, to construct the homogeneous solution to eq. Note that the constant subtracted from d/2 in the factors is always ≤ 7/3. This means that it is possible to construct such a H(d) that is free of poles and zeros at d > 14/3. Using eq. (3.10) and choosing the first form of factors does precisely that, resulting in This answer can also be obtained via the HomogeneousSolution function from the Dream package. This solution is indeed finite and free from zeroes on the stripe (6,8] as can be seen on Figure 1b.

The particular solution R(d)
With H(d) ready, the particular solution R(d) can be constructed via eq. (3.11). The numerical evaluation of the nested infinite sums in that formula is the main functionality of the Dream package, and with its help we can plot R(d) on d ∈ (6, 8] with arbitrary Because J(d) has no poles, and H −1 (d) has neither poles nor zeros on the stripe (6,8], from eq. (3.12) it follows that ω(d) must have the same pole structure as R(d).
The ansatz for the periodic function ω(d) With this information it is now possible to construct an ansatz for ω(d). Note that just as it was for VRRR integrals, applying eq. (3.16) to eq. (8.2) the asymptotic of H(d) at Im d → ±∞ can be found as With similar arguments as in Section 3, one can also show that |J(d)| and |R(d)| asymptotically behave as |Im d| α , for some α, and thus, so does |ω(d)|. For VRRR integrals ω(d) was free from poles, and together with this asymptotic that meant that ω(d) was a constant. The difference with the case of VRRV 8 is that ω(d) does have poles to compensate for the poles of R(d). To deal with them observe that if one is able to subtract from ω(d) an expression that would cancel its poles while simultaneously not altering the form of its asymptotic behavior, then the resulting difference will have no poles, and by the same argument could only be a constant.
For a function that cancels a pole at d i , but does not spoil the asymptotic behavior at Im d → ±∞, one convenient choice is a cotangent function, which in terms of eq. (3.22) is just one of the simplest functions with a single pole that is constant at both z → 0 and z → ∞.
The whole ansatz can then be constructed as a sum of a constant term and C n d i (d) for each pole of multiplicity n, Fixing the constants in the ω(d) ansatz To fix these constants numerically it is enough to insert the series for R(d i − 2ε) obtained via Dream into eq. (3.8), and for each d i from eq. (8.3) demand the cancellation of poles, Note that because J(d) is finite, it drops out from these equations.
After a i have all been fixed numerically with high enough precision, we can reconstruct their analytic values, obtaining the following: Reconstructing the rational constants here is easy enough; for the irrational ones we had to resort to using an educated guess and the Inverse Symbolic Calculator 2 . The result of inserting these values into the ansatz is plotted on Figure 1d.
With the constants from eq. (8.8) now fixed, all that is left is to use Dream to calculate the ε-series for J(4 − 2ε), and restore them in terms of MZVs. The results of this calculation for all VRRV integrals up to weight 6 are listed in Section A.3, and weight-12 results along with the full expressions in SummerTime format are available in the ancillary files, as described in Appendix D. All of these results of course match what we have calculated via direct integration in Section 4.

Solving DRR for VVRR integrals
The case of the VVRR integrals brings another complication on top of what VRRV had: not only the periodic function ω(d), but also the full normalized integral J(d) has poles now. These are ultraviolet poles coming from subdivergences. For VRRR and VRRV integrals UV poles only came from surface divergences, and thus were suppressed by the normalization factors from eq. (1.7); this is not the case for VVRR integrals.
Knowing the location and multiplicity of these poles is important because the ansatz for ω(d) from eq. (8.7) must now include them too. To determine these locations one can apply Fiesta to the same parametrization as in Section 4.3.
The way poles of J(d) complicate the calculations is by preventing the usage of eq. (8.9) when J(d i ) is divergent, because J(d i − 2ε) no longer drops out. In these cases one needs to find additional sources of information to fix the corresponding constants in the ansatz.

VVRR integrals entering Cutkosky relations alone
One source of additional information are the Cutkosky relations from Section 7. Looking closely at relations in eq. (7.7) through eq. (7.34), one can see that most of them contain no more that one VVRR integral. Only VVRR 9 with VVRR 11 , and VVRR 19 with VVRR 20 enter in pairs in eqs. (7.29) and (7.33) respectively. The rest of VVRR integrals can be directly determined from those relations, since we already have DRR solutions for all the other cuts.
Practically speaking, to obtain VVRR integrals in the same format as the other DRR solution, we only use Cutkosky relations numerically: to evaluate ω(d) in 1000 points using eq. (3.12) with 100 digits of precision. From this plot we determine the poles of ω(d), construct an ansatz in terms of the cotangent functions from eq. (8.5), fix its constants numerically via the Mathematica function FindFit, and reconstruct them analytically (more on this later).
The DRR solutions obtained this way can be double-checked by inserting their numerical values into the Cutkosky relations at multiple d with high precision, by inserting their reconstructed analytical expressions into the same relations, and also by comparing them with the series at d = 4 − 2ε calculated in Section 4. All of these checks hold.

VVRR 9 and VVRR 11
Integrals VVRR 9 and VVRR 11 enter the Cutkosky relation eq. (7.29) in a pair. Being different cuts of the same VVVV integral both have identical diagonal DRR matrix elements, M 9,9 = M 11,11 , and thus identical homogeneous solutions. It might seem like untangling them is impossible, but because they enter with different prefactors-one being 2B, and the other B * -and because each Cutkosky relation can be split into the real and imaginary parts, this turns out to be enough to uniquely fix both ω 9 (d) and ω 11 (d) just from eq. (7.29), and then proceed as explained above.

VVRR 19 and VVRR 20
In the case of VVRR 19 and VVRR 20 the Cutkosky relation from eq. (7.33) only constrains the sum of ω 19 (ν) + ω 20 (ν). For these integrals we proceed in the same way as in Section 8.1, using eq. (8.9) with d i for which J(d i ) is finite-this is enough to fix most of the constants. For the rest, we compare the ε-series for J(4 − 2ε) with the ones calculated in Section 4: using the first few terms of those series is enough to fix the remaining constants. The rest of the terms can be used to cross-check the results, together with eq. (7.33), which provides another independent check.

Restoring the analytic expressions for the ω ansatz constants
Once the constants a i,k in the ω(d) ansatz from eq. (8.7) are fixed numerically with high precision, it is desired-even though not strictly speaking necessary-to reconstruct their analytical forms.
For most integrals an educated guess and the Inverse Symbolic Calculator do the job, but VVRR 17 , VVRR 19 , and VVRR 20 proved to be a challenge.
The basic difficulty is that as we have seen in eq. (8.10) these constants contain algebraic numbers, and for example an MZV basis with rational coefficients is insufficient to reconstruct them. In the simple cases when a i,k is a product of an algebraic number and some power of π, one can try dividing its value by consecutive powers of π, and reconstructing the rest in terms of radicals. If however a i,k is a sum of terms with different π powers, or even general MZVs, one would need to split these terms beforehand.
Fortunately we have experimentally observed that the particular solutions R(d) for each of these integrals can in fact be reconstructed in terms of an MZV basis with rational coefficients. Keeping in mind that the full solution J(d) should be representable in terms of MZVs as well, from eq. (3.8) it follows that H(d) ω(d) should be too. Then, expanding around d = 4 − 2ε, where constants c k,j are rational numbers, and can be reconstructed analytically via PSLQ.
From here it is possible to evaluate the left hand side numerically using eq. (3.10) for H and eq. (8.7) for ω, reconstruct c k,j via PSLQ in the basis of MZVs, then evaluate the left hand side again symbolically, and finally solve for the constants a i,k from the ω ansatz. The radicals we have seen in eq. (8.10) will appear here as the result of the expansion of H(4 − 2ε) in an ε-series.
The practical complication here is that fixing O(20) of a i,k constants this way means working with about the same number of terms in the ε-series in eq. (8.11), and thus with higher and higher transcendentality weights, and larger and larger MZV bases, all requiring increasingly high numerical precision. This quickly becomes computationally expensive.
The trick is to exploit the observation that a i,k seem to only contain powers of π, and are free from MZVs otherwise. With this conjecture in mind, we can evaluate the left hand side of eq. (8.11) symbolically, drop all terms multiplied by MZVs that are not powers of π, then evaluate the resulting series numerically, reconstruct them via PSLQ in the basis of even powers of π, and finally return to the symbolical series, solving for a i,k in the form of where a i,k,p are algebraic numbers. Because the basis of π powers is much smaller than the general MZV basis, this procedure becomes computationally tractable.
With all ω(d) ansatz constants finally fixed, we proceed to evaluate the VVRR master integrals as ε-series around d = 4 − 2ε with 4000 digits of precision, and restore the results in terms of MZVs up to weight 12 using the basis from Appendix C. These series truncated after weight 6 are listed in Section A.2. The full weight-12 results along with the expressions in SummerTime format are available in the ancillary files, as described in Appendix D. As a final check, these reconstructed series match what we have calculated up to weight 7 via direct integration in Section 4.

Conclusions
We have presented the calculation of the previously unknown master integrals for 3-and 4-particle cuts of massless four-loop propagators. Together with the already known 2and 5-particle cuts this completes the knowledge of all such cuts. Both direct integration over the phase space and the solution of dimensional recurrence relations were used in the calculation, with the latter finally resulting in expressions that allow the numerical evaluation of the integrals as ε-series to arbitrary order with arbitrary precision via the SummerTime package.
In Appendix A we have provided analytic expressions for 3-and 4-particle cut master integrals restored via PSLQ in the basis of MZVs up to weight 6. The ancillary files described in Appendix D contain analytic results for all cut structures (including 2-and 5-particle cuts, as well as the uncut virtual loop integrals) up to weight 12, as well as the corresponding SummerTime files. We hope that these results gathered in one place will serve as a complete reference.
All the results have been cross-checked via numerical integration, by comparing values obtained via different methods, by comparing with the known results from the literature, as well as by showing consistency with Cutkosky rules.
The same methods used here are applicable to the cuts of five-loop propagators as well, although we expect that calculation to be harder for two reasons: firstly, solving IBP relations for five-loop problems is computationally much more challenging; secondly, fiveloop propagators will have multiple master integrals per sector, so we expect the appearance of coupled blocks in the DRR equations, as well as the appearance of elliptic integrals in the amplitudes. Still, this is one viable direction for further research.
Another direction to go is investigating massive propagators and their cuts, with the massless ones serving as boundary conditions for the differential equations.
Finally, our original motivation of calculating semi-inclusive cut integrals and time-like splitting functions through them now becomes feasible, with the presented cut integrals allowing us to fix the integration constants in the differential equations for the semi-inclusive master integrals. stages of the calculations; and Vsevolod Chestnov for discussions about the properties of Goncharov polylogarithms.

A Results
Here we provide the values of the master integrals for 3-and 4-particle cuts of 4-loop propagators as ε-series around d = 4 − 2ε. The normalization of the integrals is as discussed in Section 1, with prefactors B and PS n defined in eqs. (1.8) and (1.9) respectively.
For brevity, the series here are truncated after transcendentality weight 6, which is enough to cover all the ε-finite parts of the integrals. The full results up to weight 12 are available in the ancillary files on arXiv, see Appendix D for a description of those.
Note that in all our results MZVs are defined as which is a common "physicist" notation adopted by e.g. [19,54], but which has the opposite order of parameters compared to the "mathematician" notation used in [39,49].

B Table of loop integrals
Here we collect integrals used during the 1→3 amplitude calculation. Most of them are taken from [55], eq. (B.7) is from [56].
To obtain a series in ε from these integrals, the hypergeometric function p F q needs to be expanded about its parameters. This can be conveniently done using the Mathematica package HypExp [54].

*.st
SummerTime [25] files for each of the master sets. One can use these to calculate series expansion of any master set around arbitrary d, with arbitrary precision. For example, to calculate the values of the VRRV master integrals as series around d = 4 − 2ε up to order ε 2 with 30 digits of precision, use this command: