Higgs-differential cross section at NNLO in dimensional regularisation

TL;DR

This paper provides an analytic NNLO QCD calculation of the Higgs production cross section differential in Higgs momentum, including higher order dimensional regularisation terms, with numerical predictions for decay product distributions.

Contribution

It introduces a detailed analytic computation of the Higgs differential cross section at NNLO, incorporating higher order dimensional regularisation terms necessary for future N$^3$LO accuracy.

Findings

Numerical predictions for Higgs decay to gamma gamma distributions.

Analytic expressions for differential Higgs cross sections at NNLO.

Methodological advancements in handling dimensional regularisation at high orders.

Abstract

We present an analytic computation of the Higgs production cross section in the gluon fusion channel, which is differential in the components of the Higgs momentum and inclusive in the associated partonic radiation through NNLO in perturbative QCD. Our computation includes the necessary higher order terms in the dimensional regulator beyond the finite part that are required for renormalisation and collinear factorisation at N$^3$LO. We outline in detail the computational methods which we employ. We present numerical predictions for realistic final state observables, specifically distributions for the decay products of the Higgs boson in the $\gamma\gamma$ decay channel.

Figures

Equations

$$\mathrm{Proton}(P_1)+\mathrm{Proton}(P_2)\to H(p_h)+X,$$

$$P_1=\frac{\sqrt{S}}{2}\left(1,0,0,1\right),P_2=\frac{\sqrt{S}}{2}\left(1,0,0,-1\right),$$

$$i(p_1)+j(p_2)\to H(p_h)+X,$$

$$p_1=x_1{P}_1,p_2=x_2{P}_2,$$

$${\sigma }_{PP\to H+X}$$

$$p_h\equiv \left(E,p_x,p_y,p_z\right)=\left(\sqrt{p_T^2+m_h^2}\cosh Y,p_T\cos \varphi ,p_T\sin \varphi ,\sqrt{p_T^2+m_h^2}\sinh Y\right),$$

$$Y=\frac{1}{2}\log \left(\frac{E+p_z}{E-p_z}\right),p_T=\sqrt{E^2-p_z^2-m_h^2}.$$

$$\begin{gathered} {\sigma }_{PP\to H+X}\left[\mathcal{O}\right]=\sum _{i,j}{\int }_{-\infty }^{+\infty }dY{\int }_0^{\infty }d{p}_T^2{\int }_0^{2\pi }\frac{d\varphi }{2\pi }{\int }_0^{1}d{x}_{1}d{x}_2{f}_i(x_1)f_j(x_2) \\ \times \frac{d^2{\hat{\sigma }}_{ij}}{dYd{p}_T^2}(S,x_1,x_2,m_h^2,Y,p_T^2){\mathcal{J}}_{\mathcal{O}}(Y,p_T^2,\varphi ,m_h^2), \end{gathered}$$

$$s=x_1{x}_{2}S=(p_1+p_2{)}^2$$

$$\tau =\frac{m_h^2}{S},z=\frac{m_h^2}{s}=\frac{\tau }{x_1{x}_2}.$$

$$\lambda \equiv \frac{s-2p_1\cdot p_h}{s-m_h^2},x\equiv \frac{s(p_1+p_2-p_h{)}^2}{(s-2p_1\cdot p_h)(s-2p_2\cdot p_h)},$$

$$Y=\frac{1}{2}\log \left[\frac{x_1}{x_2}\frac{1-\frac{\bar{z}\bar{\lambda }}{1-\bar{z}\lambda x}}{1-\bar{z}\lambda }\right],p_T^2=s\frac{{\bar{z}}^2\lambda \bar{\lambda }\bar{x}}{1-\bar{z}x\lambda }.$$

$$\begin{gathered} {\sigma }_{PP\to H+X}\left[\mathcal{O}\right]=\tau \sum _{i,j}{\int }_{\tau }^1\frac{dz}{z}{\int }_{\frac{\tau }{z}}^1\frac{d{x}_1}{x_1}{\int }_0^{1}dx{\int }_0^{1}d\lambda {\int }_0^{2\pi }\frac{d\varphi }{2\pi } \\ \times f_i(x_1)f_j\left(\frac{\tau }{x_{1}z}\right)\frac{1}{z}\frac{d^2{\hat{\sigma }}_{ij}}{dxd\lambda }(z,x,\lambda ,m_h^2){\mathcal{J}}_{\mathcal{O}}(x_1,z,x,\lambda ,\varphi ,m_h^2). \end{gathered}$$

$${\hat{\sigma }}_{ij}\sim \sum _m\int d{\Phi }_{H+m}{\mathcal{M}}_{ij\to H+X},$$

$$d\Phi_{H+m}=\frac{d^{d}p_{h}}{(2\pi)^{d}}(2\pi)\delta_{+}(p_{h}^{2}-m_{h}^{2})\left[\prod_{i=1}^{m}\frac{d^{d}k_{i}}{(2\pi)^{d}}(2\pi)\delta_{+}(k_{i}^{2})\right](2\pi)^{d}\delta^{d}\bigg{(}p_{1}+p_{2}-p_{h}-\sum_{i=1}^{m}k_{i}\bigg{)},$$

$${\delta }_{+}(p^2-M^2)=\theta (p^0-M)\delta (p^2-M^2).$$

$$1=\int\frac{d^{d}k}{(2\pi)^{d}}(2\pi)^{d}\delta^{d}\bigg{(}k-\sum_{i=1}^{m}k_{i}\bigg{)}\int_{0}^{\infty}\frac{d\mu^{2}}{2\pi}(2\pi)\delta_{+}(k^{2}-\mu^{2}).$$

$$d{\Phi }_{H+m}={\int }_0^{\infty }\frac{d{\mu }^2}{2\pi }d{\Phi }_{HX}({\mu }^2)d{\Phi }_m({\mu }^2),$$

$$d{\Phi }_{HX}$$

$$d{\Phi }_m$$

