Contribution of third generation quarks to two-loop helicity amplitudes for W boson pair production in gluon fusion

We compute the contribution of third generation quarks (t, b) to the two-loop amplitude for on-shell W boson pair production in gluon fusion gg → WW. We present plots for the amplitude across partonic phase space as well as reference values for two kinematic points. The master integrals are efficiently evaluated by numerically solving a system of ordinary differential equations.


Introduction
Production of a W boson pair in gluon fusion, gg → W W , is a loop-induced process.Although it is expected to be strongly suppressed compared to qq → W W , there are two reasons that make it relevant.First, the large gluon flux at the LHC nearly compensates for the suppression by the strong coupling α s when compared to quark antiquark annihilation qq → W W . Second, event selection disfavours events with large longitudinal boosts that are due, primarily, to the qq → W W process [1].Hence, good understanding of the gg → W W process is needed for a reliable description of W pair production at the LHC.
The current situation is as follows.One-loop, leading order (LO) cross sections for W pair production in gluon fusion have been computed long ago [1][2][3][4].More recently, next-toleading order (NLO) QCD corrections mediated by massless quark loops were computed in refs.[5,6] using two-loop amplitudes calculated in refs.[7,8].However, these calculations ignored the contribution of the third quark generation (t, b) since top quarks cannot be treated as massless.The goal of this paper is to compute the contribution of third generation quarks to the two-loop amplitude for the gg → W W process, providing a prerequisite for improved theoretical description of W pair production at the LHC.
Massive fermion propagators that appear when top quark contributions are considered make calculations significantly more demanding compared to the massless case.Indeed, in the massive case the variety of integrals one has to consider is larger and they are more difficult to compute.
We rely on integration-by-parts (IBP) [9] identities to find linear relations between loop integrals and express the gg → W W amplitude in terms of a few (master) integrals.The Laporta algorithm [10] ensures that the system of IBP equations closes.However, multi-scale integral reductions are computationally expensive and appear to be infeasible for gg → W W with current publicly available software [11,12].
To overcome this problem, we set the mass of the top quark m t and the mass of the W boson m W to integers close to their current experimental values.Lowering the number of parameters makes the IBP reduction possible.We have used Kira [12] for the reduction as well as LiteRed [13] and Reduze [11] to find symmetry relations between integrals.
We evaluate the master integrals numerically.A widely used systematic method is that of numerical integration enabled by sector decomposition [14] which, however, is computationally expensive.Another possibility is to solve numerically a system of differential equations satisfied by the master integrals [15,16].Recently, a new method to do this was presented in ref. [17] co-authored by one of the present authors.This method is particularly suitable for problems with massive particles in the loops.We use this method to solve a system of differential equations with respect to the m 2 t variable by moving from m 2 t → −i∞ to its physical value. 1 The advantage of this method is that it allows for, essentially, arbitrary precision at low computational expense.
The remainder of the paper is organised as follows.Section 2 provides definitions of kinematic variables as well as conventions regarding the γ 5 -matrix and renormalisation of ultraviolet (UV) and infrared (IR) singularities.In section 3 we discuss Feynman diagrams involved in the calculation at one and two loops, the colour structure, and the integral reduction.A detailed presentation of the numerical approach to the evaluation of the master integrals is given in section 4. In section 5 we evaluate the gg → W W amplitude at two phase space points and present plots of the amplitude across partonic phase space.We conclude in section 6.

