Complete two-loop QCD amplitudes for tW production at hadron colliders

We calculate the complete two-loop QCD amplitudes for hadronic tW production by combining analytical and numerical techniques. The amplitudes have been ﬁrst reduced to master integrals of eight planar and seven non-planar families, which can contain at most four massive propagators. Then a rational transformation of the master integrals is found to obtain a good basis so that the dimensional parameter decouples from the kinematic variables in the denomi-nators of reduction coe ﬃ cients. The master integrals are computed by solving their di ﬀ erential equations numerically. We ﬁnd that the ﬁnite part of the two-loop squared amplitude is stable in the bulk of the phase space. After phase space integration and convolution with the parton distributions, it increases the LO cross section by about 3%.


Introduction
The top quark, first discovered at the Fermilab Tevatron [1,2], is the heaviest elementary particle in the Standard Model (SM), and may have a close relation to the electroweak symmetry breaking due to its large coupling with the Higgs boson.The single top-quark production processes deserve detailed studies because they can be used to measure the Wtb coupling structure, which may be modified by new physics.Moreover, they are often considered as important backgrounds to many new physics searches.In this paper, we focus on the tW associated production, which has the second largest rate, after the t-channel, at the large hadron collider (LHC).The inclusive and differential cross sections have been measured to an accuracy of 10% [3][4][5][6].
In order to provide precise theoretical predictions for the rate and kinematical distributions, higher-order quantum corrections are indispensable.The next-to-leading order (NLO) quantum chromodynamics (QCD) correction has been obtained for stable tW production [7][8][9][10].The correction to the process including decays was calculated in [11].The NLO QCD result was interfaced with parton shower within both the MC@NLO and POWHEG formalisms [12][13][14].The NLO electroweak correction was computed in [15].
There have been also some efforts devoted to the calculation of the corrections beyond the NLO in QCD.The approximate next-to-next-to-next-to-leading order total cross section was obtained by expanding the threshold resummation formula [16][17][18][19].The corrections induced by soft gluons have been resummed to all orders in the strong coupling [20].The N-jettiness soft function, which is one of the ingredients in a full NNLO QCD correction, has been calculated by two of the authors [21,22].Recently, we reported the analytical results of the leading color and light quark-loop contributions to the NNLO virtual corrections [23] using the two-loop master integrals obtained in [24][25][26].We have investigated the dominance of the leading color result in the one-loop squared amplitudes, of which the full color result was computed analytically in [27].However, it is ideal to have the full two-loop virtual correction, which is the aim of this paper.This is made possible due to the impressive progresses in both the analytical and numerical methods of the computation of Feynman diagrams; see, e.g., [28,29].
This paper is organised as follows.In section 2, we describe the basic setup and details in our calculation of complete two-loop bare amplitudes.We also discuss the methods to cancel the ultra-violet (UV) and infra-red (IR) divergences in the bare amplitudes.The finite part of the squared amplitude defines the hard function which is needed in a NNLO calculation.The numerical results for the NNLO hard function are presented in section 3. Finally we make the conclusion in section 4.