$${\delta }_{+}\left(q^2\right)\to {\left[\frac{1}{q^2}\right]}_{\mathrm{c}}=\frac{1}{2\pi i}\text{Disc}\left(\frac{1}{q^2}\right)=\frac{1}{2\pi i}\left(\frac{1}{q^2+i0}-\frac{1}{q^2-i0}\right).$$

$$\frac{\partial }{\partial q_{\mu }}{\left({\left[\frac{1}{q^2}\right]}_{\mathrm{c}}\right)}^{\nu }=-\nu {\left({\left[\frac{1}{q^2}\right]}_{\mathrm{c}}\right)}^{\nu +1}2q^{\mu },$$

$${\left({\left[\frac{1}{q^2}\right]}_{\mathrm{c}}\right)}^{-n}=0\text{for}n\ge 0.$$

$${\hat{\sigma }}_{ij}(S,x_1,x_2,m_h^2)\sim \sum _m{\int }_0^{\infty }\frac{d{\mu }^2}{2\pi }d{\Phi }_{HX}\times \underset{\text{Reverse Unitarity}}{\underbrace{\left[\int d{\Phi }_m{\mathcal{M}}_{ij\to H+X}\right]}}.$$

$$2k\cdot p_1=\bar{z}\lambda s,2k\cdot p_2=\bar{z}\tilde{\lambda }s.$$

$$k^2={\mu }^2=s{\bar{z}}^2\lambda \tilde{\lambda }x,\tilde{\lambda }=\frac{1-\lambda }{1-\bar{z}\lambda x}.$$

$$d{\Phi }_{HX}=\frac{s^{\frac{d}{2}-2}}{4(2\pi {)}^{d-1}}d{\Omega }_{d-2}d\lambda {\bar{z}}^{d-3}(\lambda \tilde{\lambda }{)}^{\frac{d}{2}-2}(1-x{)}^{\frac{d}{2}-2}\theta (s)\theta (\bar{z})\theta (\lambda )\theta (\bar{\lambda })\theta (\bar{x}),$$

$$\int d{\Omega }_d=\frac{2{\pi }^{\frac{d}{2}}}{\Gamma \left(\frac{d}{2}\right)}.$$

$$d{\Phi }_{H+m}=dx\frac{{\bar{z}}^2\lambda \bar{\lambda }}{(1-\bar{z}\lambda x{)}^2}d{\Phi }_{HX}(x,\lambda )d{\Phi }_m(x,\lambda ).$$

$$d{\Phi }_{H+0}=\frac{2\pi }{s}\delta (\bar{z}).$$

$$d{\Phi }_{H+1}=\frac{s^{\frac{d}{2}-2}}{4(2\pi {)}^{d-2}}{\bar{z}}^{d-3}d{\Omega }_{d-2}d\lambda (\lambda \bar{\lambda }{)}^{\frac{d}{2}-2}\theta (s)\theta (\bar{z})\theta (\lambda )\theta (\bar{\lambda }).$$

Full text

NNLO, QCD, Higgs boson, differential distributions

††preprint: CERN-TH-2017-090, SLAC-PUB-16962

1 Introduction

The discovery of the Higgs boson in 2012 at the Large Hadron Collider (LHC) by the ATLAS [1] and CMS [2] experiments founded a new era of precision Higgs physics. The impressive statistical accuracy of the experimental measurements has lead to tight constraints on the couplings of Higgs boson interactions. The Run 2 of the LHC promises a plethora of further measurements which will allow to probe the nature of electroweak symmetry breaking with an unprecedented precision.

In order to truly exploit the potential of the LHC, the excellent experimental results must be compared to equally precise theoretical predictions. The demand for high precision theory has been met in recent years with a flurry of calculations at next-to-leading (NLO) and next-to-next-to-leading-order (NNLO) in perturbative QCD. Higher order corrections are especially important for Higgs phenomenology, in part due to the slow convergence of the perturbative expansion in the strong coupling constant for the dominant mode of Higgs hadroproduction via gluon fusion.

The gluon-fusion total cross section has been computed recently through next-to-next-to-next-to-leading order (N3LO) in perturbative QCD [3, 4] in the limit of an infinite top-quark mass. At this level of precision, effects that go beyond the leading approximation of the heavy-top quark effective theory or a treatment in pure QCD, such as quark mass effects or contributions from weak boson loops, are also important. The most accurate available predictions for these effects have been combined into a single theoretical prediction for the inclusive Higgs cross section in ref. [3].

The Higgs cross section is measured in signal-rich regions of phase space which are shaped by carefully designed experimental event selections. Further restrictions on the phase space of the Higgs boson, its decay products and the associated radiation are required by the detector geometry and response. Within their defined acceptance, the experiments have superb capabilities to measure a multitude of kinematic distributions for the Higgs boson and its decay products in the years to come. In addition to the total cross section, it is therefore imperative to have precise theoretical predictions for differential cross sections.

Recently, the $pp\to H+1{\,\,\rm jet}$ fully differential cross section has been computed at NNLO [5, 6, 7, 8, 9]. Combined with the N3LO inclusive cross section, this has allowed to compute the N3LO Higgs cross section with a jet-veto [10]. This is a first example of a differential cross section in gluon-fusion at this perturbative order. If a fully differential parton-level Monte-Carlo at N3LO is achieved in the future, it will become possible to assess the efficiency of the majority of event selection criteria at the same level of accuracy in perturbative QCD as the jet-veto efficiency. To achieve this goal, one could attempt to generalise any of the available methods at NNLO (sector-decomposition [11, 12, 13, 14], slicing [15, 16, 17, 18], subtraction [19, 20, 21, 22], reweighting [23] and other methods [24]) to the next perturbative order. A generalisation of any of these methods would be a formidable task. An intermediate goal could therefore be to compute first some specific differential distributions of particular importance, which are the main ingredients for a fully differential N3LO parton-level Monte-Carlo within a slicing method.

The aim of this article is to compute the Higgs cross section fully differential in all components of the Higgs momentum and its decay products, treating the associated QCD radiation inclusively (integrating over the unrestricted phase space of all partons in the final state) at NNLO. In our computation, we maintain the full dependence of the “Higgs-differential cross section” on the dimensional regulator. This is necessary in order to enable the construction of the counter terms for ultraviolet and initial-state collinear divergences at N3LO.