Definitions
We study the contribution of third generation quarks (t, b) to the amplitude of the process keeping the exact dependence on the top quark mass m t while treating the bottom quark as massless.We only consider the case where both W bosons couple directly to the quark loop.Indeed, the single-resonant contribution of an intermediate Z can be ignored as it vanishes for on-shell W pair production [1].The process involving an intermediate Higgs boson is known [18,19].
We write the gg → W W amplitude as where a 1,2 are the colour indices of the external gluons, g W = e/ sin θ W is the weak coupling constant, and j are the polarisation vectors of external particles.We consider the CKM matrix to be an identity matrix.
We set all particles on-shell, and introduce Mandelstam variables in a standard way (2.4) These variables satisfy the relation s + t + u = 2m 2 W .We decompose the gg → W W amplitude into 38 tensor structures, The tensor structures T µν I are defined in ref. [8] and are parity-even, while S µν I are parityodd and are defined in ref. [1].Our goal is to calculate the form factors A I .
The W qq-vertex contains vector and axial parts The gg → W W amplitude can hence be written as a sum of vector-vector, axial-vector, and axial-axial contributions.The vector-vector and axial-axial terms are parity-even and can be decomposed in terms of tensor structures T µν I in eq.(2.5).On the other hand, the axialvector term is odd under parity transformations; its decomposition is possible in terms of the tensor structures S µν I in eq.(2.5).If masses of two quarks in a single generation are equal, the parity-odd contribution vanishes [1][2][3]20] and it is therefore absent in amplitudes involving only massless quark loops [7,8].
To deal with the axial part of the vertex (2.6), we employ the γ 5 -prescription of refs.[21,22] and replace γ µ γ 5 with Throughout this paper the Levi-Civita symbol ε µνρσ is defined using the convention of FORM ε 0123 = −i.Although this γ 5 -prescription is much more convenient to work with in dimensional regularisation, it violates Ward identities of the axial current.To restore the Ward identity, we have to perform a finite renormalisation [23], where J A and J A b stand for renormalised and unrenormalised axial currents respectively.

Pole structure
Throughout the calculation we employ dimensional regularisation and set the space-time dimensionality to d = 4 − 2 .The singularities ubiquitous in loop calculations appear as poles in of either ultraviolet (UV) or infrared (IR) origin.UV poles are absorbed through the introduction of renormalisation factors leading to a renormalised amplitude.Expanding the unrenormalised amplitude A b in the bare strong coupling α 0 s we have (2.9) Following [24,25] we employ a hybrid renormalisation scheme where the gluon field G µ and the top quark mass m t are in the on-shell scheme while the strong coupling constant α s is renormalised in the MS scheme.We have where µ is the renormalisation scale and S = (4π) − e γ E .The renormalisation constants are expanded in the coupling The renormalised amplitude is related to the unrenormalised one by ) b ( , m t ), (2.13) b ( , m t ) + m t Z (1)  mt C (1) b ( , m t ) where we used the fact that the tree level amplitude A (0) b vanishes.Note that while the renormalisation factors for the strong coupling and the gluon wave function are multiplicative at the level of the amplitude, the factor for the mass renormalisation is not.For this reason, the mass counterterm C (1) is calculated separately.
The relevant renormalisation factors are [26][27][28][29][30]] ) where γ g (n l + 1) = 11 6 C A − 2 3 T F (n l + 1) and n l is the number of massless fermions.The divergent contribution of the massive quark flavour in γ g (n l + 1) cancels between the gluon field and the coupling renormalisation, while the light quark contributions are absorbed into Catani's operator, as explained below.
In general, IR poles of loop amplitudes cancel against contributions of real emission processes.In the virtual amplitudes the IR singularities factorise in a universal manner [31], this allows us to write the amplitude as a sum of IR poles and a finite remainder.
Since gg → W W vanishes at tree level, the pole structure of the two-loop amplitude is particularly simple.Furthermore, the amplitude is a singlet in colour space.Hence, we can write the renormalised amplitude as where F (2) is the finite remainder.The Catani operator reads where In the following we will set In order to obtain the finite remainder of the two-loop amplitude, we require the one-loop amplitude, A (1) , expanded through O( 2).

Amplitude calculation
In this section we discuss the calculation of the amplitudes A (1) , C (1) , and A (2) .

One loop
We generate 8 one-loop diagrams using QGRAF [32] and perform colour and Dirac algebra as well as projection of the form factors shown in eq.(2.5) using FORM [33][34][35].Two triangle diagrams with g → W W transition vanish due to colour conservation.The remaining six box diagrams can be mapped to 5 independent topologies.We perform the integral reduction step using IBP identities, and express the form factors as linear combinations of 16 master integrals.The form factors for the mass counterterm amplitude C (1) are computed in a similar fashion.The one-loop amplitudes were checked against FeynArts and FeynCalc [36-39].

