Real-Virtual-Virtual contributions to the inclusive Higgs cross section at N3LO
aa r X i v : . [ h e p - ph ] N ov Preprint typeset in JHEP style - PAPER VERSION
Real-Virtual-Virtual contributions to the inclusiveHiggs cross section at N LO Falko Dulat
ETH Zurich, 8093 Zurich, SwitzerlandE-mail: [email protected]
Bernhard Mistlberger
ETH Zurich, 8093 Zurich, SwitzerlandE-mail: [email protected]
Abstract:
We present the computation of the contributions to N LO inclusive Higgsboson production due to two-loop amplitudes. Our result is a Laurent expansion in thedimensional regulator, with coefficients that are linear combinations of harmonic polylog-arithms of the ratio of the Higgs boson mass and the partonic center of mass energy. Weoutline our method of solving the differential equations for master integrals appearing inthe cross section. Solving these differential equations requires the determination of bound-ary conditions and we present a new technique for decomposing the boundary conditionsinto physical contributions. We show how these boundary conditions can be calculatedusing the method of expansion by regions.
Keywords:
Higgs, QCD, NNLO, N3LO, LHC and differential equations. . Introduction
The accurate theoretical prediction of the standard model inclusive Higgs boson crosssection is an important ingredient for the study of the Higgs boson at the Large HadronCollider. The uncertainty associated with the truncation of the perturbative series atNNLO [1] is currently estimated to be of the order of ±
10% (see for example [2]). AtN LO there exist several approximations, which, however, due to the uncertainty inherentto ad-hoc approximations are unable to significantly reduce the total theoretical uncertaintyon the cross section. With the accumulation of more data by ATLAS and CMS, theexperimental uncertainty will further decrease in the course of the next run of the LHC,making the perturbative uncertainty eventually one of the largest systematic uncertaintiesentering the extraction of coupling strengths of Higgs boson interactions.An exact computation of the inclusive Higgs cross section at N LO is expected toreduce the uncertainty to about ±
5% [3] and is therefore ultimately required for the futureof Higgs physics at the LHC. Many important steps towards this goal have already beentaken in recent years. The three-loop corrections to the gg → h amplitude have beencomputed in refs. [4]. The renormalisation of collinear and ultraviolet divergences and theassociated counterterms, as well as the partonic cross sections at lower orders have beencomputed in refs. [3,5,6]. The N LO corrections due to triple-real radiation were computedas a threshold expansion around the soft limit in ref. [7]. The two-loop soft current, whichrepresents the two-loop correction to the real emission of a single soft parton was publishedin ref. [8]. Another important contribution is the square of the one-loop correction to theemission of one parton, this was computed without approximation in ref. [9]. By combiningthese contributions together with the one-loop correction to the real emission of two softpartons the first term in the “soft” expansion of the Higgs cross section at N LO wasobtained in [10]. The one-loop corrections with two soft emissions were confirmed inref. [11]The first term in the soft expansion represents an important first step towards thecalculation of the full Higgs cross section at N LO and has been used to obtain severalapproximations of the full cross section [12]. In parallel to this work, the second terms in thethreshold expansion was made public in ref. [13] relying on results presented in this article.In ref. [13] it was shown that in order to make reliable predictions for the cross sectionthat can be used at the LHC, knowledge of the exact cross section for Higgs productionat N LO is highly desirable. It is therefore required to perform exact computations ofthe contributions that are so far only known as expansions. These contributions are thetriple-real emission contributions, the one-loop corrections to double-real emission as wellas the two-loop corrections to single-real emission.In this article we focus on the genuine two-loop corrections to the emission of a singleparton in association with the Higgs, which we denote as
RV V . This contribution is theintegration of the two-loop corrections to Higgs plus one parton over the two-particle phasespace.To conduct our calculation we employ the framework of reverse unitarity [14] to renderthe combined loop and phase space integrals appearing in the cross section accessible to– 1 –he method of differential equations [15]. We review how our differential equations canbe solved in terms of a class of special functions, so-called multiple polylogarithms [16].Furthermore, we describe our method of determining the boundary conditions that areneeded to specialise the general solutions of the differential equations. We formulate a newmethod to disentangle the different physical contributions to the boundary conditions andfacilitate their calculation by establishing a direct connection to the method of expansionby regions [17, 18].The article is organised as follows. In section 2 we introduce the general setup used tocompute the partonic cross sections. We discuss our results in section 3 and comment onthe computation of the matrix-elements in section 4. In section 5 we explain the methodof calculating our master integrals using differential equations. We provide the generalstrategy as well as a detailed description of our novel method facilitating the calculationof necessary boundary conditions. We give an example of how these boundary conditionsare calculated in section 5.5.The main contribution of this article is the calculation of the complete two-loop cor-rection to Higgs boson production at N LO and the development of a new framework fordetermining boundary condition for differential equations of master integrals.
2. The double-virtual real cross section
We consider the partonic QCD amplitudes for the production of a Higgs boson in associ-ation with one additional parton. We distinguish three different channels by their initialstate, g ( p ) + g ( p ) → g ( p ) + H ( p h ) q ( p ) + g ( p ) → q ( p ) + H ( p h ) q ( p ) + ¯ q ( p ) → g ( p ) + H ( p h ) (2.1)where q, ¯ q, g and H denote a quark, anti-quark, gluon or Higgs boson respectively with theirassociated momenta p . . . p , p h . This allows to define the following kinematic invariants s = 2 p · p + i , p = M h + i ≡ sz,s ¯ zλ = 2 p · p − i , s ¯ z ¯ λ = 2 p · p − i , (2.2)where z = M h s , ¯ z = 1 − z , ¯ λ = 1 − λ and where we explicitly indicate the small imaginarypart carried by the invariants. The partonic cross section for these processes is then givenby σ X = N X s Z d Φ X | M X | , (2.3)where X ∈ { g g → H g, q ¯ q → H g, g q → H q } . The summation sign indicates summationover final and initial state particle polarisations and colours. Throughout the paper wework in conventional dimensional regularisation with d = 4 − ǫ space-time dimensions.The process dependent factors N X containing the averaging of initial state parton colours– 2 –nd polarisations are given by N g g → H g = 14( N c − (1 − ǫ ) , N g,q → H q = 14 N c ( N c − − ǫ ) , N q ¯ q → H g = 14 N c . (2.4) N c and ( N c −
1) are the number of quark and gluon colours respectively. The phase-spacemeasure for the production of a massive Higgs boson in association with a massless partonis given by d Φ = d d p (2 π ) d δ + ( p ) d d p h (2 π ) d δ + ( p h − M h )(2 π ) d δ ( d ) ( p + p − p − p h ) , (2.5)where δ + ( p ) = (2 π ) δ ( p ) θ ( p ). Using the definitions of eq. (2.1) we can parameterise thephase-space measure as d Φ = (4 π ) − ǫ s − ǫ ¯ z − ǫ − ǫ ) dλ ( λ ¯ λ ) − ǫ θ ( λ ) θ (¯ λ ) . (2.6)In the rest of the paper we consider the mass of the top quark to be large enough for thetop quark to be decoupled from all interactions. This description can be formulated usingthe effective Lagrangian L eff = L QCD − CHG aµν G a,µν . (2.7) L QCD is the QCD Lagrangian with N f light quark flavours, H the Higgs boson field and G aµν the gluon field strength tensor. The Wilson coefficient C can be explicitly calculatedtaking into account the interactions of the top quark [19–23].We perform an expansion of the partonic scattering matrix-elements in the number ofloops | M X | = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ∞ X i =1 M ( i ) X (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , (2.8)where i runs over the number of loops. The main result of this work is the partonicscattering cross section arising due to the interference of two-loop matrix-elements withthe corresponding tree-level matrix elements and we refer to it as the double-virtual-real( RV V ) cross section. σ RV VX = N X s Z d Φ X ℜ (cid:16) M (2) X M (0) ∗ X (cid:17) (2.9)The cross section can be separated into five different contributions. σ RV VX ( z ) = X i =2 (1 − z ) − iǫ σ ( i ) RV VX ( z ) . (2.10)The individual terms σ ( i ) RV VX ( z ) do no longer contain logarithms with argument (1 − z ),i.e. they are meromorphic functions of z with at most a single pole at z = 1. Theycontain infrared and ultraviolet divergences that appear as poles in ǫ . The σ ( i ) RV VX ( z ) canbe written as Laurent series in the dimensional regulator. Each term in the series can– 3 –e expressed as a linear combination of multiple polylogarithms with rational coefficients.Multiple polylogarithms are a multivariate generalisation of the classical logarithm andpolylogarithms, ln( z ) = Z z dtt , Li n ( z ) = Z z dtt Li n − ( t ) , (2.11)with Li ( z ) = − ln(1 − z ). They can be defined analogously to eq. (2.11) via the iteratedintegral [16] G ( a , . . . , a n ; z ) = Z z dtt − a G ( a , . . . , a n ; t ) , (2.12)with G ( z ) = 1 and a i , z ∈ C . In the special case that all the a i vanish one defines G ( ~ n ; z ) = 1 n ! ln n ( z ) , (2.13)with ~a n = n z }| { ( a, . . . , a ). It is evident that the multiple polylogarithms encompass the classicalpolylogarithms as well as the harmonic polylogarithms (HPLs) G ( ~a n ; z ) = 1 n ! ln n (cid:16) − za (cid:17) , (2.14) G ( ~ n − , a ; z ) = − Li n (cid:16) za (cid:17) , (2.15) G ( a , . . . , a n ; z ) = ( − p H ( a , . . . , a n ; z ) , if a i ∈ {− , , } ∀ i, (2.16)where p is the number of elements in ( a , . . . , a n ) equal to +1. The number of elements inthe index vector ( a , . . . , a n ) is called the weight of the multiple polylogarithm.The multiple polylogarithms satisfy numerous algebraic relations, like the shuffle and stufflealgebras. A very important and useful property of the multiple polylogarithms is the factthat they satisfy a certain Hopf algebra [24, 25], which enables the algebraic derivation offunctional identities between multiple polylogarithms and is instrumental in algorithmicallyperforming iterated Feynman integrals [7, 25].
3. Results
We obtained the
RV V cross section for all partonic sub-channels completing the calculationof all two-loop contributions to the N LO Higgs boson cross section. In parallel to thiswork the authors of this article were part of a collaboration computing corrections to thethreshold expansion of the full Higgs boson cross section at N LO, which we made availablein ref. [13] using parts of this result as essential ingredients. In ref. [13] we also producedthe coefficients of the leading three threshold logarithms of the N LO Higgs boson crosssection. The result of this paper and specifically the possibility to decompose the
RV V cross section as in eq. (2.10) were key ingredients used in the derivation of the coefficientsof these logarithms. The first threshold-expansion coefficient of the
RV V cross section wasobtained in refs. [8] and agrees with the corresponding expansion coefficient of our result.Due to the length of the expressions we refrain from displaying the formulae for the
RV V cross section explicitly and make them available in a
Mathematica readable file– 4 –ogether with the arXiv submission of this paper. In this file we define the three variables sggRVV , sgqRVV and sqqbarRVV for the g g , g q and q ¯ q initial states respectively. Thesevariables contain the bare cross sections σ RV VX , σ RV VX = N X s Z d Φ X R (cid:16) M (2) X M (0) ∗ X (cid:17) , (3.1)as a Laurent expansion in the dimensional regulator. N X is given in eq. 2.4. For ourresults we separate the cross sections into contributions with a single pole at z = 1 fromthe remaining contributions that are analytic as z → σ RV VX ( z ) = X i ∈{ , , } (1 − z ) − − iǫ σ ( i ) sing X + X i =2 (1 − z ) − iǫ σ ( i ) reg X ( z ) . (3.2)For the singular terms we only expand the σ ( i ) sing X in the dimensional regulator up to order ǫ , leaving the prefactor unexpanded, while for the regular pieces we expand the product(1 − z ) − iǫ σ ( i ) reg X ( z ) up to order ǫ . We observe that the coefficients of the Laurent expansionof our cross section can be expressed as linear combinations of harmonic polylogarithmswith indices a i ∈ { , } . For convenience we have divided the expressions by a factor of κ = C α S e − ǫγ E (cid:0) s π (cid:1) − − ǫ π ( N c − . (3.3)
4. Calculation
To obtain all channels contributing to the
RV V cross section we compute the requiredtwo-loop and tree-level Feynman diagrams generated by qgraf [26]. We perform the con-traction of spinor and colour traces with custom C++ code based on the expression libraryGiNaC [27]. We work in Feynman gauge and restore gauge invariance by combining ourmatrix elements with the necessary Fadeev-Popov ghost matrix elements.Having performed all algebraic manipulations of the Feynman diagrams we arrive atour matrix-elements in terms of scalar products of loop and external momenta. Ratherthan carrying out the integration over the loop and phase-space momenta in a sequentialway, we treat all integrations on equal footing and combine them into a single integrationmeasure. d Φ = d d k (2 π ) d d d k (2 π ) d d Φ . (4.1)This combination allows us to apply the framework of reverse unitary [14]. Reverse uni-tarity exploits the duality between phase space integrals and loop integrals to treat themin a uniform way. Specifically, using Cutkosky’s rule, it is possible to express the on-shellconstraints appearing in phase space integrals through cut propagators δ + (cid:0) q (cid:1) → (cid:20) q (cid:21) c = 12 πi Disc 1 q = 12 πi (cid:20) q + i − q − i (cid:21) . (4.2)Due to the linearity of this representation, it is possible to differentiate cut propagatorswith respect to their momenta, similar to ordinary propagators. Therefore, it is possible– 5 –o derive integration-by-parts (IBP) identities [28] for phase space integrals in the sameway as for loop integrals. The fact that the cut propagator represents a delta function isimplemented by the simplifying constraint that any integral containing a cut propagatorraised to a negative power vanishes: (cid:20) q (cid:21) − nc = 0 , n ≥ . (4.3)IBP identities serve to relate different phase-space and loop integrals. The large systemof IBP identities for the integrals appearing in our calculation is solved using the Gausselimination algorithm [29], which we implemented in a private C++ code using the GiNaClibrary [27] . All integrals appearing in the cross section can be related to linear combina-tions of 72 master integrals. We discuss the methods used to solve our master integrals insection 5.
5. Calculating Master Integrals
In this section we describe the setup we used to solve our master integrals using firstorder differential equations. We start by deriving the required differential equations. Theirgeneral solution has to be constrained by fixing one boundary condition per integral. Wedescribe a novel way of facilitating the calculation of these boundary conditions. Next, wediscuss how the general solution of these differential equations can be computed. Finally,we illustrate the procedure using a simple example and demonstrate an explicit calculationof an actual
RV V boundary condition.
After integration over the final state and loop momenta eq. (4.1), the integrals are functionsof the Higgs mass M h and the partonic center of mass energy s . It is therefore convenientto define a single dimensionless ratio, z = M h s , ¯ z = 1 − z = s − M h s , (5.1)and write all master integrals as functions of this ratio, M i = M i (¯ z ) . (5.2)For brevity, we set s = 1 as the exact s dependence can be reconstructed using dimensionalanalysis.In order to evaluate the master integrals we use the method of differential equa-tions [15]. Since the master integrals are functions of a single ratio ¯ z , we can differentiatewith respect to the square of the Higgs mass, which only appears in the cut propagatorcorresponding to the Higgs on shell condition. As outlined in section 4 using the frameworkof reverse unitarity we find ∂∂ ¯ z (cid:20) p h − M h (cid:21) c = − ∂∂M h (cid:20) p h − M h (cid:21) c = (cid:20) p h − M h (cid:21) c . (5.3)– 6 –y applying this differential to our master integrals, we obtain a set of new phase spaceintegrals. Using IBP identities these integrals can again be expressed through our basis ofmaster integrals. This way we are able to express the differential of each master integralthrough the master integral itself as well as other integrals, obtaining a coupled system oflinear first order differential equations for the master integrals, ∂ ¯ z M i (¯ z ) = A ij (¯ z, ǫ ) M j (¯ z ) . (5.4)Einstein summation convention is implied. The entries of the system matrix A are ingeneral rational functions in ¯ z as well as in ǫ . The choice of basis is of course not uniqueand may be related to another one via a ¯ z and ǫ dependent transformation.We observe that the system matrix for the RV V master integrals can be written as A ij (¯ z, ǫ ) = A (0) ij (¯ z, ǫ )¯ z + A (1) ij (¯ z, ǫ )¯ z − , (5.5)where A (0) ij (¯ z, ǫ ) and A (1) ij (¯ z, ǫ ) are holomorphic functions of ¯ z . The solution to the above system of differential equations will require the specificationof one boundary condition per integral. We now present our method of obtaining theseboundary conditions.From the form of the differential equations (eqs. 5.4 and 5.5) we can see that thesolutions will have isolated singularities at ¯ z = 0 and ¯ z = 1. We can therefore pick onesingular point, in this case we choose the singularity at ¯ z = 0 for reasons we will explainsoon, and expand the differential equations around it, ∂ ¯ z ˜ M i (¯ z ) = A (0) ij (0 , ǫ )¯ z ˜ M j (¯ z ) . (5.6)Then we can formally write the limiting solution as˜ M i (¯ z ) = (cid:16) ¯ z A (0) (0 ,ǫ ) (cid:17) ij ˜ M j, . (5.7)To evaluate this matrix exponential it would be beneficial to diagonalise the matrix A (0) .However, in general, A (0) will not be diagonalisable. It is however possible to computethe Jordan decomposition of A (0) . The Jordan decomposition of any matrix A yields twomatrices R and J , so that A = RJ R − . (5.8)The matrix J is block diagonal with m blocks. Each block corresponds to an eigenvalue of A . The i th Jordan block J ( i ) corresponding to an eigenvalue λ i with geometric multiplicity n i has the dimension n i × n i , so that the n i diagonal entries each contain the eigenvalue.The elements on the first superdiagonal in each Jordan block are equal to 1. The diagonal– 7 –f J contains the eigenvalues of A , as it is the case for a diagonalised matrix. J = J (1) . . .
00 . . . 00 . . . J ( m ) , J ( i ) = λ i . . . λ i . . . λ i (5.9)The transformation matrix R consists of the generalised eigenvectors of A . In the casethat A is diagonalisable, all Jordan J i blocks have dimension n i = 1 and R consists of theeigenvectors of A , i.e. the Jordan decomposition diagonalises the matrix.To simplify the differential equations in the ¯ z → A (0) into R and J . Using the transformation matrix R we can then rotate the vector of masterintegrals ˜ M f i (¯ z ) = R ij ˜ M j , (5.10)and we find the simplified differential equation ∂ ¯ z f i (¯ z ) = J ij ( ǫ )¯ z f j (¯ z ) . (5.11)These differential equations permit the simple solutions ~f (¯ z ) = exp J (1) . . .
00 . . . 00 . . . J ( m ) log(¯ z ) ~f = e J (1) log(¯ z ) . . .
00 . . . 00 . . . e J ( m ) log(¯ z ) ~f , (5.12)with e J ( i ) log(¯ z ) = ¯ z λ i log(¯ z )1! . . . log ni − ¯ z ( n i − log ni − ¯ z ( n i − ...0 . . . . (5.13)We find therefore that the factors of ¯ z λ i , commonly referred to as integrating factors [15],appear together with log(¯ z ) raised to a power of maximally n i −
1. Note that the numberof unknown constants is still the dimension of the system of differential equations. Everyconstant f ,j is thus associated with exactly one eigenvalue λ i , while every λ i can beassociated with multiple f ,j .We can now express the limiting solutions of the original masters ˜ M i through thesolution for the f i that were just obtained.˜ M i (¯ z ) = R − ij f j (¯ z ) (5.14)We have managed to express the limiting solution of the master integrals through a linearcombination of constants f ,i , where each constant is associated with exactly one ¯ z λ i . Thisdecomposition considerably facilitates the further calculation of the unknown f ,i .– 8 –t is possible to arrive at the f ,i from a completely orthogonal point of view. Usingthe method of expansion by regions, the limiting solution of a master integral as ¯ z → regions each associated with a specific integrating factor ¯ z λ . Having identified theconstants f ,i contributing to a specific integral from the boundary decomposition, we cantherefore match the constants f ,i to the regions. Any boundary condition f ,i associatedwith an integrating factor that does not correspond to a region vanishes and therefore doesnot even require any explicit calculation of a Feynman integral.For example, analysing the RV V cross section using expansion by regions, we find onlyregions with integrating factors ¯ z a i − b i ǫ , with a i ∈ Z and b i ∈ { , , , , } . Therefore, onlyboundary conditions f ,i corresponding to λ i = a i − b i ǫ can be non-vanishing. All otherboundary constants appearing in the system are zero. Applying this boundary decomposi-tion dramatically reduces the number of boundary conditions that we needed to computefor the RV V master integrals from 72 to a mere 19.The remaining boundary conditions can be computed explicitly using expansion byregions. This step is also facilitated by the boundary decomposition, as one constant f ,i may appear in the limiting solution of more than one master integral. It is thereforereasonable to pick the simplest integrals to calculate the remaining constants.The actual computation is performed by deriving the integral representations of re-gions associated to the remaining boundary constants. This step is made especially viableby an algorithm exploiting a geometric interpretation of the parametric representation ofFeynman integrals as implemented in the code asy [18]. For a given integral and limit, asy provides a parameterisation for each region, which allows the direct expansion of theFeynman integral to obtain integral representation of the regions.In the case of the N LO Higgs production cross section a threshold expansion wasperformed and the “soft” master integrals appearing in these calculations may serve asboundary conditions for full kinematic integrals. Specifically the first term of the
RV V cross section was obtained in refs. [8] and in order to complete the full kinematic calculationwe could compare and confirm three of the boundary conditions given explicitly. However,we calculated 16 additional boundary conditions as described above. We observe thatall explicit logarithms arising from eigenvalues with geometric multiplicity larger than onevanish in the final result. We provide an explicit example of how we calculate our boundaryconditions in section 5.5.We want to stress that our algorithm of boundary decomposition can be based aroundany singular point in the differential equation. For example, we could have also calculatedthe limiting solutions for ¯ z →
1. Repeating the procedure with further singular points mayalso lead to additional constraints on the boundary conditions and therefore further reducethe number of integrals that actually have to be calculated. Furthermore, we would like topoint out that constraints on allowed eigenvalues, leading to non-vanishing f ,i can also beobtained from analyticity requirements and physical considerations [15]. In this section we discuss the method for solving differential equations of the type of– 9 –q.(5.4). In general a system of differential equations can be written as ∂ ¯ z M i (¯ z ) = A hij (¯ z, ǫ ) M j (¯ z, ǫ ) + y i (¯ z ) , (5.15)where A hij (¯ z, ǫ ) M j (¯ z, ǫ ) is the homogeneous part and y i (¯ z ) is the inhomogeneity that is zerounless a subset of master integrals has already been computed. In general, the homogeneoussolution is given by M hi (¯ z ) = (cid:16) e R d ¯ zA h (¯ z,ǫ ) (cid:17) ij M j, = H ij (¯ z, ǫ ) M j, , (5.16)where M i, is the boundary condition for master M i . Next, we need to find a particularsolution, which can depend on other master integrals. As y i (¯ z ) is known we find simply M pi (¯ z ) = H ij (¯ z, ǫ ) Z d ¯ zH − jk (¯ z, ǫ ) y k (¯ z ) , (5.17)such that the full solution can be written as M i (¯ z ) = M hi (¯ z ) + M pi (¯ z ) . (5.18)However, in general the differential equations are coupled and it is impossible to computethe matrix exponential in eq. (5.16). The desired result for our master integrals is a Laurentexpansion in the dimensional regulator. A commonly used strategy to calculate the abovematrix exponential is therefore to expand the differential equations in ǫ and decouple themorder by order.One particularly interesting version of this strategy has been proposed in ref. [30],which suggests that it is possible for Feynman integrals to find a transformation to acanonical basis such that the system takes the form ∂ ¯ z M ci (¯ z ) = ǫA cij (¯ z ) M cj (¯ z ) . (5.19)In this basis the ǫ dependence factorises completely from the system matrix. In this scenariothe inhomogeneity is zero. Furthermore, the system matrix takes the simple form A cij (¯ z ) = X k A c ( k ) ij ¯ z − ¯ z k , (5.20)where the matrices A c ( k ) have constant entries. The canonical form of the system matrix(eq. 5.20) makes the connection to multiple polylogarithms as defined in eq. (2.12) manifest.The formal solution of the canonical differential equations (eq. 5.19) can be written as M ci (¯ z ) = (cid:16) P e ǫ R d ¯ zA c (¯ z ) (cid:17) ij M cj, , (5.21)where P symbolises the path-ordered exponential and M c is a vector of boundary condi-tions. Expanding the exponential in ǫ we obtain M ci (¯ z ) = (cid:18) ǫ Z d ¯ zA cij (¯ z ) + ǫ Z d ¯ z (cid:18) A cik (¯ z ) Z d ¯ zA ckj (¯ z ) (cid:19) + . . . (cid:19) M cj, . (5.22)– 10 –t the time of our computation, no general algorithmic way to construct the transfor-mation that takes the master integrals to the canonical basis, was available . Obtainingthe canonical form is therefore a non trivial task. Some helpful insights were outlinedin [15, 32].To find a solution for differential equations the only necessary requirement is that thesystem can be sufficiently decoupled order by order in ǫ such that the matrix exponential ineq. (5.16) can be computed. While obtaining a canonical basis ensures this decoupling, thisis not the only basis that decouples the system. We choose to only transform a subsystemof 56 integrals of the complete system of 72 master integrals to the canonical basis. Theremaining 16 integrals can be easily computed using the general method.In this manner we obtain a solution for the full system depending on 72 constants ofintegration. By imposing the boundary decomposition of the limiting solution obtained inthe previous section, i.e. demanding that the full solution has the correct limit, we are ableto uniquely fix all constants. We thus calculated all 72 master integrals required for the RV V
Higgs boson cross section at N LO.
We discuss a short example, to illustrated how the methods described above proceed.Consider the integral topology T ( a , a , a ) = Z d d k (2 π ) d k − m ) a (( k + p ) − m ) a (( k + p + p ) − m ) a , (5.23)with p = p = 0 and p · p = s/
2. We choose as a basis of master integrals M = T (2 , , , M = xT (2 , , , M = ǫT (1 , , , (5.24)with x = q − m s and setting s = 1 for simplicity. With this we find the system ofdifferential equations ∂∂x ~M ( x ) = ǫ −
200 0 0 x + − x + − −
11 + x ~M ( x ) . (5.25)Next we analyse the system in the limit x →
1. This limit corresponds to the situationwhen the internal mass m is small compared to s . ∂∂x ~ ˜ M ( x ) = ǫ − x ~ ˜ M ( x ) . (5.26)Calculating the Jordan form J and the associated transformation matrix R of the systemmatrix yields J = − ǫ
10 0 − ǫ , R = − − − ǫ . (5.27) At the time of writing, an algorithm for finding a transformation to the canonical basis was publishedin [31] – 11 –he limiting solution in the Jordan basis and the original master integral basis is thus ~f ( x ) = e J log(1 − x ) ~f = f (1)0 (1 − x ) − ǫ f (2)0 + (1 − x ) − ǫ log(1 − x ) f (3)0 (1 − x ) − ǫ f (3)0 (5.28) ~ ˜ M ( x ) = − f (3)0 (1 − x ) − ǫ ǫ − f (3)0 (1 − x ) − ǫ ǫ − f (1)0 f (2)0 (1 − x ) − ǫ + f (3)0 (1 − x ) − ǫ log(1 − x ) + f (1)0 (5.29)The next step is to determine the f ( i )0 from expansion by regions. We start by determining f (3)0 . The easiest integral to compute f (3)0 is M . This integral, being a simple tadpole, isa one scale integral and as such has only one region. Fortunately, it is trivial to obtain thefull solution from the integral representation M = Z d d k (2 π ) d k − m , (5.30)and we obtain M = i (4 π ) − ǫ Γ( ǫ )(1 − x ) − ǫ (1 + x ) − ǫ . (5.31)The boundary condition is then obtained by comparing the leading term in the expansionaround x = 1 to eq. 5.29 and we have f (3)0 = − i (4 π ) − ǫ ǫ Γ( ǫ ) . (5.32)The next constant to be determined is f (1)0 . Here we choose the integral M . This integralis a massive bubble and contains the scales s and m . Analysing this integral using themethod of expansion by regions explicitly or using the code asy [18] for guidance one findsthree regions R (1)2 , R (3)2 and R (3)2 with the following scalings, R (1)2 ∝ (1 − x ) , R (2)2 ∝ (1 − x ) − ǫ and R (3)2 ∝ (1 − x ) − ǫ . (5.33)We see immediately that we do not need to compute R (3)2 as it is suppressed by one powerof (1 − x ) in comparison to the boundary conditions required. Furthermore, we knowfrom the boundary decomposition, eq. 5.29, that the region R (2)2 , proportional to (1 − x ) − ǫ corresponds to the boundary condition f (3)0 that we determined before. We therefore onlyneed to compute R (1)2 in order to obtain f (1)0 . The parametric representation of this regionis, R (1)2 = − i (4 π ) − ǫ Γ(1 + ǫ ) Z ∞ dx dx δ (1 − x − x ) x − ǫ x − − ǫ ( x + x ) − ǫ . (5.34)This integral can easily be solved in terms of beta functions and we obtain the boundarycondition f (1)0 from comparison with eq. 5.29, f (1)0 = − i (4 π ) − ǫ Γ(1 − ǫ ) Γ(1 + ǫ ) ǫ Γ(1 − ǫ ) . (5.35)– 12 –he final boundary condition f (2)0 is obtained from integral M , which has three regions R (1)3 , R (3)3 and R (3)3 with the scalings, R (1)3 ∝ (1 − x ) , R (2)3 ∝ (1 − x ) − ǫ and R (3)3 ∝ (1 − x ) − ǫ . (5.36)Here we observe a small subtlety in the computation. This integral has two regions withthe same scaling that are not suppressed relative to one another. It will therefore benecessary to compute both of them. Furthermore, the logarithm that appears in theboundary decomposition eq. 5.29, suggests that these regions will have a divergence thatis not regulated by dimensional regularisation when they are computed separately. Thisis immediately confirmed when we derive the integral representation for R (2)3 or R (3)3 . Wetherefore introduce an analytic regular ν so that the Feynman integral for M becomes M ′ = Z d d k (2 π ) d k − m )(( k + p ) − m )(( k + p + p ) − m ) ν . (5.37)Starting from this regularised integral we can perform expansion by regions as before andwe obtain the parametric representation for R (2)3 as function of ν , R (2)3 ( ν ) = ( − ν i (4 π ) − ǫ − − ǫ + ν ǫ (1 − x ) − ǫ − ν Γ(1 + ǫ + ν )Γ(1 + ν ) × Z ∞ dx dx dx δ (1 − x − x − x ) x ν ( x + x ) − ǫ + ν (cid:0) x x + ( x + x ) (cid:1) − − ǫ − ν (5.38)Performing the integrals over the parameters x i as beta functions we find R (2)3 ( ν ) = ( − ν i (4 π ) − ǫ − ǫ + ν ǫ (1 − x ) − ǫ − ν Γ( ǫ + ν ) ν Γ(1 + ν ) . (5.39)Here we see how the regulator ν regulates the divergence and we can not take the limit ν → R (3)3 with the regulator.This region has the parametric representation R (3)3 ( ν ) = ( − ν i (4 π ) − ǫ − − ǫ + ν ǫ (1 − x ) − ǫ Γ(1 + ǫ + ν )Γ(1 + ν ) × Z ∞ dx dx dx δ (1 − x − x − x )( x + x ) − ǫ + ν x ν (cid:0) ( x + x ) + 2 x x (cid:1) − − ǫ − ν . (5.40)Once again we can perform the parametric integrals in terms of beta functions and find R (3)3 ( ν ) = ( − ν i (4 π ) − ǫ − ǫ ǫ (1 − x ) − ǫ Γ( ǫ ) ν . (5.41)Also here we can see the singularity being regulated by ν . We can however combine bothregions and take the limit ν → ν → (cid:16) R (2)3 + R (3)3 (cid:17) = − i (4 π ) − ǫ − ǫ ǫ (1 − x ) − ǫ Γ( ǫ ) ( γ E + log(2) + ψ ( ǫ ) − log(1 − x )) , (5.42)– 13 –ith ψ ( x ) = d log(Γ( x )) dx and γ E = − ψ (1). Here we can see the explicit log(1 − x ) thatwas predicted by the boundary decomposition. If we compare the term proportional tolog(1 − x ) with eq. 5.29 we can confirm that it in fact corresponds to f (3)0 as predicted.The remaining boundary condition f (2) is then obtained as, f (2)0 = lim ν → (cid:16) R (2)3 + R (3)3 (cid:17)(cid:12)(cid:12)(cid:12) log(1 − x ) = − i (4 π ) − ǫ − ǫ ǫ Γ( ǫ ) ( γ E + log(2) + ψ ( ǫ )) . (5.43)With this we have determined the last remaining boundary condition. The complete systemcan now be obtained trivially by solving eq. (5.25) and demanding consistency with theabove boundary conditions. To outline our method of calculating the actual boundary conditions, we show the exampleof the double cut of the tennis court diagram, depicted in figure 1, which serves as aboundary condition to our system of differential equations. This integral was first calculatedin [8]. The momentum space representation of the integral is,
12 12
Figure 1:
The two particle cut of the tennis court which serves as a boundary condition. Z d d k (2 π ) d d d l (2 π ) d k l ( k + p ) ( l − p ) ( l + p ) ( k + p ) ( k + l + p ) − s , (5.44)with p = p − p and p = p + p − p . By using our method for decomposing theboundary conditions as outlined in section 5.2 we obtain the different boundary conditionscontributing to this integral, − ǫ + 1)(5 ǫ + 2) B ¯ z − − ǫ ǫ (2 ǫ + 1) (3 ǫ + 1) + 4(4 ǫ + 1) B ¯ z − − ǫ ǫ (2 ǫ + 1)(3 ǫ + 1) + 8(6 ǫ + 1) B ¯ z − − ǫ ǫ (3 ǫ + 1) + 12 B ¯ z − − ǫ ǫ (3 ǫ + 1)+ 9 B ¯ z − − ǫ ǫ (3 ǫ + 1) + 4(6 ǫ + 1) B ¯ z − − ǫ ǫ + (5 ǫ + 1) B ¯ z − − ǫ ǫ ( ǫ + 1)(3 ǫ + 1) − ǫ + 1) B ¯ z − − ǫ ǫ (2 ǫ + 1)(3 ǫ + 1) − ǫ + 1) B ¯ z − − ǫ ǫ (2 ǫ + 1) + B ¯ z − − ǫ ǫ (3 ǫ + 1) + B ¯ z − − ǫ (cid:18) ǫ + 14 ǫ (3 ǫ + 1) − ¯ z ǫ + 112 ǫ (3 ǫ + 1) (cid:19) . (5.45)The boundary condition that we want to determine here is B , all other boundary condi-tions can be determined independently from other integrals. To leading power in ¯ z , B isthe only boundary condition contributing. Therefore, we need to compute the region pro-portional to ¯ z − − ǫ of this integral. Using expansion by regions we can derive a momentum– 14 –pace representation or use the code asy [18] to obtain a parametric representation of therequired region: I ≡ (4 π ) − ǫ ((1 + 6 ǫ ) B (4 ǫ (1 + 3 ǫ )) (5.46)= (4 π ) − ǫ Z d Φ d d k (2 π ) d d d l (2 π ) d k l ( k − l ) (2 kp − p p )(2 lp − p p )(2 lp )( k − p ) s . Introducing Feynman parameters and transforming to projective space we obtain, I = Z d Φ Z ∞ dx dx dx dx dx dx Γ(3 + 2 ǫ ) 1 s ( x + x + x (1 + x + x )) ǫ × (cid:16) s (( x + x ) x + x ( x + x x + x + x x + x x ))+ x ( s x + s ( x + x (1 + x + x ))) (cid:17) − − ǫ . (5.47)The integration over x can be performed immediately. Using the projective transformation x → x x , the integral over x can be computed as well and we obtain, I = Z d Φ Z dx dx dx dx Γ(1 − ǫ )Γ(2 ǫ )Γ(2 + 2 ǫ )¯ z − − ǫ s ǫ s − − ǫ s − − ǫ × x − − ǫ (1 + x + x + x ) ǫ ( x + x + x (1 + x + x )) ǫ × ( x (1 + x + x ) + (1 + x )(1 + x ) x ) − − ǫ . (5.48)Next we split the second polynomial intoΓ( − z )Γ( − z − ǫ )Γ( − − ǫ ) x ǫ − z ( x + x ) z (1 + x + x ) − z +3 ǫ , (5.49)by introducing a Mellin-Barnes integral over z , such that we can perform the integral over x . After performing the projective transformations x → x x and x → x x we obtain, I = Z γ dz Γ( − ǫ )Γ(1 + 2 ǫ )Γ( − − ǫ ) Γ(2 + 3 ǫ − z )Γ( − z )Γ( − − ǫ + z )Γ( − ǫ + z ) × Z d Φ Z dx dx dx s ǫ s − − ǫ s − − ǫ x − ǫ (1 + x ) ǫ − z (1 + x ) ǫ − z (1 + x + x ) ǫ x ǫ − z × (1 + x ) ǫ (1 + x + x + (1 + x )(1 + x ) x ) − − ǫ + z , (5.50)where the contour γ is such that the poles of gamma functions with − z i in the argument(left poles) and the poles of gamma functions with + z i in the argument (right poles)are separated. Next, we introduce a second Mellin-Barnes integration to split the lastpolynomial into,Γ( − z )Γ(2 + z + 3 ǫ − z )Γ(2 + 3 ǫ − z ) x − − z − ǫ + z (1+ x ) − − z − ǫ + z (1+ x ) − − z − ǫ + z (1+ x + x ) z . (5.51)– 15 –ow we can perform the integral over x , x and x in that order and obtain I = Z d Φ s ǫ s − − ǫ s − − ǫ Γ( − ǫ )Γ(1 + 2 ǫ )Γ( − − ǫ )Γ( − − ǫ ) Z γ dz πi dz πi Γ( − z )Γ( − z ) × Γ( − − ǫ + z )Γ( − ǫ + z )Γ( − − ǫ − z )Γ( − ǫ + z )Γ(2 + 3 ǫ − z + z )(1 + 2 ǫ + z )Γ(1 + z ) × (cid:16) Γ( − ǫ )Γ(1 + z ) − Γ( − ǫ )Γ(1 + ǫ + z ) (cid:17) , (5.52)as the final Mellin-Barnes representation. Next, we need to perform the phase space inte-gral. We insert the appropriate parametrisation s = − i , s = λ, s = 1 − λ, (5.53)with λ ∈ [0 , i
0. Afterwards, we can perform the phasespace integral as a simple beta function.The contour of the Mellin-Barnes integration is defined by the requirement that itshould separate the left and right poles of the integrand. At this point, this is only satisfiedfor the integral, if ǫ is finite. Therefore, in order to be able to expand the integral in ǫ , beforethe Mellin-Barnes integration is performed, we need to analytically continue the integralto infinitesimal ǫ . This is achieved using the residue theorem, by taking the residues ofpoles that end up on the wrong side of the contour when ǫ is gradually taken to zero. Thisis automated in codes like MB and MBresolve [33]. After the analytic continuation, theintegral can be expanded in ǫ . We refrain from printing the unwieldy expansion that isobtained in this step. Afterwards, we can apply Barnes’ lemma and corollaries thereof toeliminate one of the two integrations and we are left with a one dimensional Mellin-Barnesintegral. This one dimensional integral can be easily computed by taking the residues of,e.g. the left poles of the integrand, which yields a sum representation. These sums can beperformed in terms of harmonic sums [34], which yield multiple zeta-values when evaluatedat infinity. This way we find the final result ℜ ( I ) e ǫγ E = 13 ǫ − ǫ ζ − ǫ ζ + 25716 ǫ ζ + (cid:18) ζ ζ − ζ (cid:19) + ǫ (cid:18) ζ − ζ (cid:19) + O ( ǫ ) . (5.54)Our method of solving the integrals in the Mellin-Barnes representation also provides a wayto cross check the result as the Mellin-Barnes integrals can also be evaluated numerically.A large fraction of the required boundary conditions for the RV V cross section can beobtained in a simpler fashion. For other integrals we proceed similar to the above example.
6. Conclusions
In this article we have concluded the computation of all genuine 2-loop contributions withone extra parton in the final state to the N LO Higgs boson production cross section.– 16 –his constitutes a further important step towards the computation of the full Higgs bosonproduction cross section at N LO.We have calculated numerous single real emission phase space integrals over two-loopamplitudes. We have performed these integrations using the method of differential equa-tions which we have advanced with a new method for the decomposition of boundaryconditions. This method makes the use of differential equations for inclusive phase spaceintegrals viable and makes an important connection between differential equations and themethod of expansion by regions manifest. Our method will also be a prime ingredientin the calculation of the remaining pieces of the Higgs cross section at N LO, namelythe double-real virtual piece, which requires the calculation of double-real emission phasespace integrals over one-loop amplitudes, as well as the triple-real piece, which requires thecalculation of triple-real emission phase space integrals.The results obtained in this paper have already contributed to the calculation of thenext-to-soft approximation of the Higgs production cross section, as well as to the determi-nation of the three leading threshold logarithms in general kinematics, at N LO [13] andthus have an immediate phenomenological impact.We provide our results in a separate file included with the arXiv submission of thispaper. The file contains all bare, partonic
RV V cross sections as defined in eq. (3.1) dividedby a factor given in eq. (3.3).In parallel to this calculation the
RV V
Higgs boson cross section was also computedin [35] in agreement with our result. The main difference between the two calculations isthat the authors of [35] did not employ reverse unitarity to combine the loop and phase-space integrals but rather explicitly computed phase-space integrals over two-loop ampli-tudes [36].
Acknowledgements
We are grateful to Claude Duhr, Thomas Gehrmann and Matthieu Jaquier for usefuldiscussions and comparison of closely related results prior to publication. We thank BabisAnastasiou for discussions and feedback on the manuscript. This work has been supportedby the Swiss National Science Foundation (SNF) under contract 200021-143781 and theEuropean Commission through the ERC grant “IterQCD”.
References [1] C. Anastasiou and K. Melnikov, “Higgs boson production at hadron colliders in NNLOQCD,” Nucl. Phys. B , 220 (2002), [hep-ph/0207004];R. V. Harlander and W. B. Kilgore, “Next-to-next-to-leading order Higgs production athadron colliders,” Phys. Rev. Lett. , 201801 (2002), [hep-ph/0201206];V. Ravindran, J. Smith and W. L. van Neerven, “NNLO corrections to the total cross-sectionfor Higgs boson production in hadron hadron collisions,” Nucl. Phys. B , 325 (2003),[hep-ph/0302135].[2] C. Anastasiou, S. Buehler, F. Herzog and A. Lazopoulos, JHEP , 004 (2012)[arXiv:1202.3638 [hep-ph]]. – 17 –
3] S. Buehler and A. Lazopoulos, [arXiv:1306.2223 [hep-ph]].[4] T. Gehrmann, E. W. N. Glover, T. Huber, N. Ikizlerli and C. Studerus, JHEP (2010)094 [arXiv:1004.3653 [hep-ph]], P. A. Baikov, K. G. Chetyrkin, A. V. Smirnov, V. A. Smirnovand M. Steinhauser, Phys. Rev. Lett. (2009) 212002 [arXiv:0902.3519 [hep-ph]].[5] S. Moch, J. A. M. Vermaseren and A. Vogt, Nucl. Phys. B , 101 (2004); Nucl. Phys. B , 129 (2004), C. Anastasiou, S. B¨uhler, C. Duhr and F. Herzog, JHEP , 062 (2012);M. H¨oschele, J. Hoff, A. Pak, M. Steinhauser, T. Ueda, Phys. Lett. B , 244 (2013);[6] O. V. Tarasov, A. A. Vladimirov and A. Y. .Zharkov, Phys. Lett. B , 429 (1980);S. A. Larin and J. A. M. Vermaseren, Phys. Lett. B , 334 (1993); T. van Ritbergen,J. A. M. Vermaseren and S. A. Larin, Phys. Lett. B , 379 (1997); M. Czakon, Nucl. Phys.B , 485 (2005).[7] C. Anastasiou, C. Duhr, F. Dulat and B. Mistlberger, JHEP , 003 (2013)[arXiv:1302.4379 [hep-ph]].[8] C. Duhr and T. Gehrmann, Phys. Lett. B (2013) 452 [arXiv:1309.4393 [hep-ph]], Y. Liand H. X. Zhu, JHEP (2013) 080 [arXiv:1309.4391 [hep-ph]].[9] C. Anastasiou, C. Duhr, F. Dulat, F. Herzog and B. Mistlberger, JHEP (2013) 088[arXiv:1311.1425 [hep-ph]], W. B. Kilgore, Phys. Rev. D (2014) 073008 [arXiv:1312.1296[hep-ph]].[10] C. Anastasiou, C. Duhr, F. Dulat, E. Furlan, T. Gehrmann, F. Herzog and B. Mistlberger,Phys. Lett. B (2014) 325 [arXiv:1403.4616 [hep-ph]].[11] Y. Li, A. von Manteuffel, R. M. Schabinger and H. X. Zhu, Phys. Rev. D (2014) 053006[arXiv:1404.5839 [hep-ph]].[12] M. Bonvini, R. D. Ball, S. Forte, S. Marzani and G. Ridolfi, J. Phys. G (2014) 095002[arXiv:1404.3204 [hep-ph]], M. Bonvini and S. Marzani, JHEP (2014) 007[arXiv:1405.3654 [hep-ph]], D. de Florian, J. Mazzitelli, S. Moch and A. Vogt,arXiv:1408.6277 [hep-ph].[13] C. Anastasiou, C. Duhr, F. Dulat, E. Furlan, T. Gehrmann, F. Herzog and B. Mistlberger,To appear soon.[14] C. Anastasiou, L. J. Dixon and K. Melnikov, Nucl. Phys. Proc. Suppl. , 193 (2003)[hep-ph/0211141]; C. Anastasiou, L. J. Dixon, K. Melnikov and F. Petriello, Phys. Rev. Lett. , 182002 (2003) [hep-ph/0306192]; C. Anastasiou, L. J. Dixon, K. Melnikov andF. Petriello, Phys. Rev. D , 094008 (2004) [hep-ph/0312266].[15] T. Gehrmann and E. Remiddi, Nucl. Phys. B , 485 (2000) [hep-ph/9912329];T. Gehrmann and E. Remiddi, Nucl. Phys. B (2001) 248 [hep-ph/0008287].[16] A. B. Goncharov, Math. Res. Lett. , 497 (1998) [arXiv:1105.2076 [math.AG]],A. B. Goncharov, math/0103059 [math.AG].[17] M. Beneke and V. A. Smirnov, Nucl. Phys. B (1998) 321 [hep-ph/9711391].[18] B. Jantzen, A. V. Smirnov and V. A. Smirnov, Eur. Phys. J. C (2012) 2139[arXiv:1206.0546 [hep-ph]]. A. Pak and A. Smirnov, Eur. Phys. J. C (2011) 1626[arXiv:1011.4863 [hep-ph]].[19] E. Furlan, JHEP , 115 (2011) [arXiv:1106.4024 [hep-ph]]. – 18 –
20] C. Anastasiou, R. Boughezal and E. Furlan, JHEP , 101 (2010) [arXiv:1003.4677[hep-ph]].[21] P. A. Baikov and K. G. Chetyrkin, Phys. Rev. Lett. , 061803 (2006) [hep-ph/0604194].[22] K. G. Chetyrkin, B. A. Kniehl and M. Steinhauser, Phys. Rev. Lett. , 353 (1997)[hep-ph/9705240].[23] M. Kramer, E. Laenen and M. Spira, Nucl. Phys. B , 523 (1998) [hep-ph/9611272].[24] C. Duhr, “Hopf algebras, coproducts and symbols: an application to Higgs bosonamplitudes,” JHEP , 043 (2012) [arXiv:1203.0454 [hep-ph]].[25] F. Brown, arXiv:1102.1310 [math.NT].[26] P. Nogueira, J. Comput. Phys. , 279 (1993).[27] C. W. Bauer, A. Frink and R. Kreckel, “Introduction to the GiNaC framework for symboliccomputation within the C++ programming language,” cs/0004015 [cs-sc].[28] K. G. Chetyrkin and F. V. Tkachov, Nucl. Phys. B (1981) 159; F. V. Tkachov, Phys.Lett. B (1981) 65.[29] S. Laporta, Int. J. Mod. Phys. A , 5087 (2000) [hep-ph/0102033].[30] J. M. Henn, Phys. Rev. Lett. (2013) 25, 251601 [arXiv:1304.1806 [hep-th]].[31] R. N. Lee, arXiv:1411.0911 [hep-ph].[32] M. Argeri, S. Di Vita, P. Mastrolia, E. Mirabella, J. Schlenk, U. Schubert and L. Tancredi,JHEP (2014) 082 [arXiv:1401.2979 [hep-ph]]; J. M. Henn, A. V. Smirnov andV. A. Smirnov, JHEP (2014) 088 [arXiv:1312.2588 [hep-th]]; M. Hschele, J. Hoff andT. Ueda, JHEP (2014) 116 [arXiv:1407.4049 [hep-ph]].[33] M. Czakon, Comput. Phys. Commun. , 559 (2006) [hep-ph/0511200], A. V. Smirnov andV. A. Smirnov, Eur. Phys. J. C , 445 (2009) [arXiv:0901.0386 [hep-ph]], C. Anastasiou andA. Daleo, JHEP , 031 (2006) [hep-ph/0511176].[34] J. A. M. Vermaseren, Int. J. Mod. Phys. A , 2037 (1999) [hep-ph/9806280], J. Blumleinand S. Kurth, Phys. Rev. D , 014018 (1999) [hep-ph/9810241].[35] C. Duhr, T. Gehrmann, M. Jaquier To Appear soon.[36] T. Gehrmann, M. Jaquier, E. W. N. Glover and A. Koukoutsakis, JHEP (2012) 056[arXiv:1112.3554 [hep-ph]].(2012) 056[arXiv:1112.3554 [hep-ph]].