Calculation method
We calculate the two-loop corrections to the process g(k 1 ) The Lorentz invariant Mandelstam variables are defined by where α s is the strong coupling.In this work, we do not keep the polarization information and focus only on amplitude squared.We have presented the analytic result of the one-loop squared amplitude in [27].In this paper, we calculate the interference between the two-loop and tree-level amplitudes.According to the color structure and fermion-loop contribution, it is decomposed as where N C is the number of quark colors, n l (n h ) is the total number of massless (massive) quark flavors.The coefficients of various color structures are denoted by A, B, C, D, E, F, G.In our previous paper [23], we have obtained the analytical results of the leading color A and light fermion-loop contributions E l , F l and G l .In this work, we continue to calculate the remaining contributions B, C, D, E h , F h and G h .We generate the two-loop Feynman diagrams using FeynArts [30] and perform the Dirac algebra with the help of FeynCalc [31].In our calculation, the anticommuting γ 5 scheme proposed in [32] is implemented.The traces with two γ 5 's are easy, while the traces with one γ 5 vanish because there are only three independent momenta in the tW production process.More discussion on the implementation can be found in [27].After the spin and polarization summation, all Lorentz indices are contracted, and the squared amplitude is written as a combination of a large number of scalar integrals.They are reduced to a set of basis integrals, called master integrals, making use of the integration-by-part (IBP) identities, which are automatically generated in FIRE [33].It turns out that all the scalar integrals could be expressed by 15 families of master integrals1 , which are defined according to the topologies of the Feynman integrals.Specifically, the master integrals are categorized in terms of eight planar and seven non-planar topologies, which are shown in figures 1 and 2, respectively.They have up to four massive propagators, which makes the analytical structure complicated.Name Definition Non-Planar NP1 q 2 1 , (q 2 − q 1 ) 2 , q 2 2 , (q 1 + k 1 ) 2 , (q 2 − q 1 + k 2 ) 2 , Table 1: Definitions of the master integral families in terms of the denominators D i .The above q 1 and q 2 are loop momenta while k 1 , k 2 , k 3 and k 4 are external momenta defined in Eq. ( 1).
The master integrals corresponding to the above topologies are defined as follows: where d = 4 − 2 is the spacetime dimension, γ E ≈ 0.5772 is the Euler-Mascheroni constant, and the denominators D i in each topology are listed in table 1.
The master integrals in the P1, NP1 and P2 topologies have been calculated analytically in [24,25] using the method of differential equations [34,35].There are 31,34 and 38 master integrals in the P1, NP1 and P2 topologies, respectively.And the differential equations have been transformed to the canonical form [28] after constructing a proper basis.The solutions can be expressed as multiple polylogarithms [36].The master integrals in the NP2, NP3 and NP4 topologies have also been computed after expansion in m 2 W [26].In the other integral families, however, multiple square roots of three variables appear in the differential equations.It is challenging, if not impossible, to find a transformation to rationalize all the square roots simultaneously.Another obstruction is the fact that elliptic integrals are needed in many sectors of the integral families.Despite impressive progress in the past decade [37][38][39][40][41][42][43][44][45][46][47][48][49], the approach to analytically compute the Feynman integrals depending on several elliptic curves is not as full-fledged as that for Feynman integrals which can be evaluated to multiple polylogarithms.Therefore, we choose to perform the calculation of these integrals numerically using a method that is efficient enough for phenomenology analysis.
We first identify the master integrals in each integral family, to which the other scalar integrals can be reduced.It is not required in this process that a canonical basis is chosen, though we have succeeded in obtaining such a basis in some sectors.We carry out only rational transformation of the basis and demand that the reduction coefficients have a "good" factorization property, i.e., the denominators of the coefficients can be written as N(d + M) with N independent of d and M a rational number [50,51].This is always possible because the integrals should not have poles or branch cuts depending on the spacetime dimension d.The choice of such a kind of "good" basis integrals turns out to be very important to minish the size of the coefficients after IBP reduction and to solve the differential equations efficiently.In some cases, choosing an appropriate basis avoids the cancellation of large numbers and thus makes the numerical calculation more efficient.We show the master integrals for the most complex integral families NP5, NP6 and NP7 in the appendix.
These master integrals are computed numerically.Specifically, we construct the differential equations of the master integrals with respect to s and t.And then, we calculate them at one kinematic point (s 0 , t 0 ), which is chosen in physical region, using the AMFlow [29,52] package.The results at the other phase space points can be obtained by solving the differential equations in terms of a combination of multiple series expansions.We have made use of the DESolver function implemented in the AMFlow package.The input includes the set of differential equations as well as an integration path.There may be pseudo poles on the path, which appear as poles in the differential equations but do not exist in the master integrals.In order to avoid numerical instability, we choose a contour around the pseudo poles.A typical example is shown in figure 3. The integration is performed from t n to t n+1 via a point (t n +t n+1 ) 2 + (t n+1 − t n )i in the complex plane of t, given that t s ∈ (t n , t n+1 ) is a pseudo pole on the real axis.Notice that it does not matter whether the contour goes above or below the pseudo pole.
Figure 3: A contour that was set in the numerical evaluation of the master integrals when using the AMFlow package.
We have checked our calculation of the master integrals in different methods.Firstly, we compare the results of master integrals at a phase space point (s i , t i ) obtained in two ways.One is from the direct evaluation of the AMFlow package, while the other is derived from solving the differential equation as described above.The numerical agreement between these two methods could reach more than 50 digits for the integrals of transcendental weight four 2 .Secondly, the master integrals have also been computed at some phase space points using the FIESTA [53] package that is developed with the sector decomposition method.We have found an agreement within the numerical uncertainties of FIESTA.Thirdly, the numerical results are compared against the analytic ones that are available in some families, such as the P1, P2 and NP1 topologies [24].
The two-loop amplitudes mentioned above contain both UV and IR divergences.The UV divergences cancel against the contribution from renormalization of the couplings, masses and field strength.As a result, the renormalized amplitudes can be written as where Z g,b,t represent the renormalization constants for the external states, and Z α s and Z m are renormalization factors for the strong coupling α s and the top-quark mass, respectively.Their explicit expressions can be found in [54][55][56][57].
The left IR divergences in M ren should be combined with those from real corrections to give finite predictions for the cross sections.Though we have not computed the latter explicitly in this work, the IR divergences can be reconstructed from some universal anomalous dimensions due to the understanding of the factorization structure of the amplitudes in the IR limit.We obtain a finite remainder after subtracting the IR divergences where the factor Z is known up to two-loop level for general processes in QCD [58][59][60][61][62] and to three-loop level for single top production [63].The explicit expression of this factor for tW production can be found in our previous papers [23,27].
We have made nontrivial validations at the amplitude level.First, we reproduce the leading color result which has been obtained in analytical form in [23].Second, all the IR divergences of the full color result are indeed canceled out in Eq.( 6).
We define the squared M fin as the hard function that would be used in a NNLO computation, where we have made a perturbative expansion in the second equation.Similar to Eq. ( 3), the second order of the hard function can be written as We have presented the results of H A and the terms proportional to n l in our previous paper [23].In this work, we provide the complete results of the hard function.

