# N$^3$LO Higgs and Drell-Yan production at threshold: the one-loop two-emission contribution

## Abstract

In this paper, we study phenomenologically interesting soft radiation distributions in massless QCD. Specifically, we consider the emission of two soft partons off of a pair of light-like Wilson lines, in either the fundamental or the adjoint representation, at next-to-leading order. Our results are an essential component of the next-to-next-to-next-to-leading order threshold corrections to both Higgs boson production in the gluon fusion channel and Drell-Yan lepton production. Our calculations are consistent with the recently published results for Higgs boson production. As a non-trivial cross-check on our analysis, we rederive a recent prediction for the Drell-Yan threshold cross section using a completely different strategy. Our results are compact, valid to all orders in the dimensional regularization parameter, and expressed in terms of pure functions.

^{1}

The total cross section for Higgs boson production in the gluon fusion channel and the total cross section for Drell-Yan lepton production are two of the most important Large Hadron Collider-relevant, infrared-collinear (IRC) finite observables in QCD. Due to their relatively simple kinematics, gluon fusion Higgs and Drell-Yan lepton production are benchmark processes for the application of perturbative QCD techniques. A great deal of progress along these lines has been made over the past two decades and, as a consequence, total cross sections for these processes are known to next-to-next-to-leading order (NNLO) NUPHA.B157.461 (); NUPHA.B359.283 (); NUPHA.B359.343 (); PRLTA.70.1372 (); hep-ph/9504378 (); hep-ph/9511344 (); hep-ph/9705337 (); hep-ph/0102227 (); hep-ph/0201206 (); hep-ph/0207004 (); hep-ph/0302135 (); hep-ph/0509189 (). The calculation of these IRC finite observables to one order higher in the strong coupling constant is highly desirable due to the large -factors and theoretical uncertainties inherent in the calculation of such hadron collider observables NUPHA.B359.283 (); PRLTA.70.1372 (); hep-ph/9504378 (); hep-ph/9511344 (); hep-ph/9705337 (); hep-ph/0102227 (); hep-ph/0201206 (); hep-ph/0207004 (); hep-ph/0302135 (). Unfortunately, this task is also highly non-trivial; at next-to-next-to-next-to-leading order (NLO), the current state-of-the-art, one must compute three-loop virtual corrections, single-emission, two-loop real-virtual corrections, double-emission, one-loop real-virtual corrections, triple-emission real corrections, and collinear counterterms for the parton distribution functions. All of these ingredients have severe infrared and/or collinear divergences which complicate their calculation and IRC finiteness is only achieved once all contributions have been appropriately combined.

While full fixed-order calculations at NLO remain a challenge, significant progress has been made over the last few years. For gluon fusion Higgs production, the Wilson coefficients obtained by integrating out top quark loops to construct an effective Lagrangian are known to three loops NUPHA.B510.61 (); JHEPA.0601.051 (); NUPHA.B744.121 () and the relevant ultraviolet renormalization constants have been worked out in Refs. PHLTA.B93.429 (); PHLTA.B303.334 (); PHLTA.B400.379 (); NUPHA.B710.485 (). The collinear counter terms at three loops are also known and have been obtained by the authors of Refs. JHEPA.1211.062 (); PHLTA.B721.244 (); JHEPA.1310.096 (). The three-loop virtual corrections are under control, as they can be extracted from the three-loop quark and gluon form factors computed in Refs. 0902.3519 (); 1001.2887 (); 1004.3653 (); 1010.4478 (). For finite top quark mass, an estimate of the NLO cross section for Higgs production was made in Ref. 1303.3590 () by arguing that one can exploit information encoded in the soft gluon and high energy limits to construct an accurate approximation to the full result.

In recent years, rapid progress has been made towards the analytical calculation of the relevant phase space integrals by expanding them
in the threshold limit^{2}

We describe our calculational strategy below. At a hadron collider with incoming hadrons and and center-of-mass energy , the inclusive cross section for producing a Higgs boson of mass or Drell-Yan lepton pair of invariant mass can be written as

(1) |

where or and the summation is over all partonic channels. To save space, we will often use the abbreviations and when discussing the processes that we consider in this paper. The parton distribution function of parton in hadron , , depends on the momentum fraction, , and the factorization scale, . , the partonic cross section, depends on the partonic center-of-mass energy, , the parameter defined by Eq. (NLO Higgs and Drell-Yan production at threshold: the one-loop two-emission contribution), and the renormalized strong coupling constant at renormalization scale , . In this work, we take for the sake of simplicity. In full fixed-order calculations, the dependence of partonic cross sections on and is simple but their dependence on is quite non-trivial. Throughout this paper, we shall be concerned only with a threshold expansion in the limit of the . In this regime, additional partons radiated into the final state are constrained to be soft by virtue of the phase space constraint and the partonic cross sections for both and simplify dramatically:

(2) | |||||

(3) |

