Polarized double-virtual amplitudes for heavy-quark pair production

We present the two-loop virtual amplitudes for heavy-quark pair production in light quark-antiquark annihilation and gluon fusion channels, including full spin and color dependence. We use expansions around kinematical limits and numerical integration to obtain results for the involved master integrals. From these, we determine the renormalised infrared finite remainders of the coefficients of amplitude decompositions in terms of color and spin structures. The remainders are given in form of numerical interpolation grids supported by expansions around the production threshold and the high energy limit. Finally, we provide the spin density matrix, which encodes the heavy-quark spin correlations and is sufficient for phenomenological applications. Our results are necessary for the derivation of top-quark pair production cross sections in hadron collisions in the narrow width approximation with next-to-next-to-leading order accuracy in QCD.


Introduction
The top-quark is the heaviest known particle and measurements of its properties provide important insights into the Standard Model of Particle Physics and beyond. Top-quark pair production at hadron colliders like the LHC or Tevatron is an important process for Standard Model precision measurements as well as searches for new physics. Considering the hadronic production of stable top-quark pairs, the prediction from Quantum Chromodynamics (QCD) is compelete to next-to-next-to leading order (NNLO) for the total cross section [1] and for differential distributions [2][3][4]. More recently NLO corrections from electroweak interactions [5,6] were also incorporated. A more complete modelling of pair production including decay and off-shell effects is available to NLO accuracy in QCD [7,8] in the case of the di-lepton channel and more recently also for the semi-leptonic channel [10]. These results were extended to pair production in association with a jet [11,12], which is of relevance for inclusive production at NNLO.
Corrections to the top-quark decay process are known through NNLO in QCD [13,14]. This has allowed for a partial prediction of the pair-production differential cross sections with decay modelled within the Narrow Width Approximation (NWA) [15]. The only missing piece of information is the exact contribution from NNLO production followed by LO decay of the top quarks. This requires the knowledge of polarised two-loop amplitudes for this process, which is the subject of this publication.
The evaluation of the polarised two-loop amplitudes closely follows the lines of [16]. To obtain spin and color dependence of the amplitudes, we use projection techniques, which were also successfully applied in various two-loop calculations, for instance [17]. The most demanding part of this calculation is the reduction and evaluation of involved scalar integrals. The appearing scalar integrals can be reduced to the same set of master integrals as those involved in the evaluation of the spin-summed amplitude. The evaluation of these master integrals uses a variety of analytical and numerical techniques. Exploiting the system of differential equations obeyed by these master integrals is the core idea behind these methods. Most of the physical phase space region can be accessed by solving the differential equations numerically. The regions of phase space that contain physical singularities cannot be reliably accessed using numerical integration. We perform deep power-logarithmic expansions around these singularities in order to obtain precise values for the master integrals. We provide the results in terms of an expansion around the production threshold, a high energy expansion, as well as an interpolation grid. To present and discuss some features of our results, we recast the obtained coefficients with respect to a basis in color and spin space, in terms the spin density matrix of the top quarks alone. Although our results are obtained with numerical methods, there is also progress in the analytic evaluation of master integrals for this process [18][19][20][21][22].
This paper is organised in the following way. In the next section we define the spin and color structures into which the amplitudes are decomposed. We also discuss the projection method we used to obtain the coefficients. Afterwards, we describe the methods used to obtain numerical values for the master integrals in the physical phase space region as well as the improvements we made considering the choice of the master integral basis. Next, we present and discuss the results for the obtained coefficients. We close with conclusions and outlook.
2 Structure of the amplitude

