Xin Guan ^1,,
, Xiao Liu ^1,,
, Yan-Qing Ma ^1,2,3,,
,

Corresponding author: Xin Guan, guanxin0507@pku.edu.cn

Corresponding author: Xiao Liu, xiao6@pku.edu.cn

Corresponding author: Yan-Qing Ma, yqma@pku.edu.cn

1.School of Physics and State Key Laboratory of Nuclear Physics and Technology, Peking University, Beijing 100871, China
2.Center for High Energy Physics, Peking University, Beijing 100871, China
3.Collaborative Innovation Center of Quantum Matter, Beijing 100871, China
Received Date:2020-05-07
Available Online:2020-09-01
Abstract:We reduce all the most complicated Feynman integrals in two-loop five-light-parton scattering amplitudes to basic master integrals, while other integrals can be reduced even easier. Our results are expressed as systems of linear relations in the block-triangular form, very efficient for numerical calculations. Our results are crucial for complete next-to-next-to-leading order quantum chromodynamics calculations for three-jet, photon, and/or hadron production at hadron colliders. To determine the block-triangular relations, we develop an efficient and general method, which may provide a practical solution to the bottleneck problem of reducing multiloop multiscale integrals.

HTML

--> --> -->

Owing to the good performance of the large hadron collider (LHC), we have entered the era of precision high-energy physics. Some of the most important observables are three-light particles or jet production cross-sections [1-3], which can both be used for testing the strong interaction at high energies and for determining the quantum chromodynamics (QCD) coupling constant. From the theoretical viewpoint, predictions with compatible precision are needed, which requires perturbative QCD calculations up to the next-to-next-to-leading order (NNLO). Although significant advances have been made in the past few years [4-24], a complete NNLO result is still unavailable. One of the main concurrent obstacles is computation of two-loop amplitudes.
To evaluate a two-loop five-light-parton scattering amplitude, one usually first generates an integrand, reduces all of the Feynman integrals to linear combinations of relatively simpler master integrals (MIs), and finally calculates these MIs. Because integrands can be obtained either using the unitarity method [4-9] or using the conventional Feynman diagram method, and because MIs can be calculated analytically [20-24], the bottleneck is the reduction of Feynman integrals. For example, the non-planar contribution of two-loop three-photon production at the LHC cannot be calculated, owing to the lack of such reduction for nonplanar integrals [19].
Reduction is usually achieved by integration-by-parts (IBP) identities combined with Laporta's algorithm [25-35]. Although many interesting proposals have been made recently for improving the IBP reduction [36-49], the problem of reducing multiloop multiscale integrals has not been fully solved yet. The difficulty is twofold. On the one hand, owing to the number of scales, an explicit solution of the IBP system is usually too big to be used in numerical calculations; in addition, it is very difficult to obtain [47-51]. On the other hand, although solving the IBP system numerically in a single run is feasible, one usually needs to solve it many times, for either the phase space integration or fitting analytical expressions, which is very time- and resource-consuming. For example, to reconstruct the fully analytical two-loop five-gluon all-plus helicity amplitude [17], one needs to run the numerical computation of the IBP for nearly half a million times^①. If one uses the same method for reconstructing analytical one-minus or maximal-helicity-violation amplitude, many more IBP calculation runs may be needed, which becomes prohibitive.
We note that a reduction can be obtained efficiently if a system of block-triangular relations is found, which has a small expression size and can be solved numerically very efficiently. Using our proposed series representation of Feynman integrals as input [52, 53], in Ref. [52] we described an algorithm that searched for block-triangular relations and yielded some preliminary results. Although our method developed in Ref. [52] is sufficiently good for reducing integrals with integrands having only denominators, the method is very time-consuming for physical problems that contain integrands with numerators.
In this paper, by further developing the method in Ref. [52], we propose a two-step search strategy along with a reduction scheme that is suitable for physical problems. Based on this, we successfully find out block-triangular relations to reduce integrals in two-loop five-light-parton scattering amplitudes. As expected, the relations are only 148 MB in size, and can be numerically solved hundreds of times faster than using other methods. Our work constitutes an important step towards the complete NNLO QCD calculation for three-jet, photon, or hadron production at the LHC. Because our method is efficient and general, it can be straightforwardly applied to any other process, thus providing a practical solution for the bottleneck problem of reducing Feynman integrals.

