Analytic auxiliary mass flow to compute master integrals in singular kinematics

: The computation of master integrals from their differential equations requires boundary values to be supplied by an independent method. These boundary values are often desired at singular kinematical points. We demonstrate how the auxiliary mass flow technique can be extended to compute the expansion coefficients of master integrals in a singular limit in an analytical manner, thereby providing these boundary conditions. To illustrate the application of the method, we re-compute the phase space integrals relevant to initial-final antenna functions at NNLO, now including higher-order terms in their ϵ - expansion in view of their application in third-order QCD corrections.


Introduction
Feynman integrals fulfil differential equations [1,2] in internal masses and external kinematical invariants.These differential equations can be derived at the basis of the integrand, employing integration-by-parts (IBP) reduction [3,4] to express the derivatives of the integral under consideration in terms of the integral itself plus related inhomogeneous terms.
Consequently, Feynman integrals can be computed from their differential equations by determining generic solutions that are subsequently matched to boundary values at some special values of masses and kinematics.The resulting particular solutions then produce the desired expressions for the Feynman integrals.The differential equation method [5,6] has been largely systematised and is now among the most powerful and widely used tools to compute Feynman integrals, especially in multi-loop, multi-scale settings.Variants of it are applied to obtain closed-form analytical expressions (often in terms of generalized polylogarithmic functions [7]) or numerical representations of the integrals.Fairly often, analytical and numerical aspects are combined in hybrid-type approaches.
The derivation of the differential equations for Feynman integrals can be largely automated and is potentially limited only by the complexity and size of the system of required IBP equations.Once the differential equations are established, they can at least be solved in a purely numerical manner [8] or through a step-wise series expansion algorithm [9][10][11], provided their boundary conditions are known.Some of these boundary conditions can usually be inferred from internal consistency conditions on the integrals, such as their analyticity behaviour.In most calculations of Feynman integrals with the differential equation method, these considerations are not sufficient to fix all required boundary conditions, thereby necessitating the evaluation of at least some integrals for specific kinematical limits by other methods.
Especially when fully analytical results are desired (as for example to study the transcendental structure of Feynman integrals and Feynman amplitudes), the analytical computation of the boundary conditions may become a limiting factor.
A very promising approach to compute Feynman integrals at fixed kinematical points is the Auxiliary Mass Flow Method (AMFlow, [9]), which employs differential equations in an auxiliary mass parameter to determine numerical results for Feynman integrals from high-order series expansions to any desired level of numerical precision.
In this paper, we propose a novel approach to analytically compute Feynman integrals at fixed kinematical points (both regular and exceptional), based on the idea of the AMFlow method of introducing an auxiliary mass parameter.Section 2 describes our method and highlights its essential computational steps.We illustrate the application of this method in Section 3 on the analytical computation of boundary values of infrared-divergent phase space integrals relevant to deeply inelastic scattering processes at NNLO, which we compute at their (exceptional) kinematical endpoint.Section 4 contains our conclusions and an outlook.