Our method employs reverse unitarity [25, 26, 27, 28, 29], to map the phase space integrations over partonic radiation onto their loop integral duals. We then reduce the integrals appearing in the partonic cross sections to master integrals using integration-by-part identities [30, 31] and the Laporta algorithm [32]. The master integrals corresponding to the partonic cross sections with two partons in the final state are novel and we evaluate them using two methods: the method of differential equations [33, 34, 35] and a direct integration of the angular phase space variables [36]. We are able to express all master integrals in terms of hypergeometric functions which are valid at all orders in the dimensional regulator $\epsilon=(4-d)/2$. These hypergeometric functions can be expanded to practically any order around the limit $\epsilon=0$ in terms of polylogarithms. As a result we obtain an analytic formula for all required partonic cross sections through NNLO, including the higher order terms in the $\epsilon$ expansion which are needed as input for a future N3LO calculation.

We implement our results in a computer code and use it to obtain predictions at NNLO for various differential distributions which are of interest for LHC experiments. In addition to computing the transverse momentum and rapidity distribution for the Higgs boson through NNLO, we analyse various differential properties of the production and subsequent decay of the Higgs to two photons.

This article is structured as follows. In section 2 we introduce in detail our definition of “Higgs-differential” cross sections. In section 3 we outline how we separate the integration of Higgs boson and QCD radiation degrees of freedom based on phase space factorisation and the reverse unitarity methodology. In section 4 we perform explicitly the computation of the partonic cross sections for Higgs boson production through NNLO in gluon fusion. Next, we study several key Higgs boson LHC observables in section 5 that have been obtained using a numeric code built upon our analytic results. Finally, we give our conclusions in section 6.

Note that while we only study differential Higgs boson observables in this article, our “Higgs-differential” method is not specific to Higgs boson processes. In fact, it relies solely on the fact that the Standard-Model Higgs boson is a singlet of QCD. As such, it can also be applied to compute differential distributions for other colourless final states, such as Drell-Yan or diboson production.

2 Setup for Differential Cross Sections

We consider the production of a Higgs boson in proton-proton collisions:

[TABLE]

where in parentheses we denote the momenta carried by the external particles in the process. The four-momenta of the protons in the hadronic center of mass frame are given by

[TABLE]

where $\sqrt{S}$ is the collider center of mass energy. The inclusive hadronic cross section is related to the one for the partonic processes

[TABLE]

with momenta

[TABLE]

via the factorisation formula

[TABLE]

In the above, $p_{h}^{2}=m_{h}^{2}$ is the invariant mass of the Higgs boson. The sum over $i$ and $j$ runs over all possible initial state partons. $f_{i}(x)$ are the parton distribution functions and $\hat{\sigma}_{ij}(S,x_{1},x_{2},m_{h}^{2})$ are the partonic cross sections.

In this work we are interested in deriving cross sections for observables $\mathcal{O}$ that are sensitive to the momentum of the Higgs boson $p_{h}$ and do not depend on the details of the additional radiation $X$. An observable of this type translates into a measurement function $\mathcal{J}_{\mathcal{O}}(p_{h})$ which multiplies the differential cross section giving a weight to every phase space point. A typical example for $\mathcal{O}$ is a set of kinematic cuts on the Higgs; in this case $\mathcal{J}_{\mathcal{O}}(p_{h})$ is equal to one if $p_{h}$ passes these cuts and zero otherwise, and it can be written as a product of Heaviside $\theta$-functions. More complicated, experimentally relevant observables may also be considered as long as they only depend on $p_{h}$. For example, $\mathcal{J}_{\mathcal{O}}$ may be used to weight the production cross section by appropriate Higgs boson decay matrix elements and phase spaces. We will refer to cross sections for such observables $\mathcal{O}$ that depend on the kinematic details of the production process only through the Higgs four-momentum as Higgs-differential cross sections.

One possible way to parameterise the Higgs momentum is

[TABLE]

where

[TABLE]

Here $Y$ is the rapidity of the Higgs boson, $p_{T}$ is its momentum in the plane which is transverse to the beam axis, and $\phi$ is the azimuthal angle. Due to the symmetry of scattering experiments in the azimuthal plane, partonic cross sections are always independent of $\phi$. Without loss of generality, we may therefore write Higgs-differential cross sections as

[TABLE]

where $d^{2}\hat{\sigma}_{ij}/dYdp_{T}^{2}$ is the partonic double-differential $p_{T}$ and rapidity distribution.

Results for the partonic cross sections are naturally expressed in terms of the kinematics of the partonic processes. We thus define the partonic center of mass energy

[TABLE]

and the ratios

[TABLE]

Furthermore, we introduce the two Lorentz-invariant quantities $x$ and $\lambda$

[TABLE]

and the shorthand notation $\bar{z}\equiv 1-z$, $\bar{\lambda}\equiv 1-\lambda$, $\bar{x}\equiv 1-x$. We can express $p_{T}$ and $Y$ in terms of these two new variables and the Bjorken-fractions:

[TABLE]

The Higgs-differential cross sections, when cast in terms of these variables, take the form:

[TABLE]

Note that we absorbed the Jacobian of the variable transformation from $(p_{T},Y)$ to $(x,\lambda)$ into the definition of the partonic Higgs-differential cross section. Having obtained the above formula, we are now ready to compute the cross section for any hadron collider observable which is differential in the Higgs boson momentum from a partonic Higgs-differential distribution in $x$, $\lambda$.

3 The Higgs-Differential Phase Space

In the previous section we parametrised the degrees of freedom corresponding to the momentum of the Higgs boson with variables $z$, $x$, $\lambda$ and $\phi$ defined ad hoc. In the following it will become clear why this particular parametrization of the Higgs boson degrees of freedom is beneficial to our calculation. This section is divided in three parts. In the first part we discuss how the phase space of the final state can be separated into a part that depends solely on the Higgs boson kinematics, and one that captures all dependence on the additional radiation. We will determine a set of phase space integrals in such a way that we can carry out the integration over the radiation analytically while remaining differential in the Higgs boson degrees of freedom. In the second part we discuss the computation of these integrals using reverse unitarity. Finally, we will explicitly construct the parametrisation of $p_{h}$ in terms of $z$, $x$, $\lambda$, and discuss the features of the remaining phase space.

