Integral Reduction by Unitarity Method for Two-loop Amplitudes: A Case Study

In this paper, we generalize the unitarity method to two-loop diagrams and use it to discuss the integral bases of reduction. To test out method, we focus on the four-point double-box diagram as well as its related daughter diagrams, i.e., the double-triangle diagram and the triangle-box diagram. For later two kinds of diagrams, we have given complete analytical results in general (4-2\eps)-dimension.


Introduction
Currently, the focus of high energy physics is the LHC experiment. To understand the experiment data, we need to evaluate scattering amplitudes to high accuracy level required by data. Thus for most processes, the one-loop evaluation becomes necessary. In last ten years, enormous progress has been made in the computation of one-loop scattering amplitudes(see, for example, the references [1][2][3] and citations in the papers). However, for some processes in modern colliders, such as the process gg → γγ which is an important background for searching the Higgs boson at the LHC, one-loop amplitudes do not suffice since their leading-order terms begin at one loop. Thus next-toleading order corrections require the computation of two-loop amplitudes [4][5][6]. The traditional method for amplitude calculation is through the Feynman diagram. This method is well organized and has clear physical picture. It has also been implemented into many computer programs. However, with increasing of loop level or the number of external particles, the complexity of computation increases dramatically. Thus even with the most powerful computer available, many interesting processes related to LHC experiments can not be dealt by the traditional method.
To solve the challenge, many new methods(see books [7][8][9]) have been developed, such as IBP(integrate-by-part) method [10][11][12][13][14][15][16][17][18][19](some new developments, see [20][21][22]), differential equation method [23][24][25][26][27][28][29][30], MB(Mellin-Barnes) method [31][32][33][34], etc. Among these methods, the reduction method [35][36][37][38] is one of the most useful methods. More explicitly, the reduction of an amplitude means that any amplitude A can be expanded by bases(or "master integral") as with rational coefficients c i . With this expansion, the amplitude calculation can be separated into two parts: (a) the evaluation of bases(or master integrals) at given loop order and (b) the determination of coefficients c i for a particular process. For the former part, it can be done once for all and the results can be applied to any process. Thus in the practical application, the latter part, i.e., the determination of coefficients, becomes the central focus of all calculations. Unitarity method is an ideal tool to determine coefficients [39][40][41][42][43][44][45][46][47][48][49][50][51][52][53][54]. With the expansion (1.1), if we perform unitarity cut on both sides, we will get So if both ∆A and ∆A i can be evaluated analytically, and if different ∆A i has distinguishable analytic structure(which we will call the "signature" of basis under the unitarity cut), we can compare both sides of (1.2) to determine coefficients c i , analogous to the fact that if two polynomials of x are equal, so are their coefficients of each term x n . The unitarity method has been proven to be very successful in determining coefficients for one-loop amplitudes(see reviews [55,56]). For some subsets of bases(such as box topology for one-loop and double-box topology for planar two-loop), more efficient method, the so called "generalized unitarity method"(or "maximum unitarity cut" or "leading singularity"), has been developed [43,[57][58][59][60][61][62][63][64][65][66][67]. The applicability of reduction method is based on the valid expression of expansion (1.1). Thus the determination of bases becomes the first issue. From recent study, it is realized that there are two kinds of bases: the integrand bases and the integral bases. The integrand bases are algebraically independent rational functions before performing loop integration. For one-loop, the integrand bases have been determined by OPP [68]. For two-loop or more, the computational algebraic geometry method has been proposed to determine the integrand bases [69][70][71][72][73][74][75][76][77][78][79].
In general the number of integrand bases is larger than the number of integral bases, because after loop integration, some combinations of integrand bases may vanish. For one-loop amplitudes, the difference between these two numbers is not very significant. For example, the number of integral bases is one while the number of integrand bases is seven for triangle topology of renormalizable field theories [68]. However, for two-loop amplitudes, the difference could be huge. As we will show later, for double-triangle topology, there are only several integral bases, while the number of integrand bases is about one hundred for renormalizable field theories [78]. Thus the determination of integral bases for two-loop and higher-loop becomes necessary.
Although integrand bases can be determined systematically, the determination of integral bases is far from being completely solved. It is our attempt in this paper to find an efficient method to solve the problem. Noticing that in the unitarity method, the action ∆ in (1.2) is directly acting on the integrated results, thus if the left hand side ∆A can be analytically integrated for arbitrary inputs, we can classify independent distinguishable analytic structures from these results. Each structure should correspond to one integral basis 2 . By this attempt we can determine integral bases.
In this paper, taking double-box topology and its daughter topologies as examples, we generalize unitarity method to two-loop amplitudes and try to determine the integral bases. Different from the maximal unitarity method [60], we cut only four propagators(the propagator with mixed loop momenta will not be touched). Comparing with maximal unitarity cut where solutions for loop momenta are complex number in general, our cut conditions guarantee the existence of real solutions for loop momenta, thus avoiding the affects from spurious integrations. This paper is organized as follows. In section 2 we review the one-loop unitarity method and then generalize the scheme to two-loop. For two-loop, two sub-one-loop phase space integrations should be evaluated. In section 3, we integrate the first sub-one-loop integration of triangle topology. The result is used in section 4, where integration over the second sub-one-loop of triangle topology is performed. Results obtained in this section allow us to determine integral bases for the topology A 212 . Results in section 3 is also also used in section 5, where integration over the second subone-loop of box topology is performed, and the result can be used to determine integral bases for topology A 213 . In section 6, we briefly discuss the integral bases of topology A 313 since results are well known for this topology. Finally, in section 7, a short conclusion is given.
Technical details of calculation are presented in Appendix. In Appendix A, some useful formulae for phase space integration are summarized. In Appendix B, the phase space integration is done for one-loop bubble, one-loop triangle and one-loop box topologies. In Appendix C, details of an integration for topology A 313 are discussed.

Setup
In this section, we present some general discussions about the calculation done in this paper. Firstly, we review how to do the phase space integration in unitarity method illustrated by one-loop example. Then we set up the framework in unitarity method for two-loop topologies which are the starting point of this paper.
For one-loop, the action ∆ in (1.2) is realized by putting two internal propagators on-shell. More explicitly, let us consider the following most general input 3 with massless internal propagators 4 where the inner momentum is in (4 − 2ǫ)-dimensional space and all external momenta are in pure 4D space for our regularization scheme. The unitarity cut with intermediate flowing momentum K is given by putting ℓ 2 and ( ℓ − K) 2 on-shell, and we get the expression With two delta-functions, the original (4 − 2ǫ)-dimensional integration is reduced to (2 − 2ǫ)dimensional integration. To carry out the remaining integration, we decompose ℓ as ℓ = ℓ + µ, where ℓ is the pure 4D part while µ is the (−2ǫ)-dimensional part [51], then the measure becomes Next, we split ℓ into ℓ = ℓ + zK with ℓ 2 = 0 to arrive Having the form (2.4), we can use the following well known result of spinor integration 5 [80]. Define null momentum as ℓ = tλ λ, then The most general expression for numerator will be i j (ℓ · R ij ). For each term n j=1 (ℓ · R ij ), we can construct (ℓ · R i ) n with R i = n j=1 y j R ij . Thus if we know the result for numerator (ℓ · R i ) n , we can expand it into the polynomial of y i and read out corresponding result for n j=1 (ℓ · R ij ). 4 For simplicity we consider the massless propagators, but massive propagators can be dealt similarly. 5 For one-loop, we can take either positive light cone or negative light cone, where for negative light cone, the t-integration will be 0 −∞ . For two-loop, it can happen that if we take positive light cone for ℓ 1 , then we need to take negative light cone for ℓ 2 . However, the choice of light cone only gives an overall sign and does not affect λ, λ integration. Substituting (2.5) back to (2.4), we can use remaining two delta-functions to fix t and z as After above simplification, the integral (2.2) is transformed to the following spinor form where To deal with the integral like λ|dλ λ|d λ f (λ, λ) when f (λ, λ) is a rational function, the first step is to find a function g(λ, λ) satisfying With g(λ, λ), the integration is given algebraically by the sum of residues of holomorphic pole in g(λ, λ) [46]. In Appendix B, we summarize some general results of standard one-loop integrations using above technique. It is worth to mention that for two-loop, f (λ, λ) might not be rational function. We will discuss how to deal with it later. We also want to remark that under the framework of (4 − 2ǫ)-dimensional unitarity method, coefficient of each basis will be polynomial of µ 2 (remembering the splitting ℓ = ℓ + µ). There are two ways to handle it. For the first way, one can further integrate d −2ǫ µ (µ 2 ) n to find coefficients depending on ǫ. For the second way, we just keep µ 2 , but include the dimensional shifted scalar basis [42,81], such as .

(2.10)
This is equivalent to For one-loop, dimensional shifted bases are often used. In this paper we adapt the similar strategy, i.e., keeping the µ-part and introducing the dimensional shifted bases.

Generalizing to two-loop case
In this subsection, we set up unitarity method for two-loop amplitudes, particularly for the attempt of determining integral bases. The first problem is to decide which propagators should be cut. There are three kinds of propagators: (1) propagators depending on ℓ 1 only; (2) propagators depending on ℓ 2 only; (3) propagators depending on both ℓ 1 and ℓ 2 . In principle, we can cut any propagators, but for simplicity, in this paper we will cut propagators of the first two kinds. For our choice, we cut two propagators of the first kind and two propagators of the second kind. With this arrangement, for each loop it is exactly the familiar unitarity method in one-loop case. Next we set up notation for two-loop integral. The two internal momenta are denoted as ℓ 1 , ℓ 2 in (4 − 2ǫ)-dimension, while all external momenta are in pure 4-dimension. We use n 1 , n 2 , n 12 to denote the number of each kind of propagators respectively. Then a general integrand with massless propagators 6 can be represented by 7 . (2.12) The unitarity cut action ∆ is then given by 8 With above setup, we take a well studied example [20,[60][61][62][63][64][65][66], i.e., the four-point two-loop double-box(A 313 ) integral as the target to apply the unitarity method and determine integral bases. The integrand is given by 14) and the four propagators to be cut are where K 12 + K 34 = 0. With this choice of cuts, in order to completely understand the results, we also need to consider other topologies besides double-box. The other contributions come from 6 In this paper, we consider the massless case only. For inner propagators with masses, we will leave to further projects. 7 In this paper, we use I for integrand and A for integral. 8 We have neglected some overall factors in the definition of integration since it does not matter for our discussion.

K1
K2 K3 those topologies by pinching one or more un-cut propagators of double-box, as shown in Figure  1. There are three daughter topologies I 213 , I 312 , I 303 by pinching one propagator. There are also three daughter topologies I 212 , I 302 , I 203 by pinching two propagators. Finally there is only one daughter topology I 202 by pinching three propagators. Among them, I 303 , I 203 , I 302 , I 202 are direct products of two one-loop topologies, thus their signatures are well known(see Appendix B). So in fact we need to examine two non-trivial topologies I 212 , I 213 (by symmetry I 312 is equivalent to I 213 ) together with the mother topology I 313 . Integrand of these two additional topologies are given by In the following sections, we will study I 212 , I 213 and I 313 one by one, and our basic strategy will be to integrate one loop momentum ℓ 1 first while keeping ℓ 2 arbitrarily. Then we analyze the integration of ℓ 2 based on the previous results.
3. The ℓ 1 -part integration (n 1 = 2) In this section, we do the ℓ 1 integration. Using the standard method for one-loop amplitudes(reviewed in previous section as well as in Appendix B) we get(see formula (2.12)) ∆A (a,b) where various quantities are defined as . Note that here the left cut momentum K L 1 = K 12 is the same to the right cut momentum K L 2 = K 34 up to a sign, however we keep them independently so that it is possible to formulate them to more general situations. The W 1 comes from the mixed propagator ( ℓ 1 + ℓ 2 ) 2 . Situations with non trivial topologies A 313 , A 312 , A 213 and A 212 are all included in the formula (3.1).
Let us apply our general framework to the specific case n 1 = 2. The general formula (3.1) now becomes The second line is nothing but the standard one-loop triangle integration(see Appendix B). When a = 0, the integration gives the signature of triangle part. When a ≥ 1, the integration can be decomposed into both triangle part and bubble part. We will evaluate contributions from these two parts separately.

The contribution to triangle part
The triangle signature : Based on our general formula of the standard one-loop triangle integration (B.8), the signature of the triangle part is Imposing cut conditions for ℓ 2 , i.e., δ( ℓ 2 2 − µ 2 2 ) and δ(K 2 L 2 − 2K L 2 · ℓ 2 ) we can simplify it to where we have introduced One can observe that the signature part does not depend on ℓ 2 . It is an important feature which makes ∆A (a,b) 21n 2 easier to be treated.
. (3.7) Again, using cut conditions δ( ℓ 2 2 − µ 2 2 ) and δ(K 2 L 2 − 2K L 2 · ℓ 2 ) we can do the following replacement Since all derivatives act on τ only, such replacement will not affect the result. Some algebraic manipulations shows that the coefficients of different parts are given by with To get non-zero contribution from d a dτ a (•) τ →0 , we only need to take terms with τ a power. It means that terms with τ 2 in (3.8) will always appear with terms τ 0 , therefore we can regroup Thus we can write , where F is defined as Putting all results together, the triangle part becomes To do the ℓ 2 -part integration, it is more convenient to use above form before taking the derivative over τ .

The contribution to bubble part
Again we use results given in Appendix B.
The R 3→2 [i, m] term: Using (B.10), the typical term of triangle topology to bubble is where two null momenta P 1 , P 2 are constructed as Again, to get non-zero contribution, P 1 |R 1 |P 2 ] and P 2 |R 1 |P 1 ] should always appear in pair. With a little calculations, one can see where y 1 and y 2 are defined in (3.9). Thus we can take the following replacements (3.14) After such replacements we obtain Above expression has the form λ 2 | • | λ 2 λ 2 | • | λ 2 . In order to use the results given in Appendix B, we need to rewrite them by using where we have defined 3.3 The result for n 1 = 2 after ℓ 1 -integration Collecting results from triangle part and bubble part we obtain where s, t 1 , t 2 are defined in (3.6), F in (3.10) and G 1 , G 2 in (3.17) 9 . The trick here is that instead of computing the operations d a dτ a (•) , we will firstly do the ℓ 2 -part integration. 9 Do not confuse the t 2 here with the t 2 -integration part of ℓ 2 as reviewed in (2.5).
For ℓ 2 -integration, after the t 2 -integration we are left with spinor integration given by 10 where we have defined

