Toward Three-Loop Feynman Massive Diagram Calculations

: There are many methods of searching for traces of the so-called new physics in particle physics. One of them, and the main focus of this paper, is athe study of the Z -boson decay in e + e − collisions. An improvement in the precision of calculations of the Standard Model (SM) electroweak pseudo-observables, such as scattering asymmetries, effective weak mixing angles, and decay widths, related to the Z-boson will meet severe experimental requirements at the planned e + e − colliders and will increase the chance to detect non-standard effects in experimental analysis. To reach this goal, one has to calculate the next order of perturbative SM theory, namely three-loop Feynman integrals. We discuss the complexity of the problem, as well as the methods crucial for completing three-loop calculations. We show several numerical solutions for some three-loop Feynman integrals using sector decomposition, Mellin–Barnes (MB), and differential equation methods.


Introduction
The search for non-standard effects can be conducted in various ways, from which the study of Z-boson decay in e + e − collisions is discussed in this paper. This process was essential in the LEP era, leading to the precise knowledge of essential parts of the Standard Model (SM) [1,2]. A few successors of this experiment are planned: ILC [3,4], CLIC [5,6], CEPC [7,8] and FCC-ee [9][10][11][12][13]. Up to 5 × 10 12 Z-boson decays are planned to be observed at the Z-boson resonance with the FCC-ee collider, which is about six orders of magnitude more than at LEP. This will lead to very accurate experimental measurements of the socalled Electro-Weak Pseudo-Observables (EWPOs), if the systematic experimental errors can be kept appropriately small.
This level of accuracy is sensitive to virtual beyond-the-standard-model (BSM) effects from particles with multi-TeV masses and/or very feeble interactions. The interpretations of EWPOs in terms of BSM physics requires accurate predictions of these observables within the SM. To match the precision offered by FCC-ee, at least one additional order of perturbation theory has to be included beyond what is known today. For Z-boson decays, the complete two-loop SM corrections are currently known [14][15][16], whereas three-loop corrections will be required for the FCC-ee program.
The three-loop challenge is very demanding since no universal methods exist for the calculations. The methods and corresponding tools have to be further developed so that all the contributing Feynman integrals can be efficiently calculated. A subset of leading fermionic third order corrections to electroweak and mixed EW-QCD observables was recently accomplished in [17,18]; however, these corrections did not require the evaluation of genuine three-loop integrals.

Progress in Three-Loop Calculations for Z-Boson Observables
To match the accuracy offered by the upcoming FCC-ee collider for the Z-boson observables, three-loop corrections have to be included. This is a tough challenge from a mathematical and computational point of view. At first, self-energy integrals are discussed. The goal is to calculate all of the integrals up to eight-digit accuracy where possible, aiming for 10 −5 accuracy of the whole result of the combined contributions.

Three Loop Z-Boson Self-Energy Integrals
The initial step to address three-loop calculations is to compute self-energy contributions. Those integrals are, in general, simpler than the vertices, and therefore it is a natural follow-up after completing two-loop calculations. Lacking a universal method for numerical evaluation of loop integrals, we can track the strong and weak points of the available software, which leads to further development of these programs.
In total, there are 546 integrals that have to be calculated, spanning in between five to eight propagators, zero to three inverse propagators -leading to up to rank 6 tensor integrals and zero to seven massive propagators. The next step is to include three-loop Z-boson vertices.

Three Loop Z-Boson Vertices
For the case of vertex diagrams, the following initial classifications are introduced. Bosonic diagrams refer to diagrams with no fermionic loops inside, and EW-QCD diagrams are with/without gluons inside. No tadpoles or products of lower loop diagrams are included. One can easily see a huge increase in the number of diagrams between two and three loops from around two thousands to nearly half a million in the Z → bb process.
Next, the current status of the calculations is discussed for three-loop self-energy diagrams, which is the first step in our calculations toward solving vertex integrals connected with diagrams in Table 1. All calculations are performed for the SM masses scaled to the M Z mass, namely Additionally, for Euclidean region s = −M 2 Z and s = M 2 Z in case of Minkowskian kinematics. Recently, the numerical evaluation of planar-type three-loop self-energy integrals with arbitrary masses was discussed in [19]. Table 1. The number of two and three loop vertex topologies and diagrams for the Z → e + e − process. In brackets, we give the number of topologies in which those diagrams appear. QCD stands for diagrams with gluons, and EW is the remaining part. Bosonic means no fermionic loops inside. A similar table for the Z → bb process can be found in [11].