3.1 Separation of the Phase Space

The partonic cross section at a given order in QCD perturbation theory is comprised of phase space integrals over a sum of real and virtual contributions with different multiplicities of partons in the final state. Schematically, we may write,

[TABLE]

where $X$ stands for a set of $m$ final state partons. We are going to examine the integration of a matrix element for a fixed multiplicity of partons over the phase space measure $d\Phi_{H+m}$. Due to soft and collinear singularities, such integrations are divergent in four dimensions and we regulate them in dimensional regularisation. We will elaborate on the detailed structure of the partonic matrix elements through NNLO in section 4.

The integration measure for the phase space describing the production of a Higgs boson and $m$ additional partons with outgoing momenta $k_{1},\ldots,k_{m}$ is given by

[TABLE]

where

[TABLE]

In order to compute Higgs-differential cross sections it is natural to separate the integration over the momentum of the Higgs from the integral over momenta of the final state state partons. This can be achieved by inserting unity into eq. (12) using

[TABLE]

This identity allows us to factor the $H$ plus $m$ parton phase space into the phase space of two massive particles and the phase space of $m$ massless partons:

[TABLE]

where

[TABLE]

Note that these formulae are valid for an arbitrary number of final state partons $m\geq 0$.

A key ingredient of our method is that, because of the class of observables we restricted ourselves to, the integrals over the extra QCD radiation are not constrained by the definition of the measurement. This allows us to address the integration over $d\Phi_{HX}(\mu^{2})$ and $d\Phi_{m}(\mu^{2})$ separately.

3.2 Final State Radiation and Reverse Unitarity

In the present section we address how the analytic integration over additional final state partons can be performed. In particular, we advocate that the framework of reverse unitarity [25, 26, 27, 28, 29], that has already been successfully used in the computation of the inclusive cross section [37, 38, 39, 40, 41, 42, 43, 44], provides a particularly efficient solution to this task. This framework exploits the duality between inclusive phase space integrals and loop integrals to treat them in a uniform way. Specifically, using Cutkosky’s rule [45], it is possible to express the on-shell constraints appearing in phase space integrals through cut propagators

[TABLE]

Cut propagators can be differentiated in a similar way to ordinary propagators with respect to their momenta,

[TABLE]

leading to identical integration-by-parts (IBP) identities [30, 31] for inclusive phase space integrals as for their dual loop integrals. The fact that the cut propagator represents a delta function is reflected by the simplifying constraint that any integral containing a cut propagator raised to a negative power vanishes:

[TABLE]

IBP identities serve to relate different inclusive phase space and loop integrals and express them in terms of a finite set of so-called master integrals. Once a partonic cross section is represented in terms of a linear combination of master integrals, only those integrals have to be computed by other means.

In refs. [27, 28, 29] a modification of the reverse unitarity framework was developed in order to allow for the computation of the partonic Drell-Yan and Higgs boson production cross section differential in the rapidity of the electro-weak final state boson. Here we want to take this procedure one step further in order to maintain all differential information about the Higgs boson. This can be simply achieved by refraining from including the integration over the Higgs boson momentum into the reverse unitarity framework. Exploiting again the schematic representation of our partonic cross sections eq. (11) we may write

[TABLE]

To put this into other words: the factorisation of the Higgs boson phase space and the QCD radiation phase space achieved by eq. (15) allows us to carry out the integration over the parton momenta separately from the integration over the Higgs boson momentum. This enables us to use different methods for the two phase spaces without sacrificing any information about the differential properties of the Higgs boson. In particular, we can compute the QCD radiation phase space inclusively in dimensional regularisation. This lends itself naturally to the method of reverse unitarity. It allows us to apply IBP identities to the partonic matrix elements depicted in the square brackets in eq. (21) and express them in terms of differential master integrals. These master integrals can subsequently be computed by different means such as direct integration or differential equations, as we will demonstrate below.

The resulting integrated matrix elements are given by a Laurent series in the dimensional regulator $\epsilon$. The explicit poles correspond to the regulated soft and collinear singularities of the final state radiation. In contrast to conventional methods for the computation of differential cross sections, we circumvent the problem of subtracting infrared singularities among final state partons by performing these phase space integrals analytically in dimensional regularisation. We now turn to discuss the specific parametrisation of the remaining degrees of freedom.

3.3 Parametrising the Higgs Boson Degrees of Freedom

Once the analytic integration over the partonic final state momenta has been performed, the remaining degrees of freedom are parametrised by the Higgs boson four momentum $p_{h}$, or equivalently by the collective momentum of the integrated final state partons $k$. The inclusive integration over $p_{h}$ and $k$ is plagued by infrared and collinear singularities when $k^{2}=0$ or transverse momenta vanish, $p_{h\;T}^{2}=k_{T}^{2}=0$. These kinematic limits correspond to the configuration of a single resolved real emission and to Born level kinematics, respectively. We will now parametrise the measure of eq. (16) in such a way that leftover divergences can be regulated in a straightforward manner. We shall see how the parameters $x$ and $\lambda$ introduced in section 2 are actually a simple solution of this problem.

We parametrise the scalar products of the initial state parton momenta with $k$ by

[TABLE]

We then find it convenient to exchange the variable $\tilde{\lambda}$ for $x$ where

[TABLE]

The parameter $\lambda$ then effectively measures the collective direction of the QCD radiation, while $x$ corresponds to its invariant mass. Note that these new variables range in the interval $[0,1]$. The potential leftover singularities in the $(\lambda,x)$ space are factorised and lie on the faces of the $[0,1]\times[0,1]$ square without any extra divergences as the corners are approached.

With these definitions we find after some algebra (see for example [37]) the following parametrisation for the two particle massive phase space measure,

[TABLE]

where the integral over the $(d-2)$-dimensional solid angle can be easily carried out using

[TABLE]