The integral bases of A 212 topology
With results of previous section, it is possible to discuss the integral bases of A 212 topology in this section. To do so, we need to finish the spinor integration given in (3.19) with n 2 = 2, and attempt to identify the results to the bases. We will see that there are only (dimensional shifted) scalar bases.

4.1
The λ 2 -integration for the case n 2 = 2 For the case n 2 = 2 the formula (3.19) becomes For our momentum configuration, K L 1 = −K L 2 , thus we can combine denominator together to get a simpler expression. Terms of integrand can be classified into two parts, and we evaluate them one by one.
First part: The first part can be rewritten as The second line is the standard one-loop bubble integration, thus we can use the general formulae in Appendix B.
Second part: The second part can be rewritten as The second line is again the one-loop bubble integration. After finishing the integration over λ 2 -part, we can take the derivative and the limit τ → 0, τ 1 → 0, τ → 0.

The result
Collecting all results together, we get an expression of the form where we have defined Remind from Appendix B that the signature of one-loop bubble is d −2ǫ µ(− √ 1 − u 2 ), thus the term S 202 is the signature of topology A 202 as the subscript indicates. For S 212 , since the factor can not be factorized to a form where µ 1 -part and µ 2 -part are decoupled, it can not belong to the topology A n 1 0n 2 . So it must be the signature of topology A 212 .
It is worth to mention that in the form (4.4), the dependence of a, b is completely encoded in the coefficients f 212→212 satisfying the following two conditions: (1) they are polynomials of u 1 , u 2 and s; (2) they are rational functions of external momentum K L 1 . More discussions will be given shortly after.
Having above general discussions, now we list coefficients for various a, b: Coefficients f 212→212 : Using expression given in Appendix B, the analytic results for some levels of a + b are given by (4.7) • a + b = 3: Coefficients f 212→202 : • a = 0 or b = 0: From our derivation, it can easily be seen that when a = 0 or b = 0, the coefficient must be zero, i.e., • Non-zero results:

