One-loop squared amplitudes for hadronic $tW$ production at next-to-next-to-leading order in QCD

We present the analytic results of one-loop squared amplitudes for $tW$ production at a hadron collider. The calculation is performed using the method of differential equations. After renormalization, we have checked that the infrared divergences agree with the general structure predicted by anomalous dimensions. The finite remainder contributes to the next-to-next-to-leading order hard function, one of the essential ingredients in the factorization formula of the cross section near the infrared region, which can be used in resummation of all-order soft gluon effects or a differential next-to-next-to-leading order calculation based on the phase space slicing method.


Introduction
The top quark is the heaviest elementary particle in the Standard Model (SM) and may play an essential role in the electroweak symmetry breaking due to its large coupling with the Higgs boson. Its short lifetime also allows the study on its various properties as a "bare" quark without hadronization. Meanwhile the W boson couples with many SM particles via gauge interactions and therefore could be sensitive to the new physics beyond the SM with such interactions modified. The recent precise measurement of the CDF collaboration on the W boson mass [1] has shed light on a significant tension with the other experimental observations in the assumption of the SM electroweak framework. Especially the strong correlation between the top quark and the W boson implies that the precise theoretical investigations are indispensable.
The top quark production associated with a W boson, i.e. tW -channel, as shown in the leading order (LO) Feynman diagrams in figure 1, can be used to probe the coupling between the W boson and the top quark, including the Cabibbo-Kobayashi-Maskawa matrix element V tb . The ATLAS and CMS collaborations at the CERN large hadron collider (LHC) have already observed the signal of this process based on the data accumulated at center-of-mass energies of 7, 8, 13 TeV [2][3][4][5][6][7][8][9][10]. So far, only the data recorded in 2016, corresponding to 35.9 fb −1 integrated luminosity, have been taken into account in the public analyses, and the uncertainty of the measured cross section is about 11% [11]. It is expected that the measured cross section of this process will become more precise as more data are going to be analysed in the future. To match the accuracy of the experimental measurements, the theoretical predictions must include higher-order radiative effects in the perturbative theories. Since the strong coupling is the largest among the three interactions relevant for fundamental particles, the QCD corrections are the dominant contribution. The next-to-leading order (NLO) corrections have been obtained for stable tW final state [12][13][14][15] as well as the process including the decay of the top quark and the W boson [16]. The approximate next-to-next-to-next-toleading order total cross section has been derived by expanding the threshold resummation formula [17][18][19][20]. The all-order corrections induced by soft gluons are considered in [21]. The effect of parton shower has been studied in [22][23][24].
However, the complete next-to-next-to-leading order (NNLO) QCD corrections have not been obtained yet, though a small fraction of the two-loop master integrals have recently been computed analytically [25,26]. To obtain these corrections, the double-real, real-virtual and double-virtual corrections should be calculated independently and collected together in the end to provide an infrared (IR) finite prediction. In arranging the cancellation of IR divergences, the phase space slicing method [27] is one of the promising approaches. Above the slicing cut, which can be the transverse momentum of the tW system or the N-jettiness variable [28] in this process, only the NLO divergence appears and it is appropriate to use some automatic packages to perform the calculation. In this part, it is mandatory to consider the interference with the top quark pair production with a top quark decayt → Wb in order to provide a well-defined perturbative prediction. Different schemes have been proposed to tackle this issue [29]. Below the slicing cut, a factorization of the cross section to several simpler functions holds at the leading power in the ratio of the cut over the center-of-mass energy. Expanding the factorization formula to certain orders in the strong coupling gives the required fixed-order contributions. As one of the ingredients in the factorization formula, the N-jettiness soft function has been calculated at NNLO by two of the authors [30,31]. The missing part is the NNLO hard function, which demands the calculation of one-loop squared amplitudes and the interference between two-loop and tree-level amplitudes. In this paper, we present the result of the former. This paper is organized as follows. In section 2, we show the kinematics of tW production and the notations used in the following sections. Then we describe the details in calculating the one-loop squared amplitudes in section 3, including the γ 5 scheme, renormalization procedure, and the infrared divergences cancellation. In section 4 we discuss the approximation method by expansion in m 2 W . Conclusions are made in section 5.

