Tops from Light Quarks: Full Mass Dependence at Two-Loops in QCD

I present the two-loop QCD corrections to the production of a massive quark-anti-quark pair in the massless quark-anti-quark annihilation channel. The result is obtained as a combination of a deep expansion in the mass around the high energy limit and of a numerical integration of a system of differential equations. The primary application of the outcome and developed methods is top quark pair production at the Large Hadron Collider.


Introduction
One of the most important goals of the Large Hadron Collider (LHC) is the measurement of top quark properties. This will be possible thanks to the abundant statistical samples reaching about 8 million pairs produced per year in the low luminosity phase. Besides the mass and decay parameters, the total cross section constitutes a primary observable. The experimental prospects of obtaining a measurement accuracy below ten percent for this quantity put a high demand on theoretical predictions of matching quality. At present, the known next-to-leading order corrections [1] have an error estimated from scale variation at about 12%. Soft gluon resummation [2][3][4], which has been an excellent tool for the Tevatron, and helped reduce the error to about 5%, is not safely applicable in the framework of the LHC. This is due on the one hand to the higher energy and on the other to the dominance of the gluon flux over the quark flux. Furthermore, the mentioned high statistics warrant the preparation of a Monte-Carlo generator of suitable precision, which cannot rely on resummed cross sections.
In view of these facts it is necessary to provide a result for the next-to-next-to-leading order cross section, at best in a fully differential form. This requires four separate parts at the partonic level: the two-loop virtual corrections, the one-loop squared virtual corrections, the one-loop real-virtual corrections with an additional parton in the final state, and the tree-level real corrections with two additional partons in the final state. Out of these, the second part is known from [5,6], the third from [7], where the next-to-leading order corrections to the tt+jet corrections have been derived. Unfortunately, as part of a cross section calculation for top quark pair production this result still needs subtraction terms in order to allow for integration over the full phase space. Finally, while there is no result for the real radiation, the two-loop corrections 1 have been recently evaluated in the limit of small top quark mass [14,15]. This result is applicable for highly energetic tops, for example when high p T cuts would be applied. The bulk of events comes, however, from the region much nearer to the partonic threshold.
In this paper, I present a complete result for the two-loop corrections in the quark annihilation channel valid in the whole kinematically allowed region. It has to be stressed that obtaining an amplitude expressed in analytic form through some special functions seems out of reach in the nearest future. Since the LHC will soon become operational it is necessary to resort to semi-analytic/semi-numeric methods. The method adopted here is a combination of a deep expansion in the mass around the high energy limit, which contains the power corrections to the result of [14], and of numerical integration of differential equations.
In the next section, I will first give a few definitions and then describe the power corrections. A detailed study of the numerical methods and the full result will follow in the last section of the main text.