To obtain the very much needed reduction of Feynman integrals in two-loop five-light-parton scattering amplitudes, we only need to consider integrals that originate from the four topologies shown in Fig. 1. All the other Feynman integrals are one-loop-like, and can be dealt with much easier.
Figure1. All 8-propagator families: (a) double-pentagon; (b) hexa-box; (c) penta-box; (d) hexa-triangle.

Let us consider the most complicated case, topology (a) in Fig. 1, as an example that will explain what kind of Feynman integrals do we need to reduce. There are five external momenta $ p_1,\cdots,p_5 $ flowing into the diagram, satisfying on-shell conditions $ p_i^2 = 0 $ ($ i = 1,\cdots,5 $) and momentum conservation $ \sum_{i = 1}^5 p_i = 0 $. As a result, this problem contains five independent mass scales, which can be chosen as $ \vec{s} = \{s_1,s_{2},s_{3},s_{4},s_{5}\} $ with $ s_{i}\equiv2p_i\cdot p_{i+1} $ and $ p_6\equiv p_1 $. With two loop momenta $ \ell_1 $ and $ \ell_2 $, a complete set of Lorentz scalars can be chosen as

$\begin{split} D_1 = &\ell_1^2,D_2 = (\ell_1+p_1)^2,\,D_3 = (\ell_1+p_1+p_2)^2,\,\\ D_4 =& \ell_2^2,\,D_5 = (\ell_2+p_3)^2, D_6 =(\ell_1+\ell_2+p_1+p_2+p_3)^2,\,\\ D_7 =& (\ell_1+\ell_2-p_4)^2,\,D_8 =(\ell_1+\ell_2)^2,\\ D_9 =& (\ell_2+p_1)^2,D_{10} = (\ell_2+p_2)^2,\\D_{11} =& (\ell_2+p_4)^2, \end{split}$

(1)

where the first eight are inverse propagators and the last three are introduced to make the set complete. Then, the family of integrals defined by topology (a) can be expressed as

$ I_{\vec{\nu}}(\epsilon,\vec{s}\,) = \int\frac{{\rm{d}}^{4-2\epsilon}\ell_1\, {\rm{d}}^{4-2\epsilon}\ell_2}{({\rm{i}}\pi^{2-\epsilon})^2}\,\frac{D_9^{-\nu_9}D_{10}^{-\nu_{10}}D_{11}^{-\nu_{11}}}{D_1^{\nu_1}...\,D_8^{\nu_8}}, $

(2)

where the indexes $ \nu_{1}, \cdots, \nu_8 $ are integers, $ \nu_9 $, $ \nu_{10} $ and $ \nu_{11} $ are nonpositive integers. Two integrals in this family are said to be in the same sector if the positions of their positive indexes are the same. The degree of an integral is defined by the opposite value of the summation of all of its negative indexes. Finally, we call a degree-m integral is $ \frac{m}{n} $-type if it has n positive indexes and all of these positive indexes are $ 1 $. For example, $ I_{\{1,1,1,1,1,1,1,1,-4,0,-1\}} $ is a degree-5 integral in the top sector, and it is $ \frac{5}{8} $-type.
For later convenience, we define operators $ {\hat{m}^\pm} $ (for a non-negative integer m), which generate a set of integrals in the same sector or its subsectors when acting on an integral. For any integral $ I_{\vec{\nu}} $, $ {\hat{0}^\pm}I_{\vec{\nu}} = I_{\vec{\nu}} $, $ {\widehat{m+1}^\pm}I_{\vec{\nu}} = {\hat m^\pm \hat 1^\pm}I_{\vec{\nu}} $, $ {\hat 1^-}I_{\vec{\nu}} $ generates a set of integrals with one index decreased by 1, and $ {\hat 1^+}I_{\vec{\nu}} $ generates a set of integrals with one nonzero index increased by 1. For example, we have