The Analytic Auxiliary Mass Flow Method
The analytic calculation of loop or phase space integrals in dimensional regularisation is often based on the differential equation method [1,2,5,6].In this method, differential equations in particle masses or kinematical invariants are derived at the integrand level.By matching their generic solutions to specific boundary conditions, results for the desired integrals are obtained.These boundary conditions are given typically by stating the integrals at specific values of the masses and invariants, obtained from consistency conditions or by direct evaluation.
Often there exist special limits in kinematic variables or masses, where the problem simplifies sufficiently in order to calculate the boundary conditions directly.For example, in the case of multi-loop forward scattering amplitudes relevant to deep-inelastic scattering the (kinematically unphysical) limit of vanishing parton momentum in Euclidian kinematics (x → ∞) can be implemented via a naive Taylor expansion in the momentum of the parton and leads to simple massless propagator integrals which are known up to four-loop order.This boundary condition is then used in the differential equations in x to obtain results in the physical region 0 < x ≤ 1 by analytic continuation.
For phase space integrals, the dependence on external parameters is realized through particle masses or by projecting onto a specific final state property such as the momentum fraction or rapidity of an identified final-state particle.Especially for phase space integrals, the analytic computation of a boundary condition in a specific kinematical point is often of comparable complexity to the computation of the full integral.In particular, those kinematical points that are associated with particularly simple expressions usually correspond to singular points of the phase space integral, since the external parameters often regulate the infrared behaviour of the integrals under consideration.To circumvent the need for a direct computation of boundary conditions for phase space integrals, they can be obtained in some cases indirectly by inclusive integration of the generic solution of the differential equations over a selected external parameter [12,13].This approach is however not broadly applicable, since the corresponding inclusive phase space integrals may again be of comparable complexity to the integrals under consideration.
The direct evaluation of boundary conditions for phase space integrals usually has to proceed on a process-by-process or even integral-by-integral basis, often requiring dedicated phase-space parametrizations [14] for each integral family.While this can be done for phasespace integrals with rather low multiplicity and to low orders in the dimensional regulator ϵ this procedure becomes more and more involved for higher multiplicities and expansion orders.
In the following, we describe an alternative approach to obtain analytic expressions for loop or phase space integrals in specific kinematical points, in particular in singular points.The approach is based on the idea of introducing an auxiliary scale η into the integrals under consideration, in such a way that the behaviour in a special limit, usually η → ∞, simplifies drastically.In this way one can solve the differential equations with respect to η, fix the constants of integration, and obtain the original integrals in the limit η → 0. This idea has recently been used broadly in the Auxiliary Mass Flow (AMFlow) method [9] for the numerical calculation of loop and phase space integrals.The approach of introducing an auxiliary mass was originally formulated in Ref. [15,16] and has been used subsequently in the literature to calculate in particular single-scale integrals which satisfy relatively straightforward differential equations, see e.g.Refs.[17][18][19].Its implementation in the form of the AMFlow method allows to address a much larger range of multi-scale problems in a largely automated manner, yielding results in high-precision numerical form.
The original formulation [9,16] of the method can be summarized as follows.One considers an integral family characterized by the set of denominators {D i } and the corresponding master integrals ⃗ M (ϵ, s ij ) which are functions of the space-time dimension d = 4 − 2ϵ and kinematic invariants s ij .This original integral family is modified by adding an auxiliary mass η to each propagator of the integral family (2.1) By adding the auxiliary mass η one obtains a new set of master integrals ⃗ M η (ϵ, s ij , η).Due to the increased complexity induced by the massive propagators the number of master integrals is usually larger than before.The new set of master integrals fulfils a differential equation in the auxiliary mass For these integrals, one can calculate the limit η ≫ s ij by using the large mass expansion [20] which is a special case of the expansion-by-regions method [21].Since all propagators are massive, only the region where all loop momenta l i scale proportional to the heavy mass (l i → ηl i ) contributes to the expansion.Therefore the integrals for η → ∞ reduce to fully massive tadpole integrals which are known analytically up to three-loop order [22][23][24] and to very high precision numerically up to five-loop order [25][26][27].One can then use the differential equation to compute a deep series solution around η → ∞, employing the η → ∞ behaviour to determine all constants of integration.Finally, one can transport the initial value from large values of η to the physical values of η → 0 by constructing consecutive overlapping series expansions and matching them at points where two neighbouring expansions converge.Typically, the limit η → 0 is not analytic, but also described by an asymptotic expansion, so one has to be careful when taking the limit η → 0. The drawback of this approach is that adding masses to all propagators increases the number of master integrals tremendously, which can render the reduction problem or the numerical solution of the master integrals unfeasible.In order to remedy this problem the method has been refined in Ref. [28,29].In the refined procedure, the auxiliary mass is not necessarily added to all propagators, but one chooses a smaller set of propagators to add the mass to.This allows to control the number of master integrals.The refinement comes with its own drawbacks, since the asymptotic expansion for η → ∞ becomes more involved with more than one region potentially contributing in this limit.Moreover, the limiting values at η → ∞ are not fully massive tadpoles anymore.The first point can still be handled algorithmically by the procedure of the large mass expansion, however the complexity increases because of the larger number of regions which have to be considered.The second point can be handled by applying the AMFlow method again on the η → ∞ integrals, therefore reducing the problem to fully massive tadpoles recursively [28,29].
The method has also been extended to handle phase-space integrals [30].Here the basic idea is to use reverse unitarity [31,32] to treat the phase-space integrals on the same footing as loop integrals by replacing phase-space constraints with cut propagators.When applying the AMFlow method, the auxiliary mass is not added to the cut propagators.In this way one can use the AMFlow method to recursively reduce the boundary integrals to precomputed and simple phase-space integrals.
The methods described above have been implemented in a public package AMFlow [33].Since the method increases the complexity of the integrals and therefore also of the reduction problem and the differential equations, an essential point of the implementation is that all kinematic variables s ij and the dimensional regulator ϵ are set to rational values right from the beginning.The ϵ-dependence of the result is subsequently reconstructed by fitting to an appropriate ansatz.While this is useful to speed up numerical calculations it has the drawback that the resulting integrals cannot be obtained analytically.It also restricts the application of the method to non-exceptional kinematical points, since possible singular regions of the asymptotic expansion implicitly vanish in dimensional regularization once the variable is set to the singular value and one can only recover the so-called hard region of the asymptotic expansion.However, Feynman integrals usually simplify considerably in singular limits, making them ideal candidates for the calculation of boundary conditions to the differential equation method.
To apply the auxiliary mass flow method to the computation of integrals in singu-lar points requires to retain the analytic dependence on the parameter that characterizes the singular limit.The resulting Analytic Auxiliary Mass Flow method (AAMFlow) is described in the following.
To make the discussion easier we consider an integral depending on one parameter y.The generic solution of the differential equation in y is known.To obtain the particular solution for the integral, we thus require the boundary condition, which we aim to calculate in the singular limit y → 0. After adding the auxiliary mass we have the differential equation for the master integrals (2.3) In the limit y → 0 we can use the reasoning of expansion-by-regions to find an expansion of the form where the superscript (l) enumerates the individual components of the vector of master integrals.There is only a finite number of d i to consider and d i = 0 corresponds to the so-called hard region.
The starting powers n i are specific for each region and master integral and can usually be determined beforehand or at least a lower bound can be extracted from physical arguments.We can now insert the ansatz (2.4) into the differential equation (2.3) and obtain differential equations for the expansion coefficients c Since we are only interested in the leading terms of the asymptotic expansion (i.e.j = 0) to find boundary conditions for the differential equation for general y, we only have to consider a rather small set of expansion coefficients which we denote with ⃗ c(η, ϵ).It is worth mentioning that the different branches of the asymptotic expansion, which are characterized by different d i , do not mix in the differential equation with respect to η.We therefore find simplified differential equations where the subscript i labels the factorizing branches of the differential equation.We observe that the differential equations for the coefficients ⃗ c i (η, ϵ) simplify since they do not depend on y anymore.However, some care is required in selecting those propagators that the auxiliary mass η is added to.The addition of η should not modify the singular behaviour of the integral in y → 0, since the expansions in η → 0 and y → 0 are used interchangeably and therefore have to commute.This commutativity can be validated at the level of the differential equations, which must not contain denominator factors that are finite for η → 0 or y → 0 but vanish for the simultaneous limit (η → 0 and y → 0).A targeted choice of propagators with auxiliary mass allows to avoid these pathological situations.
At this point, the expansion coefficients ⃗ c i (η, ϵ) fulfil a differential equation with respect to η and the constants of integration can be fixed in the limit η → ∞.Consequently, the y → 0 limit of the original master integrals can now be calculated numerically to high precision using the original AMFlow method.
However, the reduced complexity of the differential equation in η might allow an analytic treatment with the established techniques of multi-loop calculations.For example one can try to find a transformation to a canonical differential equation [6], i.e. a new basis ⃗ ci in which the differential equation (2.5) takes the form or to find solutions in terms of iterated integrals following the algorithms outlined in [34].
Once closed-form solutions for the ⃗ ci (η, ϵ) are found, the limit η → 0 has to be taken.As has been mentioned before, this limit is also described by an asymptotic expansion.Taking into account that each loop momentum can scale in a soft or a hard manner relative to η, we find the following ansatz for the limit η → 0 where L corresponds to the number of (un-cut or cut) loops.The physical boundary conditions that we aim for are given by the hard region (ℓ = 0).One way to obtain these physical results is therefore to expand the solution to sufficiently high orders in the dimensional regulator and use the knowledge of terms proportional to powers of log(η) to disentangle the regions.
For the two loop integrals considered in this paper, the ansatz for the η → 0 limit is explicitly given by where the expansion coefficients d (l) ij,ℓ (ϵ) still have a Laurent expansion in the regulator ϵ.We indicate their ϵ−expansion as: (2.9) At each order in ϵ the coefficients c leaving implicit the dependence of r k,m (m ≥ 0) on i, j, ℓ.We have and r k,m (m ≥ 1) representing the coefficients of the log(η) terms in the expansion.It is crucial to notice that once the integral has been solved, all these coefficients are known, and it is trivial to extract each r i,k .Our goal is to obtain the expression for the hard region, specifically each d (l),k ij,0 .We demonstrate how this can be done by looking at the first coefficient, d (l),0 ij,0 , assuming without loss of generality that the hard-region expansion starts at order O(1).We can extract d (l),0 ij,0 by examining the analytic structure of the ϵ-expansion of the c (l) ij (η, ϵ), recognizing which terms contribute to it.
We expand the ansatz in ϵ and we obtain up to O(ϵ 3 ) (2.12) The O(1) term is given by the sum of the leading coefficients of the three regions.To extract the hard-region, our strategy will consist in subtracting from the O(1) term the sum of the leading term of the η −ϵ and η −2ϵ regions, respectively d (l),0 ij,1 and d (l),0 ij,2 .These two terms can in turn be determined by the coefficients of the log(η) terms in the higher order terms of the expansion.We can see that at O(ϵ) and O(ϵ 2 ) we have that the coefficient of log(η) and of log 2 (η) are given by (2.13) This linear system of equations can be solved in terms of d and then subtracting their sum from the coefficient of order ϵ k log(η) 0 .
In the following sections, we will illustrate the application of this method on the analytical computation of phase space integrals relevant to double real radiation and to single real radiation at one loop, both in the kinematics of deeply inelastic scattering (DIS).Integrals of this type were first considered in the context of the NNLO QCD corrections to DIS structure functions [35].In the context of the development of the antenna subtraction method [36][37][38] at NNLO, the full set of phase space master integrals was computed [39] through to finite terms in ϵ, yielding the integrated NNLO antenna functions in initial-final kinematics.In view of extending the antenna subtraction method to N3LO [40][41][42], we re-compute these integrals to two more orders in their expansion in ϵ, thus providing the integrated NNLO antenna functions to order ϵ 2 .
The phase space master integrals depend on the kinematical scaling parameter z, and we compute them from their differential equations.The method outlined above is applied to find the boundary values in the singular soft limit z → 1 systematically up to transcendental weight 6, as required for the subleading terms in the ϵ-expansion of the integrated antenna functions.