Classification of integral bases
Now we need to analyze above results in order to determine the integral bases. Firstly, noticing that f 212→202 are polynomials of T 1 , T 2 , µ 1 · µ 2 , µ 2 1 , µ 2 2 as well as rational functions of external momentum K L 1 , thus we can write them more explicitly as where the tensor coefficients f The above expansion leads us to define following dimensional shifted bases and , when κ 0 = 0, we do have µ 1 · µ 2 in the numerator. Thus although there is no mixed propagator in the denominator, it contains information from the mother topology A 212 where ℓ 1 and ℓ 2 are mixed.
With above definition, we find the following reduction hinted by unitarity method 11 However, before claiming B 212 [κ 0 , κ 1 , κ 2 ] are bases of the topology A 212 studied in this paper, we need to notice that in general T i could have four independent choices in 4D, i.e., e i , i = 1, 2, 3, 4 as the momentum bases for Lorentz momenta. So if bases have non-trivial dependence of T i in the numerator, we should be careful to identify bases. This happens to topologies A 213 and A 313 . However, for the current topology A 212 , the bases B 212 [κ 0 , κ 1 , κ 2 ] are scalar bases, i.e., the numerator of bases does not depend on any external momenta T i . Now we count the number of integral bases. For pure 4D case, we can take the limit µ 2 1 , µ 2 2 , µ 1 · µ 2 → 0, thus there is only one scalar basis, with κ i = 0, i = 0, 1, 2. In [78] it is found that for planar double-triangle(i.e., the topology A 212 ), the number of integrand bases is 111 under the renormalizable conditions in pure 4D. For general (4 − 2ǫ)-dimension, if we set constraint i=0,1,2 κ i ≤ 3(i.e., the sum of the power of ℓ 1 , ℓ 2 in the numerator is less than or equal to 6) for well-behaved quantum field theories, the number of integral bases is 20.