$ \begin{split}& \hat{1}^+I_{\{1,1,1,1,1,1,1,1,-4,0,-1\}} = \{I_{\{2,1,1,1,1,1,1,1,-4,0,-1\}},\\&\quad I_{\{1,2,1,1,1,1,1,1,-4,0,-1\}}, I_{\{1,1,2,1,1,1,1,1,-4,0,-1\}},\\&\quad I_{\{1,1,1,2,1,1,1,1,-4,0,-1\}}, I_{\{1,1,1,1,2,1,1,1,-4,0,-1\}},\\ &\quad I_{\{1,1,1,1,1,2,1,1,-4,0,-1\}}, I_{\{1,1,1,1,1,1,2,1,-4,0,-1\}}, \\ &\quad I_{\{1,1,1,1,1,1,1,2,-4,0,-1\}}, I_{\{1,1,1,1,1,1,1,1,-3,0,-1\}},\\ &\quad I_{\{1,1,1,1,1,1,1,1,-4,0,0\}}\}\,, \end{split} $

(3)

and

$ \begin{split} &\hat{1}^-I_{\{1,1,1,1,1,1,1,1,-4,0,-1\}} = \{I_{\{0,1,1,1,1,1,1,1,-4,0,-1\}},\\ &\quad I_{\{1,0,1,1,1,1,1,1,-4,0,-1\}}, I_{\{1,1,0,1,1,1,1,1,-4,0,-1\}}, \\ &\quad I_{\{1,1,1,0,1,1,1,1,-4,0,-1\}}, I_{\{1,1,1,1,0,1,1,1,-4,0,-1\}},\\ &\quad I_{\{1,1,1,1,1,0,1,1,-4,0,-1\}}, I_{\{1,1,1,1,1,1,0,1,-4,0,-1\}}, \\ &\quad I_{\{1,1,1,1,1,1,1,0,-4,0,-1\}}, I_{\{1,1,1,1,1,1,1,1,-5,0,-1\}},\\ &\quad I_{\{1,1,1,1,1,1,1,1,-4,-1,-1\}},I_{\{1,1,1,1,1,1,1,1,-4,0,-2\}}\}\,. \end{split} $

(4)

We also define operators $ {\hat m}^\circleddash $, which can generate a set of integrals as a union of integrals generated by $ \{\hat m^-, \widehat{m-1}^-,\cdots,\hat 0^-\} $ when acting on an integral.
As is well-known, the most complicated^② integrals in the amplitudes are those with the highest number of propagators, i.e., $ \nu_{i} = 1\; (i = 1,\cdots,8) $, and the highest numerator degree, i.e., $ -(\nu_9+\nu_{10}+\nu_{11}) $. By studying the two-loop five-gluon scattering amplitude diagram by diagram, we find that the highest numerator degree is 5 for all integrals. Therefore, we define an integral set

$ S_{(a)} = {\hat 5}^\circleddash I_{\{1,1,1,1,1,1,1,1,0,0,0\}}, $

(5)

which contains 3914 nonzero integrals with all of the most complicated integrals in five-gluon scattering amplitude being included. Because the five-gluon scattering amplitude is sufficiently general, all of the most complicated integrals (if not all integrals) belonging to topology (a) appearing in five-light-parton scattering amplitudes are included in set $ S_{(a)} $. In fact, for the two-loop five-gluon all-plus helicity amplitude, integrals in topology (a) form a subset of $ S_{(a)} $ [5]. Therefore, for the purpose of reducing integrals in physical amplitudes, the main job for topology (a) is to reduce integrals in set $ S_{(a)} $.
For topologies (b), (c), and (d) in Fig. 1, we define sets of target integrals $ S_{(b)} $, $ S_{(c)} $ and $ S_{(d)} $, similar to $ S_{(a)} $.