PySecDec
PySecDec works very well for multiscale Feynman integrals [20,21]. Results in Euclidean kinematics obtained with various integrators that were available in pySecDec for the three-loop self-energy diagram in Figure 1 are presented in Table 2. The Monte-Carlo integrators are used, since they typically offer better accuracy for the results of high dimensional integrals compared to deterministic algorithms with the same number of calculation points. As one can see from Table 2, even in the Euclidean case, a comparison of different integrators and methods is needed for precise accuracy control.
In this and the following numerical examples, all powers of propagators are equal to 1. This and the next figures were generated using PlanarityTest [22], which is available at [23]. Table 2. Comparison of three of the integrators available in pySecDec, each evaluated using 10 7 points for the integral corresponding to the diagram shown in Figure 1 in Euclidean kinematics, s = −1. The MB result is also given; for details, see Section 3.3.

Integrator
Result Absolute Error One can also use pySecDec for the calculation of integrals below the threshold (the contour deformation option in pySecDec is not needed in this case) in Minkowskian kinematics. If we aim at 10 −5 accuracy for the combined result, in the case of O(10 3 ) three-loop self-energy diagrams, a rough estimation is that 10 −8 accuracy for single integrals is needed. Equation (2) provides a result for such an individual integral calculated with the QMC integrator and 10 7 points, s = 1 0.0255102498599 ± 6.23 × 10 −10 . ( A corresponding diagram is shown in Figure 2. Unfortunately, it appears much worse in the case of integrals with one or two internal masses, and above the threshold, even without the highest number of propagators. For the diagram with seven propagators in Figure 3, the result calculated with 10 7 points, s = 1, is Even in this case, where we have seven propagators, we already lack accuracy, and, the case of scalar integrals with eight propagators using the same number of points leads to, at best, ∼5% accuracy. These results can eventually serve as a cross-check result for other methods.
In the case of tensor integrals, for the ones with loop momenta not only in the denominator but also in the numerator, a similar discussion can be held. The complexity of tensor integrals lies in their rank and, thereby, in how many loop momenta appear in the numerator. The more of loop momenta, the more complex integral. In the case of pySecDec, the bottleneck is the available RAM memory. Memory usage grows very fast with the increased rank of tensor integrals. Most of the integrals below the threshold can be calculated even with rank 6 tensors.
However, when going above the threshold, even most of the rank 4 integrals are problematic and not always possible to calculate with the available computational resources. Concerning accuracy, there is not much difference from the scalar case. Since many tensor integrals cannot be calculated with pySecDec up to the desired accuracy, or even at all, one needs another approach. One such approach is the well-known integration by parts (IBP) method, which is further discussed.

Differential Equations and IBP Methods
In the context of the well known integration-by-parts (IBP) reductions [25] and differential equations method [26][27][28][29] for Feynman integrals, we present an improved strategy that is capable of computing an arbitrary three-loop self-energy diagram, which is needed for this project. The IBP reductions with Kira [30,31] based on the Laporta algorithm [32] are computationally demanding if the number of loops and scales increases. The three-loop self-energy diagrams require several days for a successful reduction if we keep the relevant kinematic invariants and masses generic.
If we set all masses to physical numerical values and keep only the external kinematics in symbolic form, this reduces the reduction time to few minutes. IBP identities help us to relate a large set of Feynman integrals to a smaller basis of integrals, which are typically called 'master integrals'. We may derive a linear system of differential equations for the master integrals, from which we can solve the master integrals. In order to solve the master integrals from their differential equations, we provide a numerical set of boundary values obtained through pySecDec as an expansion in the dimensional regulator = (4 − D)/2.
With the help of series expansions of the system of differential equations implemented in DiffExp [33], we may transport the boundary terms to any other point in phase-space. In our approach, we choose the kinematic point at which we compute boundary values such that the numerical evaluation of pySecDec is stable and fast, and then we transport it to the physical region using DiffExp. First, we use the method of analytic regularization that was introduced in [34][35][36][37] to change to a finite basis of Feynman integrals. The resulting integrals typically have multiple dots and shifted dimensions. They can be related to the original basis by using IBP and momentum-shift identities. Secondly, we compute our boundary terms in the Euclidean region. Both simplify the numerical evaluation of master integrals due to the absence of contour deformation and the need to resolve the divergences of Feynman integrals with pySecDec.
Thus, we obtain stable and fast numeric convergence for all master integrals with pySecDec. Using DiffExp, we may transport these boundary values to physical points in phase-space, which works by iteratively expanding and solving the differential equations in terms of series expansions. The numerical accuracy of the results in Euclidean kinematics typically propagates to the Minkowskian kinematic point, which we may check by deriving two different sets of boundary conditions in the Euclidean region, and transporting them to the same physical point, after which the error is obtained by taking the difference in the resulting expressions.
The tensor integrals are a by-product of the method described above, since any tensor integral is a sum of master integrals. We are optimistic regarding the computations of the arbitrary three-loop self-energy diagrams needed for this project by applying the aforementioned method because the IBP reductions and the PySecDec numerical evaluations simplify a great deal, and these were the previous bottlenecks of three-loop computations.
We demonstrate the application of this method with the following example ( Figure 4): with p 2 1 = s. To define the Feynman integral we follow the same conventions as in [38]. For a shorter notation, we suppress the Feynman iδ prescription. The functions in the square brackets[...] are the so called inverse propagators. The exponents a i are the propagator powers, which we treat as integer numbers. In further discussion, we study I Merc [D, {a i }, s, M 2 W ] integrals with a i > 0, i = 1, . . . , 6 and a j ≤ 0, j = 7, 8, 9. With IBP and dimension raising identities, we manage to write these integrals as a linear combination of eight master integrals in different dimensions by exploiting [39][40][41] Here, we have omitted the explicit mass scales notation. The system of differential equations with respect to the integrals in masters and the integrals in the set masters are finite in the dimensional regulator parameter expansion = (4 − D)/2 around 0. We evaluate with PySecDec these integrals in the Euclidean point with masters [3] = (−0.111853 + 0.245297i) + (0.0033199 + 0.00359974i) 2

Mellin-Barnes Method
In this part, we will show perspectives of the MB method [42] for the calculation of three-loop integrals mentioned in Section 3. As was discussed in [14], the MB method works very well for one-scale problems at the two-loop level. In going to three-loop diagrams, the situation remains the same. A natural application of the MB method is a singlescale problem and, in the case discussed in this article, these are three-loop self-energy integrals with one or two equal masses inside. Planar diagrams corresponding to the problem have at most six-dimensional MB representations. Let us show several examples of this type.
The first example is the diagram in Figure 1, and the corresponding numerical result is presented in Table 2. The MB representation for this diagram is 5-dimensional and has the following form where we abbreviated z ijk = z i + z j + z k and omitted (2πi) factors as in the AMBRE package [42,43], available at [23], used for automatic derivation of this and the following representations. The numerical result in Table 2 was obtained using the MB package [44] (http:// projects.hepforge.org/mbtools/ (accessed on 30 May 2021)) with the Cuhre integrator [24] and 10 8 integration points. Increasing the number of points by one order did not increase the integration time. Due to the compact form of the integral, the integration time remained at the same level as for the Sector Decomposition (SD) method. The accuracy is also compatible with that obtained with the SD approach.
For the next example, we chose a diagram in Figure 5 with one less propagator. The MB representation for this diagram is 4-dimensional and is shown in Equation (15).
This example is very representative. After expansion in and simplifications by the barnesroutines.m package [45], we obtained only the 1-dimensional MB representation as shown in Equation (16).
The obtained representation can be easily integrated in both Euclidian and Minkowskian kinematic regimes. In Table 3, a numerical result calculated by the MBnumerics package [14,44,[46][47][48], which is available at [23], is compared with the SD result.  Table 3. A comparison of the MB and SD methods for the integral corresponding to the diagram in Figure 5 in Minkowskian kinematics with s = 1. The last example in this section represents one of the most complicated cases for considered one-scale planar self-energy three-loop electroweak diagrams. The MB representation for a diagram in Figure 6 in a form that comes directly from the AMBRE package is 7-dimensional.

Method
It is not difficult to check that, after the transformation z 5 → z 5 − z 3 , one can simplify this representation by integrating over z 3 with the help of the so-called first Barnes lemma (details about Barnes lemmas can be found, for example, in [42]). This also can be done by barnesroutines.m after expansion. In this way, the final representation is 6dimensional.
The result obtained by the MBnumerics in Minkowskian kinematics with s = 1 is The accuracy, in this case, is not high in comparison to the 4-dimensional 2-loop vertex cases considered in [14][15][16]. Nevertheless, it can be used for crosschecks and some preliminary results. The SD result for this integral has only an accuracy of about 5-10% with the setup discussed in Section 3.2.

Reduction of Scales Using Taylor Expansion
Since, for both IBP and MB methods, it is highly desirable to have as few scales in the integral as possible, one can think of a way to reduce their number. A simple way to do that is to Taylor-expand, as discussed and applied in more general contexts in [49,50]. We apply a Taylor-expansion to the diagram depicted in Figure 7. The integrand can be expanded at e.g., M W = M Z and the expanded propagator takes the following form  Figure 7. The diagram represents an exemplary integral for which the propagator with M W has been expanded as presented in Equation (19). An artificial integral is also considered where, instead of M W , the mass of propagator is equal to M H and then expanded in M H = M Z .
As the result, one can calculate more but, in principle, simpler integrals. This can be applied recursively in order to reduce the number of scales in the integral. Below are two exemplary plots for a three-loop Z boson self-energy integral in Figures 8 and 9 that were calculated using the Taylor expansion. The results show how the integrals converge when including more terms. Calculations were made in Minkowskian kinematics, and thus for s = 1. Each point represents the value of the integral after adding another order of expansion.  The first integral, with M W inside, converges faster than the other where, instead of M W , we have M H for two reasons. First, the difference between squares of M W and M Z is smaller than for M H and M Z . Second, the difference ∆ HZ ≡ M 2 H − M 2 Z is negative, and the terms in this series are positive for even powers of ∆ HZ and negative for odd powers, which makes it an alternating series. More explicitly, for the first integral including nine terms in expansion already leads to 10 −8 accuracy, while the other one reached 10 −4 accuracy with 21 terms taken into account.

Conclusions
The SD and MB methods served well at the two-loop level of EWPOs calculations [14][15][16]. In this work, we discussed the present capabilities of these methods with suitable three-loop examples. Preliminary studies based on the present implementations of the SD and MB methods in available packages showed that these methods are insufficient for the current challenges. Then, we sketched another very promising approach based on differential equations, presenting the first three-loop examples that can push the calculations on.
The minimum set of master integrals required for the efficient evaluation of differential equations relies on catching all relations between the Feynman integrals. Therefore, we choose Kira which implements as well as PySecDec the so-called Alexey Pak algorithm to detect all relevant symmetry relations in all computations. With the method of Taylor expansion, we were able to enrich the number of relations between the Feynman integrals and, in this way, simplify the computations as well.
Since the studies at the three-loop level are dominantly numerical, it is highly desired to have a few independent methods available and, thereby, the possibility to cross-check the results. Merging all the presented methods, we hope to fulfil the theoretical requirements for future Z-boson factories.