The integral bases of A 213 topology
Encouraged by the results in previous section, in this section we determine the integral bases of A 213 topology. As it will be shown shortly after, new features will appear.

λ 2 -integration for the case n 2 = 3
For n 2 = 3 the general formula (3.19) becomes(for simplicity, we will drop "τ i → 0" from now on) Again, using K L 1 = −K L 2 we can simplify the denominator. There are also two parts we need to compute.
First part: The first part of remaining integration can be rewritten as The second line is the standard one-loop triangle integration. The one-loop triangle can be reduced to triangle part and bubble part, thus they can be interpreted as contributions from topologies A 213 and A 212 .
Second part: The second part can be written as The second line is again the standard triangle integration which contain contributions from topologies A 203 and A 202 .

Overview of results
Collecting all results together, we get an expression of the form where S 202 and S 212 have been defined in (4.5) and two new signatures are There are a few remarks for expression (5.4). Firstly it is easy to see that the signature S 203 is the direct product of signatures of one-loop bubble and one-loop triangle. Secondly there are two logarithms in the signature S 213 : one depends on both µ 1 , µ 2 and the other only depends on µ 2 . Pictorially, the first logarithm is related to the mixed propagator (ℓ 1 + ℓ 2 ) 2 while the second logarithm is related to the right hand side sub-triangle. Thirdly, all dependence of a, b are inside coefficients f while signatures are universal. However, unlike in the expression (4.4) where coefficients f are all rational functions of external momenta and polynomials of s, u 1 , u 2 , here we find that the coefficients f are in general not polynomials of s, u 1 , u 2 . In fact, factor t 2 = √ 1 − u 2 will appear in denominators. Such behavior can not be explained by dimensional shifted basis. Instead, we must regard it as the signature of new integral basis. Because of such complexity, when talking about the signature of a basis for A 213 topology, we should treat all coefficients together in a list {f 213→202 } as a single object. More explicitly we will write the expression (5.4) as The reduction of ∆A 213 is to write it as the linear combination Having above general remarks, now we present explicit results. The case b = 1: The result is Thus, at least for our choice of unitarity cuts, ∆A