The total integration measure of eq. (12) then becomes

[TABLE]

It should be noted here that our parametrisation of $\mu^{2}$ in terms of $x$ and $\lambda$ spoils the traditional factorisation of the phase space into a convolution over two independent phase spaces. Instead, our phase space is now factorised into an iterative form. First we perform the phase space integral over the QCD radiation phase space $d\Phi_{m}$, obtaining a result as a function of $x$ and $\lambda$. Afterwards we can perform the integral over the Higgs phase space $d\Phi_{HX}$ in terms of $\lambda$ and finally over the additional parameter $x$.

The expression given in eq. (26) is valid for 2 or more final state partons. The cases of zero or one final state partons need to be addressed separately as they represent limiting cases of the general parametrisation above. It is however straightforward to parametrise these cases explicitly. For zero partons we find

[TABLE]

and for one final state parton we define

[TABLE]

Now that we have derived an explicit parametrisation of the phase space for the Higgs boson it is useful to make contact again between our integration variables and the actual properties of the Higgs boson. Let $y$ be the rapidity of the Higgs boson in the partonic centre of mass frame, which is related to $Y$ by

[TABLE]

The partonic rapidity and transverse momentum of the Higgs boson are related to $\lambda$ and $x$ by

[TABLE]

It is easy to see that the parametrisation obtained here corresponds to the choice of variables introduced in eq. (9). The differential partonic cross section required to compute the Higgs-differential hadronic cross section in eq. (10) is simply obtained by performing all integrations over the final state degrees of freedom, except for the integrations over $x$ and $\lambda$.

We would like to remark that, although the phase space parametrisation in this section was derived specifically for the case of Higgs boson production for definiteness, it actually holds for any single particle colourless final state. Moreover, extending the $(x,\lambda)$ parametrisation to the case of the production of a colourless final state system composed of more particles, the separation of the phase space and reverse unitarity still work as discussed without the need of any further change.

4 Computation of the Higgs-Differential Cross Section through NNLO

In the previous section we established a framework for the computation of Higgs-differential cross sections. In this section we explicitly compute the differential cross section for the production of a Higgs boson via the gluon fusion mechanism through NNLO in the infinite top mass limit.

4.1 Partonic Cross Sections for Gluon Fusion

We compute the Higgs boson cross section in an approximation to the full Standard Model where the top quark mass is considered to be infinite and internal top quark loops can be integrated out. This leads to an effective field theory where the Higgs is directly coupled to gluons [46, 47, 48, 49] through an effective operator

[TABLE]

Here, $\mathcal{L}_{QCD}$ is the QCD Lagrangian, $h$ is the Higgs boson field, $G_{\mu\nu}$ the gluon field strength and $C^{0}$ the Wilson coefficient [50, 51, 52] that arises from matching the effective theory to the full Standard Model. The QCD Lagrangian contains $n_{f}$ massless quark fields with $n_{c}$ colours. In the following we will compute perturbative corrections in the strong coupling constant $\alpha_{S}$ through NNLO using this effective theory:

[TABLE]

where we normalised the coefficient functions $\eta_{ij}$ to the Born cross section

[TABLE]

Perturbative QCD calculations are plagued by the presence of ultraviolet and infrared divergencies that appear during the computational steps and cancel for well defined observables. We regulate these divergencies by applying the framework of dimensional regularisation in the $\overline{\mathrm{MS}}$ scheme, continuing the number of space time dimensions to $d=4-2\epsilon$. Ultraviolet finite observables are obtained by renormalising the parameters of the theory. The required redefinition of the strong coupling constant and of the Wilson coefficient are given by

[TABLE]

where $\gamma_{E}$ is the Euler-Mascheroni constant. For convenience we provide the well known factors $Z_{\alpha}(\mu^{2})$ and $Z_{C}(\mu^{2})$ in appendix A. The Wilson coefficient can be found in appendix B.

At any order in perturbation theory we distinguish 6 different initial state configurations of partons:

[TABLE]

Here $g$, $q$, $\bar{q}$ and $q^{\prime}$ represent a gluon, a quark, an anti-quark and a quark with a different flavour respectively. All other combinations of explicit (anti-)quark flavours and gluons can be obtained from the ones above. $X$ represents a specific partonic final state. The partonic coefficient function for the individual channels at order $n$ for $m$ final state partons is given by

[TABLE]

where $\sum\mathcal{M}^{(n)}_{ij\rightarrow H+X}$ is the coefficient of $\alpha_{S}^{n}$ in the coupling constant expansion of the modulus squared of all amplitudes for partons $i$ and $j$ producing the final state $H+X$, summed over polarisations. The integration measure $d\Phi_{H+m}$ was defined in eq. (12), while the initial state dependent prefactors $N_{ij}$ are given by

[TABLE]

For the differential cross section through NNLO we require amplitudes with up to two additional partons in the final state as well as contributions with up to two loops. At any given order $\alpha_{S}^{n}$, the sum of the number of loops and the number of final state partons is equal to $n$. All purely virtual partonic cross sections are identical to those required for the computation of the inclusive Higgs boson cross section that were obtain at two loops in refs. [53, 54]. At NLO we require tree level matrix elements with one additional parton in the final state (R). At NNLO we require real-virtual (RV) matrix elements with one loop and one additional parton in the final state and double-real (RR) matrix elements with two additional final state partons.

In order to compute the partonic coefficient functions we generate the necessary Feynman diagrams using QGRAF [55]. We then perform spinor, tensor and colour algebra in a private C++ code based on GiNaC [56], a code based on FORM [57] and a private Mathematica package. The resulting expressions represent the phase space and loop integrands for the partonic cross sections. Using the reverse unitarity framework discussed in section 3.2, we treat loop and phase space integration on equal footing and use IBP [30, 31] identities to express the partonic cross sections in terms of master integrals. We then compute the master integrals explicitly, as discussed in section 4.2. Inserting them into the partonic matrix elements, we obtain the partonic cross sections as a Laurent expansion in the dimensional regulator.