Kinematics and notations
The process g(k 1 )b(k 2 ) → W (k 3 )t(k 4 ) contains two massive final states with different masses. For the external particles, there are on-shell conditions k 2 The Mandelstam variables are defined as where t a 4,2 is the generator of the color SU(3) gauge group inserted between the top quark and the bottom quark. µ 1 and * µ 3 are the polarization vectors for the gluon and W boson, respectively. P L is the projection operator to the left-handed quarks.
The polarization summation over all physical polarization states for the W boson and the gluon is respectively. Here n is a light-like four-vector. Consequently, the ghost field is not needed in the external state in the loop amplitudes. In practice, we can simply use due to the Ward identity 1 . As a result, we obtain the LO amplitude squared where the factors from spin and color average in the initial states have been included. The LO hadronic cross section for tW production at the LHC is where dΦ 2 is the Lorentz invariant two-body phase space, and f i/p (x, µ f ) is the parton distribution function at the factorization scale µ f . In order to obtain the complete NNLO QCD corrections to this process, we need to obtain one-loop virtual corrections to M gb→W t up to O( 2 ), two-loop virtual corrections to M gb→W t up to O( 0 ), one-loop virtual corrections to M gb→W t+j up to O( 0 ) 2 , and tree-level M gb→W t+jj . In this paper, we focus on the first contribution, leaving the others and the combination of all contributions after phase space integration to the future work.

Bare loop amplitudes
We generate all the one-loop virtual diagrams as well as the diagrams with counterterms by FeynArts [32]. The corresponding amplitudes can be written down automatically and manipulated by FeynCalc [33]. We do not keep the polarization information at the moment and focus only on amplitude squared, i.e., we calculate the interference between the one-loop diagrams and tree-level or one-loop diagrams. All the Lorentz indices are contracted after spin sums, and therefore the loop integrals appear in the form of scalar integrals, though there may be some loop momenta in the numerators. The spin and polarization information may be useful if one wants to study the kinematic distribution of the decay products of the top quark and the W boson. This can be obtained by using the spinor-helicity formalism for the external states after decay, as in [16] for the NLO corrections, or by replacing the polarization vector of specific helicities by a linear combination of external momenta [34]. Our calculation can be readily extended to the helicity amplitudes because they share the same master integrals as the spin-summed amplitudes.
We have chosen the anticommuting γ 5 scheme following ref. [35] to calculate the traces in the squared amplitude. In this scheme, it is easy to deal with the traces with two (or more generally even number of) γ 5 matrices in our case due to the anticommutativity property. In the trace with a single (or more generally odd number of) γ 5 , since the cyclicity in traces has been sacrificed, we choose the W tb vertex in the amplitude as a reading point, i.e. the starting vertex in a trace. The non-vanishing trace with one γ 5 would contribute an anti-symmetric ε µνρσ tensor to each term in the result. Its Lorentz indices should contract only with those of the external momenta in the spin summed amplitude squared, given that there exists only one fermion line. Since only three independent momenta are involved in tW production process, the traces with one γ 5 would vanish at the end.
Then we use two different methods to reduce all the scalar integrals in the squared amplitude to a set of basis integrals, called master integrals. In the first method, the reduction is performed following the Passarino-Veltman procedure [36], and the master integrals can be found from the literature [37] or calculated by Package-X [38] up to O( 0 ). They can also be evaluated up to higher orders in using the AMFlow package [39]. In the other method, we use integration-by-parts (IBP) technique to derive relations among the scalar integrals and carry out the reduction using FIRE [40]. The master integrals are calculated using the method of differential equations [41,42]. After transforming the differential equations to a canonical form [43], the solutions can be readily obtained up to any orders in . It has been checked analytically that both methods lead to the same results for the interference of one-loop and tree-level amplitudes. For the one-loop squared amplitudes, the coefficients of −4 and −3 coincide in analytical form, while the coefficients of −2 , −1 and 0 agree numerically. We also notice that there are more master integrals using Passarino-Veltman reduction than IBP relations. The coefficients of the master integrals in the former contain no divergences, while those in the latter do. These differences indicate that the agreement between the two methods is highly non-trivial. In the following, we present more details about the IBP method.
We find that all the master integrals in the one-loop diagrams could be expressed in terms of three families of master integrals as: and We will present the details for the calculation of the first family in the text, but leave the second family to the appendix A. The third family can be related to the first using k 1 ↔ k 2 , and thus we will not present the results for the master integrals in it. For the squared one-loop amplitude, the reduction of any general integral to master integrals can be carried out as the traditional two-loop integrals. It is interesting to notice that additional Lorentz invariant combinations, such as (q 1 · q 2 ) 3 with q 2 being another loop momentum, can appear in the squared amplitude after summing the polarizations of the vector bosons. But they do not appear in the master integrals. As a consequence, we need the product of one integral in any of the above families and another in complex conjugate form.
The canonical basis of the first family is chosen to be To rationalize the roots, we define the variables x, y, z via The master integrals in eq.(3.3) satisfy the differential equation The arguments l i of this d ln form encode all the dependence of the differential equations on the kinematics, are called the alphabet, consisting of the following letters, (3.7) R i are rational matrices, given by In order to solve the above differential equations, we need to specify the boundary conditions. The integrals M 1 and M 4 are simple enough so that one can get the analytical results with full dependence on from direct computation. Since all the one-loop integrals we encounter do not have any branch cut at the kinematic points corresponding to t = 0, m W = 0, s = (m t + m W ) 2 , and ut = m 2 t m 2 W , they are regular without any singularity at these kinematic points. The boundary condition of M 2 is obtained using the fact that it is regular at t = 0. Explicitly, the coefficient of 1/t in the differential equation, which is a linear combination of some master integrals, is vanishing. Similarly, the boundary conditions of M 3 , M 5 and M 6 are derived from the regular condition at m W = 0, s = (m t + m W ) 2 and ut = m 2 t m 2 W , respectively. Since the letters are just polynomials in x, y, z, we can solve the differential equation recursively using the boundaries discussed above and write the final result in terms of the multiple polylogarithms [44], which are defined as G(x) ≡ 1 and The dimension of the vector (a 1 , a 2 , . . . , a n ) is usually called the transcendental weight of the multiple polylogarithm. We need multiple polylogarithms of at most transcendental weight four in the one-loop squared amplitudes.