Example application: phase space master integrals
In the following, we will describe in detail how the analytic auxiliary mass flow method (AAMFlow) can be applied to determine the analytic values of NNLO phase space master integrals relevant to DIS in their singular kinematical endpoint z → 1.These values are then used as boundary conditions to differential equations for these master integrals, thereby enabling their analytic computation.

Kinematics and notation
The NNLO real radiation corrections to inclusive DIS contain contributions from double real radiation (RR) processes and from the one loop virtual corrections to single real radiation (real-virtual, RV) processes.Their kinematics is described by where and p 3 is only present for the RR processes.
The initial state kinematics is fully described by the dimensionful Q 2 and a dimensionless scaling variable z: related to the Mandelstam invariant s = (q 1 + q 2 ) 2 as In the following we fix Q 2 = 1, since dimensional analysis allows us to reconstruct the dependence on Q 2 at a later stage.The integrated initial-final NNLO antenna functions (and likewise, the DIS coefficient functions) are obtained by integrating fully inclusively over the final state phase space.The RR and RV integrands can be expressed in terms of a properly chosen set of linearly independent propagator factors D j as: and The vectors v j that constitute the propagator factors D j = v 2 j are constructed from linear combinations of the q i , p i and, in the RV case, of the loop momentum k.They form a complete linearly independent basis, allowing to express any scalar product that appears in the numerator of the integrand.
The measure dΠ n for the n-particle phase-space reads To calculate phase-space integrals we employ the reverse unitarity method [31,32], which allows us to map phase-space integrals to loop integrals with forward kinematics and cut propagators, following from the Cutkosky rules [43]: where we indicate the cut of a generic propagator D i as Using reverse unitarity, the RR-integrals therefore have the following structure On the other hand, the RV-integrals have the form , α j ∈ Z.
(3.11)With this setup, both RR and RV phase space integrals are mapped to cuts of two-loop four-point integrals in forward kinematics.
The map between phase-space integrals and cut-loop integrals allows to use wellestablished techniques for the evaluation of multi-loop integrals.In particular, we reduce any RR or RV integral to a linear combination of master integrals [39] by using integrationby-parts (IBP) identities [3,4,44].The master integrals are subsequently computed from their differential equations [1, 2, 5].

