AMFlow: a Mathematica package for Feynman integrals computation via Auxiliary Mass Flow

AMFlow is a Mathematica package to numerically compute dimensionally regularized Feynman integrals via the recently proposed auxiliary mass flow method. In this framework, integrals are treated as functions of an auxiliary mass parameter and their results can be obtained by constructing and solving differential systems with respect to this parameter, in an automatic way. The usage of this package is described in detail through an explicit example of double-box family involved in two-loop $t\bar{t}$ hadroproduction.

There are many methods on the market to compute master integrals, such as: sector decomposition [22,23,24,25,26,27,28]; Mellin-Barnes representation [29,30,31,32,33,34]; difference equations [2,35]; traditional differential equations [36,37,38,39,40,41,42,43,44,45], by setting up and solving differential equations satisfied by master integrals with respect to kinematic variables s; and others [46,47,48,49,50,51,52,53,54,55]. The sector decomposition method and Mellin-Barnes representation method can be applied in principle to any integral. However, it is well known that these methods, which need to calculate multidimensional integrations directly, are very inefficient to obtain high-precision results. Difference equations and differential equations can be very efficient, but they depend on integrals reduction to set up relevant equations, which may become very nontrivial for multiloop multiscale problems. Besides, usually there is no systematic way to obtain boundary conditions for these two methods.
The auxiliary mass flow method [56,57] is also a kind of differential equations method, which calculates Feynman integrals by setting up and solving differential equations with respect to an auxiliary mass term η. This method has many advantages. First, it is systematic, because boundary conditions at η → ∞ can be obtained iteratively [57,58]. Second, as only ordinary differential equations are involved, high-precision results can be efficiently obtained [59]. Third, integrals containing linear propagators and phase-space integrations can all be calculated [60,61]. Finally, integrals reduction to set up differential equations with respective to η is usually easier than to set up differential equations with respective to s [57]. Therefore, as long as reduction tools are powerful enough to set up differential equations with respect to η, auxiliary mass flow can always provide high-precision result efficiently.

in
Ref. [57]). It is thus valuable for high-precision phenomenological studies. This paper aims to provide a public implementation of this method, including the automation of the fully iterative strategy and a high-performance numerical solver for ordinary differential equations, so that it can be more widely used for phenomenological studies.

Auxiliary mass flow
In this section we give a review of the auxiliary mass flow method, concentrating on the computation of normal loop integrals [56,57]. The extensions to compute integrals containing linear propagators or phase-space integrations can be found in Refs. [60,61].

