Systematic approximation of multi-scale Feynman integrals

An algorithm for the systematic analytical approximation of multi-scale Feynman integrals is presented. The algorithm produces algebraic expressions as functions of the kinematical parameters and mass scales appearing in the Feynman integrals, allowing for fast numerical evaluation. The results are valid in all kinematical regions, both above and below thresholds, up to in principle arbitrary orders in the dimensional regulator. The scope of the algorithm is demonstrated by presenting results for selected two-loop three-point and four-point integrals with an internal mass scale that appear in the two-loop amplitudes for Higgs+jet production.


Introduction
Feynman integrals [1] are a fundamental constituent of perturbative calculations in theoretical particle physics and many techniques have been developed to calculate them. Going beyond one loop, the calculation of multi-scale multi-loop integrals is still a limiting factor in the theoretical predictions of hard processes.
To shorten the evaluation times the results have to be algebraic in the kinematical parameters (Mandelstam invariants, external and internal particle masses). Then the evaluation at each kinematic point takes just as long as the time needed for the insertion of their numerical values. Algebraic results can be obtained if the integrands are Taylor expanded in the Feynman parameters. To ensure that the approximation is fastly-converging, each integrand must be manipulated so that it is in a form optimised for a Taylor expansion.
To obtain such results in an algorithmic way, and thus find a compromise between analytic and numerical approaches, TayInt, a program to analytically approximate loop integrals, is presented in this paper. The backbone of the method is as follows: 1. Input an integral.
2. Reduce the integral to a quasi-finite basis introduced in Refs. [58,59], such that the divergences are in the coefficient of the simplest integrals. An automated script using the libraries of the publicly available program Reduze [59][60][61] performs this.
3. For those basis integrals that are not known in analytic form, carry out a decomposition into subsectors with smoother integrands. These are obtained using the publicly available program SecDec 3 [49][50][51][52], without its contour deformation option. The subsector integrands are analytic within the integration region, but may contain integrable singularities over thresholds and at end points. 4. Use a conformal mapping to move the singularities outside of the region of integration as far away as is possible. This is done in Mathematica [62].

(a)
To produce a result valid below the kinematic thresholds, the integrand is Taylor expanded, and integrated over the Feynman parameters. This is all done in FORM [63,64].
(b) To produce a result valid above thresholds, there is a separate algorithm implemented in Mathematica. The subsectors are first mapped onto the complex half plane. The algorithm then determines which configuration to use for each sector, that is, which contour orientation to use for the multiple variable integration and how to partition the subsequent region into smaller pieces. The Taylor expansion and integration are then performed on the new integrands specified by TayInt.
In Section 2, each step of the TayInt algorithm is described in more detail. Section 3 gives an analysis of the virtues of each step of the algorithm, while applications to integrals relevant for phenomenological applications are given in Section 4. Conclusions are drawn in Section 5.