Two loops
At two loops we generate 136 diagrams using the same steps as above.To test these steps of our implementation, a subset of unreduced diagrams were crosschecked numerically against FeynArts and FeynCalc.In figure 1 we show a few representative diagrams that contribute to the two-loop gg → W W amplitude.
Since the weak bosons are not colour charged, the colour structure follows that of the gluon self-energy two-loop diagrams with a closed quark loop.These diagrams are given in table 1.All 136 diagrams contributing to gg → W W can be formed by attaching two W bosons to the quark loop and by pinching propagators.This results in four basic topologies shown in figure 2

Class Colour factor
Note that both A (2), [1] and A (2),[−1] are gauge invariant.We also observe that A (2),[−1] is finite after mass renormalisation and has no infrared poles.In order to perform an IBP reduction we express the amplitude in terms of integral families.An integral family is a set of propagators and irreducible scalar products (ISPs) that forms a basis of the linear space spanned by all scalar products containing loop momenta.For four-point kinematics in four dimensions there are 9 independent scalar products at two loops.Before the integral reduction the form factors can be written as In eq.(3.2) N T is the number of independent integral families and a T = (a 1 , . . ., a 9 ).Each a i is an integer power of the 9 independent scalar products in the integral family T .Using symmetries between diagrams, including crossing symmetry, we find a total of N T = 35 families of type (a), (b), and (d).The coefficients of integrals, c IT a T , are rational functions of s, t, m t , m W and space-time dimensionality d.
The 35 families can be further reduced to 25 two-loop irreducible families and a single one-loop squared family using IBP identities.Family definitions are given in appendix A and the corresponding topologies are shown in figure 3. We use Kira [12] to reduce all integrals I T ( a T ) appearing in the scattering amplitude to a set of 334 master integrals.
To complete the integral reduction using reasonable time and resources, several measures have been taken.First, the integrals depend on the kinematic variables s and t, the masses m t and m W , and space-time dimensionality d.Keeping all of them parametric makes the reduction formidably complicated.To overcome this problem, we set m t = 173 GeV and m W = 80 GeV in the IBP reduction, keeping only s, t, and d as parameters. 2 This simplifies the reduction tables and cuts down the run times considerably.
Second, the size of reduction tables can be reduced further by a careful choice of master integrals.Our guiding principle in choosing master integrals is absence of denominators with non-factorisable dependence on kinematic variables and space-time dimensionality d as well as avoiding denominators that lead to poles at non-integer values of d (e.g.3d − 10) [41,42].Furthermore, we aim at having simple differential equations for fast numerical evaluation.By trial and error, we find that master integrals with at most one squared propagator or a single irreducible scalar product are sufficient to satisfy the above requirements.In a few sectors integrals with two squared propagators are needed, but this appears to be an exception rather than the rule.Third, we use the select masters reduction feature implemented in Kira to project integrals onto one master integral at a time.Our final reduction tables are of the order of hundred MB each, which makes applying reduction rules to the amplitude a manageable job.