The plain method
Let us consider a dimensionally regularized Feynman integral family defined by where s is the list of all kinematic variables including Mandelstam variables and nonzero masses of particles, D = 4 − 2 is the spacetime dimension, L is the number of loops, i are loop momenta, D 1 , . . . , D K are inverse propagators, D K+1 , . . . , D N are irreducible scalar products introduced for completeness, ν 1 , . . . , ν K can be any integers, and ν K+1 , . . . , ν N can only be nonpositive integers. We next introduce an auxiliary integral family by inserting an auxiliary parameter η to each propagator of (1) Then physical results can be recovered by taking the following limit This auxiliary family, although seems to be more complicated than the original one, becomes rather simple as η approaches the infinity. This can be understood through region analysis [77,78]. More specifically, when |η| is very large, only the integration region with µ i ∼ O( √ η) can contribute, and thus every propagator can be expanded like where (ν) i ≡ Γ(ν + i)/Γ(ν) is the Pochhammer symbol. After all such kinds of expansion, what we get are combinations of equal-mass vacuum integrals, which have been intensively studied in literature [79,80,81,82,83,84]. As a result, auxiliary integrals I aux ( ν, s, , η) in the neighborhood of η = ∞ can be easily obtained and what remains is to perform analytic continuation (auxiliary mass flow ) of them to recover physical results. As auxiliary integrals can be expressed as linear combination of master integrals using integrals reduction, we only need to perform analytic continuation for master integrals, denoted by the vector I aux ( s, , η). Integrals reduction can also setup differential equations for master integrals, which look like For any fixed generic kinematic configuration s = s 0 1 , the above differential equations can be numerically solved by using series expansions, similar to numerically solving differential equations with respective to kinematic variables [85,86], which can realize the flow of η from the boundary at ∞ to physical value at i0 − .
Before describing how to solve the above differential equations, it is helpful to know some basic features of these auxiliary integrals as analytic functions of η. According to Cutkosky rules [87], integrals can be only real-valued on the real axis when η > η th , where η th is the largest threshold for the corresponding process. Thus the branch cut of the auxiliary integral can be defined as the straight line connecting η = −∞ and η = η th along the real axis, such that the Schwarz reflection principle holds everywhere except the branch cut (for real s and ). Now we can describe our strategy for analytic continuations, or solving differential equations. We first need to define a path for the analytic continuations connecting η = ∞ and η = i0 − , characterized by a list of regular points {η 0 , η 1 , . . . , η l } on which we will perform series expansions in order. A typical choice is shown in Fig. 1, where the larger (smaller) circle is defined as smallest (largest) circle centered at η = 0 that contains all singularities (no singularity) except η = ∞(η = 0). The choice of the regular points should satisfy the following rules: i) η 0 is outside of the larger circle; ii) η l is inside the smaller circle; iii) the distance between η i+1 and η i is smaller than the convergence radius of the series expansions centered at η i .
Then the flow of auxiliary mass can be divided into three main stages: i) expanding the integrals around η = ∞ and estimating at η = η 0 ; ii) expanding at η = η i and estimating at η = η i+1 for i = 0, . . . , l − 1; iii) expanding formally at η = 0 and matching at η = η l to determine the unknown coefficients in the formal asymptotic series. After these steps, we are able to take the limit η → i0 − for the expansion at η = 0 to obtain physical results.
A simple example would be helpful to explain the basic ideas of performing the expansions. For more technical details, see e.g. Refs. [56,59]. Let us consider a massless one-loop two-point integral family There is one master integral I(1, 1, ), whose result is In the aforementioned auxiliary mass flow method, we first introduce the auxiliary mass parameter to obtain Now there are two master integrals, Boundary condition for the first master integral can be computed fully analytically and the second one can be expanded near η = ∞ giving Next we define the list of regular points to perform expansions. We can read directly from the differential equations (11) that the singularities are 0, 1/4 and ∞. As a result, the list of regular points can be chosen as η 0 = −i/2, η 1 = −i/4 and η 2 = −i/8. As the first master integral in this example has been totally solved, we just consider the second one. Near η = ∞, this integral can be expanded like which is a natural generalization of its boundary condition (13). We then substitute the expansions (12) and (14) into the differential equations (11) and what comes out is a system of recurrence relations which can be used to express a n ( ) in terms of a 0 ( ), the boundary input determined by Eq. (13). Some of the results are given in the following all of which are real-valued, if is real. The expansion (14) enables us to estimate the value of I aux (1, 1, , η 0 ) through where has been used. Expansion near the regular point η = η 0 is a Taylor expansion, which looks like We again substitute this expansion along with the value of the first master integral into the differential equations (11) and obtain a system of recurrence relations, which can be used to reduce b n ( ) to b 0 ( ), the value of Then we can estimate I aux (1, 1, , η 1 ) using the expansion near η = η 0 (18) Similarly, we can expand near η = η 1 and obtain the estimation at η = η 2 At the last step, we need to consider the expansion near η = 0 and match at η = η 2 . The general form of this expansion is where the left part comes from the homogeneous equation and the right part comes from the inhomogeneous equation (sub-topology). By substituting the expansions (12) and (22) into the differential equations (11), we can obtain two sets of recurrence relations, which can be used to reduce all c n ( ) to c 0 ( ) and determine all d n ( ) respectively. For example, we have and We find there is actually only one unknown parameter, c 0 ( ), which can be determined through matching at η = η 2 . By substituting the estimation at η = η 2 (21) and the coefficients (23) and (24) into the series expansion (22), we can solve the resulting linear equation to obtain After computing these expansions, we can finally take the physical limit η → i0 − in the expansion near η = 0 (22). Note that in dimensional regularization, we have for any nonzero b. So what remains in this limit is just the leading term of the Taylor part, c 0 ( ), i.e., which agrees with the analytic result (8).