Integral families
The RR and RV phase-space integrals can be obtained as the physical three-particle and two-particle cuts of the two-loop forward scattering amplitude.The loop integrals for the latter were classified in [45], and we derive the RR and RV integral families by identifying all their topologically distinct three-particle and two-particle cuts.
All integrals can be expressed in terms of linearly independent subsets (integral families), each containing 7 elements, of the following 14 propagators: (3.12) All RR integrals must contain exactly 3 cut propagators.They can be mapped to the following four integral families: Figure 1 depicts the top sector integral (integral with three cut and four uncut propagators) of each sector, which can not have any irreducible scalar products in the numerator due to combinatorics.
The reduction of all integrals in a given family is performed using Kira [46].The list of RR master integrals, taking into account the symmetry relations between the families is where we used the notation I family list of uncut propagators , with superscripts in the list indicating propagators raised to higher powers, and I 0 corresponds to the integration of the three-  particle phase-space measure without any propagators [47]: where we choose the normalization All RV integrals must contain exactly 2 cut propagators.They can be mapped to the following three integral families, displayed in figure 2: The reduction of all integrals in a given family is again performed using Kira [46].The list of RV master integrals, taking into account the symmetry relations between the families is These bases differ from the one presented in [39].The explicit transformations are given by for the RR master integrals, and for the RV master integrals.
To compute the master integrals, we build a system of differential equations with respect to z for the master integrals and bring it into canonical form [6] using Fuchsia [48].We solve the system with iterated integrals, where we notice that the only letters appearing are the ones related to HPLs: The solution contains yet undetermined boundary conditions for all master integrals.

