Multiloop corrections for collider processes using auxiliary mass flow

With a key improvement, the auxiliary mass flow method is now able to compute Feynman integrals encountered in cutting-edge collider processes. We have successfully applied it to compute some integrals involved in two-loop electroweak corrections to $e^+e^-\to HZ$, two-loop QCD corrections to $3j$, $W/Z/H+2j$, $t\bar{t}H$ and $4j$ production at hadron colliders, and three-loop QCD corrections to $t\bar{t}$ production at hadron colliders, all of which are crucial for precision frontier in collider physics in the following decade. Our results are important building blocks and benchmarks for future studies of these processes.

Introduction -With the good performance of the Large Hadron Collider (LHC), the particle physics now enters the era of precision. By comparing high-precision experimental measurements with theoretical predictions, we can deepen our understanding of particle physics Standard Model and probe signals of new physics. With the accumulation of experimental data at the LHC and proposals of future colliders, the precision of many experiments will surpass the current theoretical predictions [1,2]. Therefore, one of the most important tasks for theorists is to compute higher-order corrections in perturbation theory to reduce theoretical uncertainties for important processes, such as two-loop electroweak corrections to HZ production at e + e − colliders, two-loop QCD corrections to 3j, W/Z/H + 2j, ttH, 4j production and three-loop QCD corrections to tt production at hadron colliders, and so on, where j means jet.
To carry out these urgently needed perturbative calculations, an important part is to calculate corresponding Feynman master integrals (MIs), which form a complete basis of general Feynman integrals. Although there have been many methods on the market, such as traditional differential equations method which sets up and solves differential equations (DEs) with respect to kinematical variables s (denoted as s-DEs) [3][4][5][6][7][8][9][10][11][12], difference equations [13,14], sector decomposition [15][16][17][18][19][20][21] and Mellin-Barnes representation [22][23][24][25][26][27], computation of MIs in these cutting-edge processes is still very challenging. For s-DEs or difference equations method, one difficulty is that setting up these equations for very complicated processes becomes hard, and another difficulty is that there is no systematic way to obtain boundary conditions. While for the sector decomposition method or Mellin-Barnes representation method, numerical integration can be too inefficient for some processes and it is hard to get high precision.
In Ref. [28], we proposed the auxiliary mass flow (AMF) method, a special case of the differential equations method. In the AMF method, MIs can be calculated by setting up and solving differential equations with respect to an auxiliary mass term η (denoted as η-DEs) with almost trivial boundary conditions at η → ∞. The method is systematic and efficient as far as η-DEs can be obtained. However, for complicated processes, the introduction of η may greatly increase the number of MIs, such that the η-DEs cannot be set up in reasonable time with current reduction techniques [13,[29][30][31][32][33][34][35][36][37][38][39][40][41][42][43][44][45][46][47][48]. As a result, the application of this method is restricted to a large extent.
In this Letter, by introducing a key improvement to reduce the number of MIs, the AMF method becomes extremely powerful so that aforementioned cutting-edge problems are now solvable (see also Refs. [49][50][51][52][53][54][55][56] for recent developments on computations of MIs in these processes). Our method and results thus provide a valuable component for current and future high-precision phenomenological studies.
Auxiliary mass flow method with an iterative strategy -Let us first give a brief review of the original AMF method proposed in Ref. [28]. For a given integral family, we introduce an auxiliary mass term to all of its propagators, which results in an auxiliary family of integrals where D = 4−2 is the spacetime dimension, D 1 , . . . , D K are inverse propagators of the original family, and D K+1 , . . . , D N are irreducible scalar products introduced for completeness. Physical integrals are defined by The introduction of η brings a nice feature: when |η| is large enough, all integrals in the auxiliary family can be expanded as linear combinations of equal-mass vacuum integrals, which have been well studied [57][58][59][60][61][62]. Then physical MIs can be obtained from MIs of the auxiliary family, denoted as I aux (η), by the flow of η from ∞ to i0 − , which can be realized by setting up and solving a system of η-DEs, Unfortunately, the introduction of η usually significantly increases the number of MIs, which makes it hard to set up the above η-DEs for very complicated processes. Our improvement is based on a simple observation that the number of MIs can be reduced if we introduce η to fewer propagators. Taking the massless two-loop doublepentagon family in Fig. 1 as an example, which contains totally 108 MIs before introducing η, we summarize relevant information in Tab. I for introducing η in several different ways. Here, the four different modes, "all", "loop", "branch" or "propagator" correspond to introducing η to all propagators, propagators forming a closed loop, propagators on the same branch, or a single propagator, respectively. It is clear that the newly introduced modes can indeed result in fewer MIs than the original "all" mode. Furthermore, we find that the computational cost to construct and solve the system of η-DEs is usually positively correlated to the number of MIs, which is as expected. As a result, new modes can provide great improvement to the AMF method.
We note that, besides the above modes, sometimes a "mass" mode can be defined, which adds η to propagators with equal nonzero mass, with the condition that this mass is independent of scales from external momenta. As the "mass" mode does not increase the number of MIs, it should be used as far as possible, e.g., Refs. [63,64]. But for the current example, there is no massive propagator to enable the use of the "mass" mode. In a general mode, however, MIs of the auxiliary family cannot be simply decomposed to linear combinations of vacuum integrals near η = ∞, because there are usually more integration regions. Fortunately, when |η| is very large, all integration regions can be systematically identified, following the general rules of region analysis [65,66] where is a linear combination of loop momenta, and κ = 0 or 1, depending on whether η is introduced to this propagator. We thus obtain vacuum integrals after expansion in this region. In the all-small region (S...S), only propagators containing η should be expanded as In this case, we get integrals in a subfamily. In mixed regions, we decompose loop momentum of each propagator as the sum of a large part L and a small part S . Then, if L = 0 or κ = 0, we can expand the propagator as Otherwise, no expansion is needed. So, after the expansion in mixed regions, the part containing large loop momenta and the part containing small loop momenta are decoupled, which results in factorized integrals. It is now clear that region analysis for the "all" mode is just a special case, where only the region (L...L) contributes.
Although the obtained integrals after expansion are more complicated than those in the "all" mode, especially the integrals obtained in all-small region, it does not bring any trouble because they have already been simpler than original integrals. Then we can further introduce η to these obtained integrals and use η-DEs to push boundary conditions to even simpler integrals. After several steps of iteration, we can always cast boundary integrals to either scaleless integrals or single-scale vacuum integrals. The former ones can be set to zero in dimensional regularization directly, and the latter ones are simple enough to calculate. Actually, we can always use the "propagator" mode to cast boundary conditions to single-mass vacuum integrals, as shown in Fig. 2, which can be related to extensively studied p-integrals [67].
Eventually, we can numerically solve the systems of η-DEs, iteratively, to obtain physical MIs. We emphasize that the step to take η → i0 − in any mode is similar to that for the "all mode" discussed in details in Ref. [28]. Using η-DEs, we can get generalized power series expansion at η = 0, and then, to obtain the limit η → i0 − , we can safely discard all nonanalytical terms in η in the spirit of dimensional regularization. A pedagogical example -As a pedagogical example, we show how to calculate MIs in the massless two-loop double-pentagon family, as shown in Fig. 1, using AMF method with the new iterative strategy. Note that MIs of this family were successfully calculated only in the last two years [50] taking advantage of the powerful canonical form of DEs [9]. The family contains five external momenta {p 1 , . . . , p 5 } flowing outside, which satisfy the on-shell conditions p 2 i = 0 and momentum conservation 5 i=1 p i = 0. Thus there are five independent kinematic variables s = {s 12 , s 23 , s 34 , s 45 , s 15 } with s ij = (p i + p j ) 2 . We denote where the first eight are inverse propagators and the last three are irreducible scalar products. As shown in Tab. I, it should be a good choice to apply the "propagator" mode to first introduce η to D 5 . Then by expanding integrands of MIs of the auxiliary family in each integration region, we find only (LLL), (SLL) and (SSS) contribute, because all other regions are scaleless. In the region (LLL), all obtained integrals after expansion can be reduced to the two-loop single-mass vacuum MI in Fig. 2 directly. While in the region (SLL), after integral reduction, only a factorized integral remains Although this integral is simple enough to calculate, we can again introduce η to its nonvacuum part to cast it to a vacuum-integral problem. This step is quite useful in practice to make the computations more automatic.
In the most complicated region (SSS), we need to handle integrals in the subfamily obtained by contracting the fifth propagator. To calculate these integrals, we again apply the "propagator" mode to further introduce η to D 3 . The region analysis is the same as before. The only nontrivial one is still the (SSS) region, where a six-propagator subfamily emerges. Then, we can repeat the above steps to simplify the expanded integrals in the (SSS) region iteratively. We show this iteration in Fig. 3, where solid lines represent propagators with η introduced. We also present the number of original MIs along with that after introducing η under each diagram. The last step is a scaleless subfamily. So, after this iteration, all integrals after expansion around η → ∞ are known to us. while results of other integrals can be found in the ancillary file. We note that results with higher precision and with higher order in expansion can be achieved easily within our method. The coefficients of n with n ≤ 0 in Eq. (10) are found to agree with analytic results obtained in Ref. [50]. The coefficients of n with n > 0 are new, and they have passed a highly nontrivial self-consistency check. On the one hand, results at another phase space point s 1 can be obtained by solving a system of one-dimension s-DEs [11,55,69,70] along the line between s 0 and s 1 , with boundary conditions at s 0 . On the other hand, results at s 1 can also be calculated directly using the AMF method. We find full agreement between results obtained in the two ways.
Applications on cutting-edge problems -Now, we use the improved AMF method to calculate MIs encountered in some cutting-edge problems. As shown in Fig. 4, integral family (a) is relevant for two-loop electroweak corrections to e + e − → HZ, families (b)-(e) are relevant for two-loop QCD corrections to W/Z/H + 2j, H + 2j,ttH and 4j production at hadron colliders, respectively; and family (f) is relevant for three-loop QCD corrections to tt production at hadron colliders. All these processes are very important in the following decade. Among them, MIs of the family (a) could in principle be calculated using the method developed in Ref. [49]; while MIs of the rest topologies may be challenging for all methods on the market and have not been calculated yet as far as we know.
To calculate MIs in families (a) and (c), we choose the "mass" mode at the first step of the iteration, because there are masses in denominators that are independent of external scales. Other families are calculated using the "propagator" mode. checked with results in the literature [55,56]. Our results are provided in the supplemental material.
On the one hand, our results can serve as highprecision boundary conditions for s-DEs if fully analytic s-DEs can be constructed, and thus MIs at any given phase space points can be achieved. On the other hand, if s-DEs are hard to obtain, MIs at any given phase space point can be calculated directly using the AMF method. Therefore, our results are important building blocks and benchmarks for various future studies.
A few comments are in order: 1) Time consumptions for setting up and solving all η-DEs are summarized in Tab. II, along with percentages at the first step in the iteration. It is clear that the first step dominates the cost in any family, and the cost of all later steps are negligible. This phenomenon implies that the "propagator" mode is usually the best choice, if a "mass" mode does not exist, because it may results in the minimal number of MIs at the first step.
2) Differential equations are solved numerically using our Mathematica package AMFlow [68]. We find the number of correct digits is in proportion to the time consumption, due to the fact that η-DEs are rational functions of η. Furthermore, it can be expected that the efficiency can be significantly improved if a low-level language, like Fortran, is used. Therefore, the number of digits is rather cheap in the AMF method, and reconstruction of analytic results via the PSLQ algorithm (see Ref. [71] and references therein) might be possible for problems with not too many mass scales, which will be investigated in the  2 916 64 1305 30 1801 63   TABLE II. Time consumption to compute various families, where "dp" denotes the double-pentagon family in Fig. 1 and "(a)"-"(f)" denote corresponding families in Fig. 4. Tsetup denotes time consumption to set up η-DEs, including steps to construct the block-triangular linear systems, T solve denotes time consumption to numerically solve the η-DEs to obtain 16-digit precision results, P1 denotes the percentage to set up and solve η-DEs at the first step in the iteration, and T s denotes time consumption to set up one-dimension s-DEs along the line between two chosen phase space points. Time consumption is counted in the unit of CPU core hours.