Iterative strategy
One interesting phenomenon from the previous example is that the number of master integrals increases after introducing η. As a result, it can be expected for much more complicated problems, the introduction of η may greatly increase the number of MIs, such that the differential equations (5) cannot be set up in reasonable time with current reduction techniques. To overcome this difficulty, in [57] we propose to apply the auxiliary mass flow method iteratively to reduce the number of master integrals, and thus the computational cost, to a reasonable level.
The key observation is that the number of master integrals can be reduced if η is introduced to fewer propagators. For example, for the two-loop fivepoint massless double-pentagon integral family with 108 master integrals shown in 2, we obtain 476 master integrals if η is introduced to all propagators ("all" mode), 319 master integrals for propagators 1-6 ("loop" mode), 233 master integrals for propagators 4-6 ("branch" mode), and the best case, 176 master integrals for the propagator 5 ("propagator" mode). For topologies where independent internal masses exist, we can do even better. We can simply treat these masses as η and thus will not introduce any extra mass scale ("mass" mode). Because "mass" and "propagator" mode introduce fewer extra number of master integrals than other modes in general, they usually perform better. However, as an expense, the boundary analysis is more complicated in general, due to more contributing integration regions as η → ∞. Following the general rules of region analysis [77,78] To obtain boundary conditions, we need to expand the integrands of master integrals in each region. Specifically, in the all-large region (L...L), each propagator should be expanded as where κ = 1 or 0, depending on whether η is introduced to this propagator or not. We thus obtain vacuum integrals 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, with propagators containing η contracted. In mixed regions, we need to 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. After the expansion, the part containing large loop momenta and the part containing small loop momenta are decoupled and we obtain factorized integrals. It turns out that usually the boundary integrals are still too complicated to evaluate directly. But this is still fine because they are already simpler than the original integrals, which means we can keep applying the previous procedure to simplify the boundary integrals until they are all known to us. For example, the double-pentagon topology can be simplified iteratively with "propagator" mode as shown in Fig. 3.
In practice, we would profit from a systematic definition of the terminal topologies. For example, we can always identify single-mass vacuum integrals as our terminals. In Ref. [58], single-mass vacuum integrals are further simplified in an iterative manner. In that way, we can simply identify the 0-loop integral (whose result is 1) as terminals, which have been proved to be more convenient.

Numerical fit
A very useful trick implemented in AMFlow is numerical fit. Consider a function f (x) which can be expanded near x = 0 as and our goal is to compute its estimation up to the k-th order with relative accuracy E n ≤ E, wheref n is the estimation of f n and E n is defined by We propose to realize this by evaluating f (x) numerically at some sample points x 0 , x 1 , . . . , x N (N ≥ k) near x = 0 and solving a system of linear equations In practice, we find two ways are useful to choose these sample points: 1. |x 0 | ∼ · · · ∼ |x N | ∼ r R, 2. x 0 , . . . , x N are distributed uniformly on the circle centered at x = 0 with radius r < R, where R is the convergence radius of the expansion (31). If one of these ways is chosen and the precision p of the samples f (x 0 ), f (x 1 ), . . . , f (x N ) is sufficiently high, then the relative accuracy of f n can be roughly estimated as It can be seen that the relative accuracy E n decreases as n increases. Thus to achieve our precision goal, we can set E k ∼ E, or equivalently This also gives a constraint about the precision p of the samples because we cannot expect a correct result if the precision of the samples is too low. So the total time consumption to obtain the estimation (32) is where t(p) is the average time needed to compute at a sample point with precision p, depending on both the nature of the problem and the numerical algorithm. Typically, in the framework of power series expansion method to solve differential equations of Feynman integrals, the dominant part of t(p) is a polynomial-like object of the number of correct digits, i.e., where α is a positive number. Therefore, we have which can be minimized by choosing and Next we can discuss how to apply this trick to the computation of Feynman integrals. To obtain numerical results of master integrals as expansions in , we can solve differential equations (5) with some numerical values of and solve a system of linear equations like Eq. (34) for each master integral 2 . We find the first way to choose sample points stated after Eq. (34) is better to use in this case, because we can always choose real values of to avoid potential complexities. Note that (41) and (42) only serve as a reference, and in practice one may need to make some adjustments to get satisfactory results.
This trick brings several benefits. First, a much simpler code structure is made possible, because in this framework all integrals are simply pure numbers rather than expansions in , which is much easier to carry out. Second, the problem of -order cancellations is totally resolved, because we never use truncated series in to express any integral throughout the calculations, which means we can always include more -orders by simply increasing the precision of the integrals. Finally, the computations at different sample points are totally independent and thus can be massively parallelized to save our waiting time.
This trick can also be applied to achieve asymptotic expansions of Feynman integrals at a given phase-space point or a given value of η. Sometimes, this becomes crucial, given the fact that there are usually many removable singularities in differential equations. With the second way to choose sample points on a circle, removable singularities inside the circle can be totally ignored.
Users can then follow the guidance outlined in README.md to install this package properly on their devices. After that, the package can be loaded by the command Get["/path/to/AMFlow.m"]; AMFlow depends on external programs to do integrals reduction. To use different reducers, one can set the following option SetReductionOptions["IBPReducer" -> reducer ]; where reducer can be any reducer whose interface with AMFlow has been built. Currently, three reducers based on Laporta's algorithm are available, including "FiniteFlow+LiteRed" [6,14], "FIRE+LiteRed" [6,17] and "Kira" [18]. Other reducers can also play their roles after users build their interfaces with AMFlow properly.