Before presenting our method for reducing two-loop five-light-parton integrals, let us first point out that for multiscale problems, expressing general integrals in terms of MIs explicitly is not preferred, even at the one-loop level. Instead, one usually sets up a system of block-triangular relations that can numerically relate all of the integrals to MIs (see [54] and references therein).
The advantage of a system of block-triangular relations over the explicit solution can be understood based on the integrals' singularities. If we express a complicated integral as a linear combination of simpler MIs, powers of Gram determinants will appear in the denominators of the coefficients of these MIs, which is necessary because only thus the linear combination of MIs can generate correct singularities of the target integral. Then, the numerators of these coefficients will have high mass dimensions and thus will have very long expressions. This difficulty can be nicely resolved using a system of block-triangular relations. Relations in each block can be very simple, but their solution can naturally generate Gram determinants in the denominator. Furthermore, correctly choosing the blocks may result in the solution involving only one Gram determinant.
Because reduction at the multiloop level is much more complicated than for the one-loop case, the above discussion implies that constructing a system of block-triangular relations may be the best way to reduce multiloop multiscale integrals. Unlike the one-loop case, where block-triangular systems can be achieved easily by analytically solving the IBP relations, block-triangular systems at the multiloop level are in general difficult to obtain.
In Ref. [52], based on our proposed series representation of Feynman integrals [52, 53] as input information, we constructed an algorithm that searched for block-triangular relations to reduce multiloop multiscale integrals. However, we found the method to be very time-consuming for physical problems, although it was efficient for reducing integrals with integrands containing only denominators. To deal with physical problems such as two-loop five-light-parton integrals, we propose here a two-step search strategy.
In the first step, we set up a system of relations that can numerically express all of the target integrals in terms of MIs. The system is allowed to be somewhat inefficient in numerical calculations; thus, the system is not required to be block-triangular. This system can be obtained either by using our series representation of Feynman integrals [52], or simply by using the well-known IBP system.
In the second step, we search for a system of block-triangular relations, which needs to be very efficient for numerical computations. The algorithm is the same as that proposed in Ref. [52] except that, instead of using our series representation of Feynman integrals, we use the numerical solution obtained in the first step as input information.
More details about the search strategy can be found in appendix.

To apply the above-proposed search strategy to physical problems, we still need to introduce the reduction scheme, which amounts to choosing target integrals and other integrals that are allowed to appear in each block. In this paper, integrals in each block are defined by operator $ {\hat m}^\circleddash $ acting on a proper integral. For example, to reduce the integrals in $ S_{(a)} $, all of the integrals are allowed to appear in the first block, and the target integrals in this block are all the 21 most complicated integrals in the top sector with degree 5. The first block enables us to express all the 21 most complicated integrals in terms of simpler integrals. Then, in the second block, we choose the most complicated integrals among the rest of the integrals as target integrals, and use operator $ {\hat m}^\circleddash $ acting on a proper integral to generate a set of integrals that covers all the target integrals. Then, the process is repeated. Eventually, any integral can be expressed in terms of simpler integrals.
Using the above method, we successfully determined systems of block-triangular relations for integrals in the four topologies in Fig. 1. The file sizes for all of these relations are acceptable, ~148 MB. To obtain these results required ~200 central processing unit (CPU) core hours to search for relations in the second step of the two-step search strategy, in addition to hundreds of CPU-core hours for generating input information by numerically solving the system obtained in the first step. Some basic information about these results is listed in Table 1.

top.	$\#$int.	$\#$MIs	$t_{\rm{search}}$/h	$t_{\rm{solve}}$/s	size/MB
(a)	$3914$	$108$	112	0.17	66
(b)	$3584$	$73$	31	0.090	40
(c)	$3458$	$61$	56	0.075	31
(d)	$2634$	$28$	8	0.035	11

Table1.Main information about the obtained reduction relations. $t_{\rm{search}}$ represents the CPU time required to search for these relations in the units of CPU-core hours. $t_{\rm{solve}}$ represents the time spent to solve these relations numerically using one CPU.

For more intuitive understanding, we show a matrix density plot for the block-triangular system of topology (a) in Fig. 2. This system contains 3914 integrals and 108 MIs, which means we need 3806 linear relations to reduce all of the target integrals. In this plot, each line represents a relation, each column corresponds to an integral, and black points represent nonzero elements in the matrix. Integrals are ordered, from the most complicated one to the simplest one, with MIs at the end of each line. The matrix is exactly block-triangular, and the largest block contains only tens of relations.
Figure2. Matrix density plot for relations of topology (a). Each row represents a relation and each column represents an integral. Black points represent nonzero elements.

Analytic expressions for all of these relations are available from the website in [57]. Technical details about our reduction scheme can be found in appendix.