The result of a = 1
In this case a non-trivial phenomenon appears, and we will show how to explain it.
Although the s 0 -part can be explained by the signature ∆A 213→213 , other coefficients are given by . where Since above four coefficients are rational functions of external momenta and polynomials of u 2 , we can claim that A (1,1) 213 is not a basis at least for our choice of unitarity cuts. There are some details we want to remark. The coefficient a 11→00 is a polynomial of T 1 and T 2 with degree one while coefficient a 11→10 is a polynomial of T 2 with degree one but rational function of T 1 . More accurately, both the denominator and the numerator of a 11→10 are polynomials of T 1 with degree one. It is against the intuition since T 1 should not appear in the denominator. However, this subtlety is resolved if one notice that the first component f 213 contains exactly the same factor [−(K 4 · T 1 )K 2 L 1 + (K 4 · K L 1 )(K L 1 · T 1 )] in its numerator, so it cancels the same factor in denominator of a 11→10 .
The case b = 2: The whole expression is too long to write down, thus we present only the general feature. Again although all coefficients contain factor 1 t 2 2 , the whole result can be expanded like the one (5.13) with coefficients as rational functions of external momenta and polynomials of s, u 1 , u 2 . Thus A (1,2) 213 is not a new integral basis.

The result of a = 2
We will encounter similar phenomenon as in the case a = 1. To get rid of tedious expressions, we will present only the main features.

Classification of integral bases
With above results, we can classify the integral bases of A 213 topology. Having shown that coefficients such as a 21→00 are polynomials of µ 1 · µ 2 , µ 2 1 , µ 2 2 and rational functions of external momenta, we can expand them, for example where the tensor coefficients a are rational functions of external momenta. This expansion leads us to define the following dimensional shifted bases Unlike the scalar basis B 212 [κ 0 , κ 1 , κ 2 ] for A 212 topology, the basis B 213;a [κ 0 , κ 1 , κ 2 ] depends on T 1 explicitly. Since T 1 is a 4-dimensional Lorentz vector, there are four independent choices and we need to clarify if different choice of T 1 gives new independent basis.
To discuss this problem we expand The momentum bases e i are constructed as follows. Using K 4 , K L 1 we can construct two null momenta , thus the momentum bases can be taken as It can be shown that T 1 = e 3,4 are spurious and the integrations are zero.
It is, in fact, equivalent to the basis B 213;0 [κ 0 , κ 1 , κ 2 ] and does not give new integral basis.
which is the true new integral basis.
Conclusion: For a = 1, the integral basis is given by B 213;1 [κ 0 , κ 1 , κ 2 ; K 4 ] The case a = 2: There are ten possible combinations ( ℓ 1 · e i )( ℓ 1 · e j ). With the explicit result we found that which is proportional to (5.24) by a factor 2K 2 L 1 . Therefore it can be reduced to the basis B 213;1 [κ 0 , κ 1 , κ 2 ], and dose not give a new integral basis.
Thus we can take either one(but only one of them) as the integral basis. We choose the combination (e i , e j ) = (e 1 , e 1 ) to be a new integral basis. For pure 4D, we just need to set µ 1 · µ 2 , µ 2 1 , µ 2 2 to zero. In this case, the factor 1 t n 2 → 1. In other words, there is only one integral bases It is useful to compare it with about 70 integrand bases found in [78] under renormalizable conditions.