The Algorithm
A generic Feynman loop integral G in an arbitrary number of dimensions D at loop level L with N propagators, wherein the propagators P j with mass m j can be raised to arbitrary powers ν j , has the momentum space representation where the q j are linear combinations of external momenta p i and loop momenta k l . The rank R of the integral is indicated by the number of loop momenta appearing in the numerator and the indices l i denote which of the L loop momenta belongs to which Lorentz index µ i . In what follows, R = 0 is taken for conciseness, although the TayInt algorithm is valid for arbitrary rank. The factor of iπ D 2 in κ in Eq. (2.2) is chosen by convention. The renormalisation scale is denoted by µ, which preserves the dimensionless nature of the coupling constant and is set to unity from here onward. The +iδ in Eq. (2.2) results from the solutions of the field equations in terms of causal Green functions.
For the calculation of unknown integrals, we rewrite every propagator in terms of Feynman parameters t j . After integrating the loop momenta, the general form of a scalar Feynman-parameterised multi-loop integral reads where the functions U and F are the first and second Symanzik polynomial, respectively, and are homogeneous in the Feynman parameters. U is positive semi-definite and F is negative semi-definite when all propagators are massless. In order to calculate loop integrals expressed in terms of Feynman parameters as rational functions of the kinematic parameters, the first universal step (U1) in the TayInt algorithm is to express the given Feynman integral G as a superposition of finite Feynman integrals G F multiplying factorised poles in . These finite integrals are either defined in a shifted number of dimensions about D = 4 − 2 , or have propagator powers greater than unity. The combination of these quasi-finite integrals which yields the original integral is found via integration-by-parts [65,66] and Lorentz invariance identities [6] and using the Laporta algorithm [67]. In practice, several scripts steer the performance of all necessary steps in the program Reduze [59,60] towards the generation of the quasi-finite basis. The divergences are restricted to the coefficients of the simpler integrals in the basis, so that the most complicated integral is always finite. Finding an optimal basis partially requires making an educated guess. The guiding principles are to express ultraviolet divergences in terms of vacuum integrals, and to relate subdivergences to sub-graphs of the original integral under consideration.
Once the original Feynman integral has a quasi-finite basis representation, the analytically unknown integrals in this basis are written in terms of their Feynman parametrisation and are then decomposed into subsectors which have smoother integrands. These subsector integrals are the building blocks of the rest of the calculation. Their improved smoothness is achieved using sector decomposition [40][41][42][43]. Thus, the second universal step (U2) in the TayInt algorithm is to perform the sector decomposition of the integrals G F in the quasi-finite basis by passing them to version 3 of the program SecDec [51]. Therein, the strategy G2, based on Ref. [68,69] and combined with the Cheng-Wu theorem [70,71] in Ref. [51], is used to yield sectors of the form where r is the number of sub sector integrals. A lj and B lj are numbers independent of the dimensional regulator . There are only N − 1 integrations in total because the first Feynman parameter t 1 is always integrated out with the δ-distribution. By construction, the deterministic algorithm results in integrands of the type where u t j and f β t j are polynomials in the Feynman parameters t j , and s 1 , s β ∈ {{q}, {m}} are kinematic invariants including masses. If the integral were not finite, the singular behaviour would now be contained entirely in the exponents A lj of Eq. (2.4).
Knowing that the integrals to be computed are finite, a sector decomposition might seem unnecessary. However, it is observed to be vital for an improved convergence of the series expansion.
The full integral can then be written in terms of its subsectors For calculations below the threshold, there still exist singular behaviours outside of the integration region. Thus the first below threshold step (BT1) in the TayInt algorithm is to maximise the distance to the nearest point of non-analyticity and so maximise the accuracy of the expansion. This is acheived by exporting the finite subsector integrands to Mathematica [62] and applying conformal mappings, For an integrand decomposed into r subsectors and containing J Feynman parameters the optimum mapping is, For the examples considered, we never mapped the final Feynman parameter in the first sector, as this parameter always appeared in the form (1 + t J ) in the denominators of the sectors. Thus, it is of no benefit to stretch the surface in that direction.
The second below threshold step (BT2) in the TayInt algorithm is to perform the Taylor expansion of the integrand, the relative simplicity of which is best suited to using FORM [63,64]. To this end, a FORM procedure for Taylor expanding functions of the form of the subsectors was written.
The third below threshold step (BT3) is to integrate the y j from y j (0) to y j (1), again in FORM. This yields results for Feynman integrals as rational functions of the kinematics valid everywhere below threshold. The precision is controlled by the order of the expansion.
However, in the kinematic region above the lowest mass threshold of a particular integral, the integrands contain discontinuities on the real axis which prevent a Taylor expansion from converging. Thus, the TayInt algorithm returns to the result of U2, specifically the subsector integrandsG F l . The Feynman +iδ prescription of Eq. (2.2) is then implemented in Mathematica. This is done by mapping the multivariate integrands of G F l onto complex half planes. The TayInt algorithm then determines the contour configuration which avoids the poles. More concretely, the over threshold part of the algorithm works as follows: 1. Implement the Feynman +iδ prescription of Eq. (2.2) by transforming the sector with j ∈ {1, . . . , J}. The mapping is chosen such that the real part stays between zero and one and the imaginary part parametrises a contour around a Landau pole.
2. Find the optimum contour configuration for each θ j with endpoints 0 and ±π, the combination of which is denoted by + or −, respectively. On this contour configuration, find the optimum variable θ * j to integrate exactly, if exact integration is possible.