Our final reduction relations have been verified numerically using an independent code $ \verb'FIRE6'$ [29] for randomly chosen phase space points, and the results of both approaches were in a good agreement.
For each given numerical point $ \epsilon $ and $ \vec{s} $, solving our reduction relations of the four families cost 0.4 s using one CPU, as is shown in Table 1. The time spent can be divided into two parts: assignment (substituting numerical $ \epsilon $ and $ \vec{s} $ into the system), which is proportional to the file size; and solving the system, which depends on both the number of relations and how these relations are coupled with each other. Because our systems are block-triangular, the time spent on the latter part is shorter. Therefore, the efficiency of numerical calculation of our reduction relations can be simply estimated by the file size.
Compared with explicit solutions, the file sizes of our reduction relations are much smaller. The file size for explicit solutions of eight-propagator integrals with degree up to 4 in topology (a), 26 integrals overall, is ~2 GB [48]; that for the explicit solutions of eight-propagator integrals with degree up to 4 in topology (b), 32 integrals overall, is ~0.8 GB [47]; and that for the solutions of all integrals in topology (c) is in excess of 20 GB for a compressed format [49]. It can be expected that our relations should be hundreds of times smaller than the complete explicit solution in terms of the file size, which results in more than a 100-fold speedup of numerical calculations, even if there is no memory deficit for storing the huge expressions of explicit solutions.
We note that the file size of trimmed IBP relations to reduce all of the integrals considered in this work is a few GB, which is also much larger than that of our reduction relations. The reason is that, although each IBP relation is simpler than ours, the IBP system involves hundreds of times more equations. Furthermore, the time spent on numerical IBP is dominated by the latter factor, because IBP relations are coupled in a complicated way. As a result, numerical IBP should be much more inefficient than our method. Through our test, numerical IBP via $ \verb"FiniteFlow"$ [35] combined with $ \verb"LiteRed"$ [34] costs about 2 min for each phase-space point, which is slower than our method by more than a 100-fold.
The above comparison reveals the advantage of our method. Numerical evaluation of explicit solutions spends too much time on assignments; while numerical IBP spends too much time on solving linear equations. Our method performs better on both parts, and therefore it is much more efficient. Similar to numerical evaluations over the field of prime numbers, our reduction relations should also be much more efficient for numerical evaluations with floating numbers, which enables phase-space integration to obtain physical cross-sections.

In this paper, we achieved the reduction of a set of integrals which covers all of the most complicated integrals in two-loop five-light-parton scattering amplitudes. Our results are expressed as systems of linear relations in the block-triangular form, which are very efficient for numerical calculations. The remaining integrals involved in amplitudes can be easily reduced using the same method, on demand. Therefore, a complete reduction of integrals in two-loop five-light-parton scattering amplitudes, which challenges all other methods, is available now. Because MIs are already known [20-24], our results provide the complete calculation of two-loop five-light-parton scattering amplitudes, and thus complete NNLO calculation of three light particles or jet-production at the LHC on the horizon.
To obtain the block-triangular relations, we developed the method in Ref. [52] by proposing a two-step search strategy along with a reduction scheme. As our newly developed method is general and efficient, other more complicated problems, like two-loop integrals for $ t\bar{t}+{\rm{jet}} $, $ t\bar{t}H $, or 4-jet hadron production, are also within reach. Our work opens the door for complete NNLO QCD calculations for production of three or more particles at the LHC.
In the current application of our method, most CPU time is allocated to solving the system obtained in the first step. Although the time spent is tolerable for the current problem, improvement may be needed for more complicated applications. There are different options. Using the method in [52], better integral sets can be explored. Another possible choice is to use trimmed IBP systems obtained by solving the Syzygy equations [44-48]. These possibilities will be addressed in future studies.
We thank K.T. Chao, F. Feng, Q.J. Jin, Z. Li, X.H. Liu, H. Luo, C. Meng, J. Usovitsch and Y. Zhang for many useful communications and discussions.

2

-->

A.1.Search strategy: step one

We take integrals originated from topology (a) in the main text as an example for explaining the details of our technique.
We want to set up a set of relations, using which we can express all integrals in $ S_{(a)} $ in terms of MIs for any given phase-space point (rational numbers for both $ \vec{s} $ and $ \epsilon $), with coefficients calculated in the finite field of a 63-bit prime number. Although the IBP method [25-35] can do this, we would like to explain in the following that our method proposed in [52] may provide a better choice.
For each given integral $ I_{\vec{\nu}} $, called a seed, there are 12 IBP relations among the integral set