Renormalization
The loop amplitude computed above contains ultra-violate (UV) and IR divergences. The UV divergences cancel with the contribution from counter-terms, arising from the renormalization of the couplings, masses and field strength.
In general, there are two methods to calculate the contribution of counter-terms. In one method, the bare Lagrangian is used to generate the Feynman rules and the loop amplitudes, and the external wave function renormalization is taken into account according to the Lehmann-Symanzik-Zimmermann (LSZ) reduction formula. We introduce the renormalization factor Z g 2 , Z b 2 , Z t 2 for the gluon field, bottom quark field and top quark field, respectively, by defining the renormalized fields as We do not consider the renormalization of the W boson field since we focus only on QCD corrections. Then the bare masses and couplings in the result should be replaced by the physical masses and couplings. The relation between the bare top quark mass m t,bare and the physical mass m t is m t,bare = Z t m m t . The strong coupling α s is renomalized in the MS scheme through the identity α bare The perturbative expansion of these renormalization factors is listed in the appendix B. As a result, the UV renormalized scattering amplitude up to two loops is where M bare is the bare amplitude and M C.T. denotes the amplitude whenever higher-order results of the renormalization factors are used. In the first line, we express the results in terms of bare quantities, while the second line contains only renormalized ones. There is not a one-to-one correspondence between the two lines, e.g., the third term in the second line could receive contributions from the second term in the first line. The other method is to use the Lagrangian in terms of renormalized fields and couplings, which includes new counter-term vertices. The top quark two-point counter-term vertex is The combination of the top quark propagator and the two-point counter-terms now becomes where we have used the notation δZ ≡ Z − 1. The omitted terms do not contribute to an NNLO calculation. Similar expressions can be obtained for the massless bottom quark by just neglecting all the mass terms. The gluon and quark three-point counterterm vertex is given by a multiplicative factor δ q At two loops, we would also need the gluon two-point and three-point counter-terms. Furthermore, no wave function renormalization factors are needed in the LSZ formula because physical fields have been employed. Therefore, the renormalized scattering amplitude up to NNLO is given by C.T. , and thus there is a one-to-one correspondence between the two lines in the above equation.
The two renormalization methods are equivalent, which can be understood by observing the cancellation of the wave function renormalization factors in eq.(3.11) and those in the combination of tree-level and counter-terms vertices 4 . As a result, only the renormalization factors associated with the external particles, as well as δZ αs and δZ t m , should be considered. In eq.(3.12), we have used a superscript to denote the order in α s for the perturbative expansion of the amplitude. Notice that M C.T. contains not only the tree-level diagrams with counter-terms at α 2 s but also the one-loop diagrams with counter-terms at α s . Since the one-loop integrals contain up to 1/ 2 poles, we need the counter-terms at O(α s ) up to O( 2 ) in M C.T. . They are also needed when we calculate the squared amplitude, e.g., M