Each of the contributions corresponding to different parton and loop multiplicities is separately infrared divergent as made manifest by explicit poles in the dimensional regulator. After summing up the contributions, some of the divergencies cancel by virtue of the KLN theorem. The remaining divergencies are absorbed by a suitable redefinition of the parton distribution functions,

[TABLE]

where the convolution indicated with $\circ$ is defined by

[TABLE]

The infrared counterterms $\Gamma_{ij}$ consist of convolutions of splitting functions $P_{ij}^{(n)}$ and can be derived from the DGLAP equation; for reference we provide its explicit form in appendix A. The remaining physical parton distribution functions $f^{R}_{i}$ are process independent and are extracted from measurements.

In inclusive calculations, it is often useful to employ the commutativity and associativity of the convolutions to rewrite them such that the infrared counter term $\Gamma$ is convoluted with the partonic cross section before convoluting the result with the bare parton distributions. The first convolution between the counter term and the partonic cross section can thus be performed analytically. The required splitting functions were obtained in refs. [58, 59] and the convolutions were performed for example in refs. [60, 61]. For the purposes of this work, it is impractical to perform these convolutions analytically, as they depend on the additional parameters $\lambda$ and $x$ due to momentum conservation. We perform the convolution of the physical parton distribution functions and the infrared counter term therefore numerically. The infrared counter terms $\Gamma$ themselves contain poles in $\epsilon$. In order to obtain the correct contributions to the finite part, it is therefore necessary to compute the lower order partonic cross sections to higher than finite power in the dimensional regulator.

To obtain a finite fixed order differential cross section we expand the product of bare parton distribution functions $f_{i}$, Wilson coefficient $C^{0}$ and partonic coefficient functions $\eta_{ij}$ and truncate the product at fixed order. After the combination of all these contributions, we complete the computation of the Higgs-differential production cross section through NNLO eq. (10). Extending the computation to one order higher in the strong coupling constant (N3LO) will require as an input the partonic cross sections computed to one order higher in the dimensional regulator. One of the main results of this work are the analytic expressions for all the necessary partonic cross sections through NNLO up to and including order $\epsilon^{1}$, as required for an N3LO computation. This result is available in electronic form in an ancillary file submitted together with this publication.

4.2 Evaluation of Master Integrals

In this section we elaborate on the computation of the master integrals that serve as building blocks for our partonic cross section. To validate our results, we follow two different strategies. The first is based on the method of differential equations [33, 34, 35] and is well established in the field of high order computations. As both strategies lead to identical results we will not discuss this method in more detail.

The second approach is based on direct evaluation of phase space and loop integrals. Here, we discuss only master integrals involving at least one final state parton as purely virtual corrections are well known [53, 54]. At NLO the only master integral is simply given by the phase space measure eq. (28).

4.2.1 Double-real master integrals

At NNLO there are contributions with two final state partons (double-real) as well as one-loop corrections to the emission of a single final state parton (real-virtual). Let us consider the double-real radiation (RR) first. We find eight master integrals $M_{i}^{RR}$, corresponding to the six diagrams in Figure 1. Propagators crossing the dashed line in the diagrams represent cut propagators, while the double line marks the massive Higgs boson propagator. The fact that we are interested in Higgs-differential master integrals means that the momentum of the Higgs boson is completely fixed, so that no integration over the degrees of freedom of the Higgs boson occurs. {fmffile}QWmasts

In accordance with the phase space definitions of eqs. (16) and (17), all RR master integrals take then the form

[TABLE]

where again we denoted with $k$ the momentum of the parton system. The generalised exponents $\alpha$ and $\beta$ and the momenta $q$ and $p$ specify the propagators that appear in the respective master integral. In particular, $p$ and $q$ are linear combinations of $p_{1}$, $p_{2}$ and $k$. The integral over $l$ stands for the integral over one of the two final state parton momenta. The result is then a function of all Lorentz invariants of the system.

The integral $I(p,q;\alpha,\beta)$ can be evaluated for generic values of the parameters in terms of Appell’s Hypergeometric function [36]:

[TABLE]

where we have defined the auxiliary quantities,

[TABLE]

The eight master integrals, shown in Figure 1, contributing to the RR partonic cross section, expressed in terms of the function $I(p,q;\alpha,\beta)$, are explicitly given by

[TABLE]

[TABLE]

where $M_{6}^{RR}$ and $M_{7}^{RR}$ are related respectively to $M_{4}^{RR}$ and $M_{3}^{RR}$ under the exchange of $p_{1}$ with $p_{2}$. The rational prefactors and the $C_{RR}$ in the above definitions serve as a normalisation for the master integrals such that their expansion in the dimensional regulator is given by a pure function of uniform transcendental weight. This is true for all but $M_{2}^{RR}$ which still contains a root of a polynomial in our variables. In particular, in the prefactor $C_{RR}$ we absorb additional factors that result from the phase space measure of the integration over the Higgs boson momentum, eq. (24),

[TABLE]

We note that the number of differential master integrals is less than the number of inclusive master integrals found in the double real cross section in ref. [25].

4.2.2 Real-virtual master integrals

In the case of the RV master integrals, the integral over the final state parton phase space becomes trivial and all master integrals are simply given by one loop integrals. We define the two well known functions,

[TABLE]

The box integral was computed for example in ref. [62]. With this we choose the following seven RV differential master integrals:

[TABLE]

The normalisation factor $C_{RV}$ again absorbs part of the integration measure of the Higgs boson momentum and is given by

[TABLE]

Note that we include the usual constants arising in the $\overline{\mathrm{MS}}$ renormalisation procedure in the definition of $C_{RR}$ and $C_{RV}$.

4.2.3 Expansion in the dimensional regulator

For further numerical evaluation, we need to expand the master integrals as a Laurent series in the dimensional regulator $\epsilon$. The Laurent expansions for the RV master integrals are well known, so we comment here only on the expansion of the RR integrals. For masters $M_{1}$, $M_{4}$ and $M_{5}$, Appell’s hypergeometric function can be reduced to Gauss’ hypergeometric function ${}_{2}F_{1}$, therefore one can obtain the $\epsilon$ expansion with well-known tools such as the HypExp package [63, 64]. For the remaining master integrals we can directly expand Appell’s hypergeometric function through weight three, starting from the following integral representation:

[TABLE]

Here we have defined $x^{\prime}=\frac{x}{x-1}$ and similarly for $y$. With this, the first integral in the second equality in eq. (4.2.3) is finite and can be expanded in $\epsilon$ before integrating. The divergence is captured by the second integral, which can be easily integrated for generic $\epsilon$. This yields the following expansion for the $F_{1}$:

[TABLE]

The above expression is real-valued for $x\in[0,1],y<0$.

This allows us to express all coefficients in the Laurent series in terms of real valued classical polylogarithms, enabling a fast and stable numerical evaluation [65].

4.3 Partonic Coefficient Functions

In the previous section we obtained analytic results for all partonic coefficient functions required for the computation of the differential Higgs boson cross section through NNLO. Moreover we computed the coefficient functions to sufficiently high order in the dimensional regulator such that the infrared subtraction term required for an N3LO computation can be constructed.

It is important to note that our coefficient functions contain single poles in the variables $\bar{z}$, $x$ and $\lambda$. These poles represent kinematic configurations where the Higgs boson degrees of freedom revert to lower multiplicity final state kinematics. The prime example is a singularity in the matrix elements when the transverse momentum of the Higgs boson tends to zero, which is the value of the transverse momentum at leading order. Specifically, these singularities are of the form

[TABLE]

where the coefficients $a_{i}$ are small integer numbers.

When integrating over the variables $\bar{z}$, $\lambda$ and $x$ we encounter singularities that lie at the boundaries of the integration range of our variables. This is the case, for example, when we are computing observables like the rapidity distribution of the Higgs boson or the inclusive cross section. In order to be able to compute such observables we will have to regulate the divergences, which we illustrate in the following example.

Consider a function $f(x)=x^{-1+a\epsilon}f_{h}(x)$, for some integer $a$ and with $f_{h}(x)$ holomorphic around $x=0$. We are interested in integrating the function over a test function $\phi(x)$ on the range $[0,1]$. In the case of our Higgs-differential cross section, the test function $\phi(x)$ corresponds to the product of the parton luminosity and the measurement function. We can explicitly subtract the divergence at $x=0$ and integrate by parts to obtain

[TABLE]

We now want to give an expression for the partonic cross section that is finite even if all inclusive integrations are performed. To this end we define in a slight abuse of notation,

[TABLE]

Here the $\delta$ distribution is to be understood as acting only on the test function and not on its coefficient in the square bracket. It is easy to see that $f_{s}(0)$ integrates to zero. We can therefore regulate the integrand $f(x)$ by subtracting $f_{s}(0)$,

[TABLE]

so that every term of its $\epsilon$ expansion can be integrated numerically.

In the case of our Higgs-differential cross sections, we need to regulate potential end-point divergences in the three remaining variables $\bar{z},x$ and $\lambda$, c.f. eq. (10). We define the distributions $\sigma_{s}$ that subtract the limits of $\sigma(\bar{z},x,\lambda)$ and label them by the kinematic limit of the cross section that they reproduce. For example $\sigma(\bar{z},0,\lambda)$ takes care of the limit of the cross section as $x$ goes to zero. After partial fractioning to avoid simultaneous singularities on both endpoints of the integral, we obtain the following decomposition,

[TABLE]

One main result of this article is the analytic computation of the partonic coefficient functions $\eta^{(k)}_{ij}(z,x,\lambda)$ as defined in eq. (32). We created finite versions of this coefficient functions in the spirit discussed above and provide them in Mathematica readable form in an ancillary file together with the arXiv submission of this article. Specifically, we provide all 18 different kinematic configurations of these coefficient functions as in eq. (4.3), each as a sum of products of rational coefficients and master integrals.

5 Numerical results for the Higgs diphoton signal

In this section, we carry out numerically the remaining integrations which are necessary in order to obtain hadronic Higgs-differential cross sections. We present results for the LHC at 13 TeV in order to showcase the types of observables which can be readily computed with our approach.

While our analytical results are original and extend the literature of NNLO Higgs-differential cross sections at $\mathcal{O}(\epsilon)$, our numerical predictions for the finite part can be compared to available computer codes. We have validated our numerical implementation against the predictions of HNNLO [17] and MCFM [66] and we have found good agreement within Monte-Carlo uncertainties.

For our numerical studies we use NNLO MMHT parton distribution functions [67] throughout, as available from LHAPDF [68]. Their default value for the strong coupling constant $\alpha_{S}(m_{Z})=0.118$ and the corresponding three-loop running are also adopted. We set the Higgs boson mass to $m_{h}=125$ GeV and neglect its width. We equate the renormalisation and factorisation scales for simplicity and we choose $\mu=m_{h}/2$ as a central value. As it is common practice, we estimate the effect of missing higher order corrections by varying $\mu$ by a factor of two around its central value.

We start by showing distributions for the inclusive production of a Higgs boson via gluon fusion. In Figure 2(a) we present the unbinned rapidity distribution of the Higgs boson. The bands correspond to the variation of the cross section at LO, NLO and NNLO in our default $\mu$ scale range. In Figure 2(b) we show the $p_{T}$ distribution of the Higgs boson. In order to have a non-vanishing transverse momentum, the Higgs boson needs to recoil against additional radiation and therefore its $p_{T}$ distribution is trivial at LO.

In addition to simple inclusive distributions for a stable Higgs boson, we can also investigate its decays. Specifically, we present differential distributions for the Higgs diphoton signal after the application of typical selection cuts for the photons. The decay of the Higgs boson to two photons allows for precise measurements of numerous properties (see e.g. [69, 70, 71, 72, 73, 74, 75]) due to its exceptionally clean experimental signature and it has played a crucial role in the discovery of the Higgs boson itself [1, 2].