Spin and color structures for virtual amplitudes
The production of heavy quark pairs at hadron-hadron colliders involves two partonic QCD processes at lowest multiplicity gg → QQ and qq → QQ . (2.1) The momenta are assigned as follows with on-shell conditions We define the following kinematic invariants where the relation s − t − u = 0 holds as a consequence of the aforementioned on-shell conditions and momentum conservation. These invariants are related to the scattering angle θ of the top quark (with respect to the beam axis in p 1 direction) and the top quark velocity β as The bare scattering amplitude can be expanded in a perturbative series in α s = g 2 s /4π and reads up to second order ) . (2.6) To facilitate the calculation of polarized virtual amplitudes, we decompose them in terms of color and spin (Lorentz) structures in the color ⊗ spin space of external particles. The color and spin decompositions of virtual amplitudes can be written as where l = number of loops, and the |C g,q i , |S g,q j on the right-hand side represent, respectively, the chosen basis structures in the color ⊗ spin space. We denote the state of the external particles by |a, b, c, d in color and |h 1 , h 2 , h 3 , h 4 in spin space, where a, b and h 1 , h 2 concern the initial state, while c, d and h 3 , h 4 the final state, in the same order as for the kinematics. With this notation, the color and spin basis structures of Eq. (2.7) can be written in full generality in the case of a gluon initial state as and similarly for a quark initial state Notice that in our work, we do not choose any specific representation of spinors or polarisation vectors. Thus, for example, h 3,4 are not necessarily helicities in the case of the heavy quarks. By providing results in terms of spin structures S i , we allow to translate the amplitudes to any particular polarisation basis. For phenomenological applications, we provide the spin density matrix, which contains all the necessary information in terms of spin vectors of the heavy quarks, see Section 2.3. The color decomposition basis of amplitudes can be chosen straightforwardly. For the gg → tt, we use the natural basis In the case of quark annihilation in the initial state, the color basis reads For each of the color structures C i , we decompose the amplitudes further in terms of spin (Lorentz) structures. To this end, we assume that all four external particles are confined to 4-dimensional space and are on-shell with physical polarization states (i.e. 4dimensional equations of motion are satisfied). Under this condition, we have in total 2 4 = 16 different physical helicity configurations both in the gg → tt and the qq → tt process. Additional symmetry properties enjoyed by the amplitudes can lead to relations which further reduce the number of linearly independent structures. Indeed, on top of the aforementioned kinematic constraints, QCD interactions are invariant with respect to parity, under which the helicity of each of the four external particles is flipped, while the color structures are left unchanged. This symmetry then reduces the linearly independent spin structures down to 8 (in 4-dimensional space), both in the gg → tt and qq → tt cases. QCD interactions are also invariant under charge conjugation. However, this operation also involves the color structure. For this reason, we did not impose C-symmetry when determining the basis of Lorentz structures for the color-stripped amplitudes. For the same reason, implications from Bose-symmetry between the two gluons are not considered at this point, but rather used as a test at the end of the calculation.
At this point, we discuss the particulars for the two specific amplitudes, since additional symmetry properties are process dependent.
Let us first consider the gg → tt case. Here, we assume that both polarisation vectors are orthogonal to both of the initial state momenta, p 1 and p 2 (see Eq. (2.2)). This reduces the number of degrees of freedom to the physical two. Spin sums should, therefore, be performed with (2.12) After stripping off the external wave functions we choose the following set of 8 Lorentz structures, with suppressed spinor indices The additional factors of m t and s are inserted such that all structures are dimensionless once multiplied with spinors (which are assumed to have mass dimension 1/2). The structures are grouped according to whether they are symmetric (S 1 to S 6 ) or anti-symmetric (S 7 and S 8 ) under the exchange of µ ↔ ν. It can be checked explicitly that each of the above Lorentz structures is mapped back to itself under the parity transformation (up to a phase factor). The Gram determinant of this set of structures is not identically zero, assuring that they are linearly independent.
In case of the qq → tt process with massless initial-state quarks and limited to QCD interactions, the massless quark line is disconnected from the massive top-quark line. Chirality conservation in QCD, therefore, implies that only half of the helicity configurations of the initial-state massless quarks are non-zero. Once this additional constraint is accounted on top of those aforementioned ones, one finds that there are only four independent helicity amplitudes left. We therefore choose the following set of 4 Lorentz structures of the form S = Γ ⊗ Γ (Γ denotes a string of γ matrices) where the left-hand side of the ⊗ symbol concerns the massless fermion line, while the right-hand side concerns the massive fermion line. Again, the linear independence of these structures can be verified via the Gram determinant. The coefficient functions of the above color and spin decomposition of amplitudes can be extracted by performing the usual projection procedure. In short, the projection of the virtual amplitude onto each of the chosen basis structures gives an equation linear in the coefficient functions. The collection of all projections onto the linearly independent, complete set of basis structures then forms an invertible linear algebraic equation system in the coefficient functions, which can be solved straightforwardly. The coefficient matrix of this linear algebraic equation system is identical to the Gram matrix of the chosen basis.
We note already at this point that in our calculation, the structures S g 6 and S q 4 have vanishing coefficients for all color structures.
Even though we performed our calculations with the basis specified in Eq. (2.10), it is possible to express the coefficient functions for the gluon channel in terms of the orthonormal basis 17) where N C = 3 is the number of colors, and 8 S , 8 A denote the symmetric and anti-symmetric octet states respectively, while 1 the singlet state. The advantage of this choice of basis is that there is no mixing between the color structures when calculating color summed amplitudes. On the other hand, with these structures the coefficient functions exhibit a simple Bose symmetry. Indeed, for the spin structures S g µν 1 , S g µν 2 , S g µν 3 , S g µν 7 , S g µν 8 the coefficients of C g 8 S and C g 1 are symmetric under the exchange cos θ → − cos θ, while the coefficient of C g 8 A is anti-symmetric under this transformation. For the spin structures S g µν 4 and S g µν 5 the situation is reversed, the coefficients of C g 8 S and C g 1 are anti-symmetric while C g 8 A have symmetric coefficients. These properties are consistent with the numerical results, and constitute a test of the calculation.