Differential equation
Having expressed the full amplitude through master integrals, we need to evaluate them.The master integrals are defined as follows where D i are denominators in one of the 26 families given in appendix A. Their topologies are shown in figure 3. Note that we absorb a factor of −i(4π) 2− e γ E per loop into the definition of master integrals.
To evaluate all 334 two-loop master integrals with massive propagators, we employ the imaginary mass method proposed in ref. [17].In the original formulation of this method, an imaginary mass term −iη is added to all propagators.The differential equation with respect to η is solved numerically starting at η → ∞ and progressing towards a physical point η = 0 + .The boundary conditions involve vacuum integrals with mass m 2 = −iη, as no physical masses or kinematic variables survive in this limit.Adding −iη to massless propagators alters the behaviour of integrals near the physical point η = 0 + and generates a singularity in the differential equation.In order to compute the result at the physical point it is then necessary to fix constants in a formal solution of the differential equation, by matching against another point within its radius of convergence.
We employ a variant of the original method, where the imaginary mass parameter −iη is introduced only to the massive propagators.Setting the mass m 2 t of massive propagators to m 2 t − iη requires no extra work as far as IBP reductions are concerned. 3Since our diagrams already have many massive propagators, the boundary conditions at η → ∞ are remarkably simple, with only 5 planar integrals and 1 nonplanar, massless 3-point integral needed to fix all master integrals at the boundary.The topologies of the boundary integrals are shown in figure 4, see appendix B for their definitions [43][44][45][46][47].
By constructing a differential equation with respect to m 2 t and choosing the boundary condition at m 2 t → −i∞, the differential equation can be used to evaluate master integrals at the physical mass, m 2 t = (173 GeV) 2 .At η → ∞, the master integrals receive contributions from 3 regions: 1.All internal momenta are comparable to m 2 t → −i∞.
2. Some internal momenta that form a closed loop are comparable to m 2 t → −i∞, while the remaining momenta are much smaller than m 2 t and are comparable to other kinematic parameters (i.e. s, t, m 2 W ).
3. All internal momenta are much smaller than m 2 t → −i∞.
In figure 5 we show a typical master integral and its regions.The boundary condition in each region can be expressed in terms of the boundary integrals given in figure 4 through IBP reduction, together with an overall scaling factor of m t .Note that if we add −iη to all propagators, only the first region contributes to the integrals.In general, the fewer propagators one takes to be infinitely massive, the more complicated boundary conditions one needs to consider.However, since we change the original integrals less, the singularity at the physical point, η = 0 + , is simpler.In fact, for the problem at hand, the differential equation is regular at the physical point and no additional complications arise.
Having fixed all boundary conditions, we are ready to solve the differential equation.We write the differential equation in terms of a dimensionless variable  and make all master integrals dimensionless using m 2 W .At each phase space point, we solve the differential equation numerically by moving along the positive imaginary axis from the boundary at x = −i∞ to the physical point x = 0.This involves three steps.
1. First, we transform the differential equation at x = −i∞ to a Fuchsian form where y = 1/x and I is the vector of master integrals.
2. Second, we use the differential equation to obtain power-logarithmic expansions of the master integrals I in terms of y in the neighbourhood of x = −i∞ or y = 0.For each master integral, we write where the numerical coefficients c ijkl are completely determined by the differential equation and the boundary conditions.The parameter N in the above equation is the desired order of the expansion; it controls the truncation error and, eventually, the precision with which I i is computed.Parameter M is the maximal power of in the series.
3. Finally, using the expansion (4.4), we move to a regular point x 0 within the radius of convergence by direct evaluation of I(y = 1/x 0 ).Then at each regular point x i along the path of integration, we Taylor expand the master integrals by expanding the equation up to order N around ) where x = x − x i .Once this is accomplished, we move on to the next point x i+1 within the radius of convergence of the new series (4.6).By repeating this expandevaluate operation, we finally arrive at the physical point x = 0. Figure 6 shows a typical situation in the complex x-plane.
There are several advantages of this method.First, it allows us to evaluate master integrals to arbitrary precision in reasonable and predictable time, which can be difficult to do using other numerical methods.The possibility to increase precision is crucial for a stable evaluation of the amplitude in quasi-singular regions, e.g.around thresholds of internal particles.Second, given an equation and a valid path of analytic continuation in the complex plane, this method produces identical results every time.This determinism makes the calculation reproducible and allows one to keep numerical errors under control.Finally, this method is fast enough for practical applications.The run time depends on the form of the equation, requested precision, depth of the -expansion, and working precision used in the calculation.However, for a result accurate to 15 digits, typical run time is about 1-10 seconds per integral, depending on the form of the equation.Evaluating all master integrals at a typical phase space point takes less than an hour on a single CPU core.
We crosschecked the evaluation of master integrals against pySecDec [48,49] and FIESTA [50].We also checked the self-consistency of the differential equations by comparing results obtained by integrating along different paths.Specifically, all master integrals have been checked against pySecDec at an unphysical phase space point up to the default precision of pySecDec (3-10 digits).Some of the master integrals have also been checked against FIESTA at another unphysical phase space point This plot shows singularities and steps of the family planar no. 1 at s = (500 GeV) 2 , t = −(300 GeV) 2 in the complex plane of x.The red crosses are the singularities of the differential equation, while the blue dots are the steps used to solve the differential equation.We approach the origin from −i∞ along the positive imaginary axis.The boundary at x = −i∞ is not shown on this plot.
In addition, evaluations at two different phase space points should be connected by a system of differential equations with respect to s and t.We pick the following two phase space points, Taking the evaluations at (s 2 , t 2 ) as boundary condition and running the equations from (s 2 , t 2 ) to (s 1 , t 1 ), we then check against a direct evaluation at (s 1 , t 1 ).We find that master integrals evaluated at (s 1 , t 1 ) in the two different ways agree up to the precision used when solving the equations (15 digits in this particular case).