We impose photon selection cuts which follow as closely as possible the diphoton analysis of ATLAS [75]. We require that the pseudorapidities of both photons satisfy $\eta_{\gamma}<2.37$, together with $\eta_{\gamma}\not\in[1.37,1.52]$, which implies that there are no photons in these regions. Furthermore, we require the photon with larger transverse momentum to satisfy $p_{T,\,\gamma_{1}}>0.35\,m_{h}$ and the one with smaller transverse momentum to have $p_{T,\,\gamma_{2}}>0.25\,m_{h}$. Since we have integrated out the associated radiation in the production of the Higgs boson, our analysis ignores photon isolation111This is an important cut for taming the diphoton background but it is a minor one for the Higgs signal process that we study here..

In Figure 3, we present the rapidity distribution of the Higgs boson after the application of the photon selection cuts described above. The bands correspond to the values of the cross section at LO, NLO and NNLO in our default $\mu$ scale range. Due to the restrictions in the coverage of the pseudorapidity of the photons, the Higgs rapidity distribution manifests some non-smooth changes. In the middle panel we show the conventional $K$-factors, i.e. we normalise the rapidity distribution at LO, NLO and NNLO for a general scale $\mu$ to the LO prediction evaluated at a fixed scale $\mu=m_{h}/2$. We observe that the relative size of QCD corrections at NLO and NNLO has a pattern which is similar to the one seen for the inclusive cross section, although the $K$-factors are not entirely uniform across bins due to the effect of the photon cuts. In the lower panel, we normalise the rapidity distribution at LO, NLO and NNLO for a general scale $\mu$ to the NLO prediction at NLO at a fixed scaled $\mu=m_{h}/2$. This shows that the relative $K$-factor from NLO to NNLO is more uniform in rapidity.

In Figure 4(a), we present the distribution of the pseudorapidity difference of the two photons. The distribution has a kinematic edge at leading order at $\Delta\eta\simeq 1.79$. Above this point it features a Sudakov shoulder: fixed order perturbative corrections are not trustworthy and resummation is required. However, the bulk of the distribution can be calculated in fixed order perturbative QCD.

In Figure 4(b), we present the $p_{T}$ distribution of the leading photon. At LO, the photon $p_{T}$ cannot exceed the value of $m_{h}/2$. In addition, the experimental selection imposes a lower $p_{T}$ value at $0.35\,m_{h}$. This severe restriction of the phase space leads to large corrections to the available perturbative results. As such, resummation would be required to obtain stable predictions in this kinematical regime. Beyond LO the phase space for larger values of $p_{T}$ opens up and the distribution becomes more well-behaved beyond 100 GeV.

These distributions are just a small selection of possible observables that can be computed in our framework. Combining the production with other decay modes is straightforward and can be used to study a number of phenomenologically relevant observables.

6 Conclusion

In this article, we presented differential distributions for Higgs boson observables at NNLO in perturbative QCD, obtained within our “Higgs-differential” framework. Our computational method is a departure from common frameworks for differential calculations in that it avoids explicit subtraction of infrared divergences at the cost of being inclusive in jet observables. This is achieved by separating the QCD radiation phase space from the phase space of the produced colour-neutral final state. By integrating the QCD radiation phase space inclusively in dimensional regularisation, soft and collinear divergences are made explicit as poles in the regulator. Furthermore, this separation enables us to employ reverse unitarity and related techniques that have been developed in inclusive calculations, in order to simplify the phase-space integrations over final-state partons which are produced in association with the Higgs boson.

The NNLO cross section has been cast in terms of few master integrals, which we compute in an arbitrary number of dimensions in terms of standard hypergeometric functions that admit expansions to arbitrarily high order in the dimensional regulator. We have presented results that go beyond the finite term in the $\epsilon$ expansion, which were unknown in the literature. These new results are essential ingredients for the construction of collinear and UV counter terms for Higgs-differential cross sections at N3LO. The master integrals encountered here do not depend on the exact nature of the colourless final state produced and could be directly re-used in a similar differential Drell-Yan calculation.

Finally, we implemented the Higgs boson cross section through NNLO in a numeric code and tested its efficiency in kinematic distributions for the Higgs boson and its decay products in the diphoton signal. The predictions for these distributions were compared against results obtained from existing Monte-Carlo generators, thus validating our approach.

The main motivation for the approach presented in this article is its extensibility to even higher orders in QCD perturbation theory. The general framework of separating the QCD radiation phase space from the Higgs phase space persists at higher orders. This completely isolates the numerical calculation of distributions from the complications of higher order QCD computations. These should only be reflected in the master integrals that appear at higher orders. In this respect, it has been encouraging that the NNLO master integrals for Higgs-differential cross sections are relatively simple to compute. Using techniques that have been honed in inclusive calculations at N3LO, it should therefore be possible to compute the required master integrals, leading to differential predictions for a hadron collider at N3LO in QCD perturbation theory.

7 Acknowledgement

We thank Babis Anastasiou for inspiring discussions and useful comments on the manuscript. We are also grateful to Simone Alioli, Claude Duhr and Achilleas Lazopoulos for fruitful discussions and to Alexander Huss for numerical comparisons and valuable exchange of views. SL, AP and CS are supported by the ETH Grant ETH-21 14-1 and the Swiss National Science Foundation (SNSF) under contracts 165772 and 160814. The work of FD is supported by the U.S. Department of Energy (DOE) under contract DE-AC02-76SF00515. BM is supported by the European Commission through the ERC grants “pertQCD” and “HICCUP”.

Appendix A Renormalization Factors

The renormalisation factors for the strong coupling constant and Wilson coefficient (34) required for a computation through N3LO [76] are given by

[TABLE]

The coefficients at the various orders in the coupling constant $\beta_{i}$ are given by the QCD beta function [77, 78, 79, 80].

The infrared counter term $\Gamma$ consists of convolutions [81] of splitting functions $P_{ij}^{(n)}$ and can be derived from the DGLAP equation. Its perturbative expansion required for an N3LO accurate calculation of the differential Higgs boson production cross section is given by

[TABLE]

Appendix B Wilson Coefficient

In the effective theory with $n_{f}$ light flavours and the top quark decoupled from the running of the strong coupling constant, the $\overline{\textrm{MS}}$-scheme Wilson coefficient reads [82, 52]

[TABLE]