Power corrections
The notation of this paper follows closely that of [14]. However, I reproduce all the necessary definitions for the convenience of the reader.
The process under scrutiny is massless quark-anti-quark annihilation into a massive quark-anti-quark pair The amplitude can be described with the help of Mandelstam variables Notice that the mass subtraction in the definition of the t and u variables was irrelevant for the results of [14], because only logarithmic terms in the mass have been retained there. The advantage of this definition lies in the symmetric range of variation where β is the velocity The results will be parameterized by two dimensionless ratios The additional scale introduced by dimensional regularization, µ, has been kept in the results as unspecified, but for the plots and numerical values reproduced in the following µ = m. The renormalization has been performed in the MS scheme with n l massless and n h massive active flavors. The necessary constants are known in the literature with sufficient precision: the strong coupling renormalization to four-loop accuracy [16,17], and the mass and field renormalization of the heavy quark in the on-shell scheme to three-loop accuracy [18][19][20]. The renormalization of the light quark field, which is non-vanishing because of the presence of closed heavy quark loops has been explicitely given in [14] with two-loop accuracy. After expanding the amplitude in the strong coupling constant up to the second order the interesting, two-loop term, contracted with the Born amplitude can be decomposed into color factors The leading behavior of the amplitude in the limit m → 0, has been derived in [14] using two different approaches. The first is based on factorization properties of QCD amplitudes [21], and exploits a relation between the massless and massive cases. Unfortunately, it does not give a handle on mass corrections or the full mass dependence. The second approach is based on a direct evaluation of occurring integrals and is an evolution of a strategy developed for Bhabha scattering [22,23]. The procedure starts with a reduction of the integrals to masters with the help of the Laporta algorithm [24], and subsequent expansion of the masters in the mass by passing through Mellin-Barnes representations [25,26]. The bulk of the work is performed by Mathematica packages MBrepresentation [27] and MB [28] together with further associated software. The derivation of the asymptotic behavior of Mellin-Barnes representations is performed recursively by closing the contours of integration and taking residues. It is obvious that arbitrarily high orders of expansion in the mass can be obtained. However, the coefficients will still be integrals requiring evaluation by summation of infinite series or by some other method. Previously, this last step has been completed with a combination of XSummer [29] and the PSLQ algorithm [30]. It has to be stressed that every next order in ǫ contains more integrals with a more complicated integrand structure. Fortunately, there exists an alternative approach based on differential equations [31,32].
Clearly, applying a derivative with respect to any invariant or the mass introduces higher powers of denominators and/or numerator structures. Furthermore, any set of integrals differing only by powers of denominators and numerators can be reduced to a smaller set of masters. In consequence one can write the following differential equation systems for the coefficients of the Laurent expansion of the master integrals The Jacobian matrices, J M and J X , have rational function elements and it is implied that any master integral is a combination of the I i (m s , x) functions The lowest power of ǫ in the sum is fixed by the singularities of the integral and cannot exceed −4 at this level of perturbation theory, whereas the highest power is defined by the coefficient in the amplitude (there are spurious poles). In practice, the deepest expansion in ǫ that occurred was down to order ǫ 3 due to the particular choice of master integrals. The differential equations Eq. 8 allow, in principle, to fix the complete functional dependence of the master integrals. Unfortunately, the functions present in the solution are not known at present, and therefore a direct integration of the system has to be postponed. Nevertheless, the differential equations in m s can be solved by means of a series expansion, with boundaries given by the small mass limit as needed for the results of [14]. Following this idea, it is possible to derive arbitrarily deep expansions of the amplitude in the mass, and thus provide the power corrections, which were out of reach of the factorization approach. Of course, the size of the intermediate expressions, combined with the available computing resources puts a natural cut-off on the highest power that can be computed. In fact, I have computed eleven terms of the series up to m 10 s . The results for the finite parts (in the ǫ expansion) of the three purely bosonic contributions A, B and C, which are also the most involved as far as the computation is concerned are shown in Fig. 1. The plots correspond to 90 degree scattering, i.e. x = 1/2. It is striking that the series do not obviously diverge, which is usually the case with this type of expansions. In fact, the series expansion for the leading color term, A, is at worst asymptotic at threshold, and can still be used to obtain an estimate accurate to a few percent at this point. On the other hand, the growth of the subleading color coefficients is indicative of the true behavior, but incorrect. As I will show in the next section, there is a true divergence due to the Coulomb singularity, which cannot be reproduced with a small mass expansion.  Figure 1: Finite parts of the purely bosonic contributions to the two-loop amplitude as expansions around the small mass limit at x = 1/2. The solid line represents the eleven terms of the series as derived for the present publication, the long dashed -ten terms of the series, the short dashed -six terms of the series, the dash-dotted -the leading behavior.
The small mass expansion is, in reality, not an expansion in m 2 /s, but rather in max(m 2 /s, −m 2 /t, −m 2 /u). For small m and at the kinematical boundary corresponding to forward scattering A similar relation holds for −m 2 /u for backward scattering. In consequence, the series will be asymptotic at best at the kinematic boundaries.
In order to study the region of convergence and usability of the series, it is necessary to specify some criteria that could be applied without reference to any external result. In fact, if the amplitude were approximated with one permille accuracy, it would be sufficient for any foreseeable applications. Customarily, the error of an expansion is estimated by the size of the last term. For geometric series, this estimate is only correct (not underestimated) if the ratio of two subsequent terms does not exceed 1/2. In the latter case, eleven terms The grey area represents the region where the series is accurate to one permille, whereas the black area the region where the leading term of the series is accurate to one percent. The dashed lines delimit the kinematically allowed region, whereas the short dashed lines inside the grey area would correspond to the convergence region derived according to Eq. 13 without the last condition.
of the series (as in the case of the present result for the amplitude) provide indeed an approximation exact to one permille. In the case of amplitudes, the series is obviously not geometric, but conditions inspired by these arguments can be imposed. Let the amplitude be written as A relatively conservative heuristic test for the one permille convergence of the result is given by the following conditions which are applied at a given (m s , x) point during a scan starting from the median x = 1/2, which has always the best convergence for a given value of m s . The last condition in Eq. 13 deserves further explanation. It turns out that for relatively small values of m s , the logarithmic terms in m s lead to slight violations of the remaining relations, but the last term of the series is still tiny. In this case, it is highly improbable that the sum of the missing terms would amount to more than a permille correction. Without the last test the region of permille convergence is, therefore, unrealistically restricted. The region resulting from the application of Eq. 13 is shown in Fig. 2. The visible discontinuities of the boundary are due precisely to the nature of the criterion (and to a lesser extent to the discretized scan). This figure shows also the region, where the leading term of the series agrees with the full result to one percent. No special criteria are needed here, of course, since eleven terms of the series expansion provide a sufficient approximation to the exact result for the purpose of determining this region.