Ultraviolet and infrared renormalisation
The chosen color ⊗ spin basis may be used in d-dimensions after extension of the spin structures by evanescent combinations. For physical applications, however, we should only need 4-dimensional quantities. Due to the presence of infrared singularities, meaningful amplitudes are only obtained after the usual ultraviolet renormalisation followed by infrared subtraction (multiplicative renormalisation). This procedure results in so-called finite remainders, which are, however, scheme dependent.
The UV renormalized amplitude reads where we used the on-shell wave function renormalisation constants Z g , Z q and Z Q . The renormalised heavy quark mass m is related to the bare mass by m 0 = Z m m. The coupling constant is renormalized in the MS scheme with n f = n l + n h active flavours As argued in [16] a decoupling of the heavy flavours from the running of α s is necessary to correctly accommodate for heavy quark mass effects in regimes where the produced heavy quarks are not very relativistic. This decoupling can be achieved by the replacement where ζ αs is the decoupling constant. The wave-function and the coupling renormalisation (including decoupling) act multiplicatively on the amplitudes and, therefore, also on the coefficients c ij . The mass renormalization counter term, on the other hand, requires an additional decomposition of the lower order amplitudes into color ⊗ spin structures. The necessary renormalisation and decoupling constants are given in the appendix A.
The UV renormalized coefficient functions still contain infrared divergences. However, the infrared structure is known in terms of lower order amplitudes [23][24][25][26][27][28][29], and can be extracted from the UV renormalised amplitude where |F n is the finite remainder amplitude, we are interested in.
It is an operator in color space and can be obtained from its renormalisation group equation where the anomalous dimension Γ is given in the appendix A. Since the Z operator acts in color-space, the terms Z (i) M (k) n have to be projected back onto the color structures to obtain the corresponding counter terms for the coefficients. The minimal definition of the renormalisation operator Z, which consists of poles in the dimensional regularisation parameter only, specifies our IR renormalisation scheme uniquely.
In [61], it was shown that the triple-color correlators of the soft anomalous dimension matrix cannot contribute to spin and color summed matrix elements. Since we keep color and spin dependence this is no longer true in our case. In fact, our calculation is the first to rely on the coefficient to correctly obtain all poles of the coefficients functions. In consequence, it constitutes the first non-trivial cross-check of this contribution to the soft anomalous dimension matrix, which was originally derived in [28].