Numerical results
We parameterize the phase space of gb → tW process using the dimensionless variables β t and cos θ in the centerof-mass frame of incoming partons, which are defined by Then the Mandelstam variables read with ).The full phase space spans over 0 ≤ β t < 1 and −1 ≤ cos θ ≤ 1 .In practice, we have generated a grid of 80 × 42 phase space points and computed the amplitude squared on this grid 3 .For simplicity, we extracted a factor e 2 g 2 s / sin 2 θ W in the numerical results of the hard functions shown in the figures.In numerical calculation, the W boson mass is set by a rational identity m 2 W /m 2 t = 3/14, which can significantly speed up the IBP reduction procedure.The renormalization scale is chosen to µ = m t .To illustrate the size of different contributions to the hard function, we show in figure 4 the results of the two-loop corrections without quark loops, the one-loop squared, the light and heavy quark loop corrections separately in 3-D diagrams as a function of β t and cos θ.One can see that these corrections are stable in a large region of β t and cos θ.The two-loop without quark loops and the one-loop squared corrections are positive while the quark-loop corrections are negative.However, they vary dramatically toward the boundary of β t = 1 or cos θ = 1.Specifically, the heavy quark-loop contributions become negative rapidly as β t → 1 and cos θ → 1, while the one-loop squared corrections are positively divergent as β t → 1.The two-loop corrections without quark loops first decrease (increase) and then increase (decrease) as β t → 1 for −1 ≤ cos θ 0.5 (0.5 cos θ ≤ 1).The light quark-loop contributions drop quickly when cos θ approaches 1 and β t is larger than about 0.6 but less than 1.
These divergence behaviors can be understood from the structure of the squared amplitude.The LO squared amplitude contains 1/(1 − β t cos θ).This fact explains the divergence near the corner of β t = 1 and cos θ = 1.The higher-order corrections develop additional collinear divergences in the form of ln(1 ± β t cos θ) and ln(1 − β 2 t ).The latter is responsible for the divergence away from the region with cos θ = 1.The explicit form of these logarithms can be predicted and resummed to all orders of the strong coupling following the method in Ref. [64].We leave this to future work.For practical phenomenology study at the LHC, the suppression of the parton distribution functions near the region of β t → 1 dominates over these large logarithms [27] and therefore it is safe to integrate over the full phase space in the calculation of the hadronic cross section.
In our previous work [23], we have provided the result of the leading color contribution.Based on the smallness of the 1/N 2 C expansion parameter, we expected that it already gives the main correction.To estimate this approximation, we show the contributions of different color structures in the NNLO hard function without quark loops (n f = n l + n h = 0) in figure 5.This is a 2-D diagram so that one can recognize the magnitude more clearly.We have drawn the curves corresponding to the parameter cos θ = ±1 and the results at other values of cos θ can be considered between them.Apart from the region around β t > 0.7, the leading color contribution is larger than the subleading color ones almost by an order of magnitude.For example where the numbers on the right hand side are aligned according to the power of N 2 C .From figure 5, the divergence behaviour of the hard function near β t = 1 is also observed unambiguously.
Lastly, but most importantly, we show the numerical result of the full NNLO hard function in figure 6.It is interesting to find that this function is flat over the large region of β t < 0.8 as a consequence of strong cancellation among different contributions.When β t becomes larger than 0.8, the NNLO hard function first decreases and then increases dramatically.It drops down only at the corner of both β t → 1 and cos θ → 1.
Notice that the hard function is not renormalization scale independent.The scale-dependent terms can be recovered by the equation where the right hand side is easy to obtain [23,27].
To have a rough estimate of the correction to the hadronic cross section at the 13 TeV LHC, we integrate over the full phase space and perform convolution with the CT14nlo parton distribution function [65] via the interface of LHAPDF [66].The factorization and renormalization scales are set at m t .We find that the full NNLO hard function provides a correction of about 3% to the LO cross section.