$\tag{A1} G_{\vec{\nu}}^{\rm{IBP}} = \{\hat{1}^+, \hat{1}^-\hat{1}^+\} I_{\vec{\nu}}. $

(6)

In addition, there are six relations owing to the Lorentz invariance [55], which can be interpreted as linear combinations of IBP relations from other seeds [56].
The above IBP relations can also be found out easily using the method proposed in [52]. To this end, we introduce a parameter $ \eta $ for all integrals in $ G_{\vec{\nu}}^{\rm{IBP}} $, and then search relations among them using input information from the series representation [52, 53]. Up to $ d_{\rm{max}} = 1 $, where $ d_{\rm{max}} $ is a half of the maximal value of the mass dimension for the coefficients of relations, we find at least 12 relations; while up to $ d_{\rm{max}} = 2 $ we find at least 12+6 relations. Because these relations are analytical in $ \eta $, we can take $ \eta\to0 $ directly and recover the aforementioned $ 12+6 $ IBP relations.
The advantage of our method in [52] is that it allows to search for relations among any set of integrals. As the simplest generalization of $ G_{\vec{\nu}}^{\rm{IBP}} $, we can define an integral set

$\tag{A2} G_{\vec{\nu}} = \{\hat{1}^+, \hat{1}^-\hat{1}^+,\hat{1}^-\} I_{\vec{\nu}}, $

(7)

and search relations among them. Up to $ d_{\rm{max}} = 2 $, there are typically 2 more relations besides 12+6 IBP relations for each seed. With more relations in hand, it is possible to select better relations to achieve a more efficient reduction. For example, our relations from all $ \frac{4}{8} $-type seeds can already reduce 15 out of all $ \frac{5}{8} $-type integrals to integrals with lower degree (these relations are available at [57]). IBP relations from these seeds cannot achieve this, because $ \frac{5}{8} $-type integrals do not show up.
One can certainly explore other integral sets for each seed, to further improve the reduction efficiency. We did not do that because the efficiency of either the IBP set (6) or the generalized set (7) is sufficient for us to deal with the problem in this work.
With integral sets in hand, we generate a system of linear equations from all seeds belonging to $ \frac{m}{n} $-type with $ 3\leqslant n \leqslant 8 $ and $ 0\leqslant m \leqslant 5 $, and use the package $ \verb"FiniteFlow"$ [35] to trim the system by removing redundant relations and solving the trimmed system numerically, which expresses all integrals in $ S_{(a)} $ as linear combinations of 108 MIs (after exploring symmetries among MIs using $\verb"LiteRed" $).
2 -->

A.2.Search strategy: step two

In this step, we search linear relations to reduce the given target integrals in $ G_1\subseteq S_{(a)} $ to simpler integrals in $ G_2\subseteq S_{(a)} $ (the reducibility can be tested numerically easily). Combining the reduction scheme that will be described in the next section, a block-triangular system can be finally obtained.
We first describe how to search linear relations among the integral set $ G: = \{I_1,\ldots,I_N\}\subseteq S_{(a)} $ of the form

$\tag{A3} \sum\limits_{i = 1}^N Q_i(\epsilon,\vec{s}\,) I_i(\epsilon,\vec{s}\,) = 0\,, $

(8)

where $ Q_i(\epsilon,\vec{s}\,) $ can be decomposed as

$\tag{A4} Q_i(\epsilon,\vec{s}\,) = \sum\limits_{\kappa = 0}^{\epsilon_{\rm{max}}} \sum\limits_{\vec{\lambda}\in\Omega_{d_i}} \tilde{Q}_i^{ \kappa\lambda_1\ldots \lambda_5} \,\epsilon^{\kappa}s_1^{\lambda_1}\cdots s_5^{\lambda_5}, $

(9)

where $ \epsilon_{\rm{max}} $ is the maximal power of $ \epsilon $ allowed to appear in the relation, $ \Omega_{d_i} = \{\vec{\lambda}\in\mathbb{N}^5|\,\lambda_1+\cdots+\lambda_5 = d_i\} $, $ d_i $ is a half of the mass dimension of $ Q_i $ which can be fixed by $ d_{\rm{max}}\equiv{\rm{max}}\{d_1,\cdots,d_N\} $, and $ \tilde{Q}_i^{ \kappa\lambda_1\ldots \lambda_5} $ are unknown rational numbers to be determined. It is crucial to point out that, for given $ \epsilon_{\rm{max}} $ and $ d_{\rm{max}} $, the number of unknowns is finite. Therefore, it can be determined by a finite number of constraints. As will be explained in the following, these unknowns can be determined by the result obtained in the first step.
Based on the system of equations in the first step, for a given numerical point $ \epsilon $ and $ \vec{s} $ every integral in G can be represented as an 108-dimensional vector, with elements being the projection onto MIs,