The integral basis of A 313 topology
In this section we turn to the topology A 313 . This topology has been extensively studied by various methods, such as IBP method [20] and maximum unitarity cut method [60][61][62][63][64][65][66], and its integral bases have been determined [20]. To determine the integral bases using our method, we need to integrate the following expression For general situation, the integration is very complicated and we postpone it to future study. In this paper, we take the following simplification. Firstly we take all out-going momenta K 2 i = 0 (i = 1, .., 4)(unlike the topologies A 212 and A 213 where K i can be massive). Secondly, based on the known results of integral bases, we focus on the specific case a = 0 and In order to make expressions compact we define some new parameters as 12 For physical unitarity cut, momentum configuration requires s 12 > 0, s 13 < 0 and s 14 < 0. So we have −1 < m < 1 , − 1 < χ < 0 (6.3) 12 It is worth to notice that s in this section is different from s in (3.6) of section 3.
by momentum conservation s 12 +s 13 +s 14 = 0. Furthermore, we define the regularization parameters γ i as where the dimensionless extra-dimensional vector ν i is defined as Under the simplification a = 0, the integration over λ 1 -part is trivial. Using (B.21) in Appendix B we can get An important feature is that the signature after λ 1 -integration depends on ℓ 2 explicitly, which is different from the signature in (3.5). Because of this, the integration over λ 2 becomes very complicated. One way to overcome is to use Thus the logarithmic part in (6.5) becomes rational function of ℓ 2 and we can use the same strategy as in previous sections. However, for the current simple situation, we can use another method. After expanding the spinor variables as

θ-integration:
The θ-dependent part of (6.5) is given by with t 2 = γ 2 −1 2 . Setting x = e iθ the integral becomes a circle contour integration with radius one . (6.10) There are three poles in total. The first one is x = 0 when b = 0 for general R 2 . The other two are roots of the quadratic polynomial in denominator It is easy to check that |x 1 x 2 | = 1. Thus one root is inside the integration contour and the other is outside. The kinematic conditions s 12 > 0, s 24 < 0, s 14 < 0 ensure that x 1 is the one inside. The residue at the pole x 1 is The case (T 2 = K 1 ): Under our simplification, we set T 2 = K 1 , thus K 1 |R 2 |K 2 ] = 0 and K 2 |R 2 |K 1 ] = 0. Because of this, there is no pole at x = 0 in (6.10). Thus after the θ-integration, (6.5) is reduced to , (6.14) in which α = γγ 1 , β = (γ 2 − 1)(γ 2 1 − 1) .
Since we use the dimensional shifted bases, the µ i part is kept and we will focus on D (6.16) We found it hard to integrate over u and get analytic results. However, in the general (4 − 2ǫ)dimensional framework, we can treat µ 2 i and µ 1 · µ 2 as small parameters and take series expansion around µ 2 i → 0. It is equivalent to taking the series expansion around γ i → 1. The details of calculation can be found in Appendix C. Up to the leading order, result for a = 0, b = 0 is given by An important check for the result (6.17) is that it has the S 3 permutation symmetry among γ 1 , γ 2 , γ. The terms ln(−χ) and Li 2 (1 + χ) do not show up for topologies A 212 and A 213 , thus they belong to the signature of A 313 . For b = 1, the result is The extra term − 1 s 2 ln −2χ γ−1 ln −2χ γ 1 −1 in (6.18) indicates that comparing to D 302 , which is the consequence of symmetry γ ↔ γ 1 in (6.15). Secondly, the appearance of term ln(−χ) is quite intriguing. There are several possible interpretations: • Under the general (4−2ǫ)-dimensional framework, D (0,2) 313 could be considered as a new integral bases.
• From the result in [20], A • In fact, ln −2χ is the result given by unitarity cut channel K 12 of one-loop massless box (K 1 , K 2 , K 3 , K 4 ) up to zero-order of (γ i − 1). It may indicate some connection with one-loop box diagram.
Finally for b = 3 we found It is obvious that D