Infrared divergences
After taking into account the contribution of counter-terms, the UV divergences disappear. But there are still IR divergences. These will cancel with the corresponding contribution of real corrections 5 in the prediction for a physical observable. This cancellation should be checked analytically before performing any reasonable numerical computation. Though the real correction is out of the scope of this paper, we can make a non-trivial check based on the understanding of the factorization property of amplitudes in the IR limit. More precisely, it has been proven that the UV renormalized amplitude can be written in the form where all the IR divergences are encoded in the factor Z and the remainder is finite as → 0 [45][46][47][48]. Following our convention, both the Z factor and the finite part have an expansion in α s , 4 An analogous equation to eq.(3.11) for the gluon propagator including counter-term vertices can be obtained in the Landau gauge. 5 At hadron colliders, renormalization of the parton distribution function is also needed. (3.14) Then it is ready to get The factor Z has been computed up to two-loop order for a general process with massive colored particles [47][48][49], and to three-loop order for single top production [50]. In the present paper, we are interested in the one-loop squared amplitudes, where the Z factor at one loop is given by The anomalous dimensions γ cusp , γ g , γ b and γ t can be found in the appendix of ref. [51]. Notice that we have chosen such a factorization scheme that Z does not have any finite part, which corresponds to the MS scheme. The presence of double poles in Z (1) requires the squared amplitudes M bare would get involved. Instead, we confine the cancellation to be numerical with high precision. In figure 2, we present our numerical results for these coefficients. In the plots, we have used the parameters β t and cos θ that are very suitable to describe the kinematics of the collision process in a finite range. More specifically, the top quark momentum is written as k µ 4 = E t (1, β t sin θ cos φ, β t sin θ sin φ, β t cos θ) in the center-of-mass frame of the incoming partons, where β t = 1 − m 2 t /E 2 t measures the velocity of the top quark and θ is the angle between the initial-state gluon and the final-state top quark.
It can be seen that the magnitudes of the coefficients span over several orders, especially growing or decreasing very fast as β t → 1. We have checked that the cancellation of these coefficients with the other parts in eq.(3.16) happens at the level of at least eight digits over almost the whole range of β t and cos θ 6 . Actually, it is not an easy task to achieve this. The main obstacle is the precise calculation of the very lengthy coefficients of the master integrals in the squared amplitudes, which contain high powers of s, t and u variables. And they are represented in terms of β t and cos θ as In our calculation, for each phase space point (β t , cos θ), we first compute the above s, t and u variables in high precision floating number, e.g., by keeping over 20 digits, and then transform them to rational numbers before inserting them to the expressions of the coefficients of the master integrals. We have observed a strong cancellation among different terms in each coefficient for small or large β t . The cancellation of all the divergences in eq.(3.16), either analytically or numerically, is an important check of the whole calculation.
At last, we turn to the finite part. Figure 3 shows the numerical result of M (1) fin 2 as a function of β t and cos θ. It can be seen that the finite one-loop squared amplitude increases remarkably when β t becomes large. The reason is that the mass of the top quark becomes vanishing in the limit β t → 1 7 . New collinear divergences, which have been regularized by m t and appeared as logarithms of 1 − β 2 t , will appear. From the right plot of figure 3, we observe that the finite part is insensitive to cos θ when β t is small. But the dependence of the finite part on cos θ becomes stronger with the increasing of β t . 6 We do not cover the phase space points with βt = 1, where the coefficients become divergent due to new collinear singularities. But in practice, we will not reach such a point since the top quark mass is not vanishing (or the collider energy is limited). 7 Or one may think that it is the region where the collider energy goes to infinity. This does not affect the discussion on the divergences. For comparison, we present the LO squared amplitude in the left plot in figure 4. We see that it is divergent only in the limit β t → 1 and cos θ → 1. This is the region when the tchannel propagator in figure 1 becomes on-shell. Notice that the two kinds of divergences in the limit β t → 1 are different, because the collinear divergences discussed above can appear for any values of cos θ. In practice, both of them are not harmful to the hadronic cross section due to the suppression of parton distribution function in the large β t region. After phase space integration, the LO partonic cross section is proportional to ln(1 − β t ), but the hadronic cross section is vanishing as β t → 1, as shown in the right plot of figure 4. This means that the suppression of the parton distribution function surpasses the enhancement of the partonic cross section. This is also the case at higher orders, since only ln n (1 − β t ) can appear in the partonic cross section. From the right plot of figure 4, we also observe that the suppression becomes weaker with the increasing of the collider energy.
To see the impact of the one-loop amplitude squared, we show the ratio R of |M (1) fin | 2 to |M (0) fin | 2 in figure 5. We find that the ratio R changes slowly away from the region of β t → 1 or cos θ → 1 and grows dramatically for large β t . For the phase space with β t < 0.9, which is the dominant region contributing to a real collider process, the corrections from one-loop amplitude squared are at the percent level, if the strong coupling is taken to be around 0.1.