Determine the optimum n-fold partitioning
with h n,j = ±π and l 1,j = 0, to use for the Taylor expansion of each sector integrand.
4. Perform the Taylor series expansion up to order p in each sector (2.12) and estimate the uncertainty of the full result by comparing the relative size of the contribution of the p-th order to the full Taylor series expansion, adding the contribution from each p-th order expanded sector in quadrature, .
Because of the use of exact one-fold integrations where possible, and the partitioning of the surface, the over threshold algorithm combines algebrauc and analytic manipulations, requiring flexibility. Therefore, it is implemented in Mathematica.
To elaborate, the first over threshold step (OT1) is to transform the Feynman parameters of the r subsectors, according to the rule, t j → 1 2 + 1 2 exp (iθ j ), and also generate a representative sample of the kinematic region in which the results are to be valid. This is a nested list of values for the β scales in the integral at γ points in the kinematic region within which we desire results, K = {{s 1 , . . . , s β } 1 , . . . , {s 1 , . . . , s β } γ } = {K 1 , . . . , K γ }. After that, the second over threshold step (OT2) uses the mean absolute value of the θ j derivatives (MAD:m l ), with the kinematic invariants set to the sample values, . (2.14) Note that the mean is also taken over the kinematic sample and the points for the θ j inserted along the surface, Θ As this is a J-fold surface, the third over threshold step (OT3) is to perform all possible one-fold integrations in the θ j exactly, without using an integrand expansion. If this is not possible, TayInt uses the Θ o(1),...,o(J) from OT2.
Next, the fourth over threshold step (OT4) is to determine the optimum post-integration surface, Θ o(1),...,o(J−1) . To achieve this, all the one fold exact integrations are performed and all resultant J − 1 variable integrands examined, using the two-step surface scanning process with the MAD. The Θ o(1),...,o(J−1) which minimises the MAD is selected. If the mean absolute derivatives of each of the possible J − 1 surfaces are, within a relative difference of 0.01, equally smooth, then we stop the process. This is because all the possible surfaces are then too similar and the potential for further optimisation is negligible. Provided there is another J variable contour which has a MAD within a relative difference of 0.01, the second best surface is taken and the procedure starts again. The MADs of each surface, with and without exact integration, are finally compared, and the surface with the overall minimum MAD is selected. In the vast majority of cases the post-integration surfaces are chosen, if an exact integration is possible.
The fifth above threshold step (OT5) determines the optimal way to partition the surface into sections P j within which the expansions are carried out. Using the MADs the relative size of the fluctuations in each section of the surface is calculated. A suitable partitioning is then chosen, i.e., the more fluctuations, the smaller that section. If no exact integration can be performed, TayInt uses two more orders and twice as many partitions to maintain the same degree of accuracy.
In summary, the over threshold part of the algorithm makes a complex mapping in several variables, determines the optimum contour configuration, variable to integrate exactly and partitions for the integrand expansion. The resulting integrands of each sector, G F l (θ 1 , . . . , θ J ) are then expanded and integrated in the sixth above threshold step (OT6), where e j are the expansion points, m j the order of the expansion and n j the number of partitions, in each parameter θ j . Adding up the results from each sector generates an expression for the full finite Feynman integral in terms of rational functions of the kinematic scales that are valid everywhere above thresholds in each invariant. Thus, systematic approximations are obtained for Feynman integrals with full kinematic dependence, valid in all kinematic regions, above and below mass thresholds. The precision of the approximation is controlled by the order of the expansion and the resolution of the partitioning. Finally, to estimate the uncertainty in the TayInt calculation, the truncation error in the Taylor expansion is calculated. To do this, the highest order contribution of all T F l , where s j = m j , are considered, summed in quadrature and divided by the full result, T F , where s j runs from 0 to m j . The resulting uncertainties for each sector are then added in quadrature to produce the final uncertainty estimate in the result. Note that the uncertainty will be overestimated in the vicinity of any kinematical point for which one of the sectors evaluates to a numerical zero by TayInt, meaning that theG F l (θ 1 , . . . , θ J ) is oscillatory. Nevertheless it always constitutes an overestimation of the uncertainty when less reliable. Moreover, the uncertainty estimate is always highly conservative as the pth order of the Taylor expansion is used to estimate the truncation errors in the results calculated using an expansion up to order p, rather than the order p + 1. U1: reduce the Feynman Integral to a quasi finite basis U2: perform a sector decomposition on the finite integrals in the basis below threshold above threshold A summary of all steps of the TayInt algorithm is given in Tab. 1.