Conclusion
In this paper we applied the unitarity method to two-loop diagrams to determine their integral bases. Two propagators for each loop are cut while mixed propagators are untouched. Integrations for the reduced phase space have been done in the spinor form analytically. Based on these results, analytical structures have been identified and integral bases have been determined.
To demonstrate, we applied our method to investigate the double-box topology and its daughters, with appropriate choice of cut momenta and kinematic region. For the A 212 topology, we found that there is only one scalar basis for the pure 4D case, while for general (4 − 2ǫ)-dimension, if we use the dimensional shifted bases, there are 20 scalar bases under good renormalizability conditions. For the A 213 topology, there is also only one scalar basis for the pure 4D case, but for the (4 − 2ǫ)-dimension, scalar bases are not enough even considering the dimensional shifted bases. We found that there are 48 dimensional-shifted integral bases for renoramalizable theories. For the A 313 topology, it is difficult to get an exact expression for general (4 − 2ǫ)-dimension case. Thus we only considered a specific case A (0,b) 313 with T 2 = K 1 . We presented results to the zeroth-order and found three bases for general (4 − 2ǫ)-dimension if we do not allow coefficients depending on ǫ.
Based on the method demonstrated in this paper, several possible directions can be done in the future. Firstly, for the A 313 topology, the exact result for the specific case a = 0 is still missing.
The general value of a should also be considered. Secondly, topologies discussed in this paper are not the most general cases. The most general configurations are those that each vertex has external momenta attached as well as massive propagators. Results of these more general cases are necessary. Thirdly, to obtain a complete set of integral bases, we need to investigate other topologies classified in [78]. Finally, besides determining the integral bases, the unitarity method is also powerful for finding rational coefficients of bases in the reduction. We expect that, after the complete set of integral bases being obtained, such method can be useful for practical two-loop calculations 13 .

A. Some useful formulae
In this section, we present some useful formulae appearing in various calculations in the paper.

Total derivative:
For the holomorphic anomaly method, it is important to write an expression into the total derivative form. Here we list results for two typical inputs: Pole of ℓ|QK|ℓ : In the calculation, we will meet pole of the form ℓ|QK|ℓ frequently. It contains two poles and we need to separate them. If both Q, K are massless we can write it as ℓ|Q [Q|K] K|ℓ . If at least one of them is massive, for example K, we can construct two massless momenta as P i = Q + x i K, i = 1, 2, where Using this we have Residue of high order pole: Poles we met are often not single poles. To read out residues of poles with high order we can do as follows. Using the expression It is very important to emphasize that the |ℓ] part has been set to |η], while the |ℓ is replaced by (|η − τ |s ) before taking the derivative. Evaluation of P 1 |R|P 2 ] P 2 |S|P 1 ]: We often encounter expression P 1 |R|P 2 ] P 2 |S|P 1 ], which can be evaluated as where ǫ(P 1 RP 2 S) denotes ǫ µνρσ P µ 1 R ν P ρ 2 S σ . To evaluate ǫ(P 1 RP 2 S) 2 , a simple way is to consider P 1 |R|P 2 ] P 2 |S|P 1 ] P 1 |S|P 2 ] P 2 |R|P 1 ].

B. Standard one-loop integrations
In this section we list some standard one-loop results. We focus on the following standard integral [52] which is the integration in (2.7). In our application, we only need cases n = 2, 3, 4.

B.1 The bubble integration
When n = 2, we have ∆A where (A.1) has been used. For the pole λ|RK|λ , we use the construction given in Appendix A to read out two poles λ|P 1 and λ|P 2 with P i = R + x i K(see (A.3)). For the first pole |λ = |P 1 , the residue is while for the second pole |λ = |P 2 , the residue is Putting them together we obtain Let us give a few examples: The case a = 0 gives the analytic signature S bub = (− √ 1 − u) for the one-loop scalar bubble basis. For cases a = 1, 2, the part inside the curly bracket is indeed polynomial of u.

B.2 The triangle integration
For the case n = 3, we can split the integrand as follows After the splitting, the first term inside the big bracket produces the signature of triangle, while the second term produces the signature of bubble. Thus we have the following two standard integrations.

B.2.1 Triangle-to-triangle part
For the first term, writing into total derivative we have The pole is given by factor λ|QK|λ a+1 . Using results in Appendix A, for the pole η = P 1 with auxiliary spinor s = P 2 the residue is .
For the pole η = P 2 with auxiliary spinor s = P 1 the residue is .
One can observe that the derivative part is in fact the same for both contributions after taking the limit τ → 0. Thus the sum of two contributions is .
After some manipulation, we finally have where S tri is the signature of triangle and C (a) 3→3 is the corresponding coefficient: The a = 0 case gives the result for standard scalar triangle and other a's, give the corresponding coefficients under the reduction. One can verify that the coefficients are indeed rational functions.

B.2.2 Triangle-to-bubble part
The typical term in (B.5) for triangle-to-bubble part is The residue of pole λ|QK|λ i+1 can be read out as in previous subsubsection and we get . (B.10) To get a Lorentz contracted form, we need to use the following key fact: to have non-zero contribution, factors P 1 |R|P 2 ] and P 2 |R|P 1 ] should always appear in pair. Thus we can transfer (B.10) to Thus the contribution for the triangle-to-bubble part is given by Putting two parts together, we get 3→2 . (B.14)