Numerical solution
Since the series expansion does not satisfactorily approximate the result over the whole range of variation of kinematic parameters, it is only natural to try to solve the differential equations for the master integrals numerically and thus obtain a purely numerical description of the amplitude. This idea has originally been put forward in [33] for master integrals corresponding to the two-loop sunrise graph without, however, relating to a concrete physical problem. In fact, to the best of my knowledge, the only applications to physical processes have been attempted in [34,35]. The problem at hand pushes the difficulty level substantially further, and requires, therefore, a careful assessment of feasibility. Before I show that high precision numerical results may indeed be obtained at all relevant points of phase space, there are several issues that have to be clarified, as it can hardly be overstressed that the right choice of numerical algorithms will make the difference between success and failure.

Implementation
At first, it is necessary to determine the type of differential equation system under consideration. It turns out that different methods are available for stiff and non-stiff problems. The distinction between the two is somewhat fuzzy in the professional literature, but there is agreement that stiff problems involve exponentially decaying solutions. The criterion used in practice is the existence of large negative eigenvalues of the Jacobian. In the present case it is, however, easier to use heuristic arguments instead of performing a numerical analysis.
In fact, experience accumulated in numerous higher order calculations shows that exponentially decaying components would be rather unusual. I will, therefore, assume without further consideration that the system is non-stiff. In this case, there is still a large number of algorithms available. However, because the solutions, i.e. the master integrals, must be very smooth (we remain above all thresholds) and high precision will be requested, a variable coefficient multistep method [36] is expected to be most efficient [37]. Indeed, this kind of methods is based on polynomial interpolation/extrapolation with polynomials of order up to 12, which is a guarantee of very fast convergence under the assumption of suitable smoothness (if higher order derivatives were large, the errors would obviously grow severely with the order). The next choice concerns the treatment of master integrals, as they can be considered to be either complex or two-component real functions (after decomposition into real and imaginary parts). Clearly, working with complex functions reduces the size of the Jacobian, which may give substantial improvements in the running time. Furthermore, the real function approach suffers from the fact that the imaginary parts of many integrals vanish for real arguments, which poses problems as far as error control is concerned. Indeed, in such cases it is impossible to use a relative error criterion and absolute errors must be used. The size of the latter can only be determined in conjunction with the final amplitude in order to have control over the precision of the outcome. There seems to be only one, albeit very strong, argument in favor of real functions, namely that most of the available software works with real numbers exclusively. A quick glance at the literature of the subject shows, that writing a code from scratch should be avoided. Fortunately, one of the most advanced software packages implementing the variable coefficient multistep method [38] has recently been translated to complex arguments.
Closely connected to the choice made above is the problem of error control. Customarily, working with complex functions implies that the error is given by the modulus of the difference between the exact value and the approximation. In the present case, however, we are only interested in the real part of the amplitude, hence the imaginary components of the master integrals will be discarded. In consequence, we need to control the error of the real part and not that of the modulus. Fortunately, unless the imaginary part is much larger than the real, the two are not much different. It turns out that in the present calculation, only about 6% of the evaluated phase space points involved an integral, for which the imaginary part was more than 10 4 larger than the real part. Therefore, for simplicity reasons, the traditional error estimate has been used in the following. Notice also that the imaginary parts could have been discarded from the start, since the system of differential equations is linear and we do not cross any singularities.
After settling the implementation questions, it is necessary to decide on the position of the boundary conditions. In the original publication [33], it has been proposed to start from a threshold or a pseudo-threshold, since the values of the integrals at these points were known. This approach has the drawback that these points are at the same time singularities of the differential equations, which requires slight modifications of the algorithm and leads invariably to a substantial loss of precision when evolving further from the boundary. Here, I use the series expansion of the previous section to compute the values of the integrals to very high precision. In fact at the relative error estimated by the size of the last term (very conservative) does never exceed 10 −18 . The second condition in Eq. 14 deserves explanation, because the median point x = 1/2 would have led to better convergence. Unfortunately, we will see later that it is also a singular point of the system and thus cannot be used. The choice taken results in the loss of about two digits and is compensated by a twice smaller value of m s . Let me now turn to the discussion of the contour of integration. It is clear [33] that evolving along the real axis should be avoided, in order not to stumble on the singularities. Since we are working with complex functions anyway, it is not a problem to deform the  contour into the complex plane. In general, this deformation is only restricted by causality 2 . In the present case, however, we always remain above thresholds, which means that the Riemann sheet has already been chosen when fixing the boundary conditions, and any curve will be appropriate. To take full advantage of the multistep algorithm for integration, which depends on several previous values of the system, it would be desirable to use a single smooth curve to reach a given point starting from the boundary. Unfortunately, we need an integration in two variables, and experimentation has shown that it is more efficient to perform two separate evolutions. For each of these, I will use an elliptic contour, with a user specified eccentricity as depicted in Fig. 3. The latter freedom allows for a relatively easy estimate of the final global error, by computing the desired amplitude with several different contours. It is also interesting to note that in practice the convergence turns out to be faster for more circular contours.
A rather unpleasant feature of the system of differential equations at hand is the presence of a number of singularities in the Jacobians. They are summarized in Tab. 1. It is crucial that aside from thresholds the solution is regular at these points. Therefore, in the case of the four different singularities, which occur inside the kinematically allowed region it is sufficient to resort to interpolation. The problem is more acute for the singularities at the boundary (forward/backward scattering). These approach the true branching points at the thresholds when m → 0, and therefore interpolation is only efficient for moderate values of m, when the necessary points outside the kinematically allowed region are not trapped too close to the branching points. Otherwise, extrapolation is necessary. Notice that the concept of "dangerous distance to a singularity" requires specification. In fact, it is defined by the available numerical precision and the strength of the singularity (power of the singular polynomial in the denominator of a coefficient). Finally, one has to remember that the solution for the master integrals will be input into the expression for the am-Jacobian singularity branching allowed interpretation m s = 0 yes collinear singularity m s = 1/4 yes s-channel threshold m s = −1/4 x = 0 yes t-channel threshold x = 1 yes u-channel threshold x = 1/2 yes perpendicular scattering Table 1: Complete list of singularities of the Jacobians, J M and J X , of the system of differential equations. The singularities occur in both systems of differential equations in m s and x apart from the point m s = −1/4, which is present only in the differential equations in m s . The table indicates in addition the presence of a branching point at a given singularity (a blank entry denotes a regular point of the solution), and specifies whether a singularity occurs within the kinematically allowed region (a blank entry denotes a point outside the allowed region).
plitude, which may lead to further cancellations, and hence instabilities. For the present calculation, I have adopted extended precision (quadruple) in the numerics to overcome this problem.