Input
First, we should use the function AMFlowInfo to define globally used objects during the computation, like AMFlowInfo[key ] = obj ; where key should be a string pre-defined in AMFlow and obj should be the corresponding object. We list most frequently used pre-defined strings and the meaning of their corresponding objects below: "Family" -the name of the integral family; "Loop" -a list of all loop momenta; "Leg" -a list of all external momenta; "Conservation" -a list of replacement rules for momentum conservation; "Replacement" -a list of complete replacement rules for scalar products among external legs; "Propagator" -a list of complete inverse propagators; "Numeric" -a list of replacement rules indicating the numerical kinematics where to perform the computation;

Automatic computation
AMFlow provides a function named SolveIntegrals to perform automatic computations of Feynman integrals. A general usage of this function should be like auto = SolveIntegrals[target, goal, epsorder]; where target is a list of target integrals, goal represents the precision goal and epsorder means the length of expansion in the final expansions, i.e., starting from −2L and ending at −2L+order with L the number of loops. This function will first reduce the target integrals to master integrals and then compute master integrals using auxiliary mass flow. The output auto is a list of replacement rules from integrals to their values.

Manual computation
Although SolveIntegrals is designed for most general purposes, there could be some extreme cases where this function may not be able to produce satisfactory results. So we introduce a more involved way to compute integrals in this section.
We first use the function GenerateNumericalConfig to regenerate the parameters for numerical evaluation suggested by SolveIntegrals {epslist, workingpre, xorder} = GenerateNumericalConfig[ goal, epsorder]; where goal and epsorder have been defined in the previous section. The output is a triblet: epslist is a list of suggested sample points of , workingpre is the suggested working precision and xorder is the suggested truncated order of the power series expansions. In principle, if these suggested parameters are used, we will obtain exactly the same results as SolveIntegrals. So, when SolveIntegrals fails to generate satisfactory results, users can define their own epslist, workingpre and xorder.
We then tell the program our preferred parameters by can be set to −4 in this example. The output exp is just the list of expansions in eps for target integrals.

Other functions
There are also other useful functions in AMFlow. Here we just give a brief summary. For more details, users can investigate corresponding examples provided in the folder examples.
See aotumatic_phasespace and feynman_prescription. 3. Computation of asymptotic expansions using the differential equations solver provided in AMFlow, either by traditional matching or numerical fit introduced in section 3. See differential_equation_solver.

Computation of integrals with complex kinematic parameters.
See complex_kinematics.

Computation of integrals in arbitrary space-time dimension.
See spacetime_dimension.

Summary of options
AMFlow allows users to set global options through SetAMFOptions, SetReductionOptions and SetReducerOptions. Here we list and describe the most frequently used options. For other options, we refer the users to the file options_summary.

Option and Default
Description SetAMFOptions "D0"→4 A rational number D 0 such that the integrals will be computed with D = D 0 − 2 . "WorkingPre"→100 Working precision when performing numerical computations, including solving differential equations and fitting. "XOrder"→100 Truncated order of expansions when solving differential equations. SetReductionOptions "IBPReducer"→ "FiniteFlow+LiteRed" Integration-by-parts reducer. Available reducers include "FiniteFlow+LiteRed", "Kira" and "FIRE+LiteRed". "BlackBoxRank"→3 Suggested maximal rank of seed integrals when constructing IBP systems. But if the maximal rank of target integrals s is larger than the value of this option, then the maximal rank of seed integrals will be adjusted to s internally. "BlackBoxDot"→0 Suggested maximal dot of seed integrals when constructing IBP systems. But if the maximal dot of target integrals r is larger than the value of this option, then the maximal dot of seed integrals will be adjusted to r internally. SetReducerOptions ("FiniteFlow+LiteRed" or "FIRE+LiteRed" used) "EMSymmetry"→False A parameter indicating whether symmetries among external legs should be exploited when preparing the topology using LiteRed. SetReducerOptions ("Kira" used) "IntegralOrder"→5 A positive integer ranging from 1 to 8 specifying the integral ordering for Kira. For more details, see Ref. [88]. "ReductionMode"→ "Kira" Reduction mode for Kira. Available modes include "Kira", "FireFly", "Mixed" and "NoFactorScan". See Ref. [18] for more details.

Summary and outlook
In this paper, the Mathematica package AMFlow is presented together with some explicit examples. We have highlighted the numerical fit strategy, which can overcome many difficulties when numerically solving differential equations.
The differential equations solver provided in AMFlow is of high performance and very suitable for high precision computations. In later version of AMFlow, we will provide some functions for users to access this solver in a more convenient way.
With the auxiliary mass flow method, integral reduction will be the only input for calculating Feynman integrals [58]. In the near future, a public implementation of the reduction method developed in Refs. [11,12] will be available, which can typically reduce the time consumption by 2 orders of magnitude comparing with other methods on the market. With this powerful reduction package and AMFlow, complicated integrals such as those in Ref. [57], can be computed automatically.