Conclusion
In this paper, we present the calculation of complete two-loop amplitudes for hadronic tW production.The master integrals contain up to four massive propagators, and the corresponding differential equations involve multiple square roots that can not be rationalized simultaneously.Moreover, some of the master integrals rely on several elliptic curves, which poses a challenge to analytical calculation.We first choose such a good integral basis using only rational transformation that the dimensional parameter is decoupled from kinematic variables in the denominators of the coefficients in the differential equations as well as the reduction coefficients in the amplitudes.And then we calculate the boundary conditions and solve the differential equations we derived numerically with the AMFlow package.We have implemented a contour in the complex plane of the integration variable to maintain numerical stability in the case of pseudo poles.We have made nontrivial checks on our calculation for both the master integrals and the full amplitudes.In particular, the total divergences arising from many Feynman diagrams agree with the ones predicted from universal anomalous dimensions.The finite remainder contributes to the hard function that would be used in a NNLO computation.The NNLO hard function is stable when the top-quark velocity β t is less than 0.8.As β t increases, the hard function changes dramatically due to the large logarithmic enhancement from ln(1 − β t ).After phase space integration and convolution with the parton distribution function, the NNLO hard function increases the LO cross section by about 3%.

Figure 4 :
Figure 4: Different contributions to the NNLO hard function as a function of β t and cos θ.Note that the 3-D diagrams have been rotated to an appropriate angle for visualization convenience.The n l (n h ) contribution denotes the light (heavy) quark loop correction.

Figure 5 :
Figure 5: Contributions of different color structures in the NNLO hard function without quark loops normalized by the LO result.The red, blue, cyan, and dark blue lines correspond to H A N 4 C , H B N 2 C , 50 × H C and 1000 × H D /N 2 C , respectively.The solid and dashed lines represent the results of cos θ = −1 and cos θ = 1, respectively.The lower panel shows the ratio of the leading color contribution to the hard function without quark loops.

Figure 6 :
Figure 6: NNLO hard function normalized by the LO result as a function of β t and cos θ.
Non-planar integral topologies for gb → tW.The black and red thick lines represent the top quark and the W boson, respectively.The black thin lines denote massless particles.For the NP2 and NP4 topologies, the symmetric diagrams (k 1 ↔ k 2 ) are needed in the amplitude.