For simplicity, we keep the running coupling dependence of the functions introduced on the right-hand side of Eqs. (2) and (3) implicit.
In dimensions, and are given respectively
by^{3}

(4) |

where is the effective coupling of the Higgs boson to gluons in the limit of infinite top quark mass Shifman:1979eb (),

(5) |

is the number of colors, is the fine structure constant at renormalization scale , and is the electric charge of quark . It is worth pointing out that has a perturbative expansion in the renormalized strong coupling constant of its own (see Ref. NUPHA.B744.121 () for explicit expressions in a slightly different normalization) and, throughout this work, we employ the conventional dimensional regularization scheme. The coefficient functions and contain singular distributions of the form

(6) |

where the plus distribution acts on functions regular in the limit as

(7) |

The functional form of the coefficient functions can be obtained by taking the limit of the full partonic cross section. Alternatively, they can be directly calculated in the soft limit of QCD, using a factorization formula derived from the soft-collinear effective theory of QCD hep-ph/0005275 (); hep-ph/0011336 (); hep-ph/0109045 (); hep-ph/0206152 (). At our chosen scale, the factorize into products of the form hep-ph/0605068 (); 0710.0680 (); 0808.3008 ()

(8) |

where the are hard functions encoding the virtual corrections at threshold and the are renormalized soft functions encoding the real radiative corrections at threshold. At first, it might seem peculiar that we have not taken the soft functions in Eq. (8) to be a functions of , where is twice the energy of the soft QCD radiation in the final state. However, in threshold kinematics, we have

(9) |

and we see that the form of Eq. (8) follows naturally from Eq. (9).

In general, hard functions in soft-collinear effective theory are complex squares of hard Wilson coefficients obtained by matching QCD onto soft-collinear effective theory. The perturbative expansions of the Wilson coefficients relevant to Higgs and Drell-Yan production are known to the order required for an NLO calculation of the . They can be extracted from the three-loop quark and gluon form factors 0902.3519 (); 1001.2887 (); 1004.3653 (); 1010.4478 (). In fact, they are written down explicitly in Ref. 1004.3653 (), Eqs. (7.4) - (7.9). The actual hard coefficients that we need for our analysis are derived by taking the complex square of Eq. (7.3) in that work, setting , and then expanding in to third order.

The bare soft functions can be written as squares of time-ordered matrix elements of pairs of semi-infinite Wilson line operators,

(10) |

As mentioned above, the soft functions defined in Eq. (10) depend on the ratio of twice the energy of the soft QCD radiation to (i.e. on ), as well as the bare strong coupling constant . is a conventional normalization constant,

(11) |

depending on whether the soft Wilson lines are in the adjoint or the fundamental representation of . The summation in Eq. (10) is over all possible soft parton final states, . The operator acts on the final state according to

(12) |

where is the energy of the soft radiation in final state . The Wilson line operators, and , are respectively defined as in-coming, path-ordered and anti-path-ordered exponentials hep-ph/0412110 (),

(13) |

In the above, , where the are either adjoint (for ) or fundamental (for ) matrices. As usual, and are light-like vectors whose space-like components are back-to-back and determine the beam axis. For a generic four-vector, , we have

(14) |

given the usual definitions and .

In the perturbative regime, we can calculate the bare soft functions order-by-order in . They admit expansions of the form

(15) |

In Eq. (15), is the usual factor,

(16) |

where is the Euler constant and is the ’t Hooft scale introduced from continuing the space-time dimension to dimension. After renormalization in the scheme, it will be replaced by the renormalization scale, . The factor can be expanded in in terms of the plus distributions introduced above (Eqs. (6) and (7)),

(17) |

The expansion coefficients of the bare soft functions, , are functions of and only. In the case of , the one- and two-loop results are well-known and
were calculated to all orders in in Ref. hep-ph/9808389 ().
The structure of the eikonal approximation is such that, through two-loop order, the color degrees of freedom completely
factorize from the part of the bare expansion coefficients which encodes the non-trivial soft dynamics; by simply replacing
everywhere with , the results reported in hep-ph/9808389 ()
can be carried over in a straightforward manner to the case of as well^{4}^{5}

As mentioned briefly in the introduction and in the last paragraph, there are several types of real radiative corrections in the threshold limit at three-loop order involving the emission of one or more soft partons. To obtain
complete results for the , one must compute soft single-emission, two-loop 1309.4391 (); 1309.4393 () and one-loop squared real-virtual corrections PHRVA.D60.116001 (); hep-ph/0007142 (); JHEPA.1312.088 (); 1312.1296 (),
soft triple-emission real corrections 1302.4379 (), and soft double-emission, one-loop real-virtual corrections, which are not yet publicly available^{6}