$\tag{A5} I_i = \{C_{i,1}, \cdots, C_{i,108}\}\,,\quad i = 1,\cdots,N\,. $

(10)

By inserting these numerical vectors into Eq. (8), we obtain a vector equation, which results in at most 108 independent constraints over the unknowns. By repeating the above procedure many times (at most several thousand in this work), a sufficient number of constraints can be obtained, for determining all of the unknowns. As the above values are actually calculated in the finite field of a given prime number, we still need to repeat the procedure for several different prime numbers (at most 15 in this work) and use the Chinese remainder theorem to reconstruct the real results of the unknowns. Finally, linear relations with given $ d_{\rm{max}} $ and $ \epsilon_{\rm{max}} $ are obtained.
To reduce $ G_1 $ to $ G_2 $, we just set $ G: = G_1\cup G_2 $ and search relations among G with different values of $ d_{\rm{max}} $ and $ \epsilon_{\rm{max}} $. For the purpose of the current work, we find it is sufficient to fix $ \epsilon_{\rm{max}} = 3 $. To find out simple relations, we follow the algorithm proposed in [52] by starting the search procedure with $ d_{\rm{max}} = 0 $ and increasing $ d_{\rm{max}} $ by 1 each time, until enough relations are obtained to reduce $ G_1 $ to $ G_2 $.
2 -->

A.3.Reduction scheme

The reduction scheme determines which integrals should be involved in each block. We generate the integrals through previously defined operator $ {\hat m}^\circleddash $ acting on properly chosen integrals.
For example, in the first block for topology (a), we need to reduce the most complicated $ \frac{5}{8} $-type integrals. To this end, we set $ G: = S_{(a)} = {\hat 5}^\circleddash I_{\{1,1,1,1,1,1,1,1,0,0,0\}} $ with $ G_1 $ chosen as all 21 $ \frac{5}{8} $-type integrals. We indeed find out 21 independent relations, which can reduce all $ \frac{5}{8} $-type integrals to simpler integrals. The most complicated relation corresponds to $ d_{\rm{max}} = 7 $, which means that the coefficients of $ \frac{5}{8} $-type integrals are degree-2 polynomials in $ \vec{s} $. We then reduce $ \frac{4}{8} $-type integrals, which can be realized by setting $ G = {\hat 4}^\circleddash I_{\{1,1,1,1,1,1,1,1,0,0,0\}} $ with $ G_1 $ chosen as all 15 $ \frac{4}{8} $-type integrals. To reduce the rest of the top-sector integrals, we set $ G = {\hat 3}^\circleddash I_{\{1,1,1,1,1,1,1,1,0,0,0\}} $ with $ G_1 $ chosen as 11 top-sector integrals that are not MIs.
After reducing the top-sector integrals, we still need to reduce the integrals in the subsectors. For example, for the seven-propagator sector $ I_{\{1,1,1,1,1,1,1,0,0,0,0\}} $, whose most complicated integrals in $ S_{(a)} $ are of $ \frac{4}{7} $-type, we set $ G = {\hat 4}^\circleddash I_{\{1,1,1,1,1,1,1,0,0,0,0\}} $ with $ G_1 $ chosen as all 35 $ \frac{4}{7} $-type integrals in this sector.
Based on the above scheme, we obtain 3801 reduction relations. By introducing additional 5 symmetry relations among MIs, we have 3806 relations in total that can express 3914 integrals in $ S_{(a)} $ as linear combinations of 108 MIs.
We note that there is a way to further reduce the block size that has not been applied in this work. For example, by setting $ G = {\hat 3}^\circleddash I_{\{0,1,1,1,1,1,1,1,-1,0,0\}} $, we can generate smaller-size blocks to reduce a part of $ \frac{4}{7} $-type integrals.