future.
3) Integrals reduction is needed to set up DEs. We utilize the method developed in Refs. [38,39] to construct block-triangular linear systems among Feynman integrals, which are found to be much more efficient than traditional integration-by-parts systems constructed by LiteRed [33] and solved by FiniteFlow [41]. To set up η-DEs in this Letter, the block-triangular systems can roughly speed up by 2 orders of magnitude. 4) Time consumptions for constructing systems of onedimension s-DEs of physical MIs are also summarized in Tab. II. It looks first surprising that they are usually larger than the costs to set up η-DEs, although η-DEs may involve more MIs. This is understandable because the derivative of Feynman integrals with respect to η only increase the degree of denominators by 1, but that with respect to s increases not only the degree of denominators but also the degree of numerators, which results in a much larger reduction system that is much harder to solve.
Summary and outlook -We find that, by introducing the auxiliary mass η to fewer propagators at a time, the AMF method has a much better performance, and meanwhile the systematics of AMF method is untouched. We show that the AMF method can provide high-precision boundary conditions for traditional DEs as well as provide a highly nontrivial self-consistency check. Furthermore, an interesting observation is that it is usually easier to set up η-DEs comparing with s-DEs. Therefore, the AMF method may still work even if a problem is so difficult that s-DEs becomes hard to set up.
Equipped with the improved AMF method, we have successfully calculated some very difficult Feynman integrals encountered in two-loop electroweak corrections to e + e − → HZ, two-loop QCD corrections to 3j, W/Z/H + 2j, ttH and 4j production at hadron colliders, and threeloop QCD corrections to tt production at hadron colliders. As all these processes are crucial in the following decade, our results are important building blocks and benchmarks for future studies.
Our method has been implemented in our package AMFlow [68]. In the near future, with visible developments of integral reduction techniques, we expect that master integrals of all next-generation processes can be calculable using the AMF method.