Required boundary conditions
A priori, one boundary condition per master integral must be determined from a calculation of the integral at some fixed value of z.However, for many of the master integrals, this boundary condition can be inferred the self-consistency of the solution of the differential equation in the limit z → 1.
Expanding the differential equations around z = 1, one finds that all RR master integrals must behave as where n i is a constant that can be determined for each integral and depends on the soft scaling of its propagators.In the case of the RV master integrals the loop momentum is not constrained to the soft scaling.Therefore, we find that in the limit z → 1 they can develop two possible branches For the RR master integrals we can factor out the uniform scaling behaviour and impose the cancellation of all the terms proportional to in the limit z → 1.This yields a system of equations for the boundary terms at z = 1, which leads to a substantial reduction of the unknown boundary terms.After solving this system, we find that the only RR integrals for which we need to compute the boundaries are We note that only two of them are top-sector integrals, from integral families B and C.
For the RV masters we separate the discussion into the two branches discussed in Eq. (3.23).Since the determination of the leading behaviour is more involved in the RV case, we only divide out the ϵ-dependent scaling in the limit z → 1 for the individual branches and demand the cancellation of all logarithms in 1 − z.For the first branch, we find that we have to compute the following boundaries and for the second branch We calculate the boundary conditions for (3.25) in the soft limit z → 1, by parametrizing the integrals as follows The boundary conditions to the differential equations are fully contained in the leading coefficients c 0 (ϵ), d 0 (ϵ) or e 0 (ϵ) for each integral.These will be computed below by employing the AAMFlow method.3.2 Boundary terms for RR master integrals