Efficiency and numerical stability
In order to illustrate the efficiency of the approach, I show in Tab. 2 different timings and other related informations for the complete solution at the point m s = .2 and x = .45. Since the evolution is performed separately first in m s and then in x, I require two more digits of local precision in the first step, so that the estimate of the error will be dominated by the second. Note that in the first evolution the x value is fixed from the very beginning as specified in the boundary condition Eq. 14. This value has been hardcoded in the Jacobian, which not only leads to a much faster evaluation, but also to less severe numerical instabilities in the coefficients themselves. The latter feature is actually crucial in reaching the higher precision in the first evolution. Returning to the error estimate, it has to be   Table 2: Timing and efficiency information for the numerical integration of the system of differential equations to the point m s = .2, x = .45. The numbers have been obtained on a 2GHz Intel Core 2 Duo system, after compilation with the Intel Fortran compiler. Quadruple precision is an option of the compiler.
stressed that in numerical integration of systems of differential equations, the precision is specified locally, i.e. one requires that the error of the approximation does not exceed a certain value at every step. Therefore, the global relative error is estimated by the product of the number of steps taken and the requested local error. In practice, the number of steps in the m s evolution never exceeded ten thousand, and therefore the final precision for a local relative error of 10 −20 (in quadruple precision) was roughly sixteen digits. The number of steps needed in the second evolution is usually smaller, because the starting point is far from any singularities. In consequence for the requested local error of 10 −18 , the final global error should not exceed about 10 −15 . There is one additional source of error connected to roundoff and numerical cancellations in the coefficients. Its presence is visible in the table, when comparing the solution in double precision and quadruple precision for the same requested local relative error. In the case of the x evolution of the full system, the number of steps is larger in double precision, precisely because the error estimates are not satisfied due to random roundoff errors. The software tries to reduce the step size until the error estimate satisfies the bound, which eventually happens because the random variations around the true value must, sooner or later, turn close to it. Let me finally comment on the running time. Clearly, if only the leading color coefficient were needed with moderate precision, a value at a single phase space point could be obtained within less than a second. For the full color structure in quadruple precision, as much as fifteen minutes are needed. Although this does not allow for direct implementation in a  Monte-Carlo generator, the functions are smooth enough to be interpolated starting from a grid of values. The efficiency is by far sufficient to obtain dense grids. Therefore, a grid of numerical values for moderate m s together with the series expansion of the previous section for small m s is a complete solution to the problem of evaluation of the two-loop amplitudes for the production of a heavy quark-anti-quark pair in massless quark-antiquark annihilation. In Tab. 3, I give the values of all the color coefficients with ten digits precision at the point m s = .2, x = .45, which is well outside the region of convergence of the series expansion. The actual precision of the result estimated by the variation with respect to the change of the integration contour in x was roughly fourteen digits.
In the course of preparation of the present work, I derived the necessary grids mentioned above. For simplicity the singularities at the kinematic boundaries (forward/backward scattering) have been avoided by keeping a distance of 10 −3 in x, which is unnoticeable on the plots, but can be improved upon in the future. In this respect, the actual needs will only be apparent once the virtual corrections will be combined with real radiation. The range of variation in m s was chosen such that the distance to the threshold parameterized by was at least 10 −3 . This is safely sufficient for any practical applications. The plots of the finite parts of the purely bosonic corrections can be found in Fig. 4, whereas those for the contributions containing a single closed fermionic loop in Fig. 5. The interesting feature is the large variation of the subleading color contributions, when nearing the threshold. Of course, this is due to the Coulomb singularity. In particular, the C 0 coefficient behaves like 1/β 2 , which cannot be compensated by phase space integration. This leads to a true divergence, to be taken care of by resummation in the complete analysis. Finally, let me point two possible improvements of the implementation. The first one is related to the numerical precision. Clearly, different choices of the basis functions may be used to milden the strength of the singularities. For example, for two functions, which suffer from large cancellations, one could introduce a mixture such that one of the functions is small, while the other is large, but contributes with a small coefficient. The second possible improvement is at a lower level and concerns the time of evaluation. Since the computation is done in quadruple precision, one could try different libraries instead of the compiler's built-in routines. In fact, the implementation of [39] has proven to be about three times faster on some problems. Of course, this is of lesser importance than precision, since even times of the order of five minutes per point would still not allow to perform the integration in real time within the framework of a Monte Carlo generator.

Conclusions
In this paper, I have presented a complete solution to the problem of evaluation of the two-loop amplitude for heavy quark production in light quark annihilation. The result, a numeric approximation obtained by a combination of a small mass expansion and integration of a system of differential equations, is satisfactory from the point of view of applications. Aiming at a complete description of the top quark pair production cross section at the next-to-next-to-leading order at the LHC, the next steps to perform will be the evaluation of the two-loop amplitude in gluon fusion and the computation of the real radiation contributions. The mixed real-virtual contributions are in principle known, but have to be supplied with suitable subtraction terms in order to allow for integration over the full phase space.
The result for the series expansion of the amplitude, as well as the grid of values obtained by numerical integration of differential equations (with all numbers rounded at 5 digits) are available in the form of Mathematica files, qqQQ2Lseries.m and qqQQ2Lnumeric.m respectively, attached to the source of the paper on the preprint server http://arXiv.org, but can also be obtained from the author upon request.