Numerical evaluation
We parametrise the phase space of the W bosons using the angle θ between p 1 and p 3 in the centre of mass frame, see figure 7, and the relative velocity of the W boson pair, β.
These quantities are related to Mandelstam invariants through the following equations (5.1) Figure 7: Illustration of the scattering angle θ between incoming momentum p 1 , parallel to the z-axis, and outgoing momentum p 3 .
For further references, we present helicity amplitudes evaluated at two phase-space points P2: This corresponds to √ s ≈ 185 GeV and θ ≈ 102 • and √ s ≈ 367 GeV and θ ≈ 37 • for P1 and P2 respectively.
To evaluate the tensor structures in (2.5) we construct polarisation vectors for the two gluons using spinor-helicity formalism.The polarisation vectors are given by µ Keeping in mind that, eventually, we will be interested in the decay of the W bosons into leptons, we express the polarisation vectors of on-shell W boson states through a current that describes W − → eν and W + → eν transitions The massless momenta are constructed by flattening the massive momenta Table 2: Massless incoming and outgoing momenta in units of GeV in the centre of mass frame for the phase space points defined in eq. ( 5.2) and eq.( 5.

3).
A Table 3: Evaluation of the two independent helicity amplitudes at one loop for the phase space points defined in eq. ( 5.2) and eq.( 5.3).
We label the helicity amplitudes by the helicities of the two incoming gluons and two of the out-going leptons λ 1 λ 2 λ 5 λ 7 , where λ i = L, R. W bosons are left-handed and there are four helicity configurations for the gluons.Two operations allow us to establish relations between these configurations.First, we can flip all helicities simultaneously by complex conjugation of the polarisation vectors.Second, we can flip the helicities of the W bosons only by swapping the momenta (p 5 ↔ p 6 and p 7 ↔ p 8 ) in the currents of eq.(5.6).Hence, only two helicity configurations are independent, we choose to present LLLL and LRLL.
One-loop helicity amplitudes for the two phase space points defined in eq. ( 5.2) and eq.( 5.3) are given in table 3.
We present -expansions of the two-loop amplitudes for leading and sub-leading colour in tables 4 and 5 respectively.Comparison with the predicted structure of IR poles is also shown.The renormalisation scale µ is set to 2m W .We note that the renormalised two-loop amplitudes (2.14) receive a trivial contribution from the counter term amplitude.It is independent of N C and we do not include it here.

Two loops LLLL
, [1] /A -Table 4: Evaluation of the two-loop helicity amplitudes, LLLL and LRLL, for the phase space points defined in eqs.(5.2) and (5.3) for the leading colour contribution.We normalise by the one-loop amplitude A (1) | =0 and show the infrared pole structure for comparison.We set the renormalisation scale µ = 2m W .
Two loops 1.48368538287541 + 1.38326340829964i Table 5: Evaluation of the two-loop helicity amplitudes for the phase space points defined in eqs.For the one-loop amplitudes we construct a uniform, dense 99 by 99 grid in terms of the variables β and cos θ defined in (5.1) with step sizes of 0.01 and 0.02 in the ranges [0.01, 0.99] and [−0.98, 0.98] respectively.The absolute value of the two independent helicity amplitudes are plotted in figure 8.We stress that the helicity amplitudes presented here depend on the polarisation vectors of the on-shell W bosons, see eqs.(5.7) and (5.8).To avoid this one can project onto helicity dependent form factors defined in refs.[1,7,8].
For the two loop amplitude we use a sparse grid for the bulk of phase space, 0.1 ≤ β < 0.8.The step size in β is 0.1 and cos θ ranges from −0.8 to 0.8 in steps of 0.2 with an additional border at cos θ = ±0.96.For the production threshold, 0.01 ≤ β < 0.1 we use a step size of 0.01 and same resolution for cos θ as in the bulk region.In the high-energy region 0.8 ≤ β ≤ 0.99 we also use the step size of 0.01 for β, but increase resolution in cos θ with a step size of 0.04 in the range from −0.96 to 0.96.
In total 1156 points have been computed to produce plots for the two-loop helicity amplitudes.In figure 9 and 10 we plot the interference of the finite remainder with the   (5.9) For illustration purposes we only show the vector-vector and axial-axial part of the amplitudes in these plots.