Set-up of the analytic AMFlow
Following the prescription presented in Sec. 2, we create new auxiliary integral families for the original families B and C by adding an auxiliary mass to certain propagators.We define the following auxiliary families, illustrated by their corner integrals in Figure 3: where we indicate the propagators with an auxiliary mass as We separately reduce each family to master integrals using Kira, without mapping the two via symmetry relations.We obtain 8 master integrals for family B aux and 5 for family C aux .These families are the only ones we need to calculate for the boundaries, since I A 1,7 is mapped via symmetry relations to I B 1,2 .For each auxiliary family, we build a system of differential equations with respect to the new variables u = 1/η 2 and y = 1 − z, as these are more convenient for fixing the constants of integration in the large mass limit η 2 → ∞ (u → 0).
The physical value of the master integrals around z = 1 is then obtained by solving the differential equations in u in terms of iterated integrals, retaining only the leading term in small y, followed by taking the limit of vanishing auxiliary mass η → 0 (u → ∞).This implies performing a change of variables from u back to η 2 , which includes an analytic continuation of the iterated integrals.

Family B
We set up a differential equation with respect to u = 1/η 2 for the integrals in family B aux .The master integrals are: Performing the same scaling analysis we did in Sect.3.1.3,we find the following behaviour in the limit y → 0 We need the leading term of O(y −1−2ϵ ) of the top sector to calculate the boundary condition for the physical integral I B 1,2,3,8 .Therefore we only need to solve the differential equation associated to the top sector up to this order in y.It turns out that at this order, the inhomogeneous term of the differential equation for B aux 1,2,3,8 receives contributions only from the integrals B 0 , B aux 1 , B aux 1,2 and B aux 2,3 .The differential equations for these integrals do however involve all other masters of the B aux topology.We give details on the calculation in the following.
The integrals B aux 1 and B aux 3/1 form a coupled system of differential equations: We decouple it into a second-order differential equation for B aux 1 , and an algebraic relation for B aux 3/1 .
solve the second order differential equation by inserting a symbolic ansatz: and solving for the c ij (u), up to integration constants.The large-mass limit of B aux 1 is given by and allows us to fix all the remaining constants of integration.The master integrals B aux 1,2 , B aux 3/1,2 , B aux 3 2 /1,2 fulfil a coupled 3 × 3 system of inhomogeous differential equations.We decouple this system to get a third order differential equation for B aux 1,2 using OreSys [49] (which is based on Sigma [50,51]).We solve this differential equation using HarmonicSums [52][53][54][55][56][57][58][59][60][61][62][63][64][65] and fix its constants of integration with the largemass expansion of B aux 1,2 The remaining integral B aux 2,3 fulfils a simple first-order inhomogeneous differential equation: We insert a symbolic ansatz for B aux 2,3 that respects its scaling in y (B aux 2,3 ∼ y 0−2ϵ ): and we solve the differential equations for the c i (u) order by order in ϵ.The remaining constants of integration are fixed from the large-mass behaviour We can now consider the differential equation for the top sector.From the scaling in y of the physical top sector, we know that it is proportional to y −1−2ϵ and therefore we require only that order in the expansion of its differential equation in order to determine the leading term of its y expansion.Taking into account the y-scaling of all subsector master integrals we then expand the top-sector differential equation for y → 0 and keeping only the term of O(y −1−2ϵ ): Then, we insert a symbolic ϵ-expansion for each integral and build a differential equation system for each ϵ-order into which the solutions for all the subsectors are inserted.The resulting differential equations are solved using HarmonicSums.The large-mass expansion of the top sector is given by which is a subsector integral not reducible to the phase-space.However, for the purpose of fixing the boundary values it is sufficient to impose that the solution scales as u 2 .
After fixing the constants of integration we find The result has been obtained up to and including terms of O ϵ 5 .The result is given by iterated integrals over the alphabet In the final step, the analytic continuation to η = 1/u → 0 has to be performed.The analytic continuation and subsequent expansion in η can be performed with routines implemented in HarmonicSums, e.g.we find Due to the analytic continuation, we are left with constants given by iterated integrals over the alphabet As long as only the first two letters of the alphabet in Eq. (3.49) contribute, the constants are well known and reduction tables up to high weight have been computed [54, Once the other letters of the alphabet contribute, we do not have such a complete reduction of constants.A practical way to deal with the new kind of constants appearing is to evaluate them with high precision utilizing for example ginac [7] and reconstruct the analytic structure using the PSLQ algorithm [67].This procedure will succeed since we know that for the physical branch in the limit z → 1 only multiple ζ values contribute.
Alternatively, the constants can also be reduced analytically by replacing the integer 2 in the alphabet (3.49) by a parameter t and performing a subsequent fibration to HPLs of argument t.The resulting HPLs at argument t = 2 are then evaluated in terms of known transcendental constants.
We find for the leading expansion term in the limit η → 0 where we again only give explicit results up to O(ϵ 0 ).We extract the hard region according to Eq. (2.14).As a first system of equations we find (where we suppress the global factor of S Γ π): with the solution It turns out that the branch with n = 2 does not contribute in general for this integral.
Proceeding like this also for the higher orders in ϵ, we finally obtain the physical boundary condition We notice that the boundary is of uniform transcendental weight.

Family C
For the auxiliary family C we find the following master integrals: Performing the same scaling analysis we did in Sect.3.1.3,we find the following behaviour in the limit y → 0 The differential equations in auxiliary family C are not coupled and we can simply solve them sequentially in a bottom-up approach.Starting from C aux 1 , for each subsector integral we solve the differential equation with respect to u and fix the constants of integration in the large mass limit.We notice that there is one subsector integral, C aux 2,6 , for which the large mass limit does not correspond to the phase space.Its constants of integration are however easily fixed by imposing the correct scaling in u of the solution.
Having the solution for all the subsectors, we can turn to the differential equation for the top sector C aux 1,2,6,9 , which is the one we actually need for the physical boundary conditions.Its differential equation reads To have a smooth workflow, we first substitute a symbolic ansatz for each master integral (except for I 0 ), where the leading order in y is determined in eq.(3.53), and expand the obtained system up to the y-order needed to fix the leading order of the top sector, which is of O(y −1−2ϵ ).We then insert the actual expression of the coefficients.The differential equation is solved using HarmonicSums and the algorithm presented in [34].The large-mass expansion for C aux 1,2,6,9 reads Also in this case all constants of integration could be fixed by requiring that the terms of order u i , i ≤ 1 vanish, without explicitly calculating C aux 6,9 .Following the same steps presented in Sect.3.2.2,we analytically continue our solution to the variable η 2 and take limit η 2 → 0. Extracting the hard-region of this limit allows us to calculate the first coefficient of the y expansion of I C 1,2,6,9 : thereby providing the boundary condition in the limit z → 1 for I C 1,2,6,9 .

Boundary terms for RV master integrals
The procedure outlined above for the boundary conditions to the RR master integrals can also be applied in the case of the RV master integrals.However, since masses are now also added to loop-momentum propagators, the asymptotic expansion for η → ∞ gets more involved.
We illustrate this here on the example of fixing the boundary conditions for the topsector RV master integral I G,(2) 1,2,8,9,10 .It is noted that the other RV master integrals I E 4,5 , I E 5,7 and I F 4,12 all correspond to one-loop bubble graphs, whose boundary conditions are obtained in a straightforward manner by direct integration.This example allows us to document the more complicated structures arising when applying our workflow to loop integrals.

Set-up of the analytic AMFlow
We add the auxiliary mass only to propagator D 2 , as depicted in Fig. 4.This minimal choice allows to limit the growth in the number of master integrals and makes the y and η 2 expansion independent.In the following, we indicate this auxiliary family as G aux : We find nine master integrals: We know that these integrals contain two branches in the variable y, given by the different scalings of the loop momentum k.We therefore can write each master integral as the sum of two contributions: We derive the differential equation for the master integrals with respect to u By inserting the ansatz (3.59), we split the system into two independent systems, one for each branch in the asymptotic expansion: We can now exploit the independence of the y expansion with respect to η 2 to simplify the differential equation for each branch.In fact, due to this independence, also the large mass limit (η 2 → ∞) of the master integrals presents an analogous structure of the form (3.59) and has the advantage to be more easily calculable.In particular, if one of the two terms G i , G i is equal to zero, then it is also zero for arbitrary values of η 2 , by the independence of the two expansions.This allows to nullify specific entries in M for each branch of y.As a result, our system now reads subsector integral G aux 8,9 is independent of the auxiliary mass.In particular, it consists of the integration over the 2 → 2 particle phase space of a massless bubble.We can easily calculate it in closed form: We notice that the only y branch contributing is y −2ϵ , so G aux 8,9 = G aux,(2) 8,9 . The subsector integral G aux,(2) 2,8,9 , satisfies the following differential equation that can be easily solved by inserting G aux 8,9 (3.63).The integration constants are fixed in the limit of large auxiliary mass.
As a general feature, when applying our workflow to one-loop integrals, we can distinguish two regions, where the loop momentum k scales either as large (L) as the auxiliary mass, k ∼ O(η), or small (S) with respect to it, k ∼ O(1).The sum of the two regions gives the full large-mass expansion: In this particular case, since we are interested in the y −2ϵ branch of the large-mass expansion, it is sufficient to examine only the small (S) region.The reason is the following: the integration over the 2 → 2 particle phase space produces a factor y −ϵ .In the large momentum region, all integrals over the 2 → 2 particle phase space contain massive tadpoles, which do not present any y dependence, and therefore can not give any extra factor of y −ϵ .For those integrals where it is non-vanishing, the large-mass expansion for the y −2ϵ branch must therefore arise from the small (S) region.
We find the following u → 0 behaviour for G

Conclusions
The auxiliary mass flow method [9] has proven to be an immensely useful tool to determine loop and phase space master integrals for specific values of masses and kinematical parameters.It relies on differential equations in an auxiliary mass parameter, which are solved through high-precision series expansions, yielding its answers in purely numerical form that can in principle be obtained to a high number of desired digits.The method is often used in combination with the differential equation method [1,2,5,6], where it serves to provide the boundary conditions that are required to obtain the particular solutions for master integrals from the generic solutions of the differential equations.
In this work, we took up the idea of auxiliary mass flow to devise an analytic procedure to compute master integrals at specific kinematical points.The resulting Analytic Auxiliary Mass Flow (AAMFlow) method enables us to retain full control on the asymptotic structure of the solution in the physical limit, where the method is used to determine boundary conditions to differential equations.This feature allows in particular to compute boundary conditions at exceptional kinematical points, where master integrals often assume a particularly simple structure in terms of transcendental constants.
To illustrate the new method, we computed phase space master integrals relevant to the NNLO corrections to deeply inelastic coefficient functions and to initial-final antenna functions to transcendental weight 6.Our results agree with previous derivations up to weight 4, and the novel terms will become relevant to extend the antenna subtraction method to N3LO.The analytic results are available in computer readable form as ancillary file to this publication.
By providing analytic expressions and allowing to address exceptional kinematical points, the AAMFlow method could in particular become relevant for few-scale loop or phase space master integrals, whose differential equations are amenable to the derivation of an analytical generic solution.By also providing the boundary terms in analytical form, the method will thus enable closed-form expressions for master integrals and amplitudes, which are ideally suited for mathematical investigations of their properties.
(l),0 ij,0 easily.The outlined procedure can be easily extended to fix any coefficient d (l),k ij,0 by first fixing the d (l),k ij,1 and d (l),k ij,2 by solving the system

Figure 1 :
Figure 1: Top sectors of integral families A,B,C,D.

Figure 2 :
Figure 2: Top sectors of integral families E,F,G.
Corner integral of B aux .(b) Corner integral of C aux .

Figure 3 :
Figure 3: Corner integrals of the auxiliary families.