Spin density matrix
For illustration of our results, we choose to recast the amplitude into a more convenient form. In general, since each coefficient has a real and imaginary part, our calculation yields 54 real functions. However, not all of them enter independently into physical predictions. Therefore, we also evaluate the spin density matrix, which contains all the necessary information on the top-quark spin dependence and is sufficient for phenomenological applications.
The spins of the top-quarks in their rest frame can be described by two normalised spin 3-vectorsŝ t andŝt. They correspond to two four-vectors s t and st in the center of mass frame which have the properties The vectors s t and st enter the matrix element through the insertion of the spin projectors u(p 3 , s t )ū(p 3 , s t ) = / p 3 + m 1 Since we work with finite remainders without any divergences, the presence of the γ 5 matrix does not constitute any complication. Indeed, the spin density matrix is simply evaluated in 4 dimensions. The two-loop contribution to the spin-density matrix for both partonic processes can be decomposed as These functions are related to the 2-loop components of the spin density matrix R q,g as defined in [51] through The coefficients of the occurring structures are functions of cos θ and β only. In pure QCD C,P and CP invariance hold and imply that B t = Bt = B as well as D 1 = D 2 = D [51] for both channels. Therefore we are left with R 2−loop q,g = A q,g + (B) q,g µναβ p 1µ p 2ν p 3α s tβ + µναβ p 1µ p 2ν p 3α st β + (C) q,g (s t · st) + (D) q,g (p 1 · s t )(p 1 · st) + (p 2 · s t )(p 2 · st) In the gluon case we have an additional bose-symmetry which implies that the functions A g , C g , D g are symmetric in cos θ and that B g has to be an antisymmetric function in cos θ. It also implies the relation E 12g (cos θ) = E 21g (− cos θ).

Scalar integrals
The coefficients functions c ij are given by linear combinations of a large number of scalar integrals with rational coefficients in s, t, m 2 and . These scalar integrals are expressed through linear combinations of master integrals using an Integration-by-Parts (IBP) reduction. We can rewrite the coefficients, as well as the master integrals, in terms of dimensionless variables m s = m 2 /s and x = t/s. From the IBP relations we can obtain a system of differential equations for the master integrals x where A (ms) and A (x) are matrices whose elements are rational functions in m s , x and . We do not choose the same set of master integrals as in the spin summed calculation. There is an enhancement of the matrix elements at high energies and small/large scattering angles resulting from diagrams of a t/u-channel type as indicated in Fig. 1. These enhancements require numerically very stable results for the master integrals in these phase space regions. For this reason, we decided to try a basis for a subset of the integrals, which corresponds to the -form of the differential equation (see next section). Our hope was that in this basis the numerical evaluation will become more stable. Ultimately, this turned out not to be the case. We stress, nevertheless, that the results obtained in the old and new bases for the spin summed amplitudes agree to several digits, within the accuracy of the calculation.

Canonicalization
With the hope to achieve a better stability when numerically solving the differential equations for master integrals involved in the two-loop gg → tt process, we choose to put the equation system partially into the -form [30], where the right-hand side of the differential equation system is proportional to = d−4 2 and the singularities are only simple poles in the kinematic variables. Algorithmic approaches have been devised to arrive at the -form for a given differential equation system for master-integrals in a single variable [31,32]. They have been implemented in Fushsia [34] and Epsilon [35] which are publicly available. In the case of multiple variables there also exists an algorithm, presented in Refs. [36,37] and implemented in a program called Canonica. It is well known that, for a given set of master integrals, an -form is not always achievable by a rational transformation of the integral basis. It is also not always possible [32,38,39] even with more general transformations. In particular, Ref. [32] provides a strict criterion for the existence of an -form in the case of master integrals of a single variable. In particular, the 4-dimensional homogeneous part of the differential equation system typically corresponds to high-order Picard-Fuchs differential equations that do not factorize completely [39]. The simplest counter example is given by the differential equations of the master integrals of the two-loop sunset diagrams with identical masses [40,41], where solutions involve elliptic integrals.
The topology of the two-loop sunset diagrams with equal masses appears in the IBPreduction of the master integrals of the two-loop gg → tt diagrams. Thus, it is not a surprise that the full system of differential equations of the 422 master integrals involved cannot be completely put into the -form. In addition, a considerable amount of sectors require individual coordinate transformations in order to arrive at their respective -forms by rational transformations (in the new variables). Since we would like to numerically integrate the full differential equation system of all master integrals in one go, we are then forced to divide the master integrals into two subsets: those that can be directly put into the -form via a rational transformation in the original variables and those that cannot. The second subset essentially consists of master integrals fulfilling any of the following three conditions: 1) their expressions involve elliptic integrals; 2) coordinate transformations are required in order to reach their -form; 3) their derivatives involve any one of the aforementioned two kinds of master integrals.
Under such tight selection criteria, there are only 65 master integrals that can be directly transformed into the basis observing the -form (in the original variables). They are identified and then subsequently moved to the front of the differential equation system of the 422 master integrals, without spoiling the block-wise triangular structure of the differential equation system. The numerical evaluation of the complicated master integrals are expected to benefit from the -form of these 65 master integrals. The differential equation system of these 65 master integrals in question involves more than one variable 1 , and we employ the package CANONICA [37] to find the rational transformation matrix needed for obtaining the -form. A few modifications of the program were made in order to tackle this 65-by-65 system with less time consumption.
As a side remark, we would like to briefly mention the following point. Due to the existence of remnant rational transformations that preserve the -form of a differential equation system, the new basis integrals defined by the rational transformation matrix returned by the package CANONICA [37] are not guaranteed to be of uniform weight [30,38]. Upon a closer examination, we find that in general not all of these remnant rational transformations respecting the -form (if it exists) are of weight 0 according to the counting rules laid in [30,38].
To be more specific about this, we find this remnant freedom to be the following. Under any rational transformation mixed with any coordinate-transformation under the condition of keeping the resulting differential system still rational in the new variables, the remnant rational transformations that preserve the -form of the differential equation system (assuming it exists and is represented by dÃ), read whereĈ can be any invertible constant matrix of rational numbers of the same dimension as the -form coefficient-matrix dÃ.T I ( ) is independent of kinematics but possibly possesses some non-trivial -dependence. It can be any element from the invariance symmetry group of the coefficient-matrix dÃ with matrix elements being Laurent-polynomials in with rational numerical coefficients. In other words,T I ( ) is a matrix that is invertible and commutes with dÃ with all its matrix-elements living in the ring of over rational numbers. Owing to its invertibility and the fact that it consists of rational numbers only, the matrix C always preserves the uniform weight feature of the vector on which it acts (if this feature is there in the first place). However, some ofT I ( ) with non-trivial -dependence can turn a list of uniform weight master-integrals into a list of non-uniform-weight integrals, and vice versa, even though both perfectly observe -form differential equations. Since there is no reference to the concrete boundary conditions of the differential equations in the process of finding the rational transformation done by the package CANONICA, it is therefore not guaranteed that the solutions of the -form differential equations thus obtained are of uniform weight.