Conclusions
In this paper we computed the contribution of the third generation quarks to the twoloop helicity amplitudes for W boson pair production in gluon fusion.We use projection operators to obtain form factors that can be calculated using integration-by-parts integral reduction.To overcome the computational bottleneck of the reduction step, we fix the masses of the top quark and the W bosons to integer numbers close to their experimentally determined values.The master integrals are evaluated numerically by solving a system of differential equations with respect to the top mass parameter.This approach allows for arbitrary precision and is especially efficient for processes involving massive internal particles.
The present calculation opens up the possibility to include the contribution of third generation quarks into NLO QCD corrections to the cross section of W pair production in gluon fusion.More generally, this method can be used for numerical calculations of many loop amplitudes with massive particles.As higher order virtual corrections to many processes involving internal masses are currently beyond the reach of analytic methods, this approach represents an alternative way forward.

B Boundary condition of the differential equation
The explicit expressions for the boundary integrals in figure 4 are listed below [43][44][45][46][47], where q 2 corresponds to the four-momentum squared of the external legs.Note that I 1 and I 3 are one-loop integrals, thus they enter the boundary condition through the products among themselves. )

Figure 1 :
Figure 1: Representative two-loop Feynman diagrams for gg → W W .

Figure 2 :
Figure 2: The four basic topologies.The internal lines can be both massive and massless.The first three (a)-(c) are planar and can at most have 5 massive internal lines, while the nonplanar topology (d) can at most have 4 massive internal lines.

Figure 3 :
Figure 3: Topologies of integral families.Solid and dashed lines correspond to massive and massless particles respectively.Internal massive particles have mass m t while external massive particles have mass m W .All nine planar topologies and nonplanar no. 2 and 3 can be crossed (p 1 ↔ p 2 ) giving a total of 26 topologies.

(a) I 1 (b) I 2 (c) I 3 (d) I 4 (e) I 5 (f) I 6 Figure 4 :
Figure 4: Topologies of boundary integrals.Solid and dashed lines correspond to massive and massless particles respectively.See appendix B for their explicit expressions.

Figure 5 :
Figure 5: A typical master integral and its leading regions.Solid and dashed lines correspond to massive and massless particles respectively.Internal massive particles have mass m t , while external massive particles have mass m W .

50 Figure 6 :
Figure6: A typical path to solve the differential equation at a certain phase space point.This plot shows singularities and steps of the family planar no. 1 at s = (500 GeV) 2 , t = −(300 GeV) 2 in the complex plane of x.The red crosses are the singularities of the differential equation, while the blue dots are the steps used to solve the differential equation.We approach the origin from −i∞ along the positive imaginary axis.The boundary at x = −i∞ is not shown on this plot.

Figure 8 :
Figure 8: Absolute value of the vector-vector plus axial-axial part of the one-loop helicity amplitudes.

Figure 9 :
Figure9: Finite remainder of the vector-vector and axial-axial part of the two-loop, leading colour (N C ), helicity amplitudes interfered and normalised by the leading order amplitude, see eq. (5.9).We set the renormalisation scale µ = 2m W .

Figure 10 :
Figure 10: Finite remainder of the vector-vector and axial-axial part of two-loop, subleading in colour (1/N C ), helicity amplitudes interfered and normalised by the leading order amplitude, see eq. (5.9).We set the renormalisation scale µ = 2m W .

Table 1 :
The classification of diagrams in colour structures follows that of the two-loop gluon self-energy diagrams involving a closed fermion loop.This motivates splitting the amplitude into leading and sub-leading colour.nonplanarcontributions.We find 33 non-vanishing diagrams in class L (of which 17 are nonplanar), 20 in class S and 40 in LS.This classification motivates splitting the amplitude into leading (N C ) and sub-leading (1/N C ) colour contributions, )