W
The complexity of the loop amplitude mounts up prominently as the number of the independent scales increases. There are four independent scales in the process of tW production, i.e., m 2 t , m 2 W , s, t. It is beneficial to computation if one can expand the amplitude in one scale. For the process we are interested in, setting m W = 0 does not bring any new IR divergences over the full phase space due to the massive top quark propagator. Therefore, it is feasible to take the Taylor series expansion of the amplitude with respect to m 2 W . After setting m W = 0, one only needs to deal with integrals depending on three variables, i.e., s, t, m 2 t . As a consequence, both the IBP reduction and calculation of master integrals are simplified considerably. The square roots in eq.(3.3) will disappear, and thus the one-loop results become easier. Furthermore, at two-loop level, small mass expansion can make the calculation of some master integrals that would contain elliptic functions when keeping finite m W possible in terms of multiple polylogarithms. A similar method has been adopted in the analytic calculation of the two-loop gg → HH amplitude [52,53].
Given that the analytic result is still out of reach, we want to apply the method of Taylor expansion in m 2 W to compute the approximate two-loop amplitude of tW production. In this section, we demonstrate the validity of this method using the exact one-loop result we have obtained. As shown in eq.(2.2), the LO scattering amplitude contains the denominator s and u − m 2 t . One reasonable choice is to keep the invariant s and u and replace t → −s − u + m 2 t + m 2 W before expanding the amplitudes in a series of m 2 W . Regarding the calculation of master integrals, we explicitly show the expansion of I i n 1 ,n 2 ,n 3 ,n 4 defined in eq. .
(4.1) -13 - The differential operator of m 2 W should be written in a form that can directly act on a loop integral, with the coefficients Notice that these c i coefficients should be expanded further in a series of m 2 W before taking the limit m 2 W → 0 in eq.(4.1). Actually, this differential operators can act not only on the master integrals, but also on the whole squared amplitude. However, we do not apply them to the variable m 2 W in the summation of the polarizations of the W boson in eq.(2.3). This denominator arises because the longitudinal polarization state of the W boson has a different coupling to quarks from the transverse ones. If we separate the contribution from different polarizations and express the coupling of the longitudinal polarization state and quarks in terms of the top quark Yukawa coupling, there is no such m 2 W factor in the denominator, according to the Goldstone boson equivalence theorem [54,55]. This property can also be seen if we use Feynman-'t Hooft gauge to sum the polarizations of the W boson and add the contribution of the Goldstone boson production diagram, which does not involve a polarization sum. After applying the differential operator to master integrals, each term in the series of m 2 W can be expressed as a linear combination of scalar integrals. Then we reduce them to a set of master integrals with m W = 0. The advantages of the expansion come from two aspects. First, the number of these master integrals is less than that before expansion. Second, the analytic computation is easier to perform. For example, the square root in eq.(3.3) does not exist at m W = 0. Further, the set of master integrals is the same at all orders in the expansion, therefore no additional computation cost needs to be paid if we want to obtain the result of even higher orders in the expansion of m 2 W in principle. Figure 6 demonstrates the accuracy of the m W expansion, taking the real part of I 1 1,1,1,1 as an example. We compare the expanded results for the coefficients of i with i = −2, · · · , 2 to the exact ones. As more higher orders in m 2 W are included, the expansion approximation becomes better. The difference between the expansion up to O(m 6 W ) and the exact values are at the permille level or better. We also find that expansion series converge faster in the larger β t region. This is reasonable because the expansion parameter m 2 W /s is generally smaller in this region. In the plots, there are some regions where the curves seem to be divergent. We have checked the absolute values of both the expanded and exact results, and found that they are almost vanishing in these regions.
Then we consider the expansion of the squared amplitudes. As mentioned above, the differential operator can be applied to the squared amplitude. This can be done just after creating the amplitude from the Feynman diagrams. Firstly, we replace k 4 by k 1 + k 2 − k 3 everywhere. Secondly, we contract the Lorentz indices in the spin summed squared amplitude. Thirdly, we apply the differential operator in eq.(4.2) to the squared amplitude 8 , and reduce all the scalar integrals to master integrals with k 2 3 = 0. However, we will use another method. Since we have expressed the squared amplitude as a combination of the master integrals and their coefficients, we only have to expand the coefficients in m 2 W , given that the master integrals have been expanded as discussed above. Figure 7 shows |M (1) bare | 2 and their expansions in m 2 W with cos θ = 0.5 as a function of β t . It can be seen that the differences between the exact results and the expanded ones up to O(m 4 W ) are already notably small. The corrections of order m 6 W are negligible. We have also tried to expand only all the master integrals but keep the exact form of the coefficients of the master integrals. One may expect that this approach would have a better approximation to the exact result. On the contrary, we find that the expansion is worse, sometimes even seems divergent in small β t regions. As mentioned before, the coefficients contain a lot of terms, in which the powers of the kinematic variables are very high. For specific values of the kinematic variables, the cancellation among these terms happens over about ten digits. In the expansion in m 2 W , this cancellation holds at each order. If the coefficients contain higher orders in m 2 W than the master integrals, the cancellation is not complete. As a result, the deviation between the expanded results and the exact ones may be large.

Conclusions
The precise prediction of tW associated production at a hadron collider is important to study the properties of the top quark and the W boson. We have computed the analytic results for the one-loop squared amplitudes for this process using canonical differential equations. The renormalized amplitude squared has up to −4 poles, which have been checked against the general infrared structures predicted by anomalous dimensions. The finite part gives rise to about a few percent corrections compared to the corresponding LO results. We also investigate how to obtain approximated results using the method of expansion in m 2 W , and find that the expanded results up to O(m 4 W ) are already accurate enough to present the exact results. It is promising to apply this expansion method in calculating the interference of two-loop and tree level amplitudes.

B Perturbative expansion of renormalization factors
The strong coupling α s is renomalized in the MS scheme. The renormalization factor is given by [56,57] Z αs = 1 − α s 4π β 0 + α s 4π In QCD, we have C A = 3, C F = 4/3, T F = 1/2. Here n f is the total number of quark flavors, including both massless and massive quarks. As a result, the running strong coupling satisfies the renomalization equation