Master integral evaluation
We subdivide the physical phase space region into three regions, the high energy limit where m s → 0, the threshold region β → 0, and the "bulk" which describes the rest of the phase space region.
For the numerical integration of the new set of master integrals, high precision boundaries are needed. They are obtained from the power-logarithmic expansion in the high energy limit from the original set of master integrals. The first few terms of those expansions were obtained with Mellin-Barns techniques using the MB package [46]. These expansions are exact in t, where in some cases the differential equations were used to get the exact behaviour from the limit t = 0. In this double limit m 2 → 0 and t → 0 the integrals were evaluated numerically with very high precision and then resummed with PSLQ algorithm [47] or XSummer [48]. These first terms were used to derive deep expansions in m s by using the available differential equations. These deep expansions were subsequently used to compute high precision boundaries for the numerical integration.
Starting from the numerical results obtained from the deep power logarithmic expansion, we perform a numerical integration along contours in the complex plane. In our programs we incorporate software from [49] for solving the differential equations and [50] to handle higher precision numbers. The endpoints of the contours define an interpolation grid. In the region which is accessible with this method, i.e. where the coefficient functions are not too singular, we evaluate the amplitudes using this interpolation grid. The sampling points are the same as in [1,[52][53][54] with an extension to higher values of β.
In the limit β → 0 some master integrals show singular behavior and are difficult to obtain with the method of numerical integration. We perform a deep power-logarithmic expansion of the master integrals in β by again exploiting the differential equations up to O β 50 and O ln 10 β . This expansion is done for several fixed angles cos θ with unknown boundary conditions, which are finally determined by matching to the results obtained by numerical integration. In this publication, we provide results for the finite remainders of all coefficient functions. They are given in the form of an interpolation grid as well as kinematic expansions in the high energy limit and near the production threshold. The decomposition into the At tree-level we find the following non-vanishing coefficients in case of gluons to be non-vanishing. We find that in both cases, quark and gluon initial state, one spin structure has vanishing coefficients for every color structure at one and two-loops. Indeed, in the case of gluons the coefficients c i6 vanish, as do the coefficients c i4 in the case of quarks. All other spin structures have non-vanishing coefficients.
The high energy limit of the coefficients was calculated as an analytic power-logarithmic expansion in m s = mt s up to O m 4 s , using the boundary expressions for the master integrals. This expansion assumes that t, u m 2 t and is therefore not valid in the region of high-energy forward/backward scattering. The results were cross-checked against the spin summed amplitude.
We want to mention that the depth of the expansion does not translate easily to the expansion depth of the square summed or spin correlated matrix element since there is a non trivial dependence on m s (or β in case of the threshold expansion) hidden in the spin structures themselves.
The "bulk" region is parameterized on a grid which is specified by equally spaced points in β and two additional points close to the high energy boundary. The point β 80 = 0.999 is sufficient for LHC with 8 TeV center-of-mass energy, which was the extent of the interpolation grid of the spin summed calculation. Here, we extend the grid to β 81 = 0.9997 which corresponds to a center of mass energy of 14 TeV, for contemporary applications. For cos θ we choose 42 points obtained from where we chose the x i as the 21 points obtained from the Gauss-Kronrod integration rule of degree 10. Values for β < 0.1 were obtained from the threshold expansion of the master integrals. The dependence on the number of light fermions is kept. In order to illustrate our results we plot the coefficients of the spin density matrix. We introduce the following normalization factors, which were also used for the presentation of the results in [16] Figure 4. The difference between the threshold expansion for coefficient A g up to β n with n = 0, 2, 4, 6 and results from numerical integration for a fixed angle θ. and define is the spinsummed and averaged two-loop finite remainder and was checked against the result from [16]. The function B n l =5 q,g describes the transverse polarization of the top quarks resulting from absorptive parts of the amplitude. At tree-level, this coefficient vanishes due to the absence of complex couplings in QCD. At higher orders, the non-vanishing imaginary part of the virtual amplitudes yields non-zero coefficients. The remaining functions encode the spin correlations between the top and anti-top quark since their structures contain both spin vectors. In the gluon channel, all expected symmetry properties under cos θ → − cos θ of the coefficient functions are clearly fulfilled.
The threshold region is covered by points obtained from the deep power-logarithmic expansions of the master integrals. In addition we perform a power-log expansion in β for all coefficients up to β 2 . This is done for different but fixed scattering angles cos θ c ij (β, cos θ n ) = 2 k=−2 2 l=0c ij,kl,n β k ln l β The dependence on θ was recovered by performing a fit for each set {c ij,kl,n } n to a polynomialc ij,kl = 2+k n=0 a n cos n θ separately for the real and imaginary part. The results are also available in electronic format together with this paper. This expansion was used to (d) Figure 5. The difference between the threshold expansion for coefficient B g , C g , D g , E 12g up to β n with n = 0, 2, 4, 6 and results from numerical integration for a fixed angle θ.
determine the corresponding coefficients of the spin density matrix as well. Up to O β 0 we reproduce the analytic result obtained for the spin-summed case in [16].
To study the quality and convergence of the expansion, we also calculated the density matrix for a fixed angle (chosen to be the point x 9 ) up to order O β 6 . In Figs. 4 and 5 we compare this expansion against the results obtained from the interpolation grid. We show the difference 2 X n l =5 diff (β, x 9 ) = X n l =5 thres (β, x 9 ) − X n l =5 grid (β, x 9 ) , (4.10) with X ∈ {A g , B g , C g , D g , E 12g } for different expansion depths of X n l =5 thres (β, x 9 ). The series seems to converge nicely and, if expanded up to O β 6 , provides a reasonable description of the amplitude in the region β < 0.3.

Conclusions and Outlook
We presented the decomposition of the heavy-quark pair production amplitude in terms of spin and color structures at two-loop level. We provide results in terms of interpolation grids and kinematic expansions for these coefficients. As a first application we calculated the spin-density matrix for top-quark pairs. With this work we provide the missing piece needed for the calculation of on-shell top-quark pair production and decay at NNLO in QCD including spin-correlation effects in the narrow width approximation. We improved the numerical results obtained for the involved master integrals by changing to a partly canonical basis. The incorporation of these amplitudes in a full-fledged calculation of top-quark pair production and decay is work in progress.
The full set of results of this paper is available at: https://git.rwth-aachen.de/mczakon/PolarizedTTNNLO.