To generate the full set of cut eikonal Feynman diagrams, we use QGRAF Nogueira:1991ex (). Representative cut diagrams are shown in Fig. 1. In order to test the correctness of our squared matrix elements while computing them, we found it convenient to use a general gauge for our internal gluons and lightcone gauge for our final state gluons. As we explain below, the fact that our problem supplies two natural light-like reference vectors, and , allows for a very stringent consistency check on our construction of the integrand. Once the diagrams are generated, they are processed by an in-house FORM math-ph/0010025 () code which dresses the cut eikonal diagrams with Feynman rules and performs all relevant Dirac and color algebra. After that, we use an in-house Maple code to express all cut eikonal diagrams as a linear combination of Feynman integrals belonging to one of 28 integral families introduced for the eventual purpose of performing an integration by parts reduction PHLTA.B100.65 (); NUPHA.B192.159 () on the integrand using Laporta’s algorithm hep-ph/0102033 (). Mapping Feynman integrals to specific integral families is not entirely straightforward in this case and requires the application of numerous partial fraction identities to the raw integrand produced by our FORM code. This complication is due in part to the constraint on the kinematics imposed by the delta function in the operator definition of the soft functions (see Eq. (10)). We perform the integration by parts reduction using the development version of Reduze 2 1201.4330 () which, among other new features, supports the reduction of phase space integrals. The idea behind this was worked out some time ago and is usually called the reverse unitarity method hep-ph/0207004 (); hep-ph/0306192 (); hep-ph/0312266 (). The key insight is that, for the purpose of integral reduction, one can use the relation

(19) |

to replace delta function constraints with propagator denominators. After integration by parts reduction, the result is a rather simple linear combination of nine master integrals.

To present our result in a compact fashion, let us first introduce some additional notation. Besides the Casimir invariants and , already defined above, and the number of massless flavors, , our result for will depend on an additional Casimir invariant which itself depends on the index :

(20) |

Finally, if we make the definitions

(21) |

the nine master integrals mentioned above are given by

(22) | |||||

(23) | |||||

(24) | |||||

(25) | |||||

(27) | |||||

(28) | |||||

(29) | |||||

(30) | |||||

While integrals through can be simply expressed in terms of gamma functions and generalized hypergeometric functions, is more complicated. In particular, it involves the Kampé de Fériet function (see Ref. Exton (), Eq. (1.3.2.1) for our Kampé de Fériet function conventions)

(31) |

A few words about our derivation of the master integrals are in order. Although, for the most part, we wish to defer the discussion of the novel method used to do the integrals
to a more technical future publication, let us briefly describe some consistency checks on Eqs. (22)-(30) that we found useful. For all integrals, a good first step is to integrate out the virtual momentum
since, in all cases, one obtains simple expressions built out of gamma functions and Gauss hypergeometric functions.
The numerical sector decomposition hep-ph/0004013 (); 0709.4092 () code FIESTA 3 1312.3186 () was successfully used to check our virtual integrations for sign errors.^{7}^{8}^{9}

The master integral is somewhat more subtle than -. In an attempt to expand the second line of Eq. (30) in , we encountered the single-parameter integral

(32) |

Historically, such integrals were considered pathological due to the fact that, naïvely, they do not appear to be sector-decomposable. The problem is that the integrand of (32) has small asymptotics of the too-singular form . Integrals such as (32) have appeared in other QCD computations. For example, in Ref. hep-ph/0404293 (), a dedicated unitarity-based sewing method was used to avoid classes of integrals with similar power-law singularities. In fact, we found that we were able to treat integrals like (32) in a very direct way. First, observe that the above integral converges for all complex such that . In spite of the fact that the integral does not converge in a neighborhood of , one can expand (32) in a Laurent series about anyway by carefully applying the principle of analytical continuation. To understand this, begin by considering a fictitious, alternative representation of our result, Eq. (33), written in terms of integrals free of power-law divergences which converge in a neighborhood of . This hypothetical representation should always exist if all of the divergences in the calculation are regulated by . Using the fact that all Feynman integrals in dimensional regularization are analytic functions of , one can perform an analytical continuation to a value of which happens to simultaneously lie in the region of convergence of (32) and through . Now, by performing an integration by parts reduction on the fictitious representation, one can recover the precise form of Eq. (33), written in terms of our preferred basis of master integrals. Due to the fact that, by assumption, the point in the complex plane to which the original analytical continuation was performed lies within the domain of convergence of all of our master integrals, one can derive the closed formulas for - given above (the second lines of Eqs. (22)-(30)) in this region of the complex plane. Finally, one can analytically continue back to the point and derive Laurent expansions for all of the masters integrals using standard techniques. Carrying out this program explicitly was necessary to Laurent expand about and our analysis will be presented in a future publication which features an in-depth discussion of the master integrals.

Our choice of integral basis also requires some explanation. As is clear from the Taylor series expansions of our master integrals, the definitions we have made are such that all of the are pure functions in the sense of reference 1304.1806 (). As we shall see, the all-orders-in- result for assumes a particularly simple form when expressed in terms of the integral basis defined above. After integration by parts reduction, it can be written as

(33) | |||||