B.3 The box integration
The box integration is given by Now we can evaluate various parts one by one.

B.3.1 The box-to-box part
This part comes from pole λ|Q 1 Q 2 |λ in (B.16). Using Q 1 + x i Q 2 to construct two null momenta P i , we get the residue which can be written as

(B.18)
Box part: The first term in (B.18) produces the signature of box as well as the coefficient It can be shown that there is a recursion relation where Triangle part: The second term in (B.18) produces the signature of triangle. Using the second term in (B.18) can be rewritten as We will combine (B.24) with results in the next subsubsection to produce the complete triangle part.

B.3.2 The box-to-triangle part
Since Q 1 and Q 2 are symmetric, we will focus on the triangle constructed by K, Q 1 . The contribution comes from the first term of (B.16). This term contains two kinds of poles: λ|Q 1 Q 2 |λ and λ|KQ|λ . The contribution of pole λ|Q 1 Q 2 |λ has been evaluated in previous subsubsection. For the second pole, after writing it into total derivative, it is 1 λ|KQ|λ a . Using Q 1 , K to construct two null momenta P 1 , P 2 , the residue is given by two parts. The first part contains ln(x 1 x 2 ) (which is nothing but ln It cancels the first term of (B.24), since the sum of all residues of a holomorphic function is zero 14 . The second part contains ln(x 1 /x 2 ) which is the signature of triangle. The contribution can be written as To write the spinor form to the Lorentz contracted form, we can take similar manipulation as the one from (B.10) to (B.11). The result is where we have defined It is worth to mention that T 8 appears as T 2 8 , thus the Levi-Civita symbol has been removed.
14 It is worth to notice that by power counting, infinity does not contribute residue.

B.3.3 The box-to-bubble part
Having finished the computation of (B.16), we turn to the (B.17). The total result can be expressed as where the typical term is There are three poles for this part: λ|Q 1 Q 2 |λ , λ|Q 1 K|λ and λ|Q 2 K|λ . The contribution from λ|Q 1 Q 2 |λ is zero when summing up the two lines in (B.17). For the remaining two poles, because of the symmetry Q 1 ↔ Q 2 , we will focus on R 4→2 (Q 1 )[i, m] only. We use Q 1 , K to construct two null momenta P 1 , P 2 and get residue . (B.31) We can rewrite the expression to the following Lorentz contracted form One can verify that T 8 will appear as T 2 8 after summing R 4→2 (Q 1 )[i, a] and R 4→2 (Q 2 )[i, a], thus the Levi-Civita symbol does not appear in the final result.

C. The integration for topology A 313
It is hard to get the explicit result for (6.15). In this Appendix we develop a method to find approximate expressions. Technically the case b = 0 is the most complicated one, while the b ≥ 1 cases can be reduced to the case b = 0 plus some simple integration. Before working out the integration case by case, we give two explicit integrations where we have used the conditions −1 < m < 1 and γ i ≥ 1. These two results are useful for our further discussion.
For the limit a → 0, it is easy to see that only in the case n = 0 it is divergent and we have We can use parameterizing method to sum up above awesome form. If we define then G(x) satisfies the differential equation Using the same method, we can find the differential equation for G(0, x). After some variable replacement we get G(0) = G(0, 1) = 2 0 ρ g(x) dx .
(C. 16) Combining (C.14) with (C.16) and exchanging the integration and the summation finally the finite term can be written as Now we focus on the indefinite integration. If we change the integration variable as cosh y ≡ α + x β , cosh y ± ≡ α ± 1 β , cosh y 0 ≡ α − mγ 2 β , (C. 19) 313 | f inite is given by F (y + ) + F (y − ) − 2F (y 0 ) up to an overall factor. Above calculations are done for pure 4D limit of γ 2 . After taking the pure 4D limit of γ 1 , γ we finally reach after combining with the divergent term (C.9). It means that as a function of γ 2 , above combination is finite under the limit γ 2 → 1. Thus to our zero-order(i.e., the pure 4D case), we can just set γ 2 = 1 before doing the integration and get For the case b=2,3 we will not show the computation details again. The main point is that in the first step, we choose a proper linear combination of D Then we can prove that these combinations are convergent at γ 2 → 1 just as in the case b = 1. Thus we can take γ 2 = 1 before integrating those combinations. The second step is to integrate these dramatically simplified integrands. In this step the most efficient way is to use the variable replacement in (C.19) and we can integrate quickly. For example in the case b = 2, combination D (0,2) 313 − χ sD (0,1) 313 is given by 35) and the corresponding expressions for b = 3 are even longer. To find the approximate results in the pure 4D, in the last step we take limit γ, γ 1 → 1. Carrying out these steps, finally we get