Discourse on the Method
To facilitate a deeper understanding of the method, the different steps are illustrated for the integrals in Fig. 1, their characteristics being listed in Tab. 2. The finite sunrise S14 01220 [72] and triangle graph T41 [72] serve as examples to illustrate the TayInt algorithm. The integrals I10 [73], I21 [73] and I39 appear in the two-loop amplitudes [31,56,74] for Higgsplus-jet production in gluon fusion, mediated through a massive top quark loop. They demonstrate the applicability of TayInt to complicated multi-loop multi-scale integrals. The number of Feynman parameters quoted in Tab. 2 is counted after performing the integration of the δ-distribution of Eq. (2.3). In what follows and where given, the powers of each propagator are denoted by superscripts, i.e. X ijk . The kinematic invariants s, and u are defined as s = (p 1 + p 2 ) 2 and u = ( Figure 1. The two-loop finite sunrise S14 01220 and triangle T41, I10, I21 graphs for which analytical results are available [72,73]. The box-type integral I39, (e), is so far unknown analytically. Dashed lines indicate massless, solid internal lines massive, and dots squared propagators. Solid external lines denote, where indicated, massive and else off-shell particles.
As an example of the step U1 of the TayInt algorithm, the divergent sunrise integral S14 01110 [72] is written in terms of the finite integrals S14 01220 , S14 01320 and the tadpole S6 30300 , where the poles in can be seen as the (−4 + D) −1 terms.
As an example of the TayInt step U2, the O( 0 ) coefficient of the first subsector of S14 01220 reads S14 01220 and has two Feynman parameters, t 1 , t 2 and two scales, p 2 and m 2 . The O( 0 ) coefficient of the first and second subsector of the I10 [73] integral read They have three Feynman parameters and three kinematic scales, m 1 , m 2 and u.
The most precise way of computing the resulting subsector integrals was investigated by comparing the results obtained using various ways of expanding integrands to the literature for analytically known Feynman integrals. In particular, a comparison of expansions into Taylor series, geometric series, reverse Padé approximations, Chebyshev and Gegenbauer polynomials exposed the Taylor expansion as the most accurate given a variety of test cases. In Fig. 2, the relative difference between an order five Taylor expansion of, and the actual S14 01220 1 subsector integrand is plotted. While the expansion is rather accurate around the expansion point (t 1 , t 2 ) = ( 1 2 , 1 2 ), the differences between an ordinary Taylor expansion and the actual integrand can become large close to the edges of the integration region. In the case of S14 01220 1 , these amount to roughly 1%. The Taylor expansion turns the rational function R(t j ) into a polynomial P (t j ) which can be integrated analytically. Thus, it is necessary to check that accurate results for integrals can still be obtained after expanding the integrand. To this end, the example integral T41 was calculated by Taylor expanding the integrand to sixth order. In Fig. 3(a), the result is compared to the exact result of Ref. [72], with the truncation error of the Taylor expansion shown in the lower half of the plot. The Taylor-obtained result is plotted as a solid blue line in contrast to the literature result in a red dot dash line. The order six truncation error is plotted below using a lilac band. It can be observed that the combination of expansion and integration is effective for calculating Feynman integrals via  their subsectors for most values of For x ≈ 100, the discrepancy between approximated and exact result roughly reaches an unacceptable 9%.
To improve on the quality of the approximation, methods to maximise the distance to the nearest point of non-analyticity were investigated. The underlying rationale is that the convergence radius of a Taylor expansion is limited by the distance from the expansion point to the nearest point of non-analyticity. The latter are found in the region t j ∈ [0, −∞] for the subsectors. The use of conformal mappings in the t j , as detailed in Eqs.  integral I10 1 with two Feynman parameters set to zero, is shown. It has a point of non-analyticity outside the integration region, at t 0 = −1.
Applying the conformal mapping, on the integrand, as detailed in Eq. (2.9), the region of integration changes to y 0 ∈ [−1, − 1 2 ] and the non-analytic point is stretched to infinity, as shown in the lower plot of Fig. 4. The plot insets zoom into the actual region of integration.
The quantitative improvement can be seen by comparing Figs. 3(a) and 3(b). In Fig. 3(a), the order-six Taylor expanded result for the O( 0 ) coefficient of the T41 integral is compared to the exact result known from the literature [32]. In Fig. 3(b), a conformal mapping is applied to the T41 integrand before performing a Taylor expansion up to order six. For x ≈ 100, the discrepancy between approximated and exact result decreases to less than 3% when using a conformal mapping.
To view the effect of applying a conformal mapping to a more complicated example, the ratio of the result computed with an integrand Taylor expansion up to sixth order and the result computed numerically using SecDec is plotted for the O( 0 ) coefficient of the full integral I10, see Fig. 5. Error bars on the numerical results from SecDec are not plotted due to the high requested numerical accuracy of 10 −8 . The ratio is plotted over a kinematic range below threshold, the result with a conformal mapping shown in green and the one without mapping in blue. The average ratio between SecDec and the TayInt result is 1.00043 with the conformal mapping, and 1.00134 without it. Using the conformal mapping therefore increases the precision by more than a factor of three.  No further steps are required to calculate Feynman integrals below threshold. However, above threshold there are integrable singularities in the subsectors, within the region of integration, which renders all steps beyond the sector decomposition moot. By transforming to the complex plane (OT1), these singularities can be avoided. In Fig. 6 a slice of I10 2 is plotted in t 0 without, and in θ 0 with a complex mapping, respectively. In Fig. 6(a), the integrand without an analytical continuation to the complex plane contains a series of threshold singularities. The ridges are cut for better comparison with Fig. 6(b). The complex integrand in Fig. 6(b) shows a smooth behaviour everywhere in the integration region. This demonstrates how integrable poles in the physical region can be avoided and confirms that an ordinary series expansion of the integrand in the Feynman parameters t j is not sufficient to reproduce the actual result.
OT1 yields a version of the subsectors which can be expanded and integrated. How- ever, the optimum contour configuration from the possible Θ o(1),...o(J) surfaces must be determined. So must the optimum variable to integrate exactly, θ * j , if exact integration is possible. This is done in steps OT2-4. In Fig. 7, the integrand of I10 2 is plotted, with one integration performed exactly. Figs. 7(a) and 7(b) show that either choosing the wrong variable to integrate out exactly or choosing the wrong contour, leads to instabilities in the integrand. The optimal result is achieved when both are determined by TayInt, see Fig. 7(c). After the analytic continuation of the subsector integrands is optimised, large gradients along the edges of the complex surfaces, see e.g. the integrand surface along θ 0 = π of Fig. 7(c), are addressed in OT5. To maximise the precision of the result, the surfaces are partitioned and expansions performed around the central value of each partition. The Taylor expanded sections are then integrated and combined to yield a result for the entire sector. The rationale behind the partitioning is demonstrated in Fig. 8, where the one-dimensional integrand of Eq. (3.6) is shown for kinematic invariants set to arbitrarily chosen over threshold values m 2 = 781.249 GeV and m 1 = 173 GeV. In the upper plot, a discontinuity along the real line arising from the missing implementation of an analytic continuation of the integrand into the complex plane can be observed. To remedy this problem, the integrand is analytically continued into the complex plane, as described in Eq. (2.10) and demonstrated by showing the transformed real and imaginary part of the integrand in the lower plot of Fig. 8, in green and blue, respectively. Even though the integrand is now suited to a Taylor expansion, the gradient of the integrand is large for θ 0 → 0. Thus, the integration region is split according to Eq. (2.11), with the new integration boundaries (l, h) k,j marked by black lines in the bottom plot of Fig. 8. The expansion and integration is performed in each partition individually. The algorithm splits the integral such that within each partition the gradient is small. Hence, the convergence of the Taylor expansion within each integral piece is faster than that of an expansion of the whole integrand.
For generating results valid above threshold, one might ask why an ordinary Taylor   instructive to take a closer look at the first subsector of the integral I39, (3.8) All subsectors of I39 have four Feynman parameters and four scales, s, u, m 1 and m 2 .
Because of the complexity of these sectors after performing the first integration exactly, an expansion beyond (typically) tenth order in the Taylor series is not possible beyond O( 1 ) over threshold, due to the size of the algebraic expressions that are generated at intermediate stages. However, the result can still be made as precise as required by using more partitions. It is important to notice that the increase in the number of partitions allows the circumvention of memory bottlenecks. The expressions for each partition can be truncated at a lower order in the Taylor series than the full expression, owing to the smaller distance from the expansion point. Increasing the number of partitions does increase the algebraic computation time required to obtain the series expansion, which can however be parallelised trivially. Once the result is computed, an instant evaluation of arbitrary phase space points is possible. Thus, an increased partitioning enables the result to meet a target precision. For example, the I10 sectors have three Feynman parameters and three scales, and the Taylor series can be computed to beyond tenth order. Doubling the number of partitions at order six reduces the error in the real part by a factor of 11.85 and that of the imaginary part by a factor of 7.25.

Application to three-scale two-loop four-point integrals
To illustrate the power of TayInt, explicit results for the integrals I10 and I39 are presented below and above the threshold, and for different orders in , respectively. In Fig. 10, results for the finite I10 integral below threshold and up to O( 2 ) are shown. In the upper half of the plots, the approximated TayInt result is shown as a blue solid line, overlaid with results generated with the program SecDec, depicted as orange crosses. The only, though hardly noticeable, deviation can be seen directly on and around the threshold, at u = 4 m 2 1 , where the difference of the TayInt with respect to the exact result reaches at most 0.7%. In the lower half of the plots, the uncertainty band of TayInt is shown in lilac and can be compared to the uncertainties coming from SecDec shown in green. The SecDec results were computed using default numerical integration parameters and the integrator Vegas [75], asking for a relative accuracy of 10 −3 . The relative accuracy is adapted to the accuracy of the TayInt results. The same colour coding is used for all subsequent plots of results.
There is no appreciable precision loss as the order in increases. For I10 the mean SecDec TayInt ratios are 1.0006, 1.00039, 1.00021 at 0 , 1 and 2 , respectively. In Fig. 11, results for the finite I39 integral up to O( 2 ) are plotted in s below threshold, asking for a relative SecDec accuracy of 10 −3 . Again, the maximal deviation between the  difference reaches 0.1%, rounded up. There is no appreciable precision loss as the order in increases. In fact, quite the contrary, the mean SecDec TayInt ratios are 1.0003, 1.00017, 1.00009 at 0 , 1 , 2 respectively. In Fig. 12 the TayInt result for I10 at O( 0 ), obtained with a sixth order, four-fold partitioned, Taylor expansion, is plotted and compared to results from the program SecDec. The SecDec results were computed using default numerical integration parameters and the integrator Vegas, asking for a relative accuracy of 10 −4 . The plot shows the dependence on the scale u in the threshold region around u = 4 m 2 1 ∼ 120000 GeV 2 and above threshold. By referring to Fig. 10(a), a smooth transition from the below threshold expansion to the above threshold expansion can be observed. The size of the error directly on threshold is rooted in the fact that it displays a Landau singularity, a physical discontinuity, for which a Taylor series expansion has to break down by construction. Nonetheless, even on threshold, the relative accuracy of the TayInt result only drops to 10 −2 . To generate the SecDec results for Figs. 12 and 13 a relative accuracy of 10 −4 , and the integrator Vegas were chosen. Fig. 13 is a zoom into the region of larger u values. The TayInt truncation errors in the lower half of the plot are shown in yellow, using a four-fold partitioning, and in lilac, using an eight-fold partitioning. A strong increase in accuracy can be observed when the integrand is partitioned more often. More specifically, the relative truncation errors decrease from O(10 −4 ) with 4 partitions, to O(10 −5 ) when 8 partitions are used.  Close to threshold, there can be occasional rapid changes at the endpoints of integration regions which lead to larger truncation errors, for example increases of 10 % are observed for the O( 0 ) coefficient of I10. Beyond u = 30 m 2 1 in the over threshold region, the points of rapid fluctuation move closer to the boundary of the integration region, however do not enter it. This also leads to a small reduction in the precision of the TayInt results, however this loss is at the 0.01% level for the O( 0 ) coefficient of I10. This is illustrated in Fig. 14 by potting the first I10 sector near threshold, reasonably over threshold, and very far over the threshold in u.
In Fig. 15  integrands cannot be integrated exactly. However, the same level of accuracy with respect to the lower orders in is achieved by increasing the number of partitions used in OT5, with respect to the lower orders. To make the behaviour of the uncertainty estimates in the near threshold region clearer, more points are plotted in the u ∈ [4m 2 1 , 6m 2 1 ] region. It is important to observe that the O( 0 ) and O( 2 ) results for I10 of Fig. 12 and 16, each display a region (well above threshold) with an apparent deterioration of the TayInt error. This is not a result of the Taylor expansion having a smaller radius of convergence, as is the case near or well above threshold. Rather, this is a parametric effect that arises from two of the subsector integrands exhibiting oscillatory behaviour around zero in the vicinity of these kinematical points. These result in enhanced numerical cancellations among consecutive orders of the expansion. As a consequence, the uncertainty is grossly overestimated by the truncation error.
In Fig. 17 the TayInt approach is applied to the O( 0 ) coefficient of the integral I39, with more propagators and scales, without a loss in accuracy compared to the simpler examples discussed above. For the calculation an eight-fold partitioning and a fourth order Taylor expansion was used. The TayInt results agree well with the SecDec results, and are stable over the whole kinematic region, given the truncation error only ever varies by 0.01%. For the SecDec results a relative accuracy of 10 −3 was chosen. In Fig. 18  integration could be performed and a fourth order Taylor expansion with six partitions was used, and the TayInt results are consistent with the SecDec results across their full kinematic range. Relative to the complexity of the integrand, a small number of partitions was used, however no further precision was needed to obtain agreement between the TayInt result within the associated uncertainty. To put this into context, the TayInt result, 0.000934066 + 0.000126179i, at the ninth kinematic point at which SecDec is evaluated, the result of which is 0.000933183 + 0.000127422i, has an associated absolute uncertainty of 3.5495 · 10 −6 + 3.48459 · 10 −6 i. Thus there is agreement between TayInt and SecDec within the TayInt uncertainty. For the sake of comparison, the SecDec results for Fig. 18 were obtained asking for a relative accuracy of 10 −3 . In Fig. 19, the O( 2 ) coefficient of the integral I39 is calculated with TayInt. No exact integration could be performed and a fourth order Taylor expansion with eight partitions was used, and the algebraic TayInt result coincides with the SecDec results at the set of points considered above the 4m 2 1 threshold. The O( 2 ) SecDec results were obtained at a requested relative accuracy of 10 −3 . It is already apparent from Fig. 11 that the accuracy of the results obtained for I39 below s = 4m 2 1 is independent of the order in , and this is now shown to be true over the 4m 2 1 threshold as well for I39, an integral for which there is, thus far, no analytic result. This demonstrates the predictive power of TayInt, as the accuracy of its results is independent of the number of scales, the number of propagators, the kinematic  region or the order. To reinforce this, the mean difference between the TayInt and SecDec results for I10, I21 and I39 over threshold, ∆, is tabulated in Tab. 4. All TayInt results are based on a sixth-order, four partition, series expansion. The results for I10 are 792.8 MB in size, while for I39 they total 1.3488 GB. The file sizes refer to their unsimplified version. Using the plain Simplify command in Mathematica the size of the results can already be reduced by a factor of 1.33. Furthermore, removal of floating point zeros, 0., that appear in the results files leads to a reduction in size by a factor of 2.6 from the sizes quoted here. With the planned improved automation of TayInt, the simplification of the result files will also be addressed. One of the most attractive features is that the precision of the TayInt results is independent of the order. Nevertheless the computation time, both for the configuration determination and the actual calculation, increases significantly when going to higher orders in .
In Fig. 20 the truncation errors obtained using 8 and 16 partitions are plotted for the integral I39 at O( 0 ) in yellow and purple bands, respectively. The y-axis is truncated so that the larger ratios near threshold, due to the TayInt and SecDec results being consistent with zero, cannot be seen. However, given the size of the SecDec errors in the imaginary part, an inset provides a closer look at the TayInt error bands. The TayInt result is based on a fourth-order series expansion. The SecDec results have a relative accuracy of 10 −3 . The TayInt uncertainties decrease by an order of magnitude when the number of partitions is doubled, replicating the effect seen for the I10 integral. The average SecDec evaluation time at a phase-space point increases by a factor of 1.6 for each order of magnitude raise in relative precision, while the evaluation using the TayInt result is always instantaneous.

Conclusion
TayInt is a new algorithm which calculates systematic approximations to Feynman integrals algebraically in the kinematic invariants such that the results have validity in all kinematic regions and can be made arbitrarily precise.
The algorithm takes the propagators as input and works with subsector integrals generated by the program SecDec. The actual integration is facilitated via a Taylor expansion in the integration parameters. The accuracy is bolstered by conformal mappings and partitioning of the integrand before performing the Taylor expansion. The validity over threshold is ensured by performing a variable transformation which implements the correct analytical continuation of the integrand into the complex plane. Results can be obtained to higher orders in the dimensional regulator , both above and below mass thresholds.
We demonstrated the application of TayInt on two-loop three-point and four-point integrals with an internal mass, and used these examples to illustrate several features and virtues of the TayInt algorithm.
The TayInt algorithm expands upon the SecDec framework, and is in principle applicable to Feynman integrals with any number of loops and any number of kinematical scales. Its practical application could be limited by the algebraic complexity of intermediate expressions and final expansions. Based on the proof-of-principle applications considered here, we are confident that our expansion algorithm can be applied to many two-loop and three-loop problems of high phenomenological interest, where closed analytical expressions can not be obtained.