Double Higgs boson production at NLO: combining the exact numerical result and high-energy expansion

TL;DR

This paper enhances predictions for Higgs boson pair production at NLO by combining exact numerical results with high-energy expansions, improving accuracy for boosted Higgs studies at current and future colliders.

Contribution

It introduces a method to merge exact NLO calculations with high-energy expansions using Padé approximations, extending the validity of virtual corrections across phase space.

Findings

Extended the virtual matrix element grid to higher transverse momenta.

Provided improved NLO predictions for 14 TeV and 100 TeV colliders.

Made the updated computational grid publicly available.

Abstract

We consider the next-to-leading order QCD corrections to Higgs boson pair production, using our recent calculation of the form factors in the high-energy limit. We compute the virtual corrections to the partonic cross section, applying Pad\'e approximations to extend the range of validity of the high-energy expansion. This enables us to compare to the exact numerical calculation in a significant part of the phase space and allows us to extend the virtual matrix element grid, based on the exact numerical calculation, to larger values of the (partonic) transverse momentum of the Higgs boson, which is important for boosted Higgs studies. Improved predictions for hadron colliders with centre-of-mass energies of $14\ \mathrm{TeV}$ and $100\ \mathrm{TeV}$ are presented. The updated grid is made publicly available.

Equations

$${\mathcal{M}}^{ab}$$

$$A_1^{\mu \nu }$$

$$A_2^{\mu \nu }$$

$$q_{ij}$$

$$X_0$$

$$F$$

$$F_1^{(0)}$$

$$F_2^{(0)}$$

$$F_1^{(1)}$$

$$F_2^{(1)}$$

$${\alpha }_s^{(6)}(\mu )$$

$$A_{\nu }^{(6)}(\mu )$$

$$F^{\mathrm{fin},(1)}=F^{(1),\mathrm{IR}}-K_g^{(1)}{F}^{(0)},$$

$$K_g^{(1)}$$

$$F^{\mathrm{fin},(1)}$$

$$t$$

$${\tilde{\mathcal{V}}}_{\text{fin}}$$

$$C$$

$${\mathcal{V}}_{\text{fin}}$$

$$\left(\beta (s_i),\cos ({\theta }_i),{\mathcal{V}}_{\mathrm{fin}}^{\mathrm{grid}}\left(P_i\right),\pm {\delta }_i\right)$$

$$\beta (s)=\sqrt{1-\frac{4m_h^2}{s}}\mathrm{and}\cos (\theta )=\frac{s+2t-2m_h^2}{s\beta (s)}.$$

$$p_T^2=\frac{tu-m_h^4}{s}\mathrm{and}u=2m_h^2-s-t$$

$${\mathcal{V}}_{\mathrm{fin}}^N$$

$${\omega }_i$$

$$\alpha$$

$${\delta }_{\alpha }$$

Full text

1 Introduction

A primary goal of the Large Hadron Collider (LHC) and future colliders is the exploration of the electroweak symmetry breaking sector of the Standard Model (SM). It is important to find out whether the Higgs potential is indeed of the form suggested by the Standard Model, where the trilinear and quartic Higgs boson self-couplings are completely determined by the Higgs boson mass and its vacuum expectation value. As a deviation of the trilinear coupling from its SM value would be a clear sign of physics beyond the Standard Model, it is crucial to have precise predictions for processes which allow the measurement of this coupling.

An important process in this regard is Higgs boson pair production via gluon fusion. While the cross section is about a factor of one thousand smaller than that of single Higgs boson production, it is nevertheless the Higgs boson pair production channel with the largest cross section. It also has the interesting feature that there is a delicate cancellation between triangle-type diagrams, containing the trilinear Higgs boson coupling $\lambda$, and box-type diagrams, containing only Yukawa couplings, such that deviations of the trilinear coupling from the SM value can lead to distinct features in observables such as the Higgs boson pair invariant mass ($m_{hh}$) distribution.

Measurements of double Higgs boson production in gluon fusion at the LHC already have led to constraints on the ratio $\kappa_{\lambda}=\lambda_{\rm{BSM}}/\lambda_{\rm{SM}}$ [1, 2, 3, 4], where currently $-5.0\leq\kappa_{\lambda}\leq 12.0$ at 95% confidence level [3] is the most stringent bound derived from Higgs boson pair production measurements. The tightest bounds typically result from the combination of various Higgs boson decay channels. Among these, an important channel is the $b\bar{b}b\bar{b}$ channel [5, 6] due to the large branching ratio of $H\to b\bar{b}$. Reconstructing the Higgs bosons from boosted jets is promising, not only in view of a potential 27 TeV or 100 TeV collider, but also at the HL-LHC [7]. However, such an analysis also requires precise predictions in the high-$p_{T}$, or large-$m_{hh}$ regime of the Higgs bosons, reaching $m_{hh}$ values of about 3 TeV at the LHC, which is a region where high-energy expansions are fully justified. In this paper we will combine the high energy expansion of Refs. [8, 9, 10] with the full NLO calculation of Refs. [11, 12, 13] to arrive at predictions which combine the virtues of both approaches in the kinematic ranges where they work best.

The leading-order (LO) contribution to Higgs boson pair production in gluon fusion already contains one loop, therefore the calculation of higher-order corrections is a complex task. While the LO calculation was performed some time ago [14, 15, 16], next-to-leading order (NLO) results with full top quark mass dependence became available only recently [11, 12, 17], based on a numerical evaluation of the corresponding two-loop integrals.

Analytic higher-order results are known in various approximations. NLO QCD corrections in the $m_{t}\to\infty$ limit, sometimes also called “Heavy Top Limit (HTL)”, or “Higgs Effective Field Theory (HEFT)” approximation, have been calculated in Ref. [18] using the so-called “Born-improved HTL” approximation, which involves rescaling the NLO results in the $m_{t}\to\infty$ limit by a factor $B_{\rm FT}/B_{\rm HTL}$, where $B_{\rm FT}$ denotes the squared LO matrix element in the full theory. In Ref. [19] an approximation called “FT${}_{\small\mbox{approx}}$” was introduced, which contains the real radiation matrix elements with full top quark mass dependence, while the virtual part is calculated in the Born-improved HTL approximation.

The next-to-next-to-leading order (NNLO) QCD corrections in the $m_{t}\to\infty$ limit have been computed in Refs. [20, 21, 22, 23]. The HTL results have also been improved in various ways: the virtual corrections have been supplemented by an expansion in $1/m_{t}^{2}$ in [24, 25] up to order $1/m_{t}^{12}$ at NLO and $1/m_{t}^{4}$ at NNLO. Real radiation corrections, which involve three closed top quark loops have been expanded up to $1/m_{t}^{8}$ in Ref. [26]. Soft gluon resummation has been performed at NNLO+NNLL level in [27]. In Ref. [28], the NNLO calculation in the HTL of Ref. [23] has been combined with results including the full top quark mass dependence at NLO as well as in the $2\to 4$ matrix elements present in the NNLO real radiation. The latter results have been supplemented by soft gluon resummation in Ref. [29]. Analytic approximations for the top quark mass dependence of the two-loop amplitudes entering $gg\to HH$ at NLO have also been studied in the high-energy limit [8, 9, 10], around the top pair threshold expansion combined with large mass expansion [30], and for small Higgs boson transverse momentum [31].

The full NLO calculation of Refs. [11, 12] has been combined [13, 32, 33] with parton showers within the POWHEG-BOX-V2 [34, 35, 36] and MG5_aMC@NLO [37, 38] frameworks as well as within Sherpa [39]. Ref. [33] contains a discussion of showered results for non-SM values of the trilinear Higgs coupling, as well as a comparison of Pythia8.2 [40] and Herwig7.1 [41] showers in combination with Powheg.

The purpose of this paper is to provide results for the process $gg\to HH$ at NLO which are valid and accurate in the low-, medium- and high-energy regimes. This is achieved by combining the high-energy expansion, computed in Refs. [8, 9, 10], with the existing grid of the exact NLO result [13, 42], such that the finite part of the virtual amplitude can be evaluated at any phase space point without having to do costly two-loop numerical integrations. Previously, the grid of the exact NLO result was constructed based only on unweighted events, which are sparse in the high-energy region, and the grid was therefore statistically limited in the high-energy region. Extending the grid to higher energies using the exact NLO result would require the costly evaluation of additional phase-space points in a regime where the numerical convergence of the two-loop integrals can be poor. Instead, by combining the existing grid with analytic results obtained through a high-energy expansion, after a careful assessment of the regions in which the latter leads to an improvement, we are able to present results with small uncertainties over the full kinematic range. This improvement is particularly relevant for highly boosted Higgs bosons, for which the previous grid was unreliable. Parton shower Monte Carlo programs based on the new grid, presented here, can reliably be used to make predictions in an extended kinematic range.

The remainder of the paper is structured as follows. In the next section we introduce our notation and in Section 3 we describe our approach to obtain Padé approximations for the NLO virtual corrections based on the high-energy expansion of the form factors. This approach is validated in Section 4 at the level of the master integrals. In Section 5 we present numerical results for the virtual corrections and, in Section 6, we study their impact on the transverse momentum and invariant mass distributions. We conclude in Section 7.

2 Notation and conventions

The analysis we perform in this paper is based on the results for the form factors obtained in Refs. [8, 9]. Let us briefly repeat the notation and conventions introduced in these references.

The amplitude for the process $g(q_{1})g(q_{2})\to H(q_{3})H(q_{4})$, with all momenta $q_{i}$ defined to be incoming, can be decomposed into two Lorentz structures

[TABLE]

where $a$ and $b$ are adjoint colour indices, $s=(q_{1}+q_{2})^{2}$ is the squared partonic centre-of-mass energy and the two Lorentz structures are given by

[TABLE]

with

[TABLE]

where $s,{t=(q_{1}+q_{3})^{2}}$ and ${u=(q_{2}+q_{3})^{2}}$ are Mandelstam variables which fulfill $s+t+u=2m_{h}^{2}$, $T_{F}=1/2$, $G_{F}$ is Fermi’s constant and $\alpha_{s}(\mu)$ is the strong coupling constant evaluated at the renormalization scale $\mu$.

We define the expansion in $\alpha_{s}$ of the form factors as

[TABLE]

and decompose the functions $F_{1}$ and $F_{2}$ introduced in Eq. (1) into “triangle” and “box” form factors. We thus cast the one- and two-loop corrections in the form

[TABLE]

$F^{(1)}_{\rm dt1}$ and $F^{(1)}_{\rm dt2}$ denote the contribution from one-particle reducible diagrams such as the one shown in Fig. 1(f). In Ref. [9] this contribution has not been considered since the full top quark mass dependence is available from Eqs. (24), (25) and (26) of Ref. [43].

At this point a comment on the definition of $\alpha_{s}$ is in order. In Ref. [9] $\alpha_{s}$ has been defined with six active flavours which is an appropriate choice for the high-energy limit. In this paper, we compare to Ref. [12] where a five-flavour $\alpha_{s}$ has been used. Thus, we have to transform $\alpha_{s}$ and the gluon wave function from the six-flavour to the five-flavour theory using the relations

[TABLE]

where $A_{\nu}$ is the gluon wave function. As can be seen from these expressions the additional terms cancel because the number of external gluon fields equals the number of strong couplings $g_{s}$ in the Born amplitude, such that the resulting analytic expressions remain identical.

After the renormalization of the ultraviolet divergences, the form factors are still infrared divergent. Infrared finite results can be obtained by making a suitable subtraction [44]

[TABLE]

where $K_{g}^{(1)}$ is given by

[TABLE]

$\gamma_{E}$ is Euler’s constant and $C_{A}=3$ is a colour factor. After the decoupling of the top quark we have $n_{l}=5$ as the number of active flavours.111In Ref. [9] infrared subtraction has been performed in QCD with six active flavours. Note that the choice of $K_{g}^{(1)}$ is not unique. For example, finite form factors are also obtained if the $\mu$-dependent factor multiplies only the $1/\epsilon^{2}$ term inside the square brackets of Eq. (9), and not the $1/\epsilon$ term. The resulting form factors differ by terms proportional to $\log(\mu^{2}/(-s-i\delta))$. For the definition of $K_{g}^{(1)}$ in Eq. (9) it is convenient to introduce

[TABLE]

where $\tilde{F}^{(i)}=F^{{\rm fin},(i)}(\mu^{2}=-s)$, and $\beta_{0}=11C_{A}/12-T_{F}n_{l}/3$.

In Ref. [9] we express the analytic results for the form factors in terms of $m_{t},s,t$ and $m_{h}$. Note that our two-loop expressions are Taylor expanded222Let us stress that only the form factors are expanded in $m_{h}$ and the factor $3m_{h}^{2}/(s-m_{h}^{2})$ in Eq. (5) is kept exact. in $m_{h}$ including terms up to $m_{h}^{2}$. In Section 3 we use the (infrared-finite) form factors to construct the virtual corrections for the process $gg\to HH$. We adapt the notation of Ref. [12] and express our results in terms of the variables $m_{t},s,m_{h}$ and $p_{T}^{2}$. This is achieved using the relation

[TABLE]

and a subsequent re-expansion of the two-loop form factors in $m_{h}$ up to order $m_{h}^{2}$. We use the exact expression for the one-loop corrections [15, 16] and thus no expansion is necessary.

3 Padé improved virtual corrections

We adapt the notation of Ref. [13] and define (see also Eq. (4.1) of Ref. [30])

[TABLE]

with

[TABLE]

Here $\alpha_{s}$ corresponds to the five-flavour strong coupling constant. Furthermore, we introduce

[TABLE]

In Ref. [42] a grid of 3398 phase-space points is provided in the file Virt_full_noas.grid where the result for the phase-space point $P_{i}=(s_{i},t_{i})$ is given in the format

[TABLE]

with the coordinates $\beta(s)$ and $\cos(\theta)$ given by (see page 4 of Ref. [13])

[TABLE]

We use Eq. (11) together with

[TABLE]

in order to switch to the coordinates $({\sqrt{s}},p_{T})$ in the following.

For the numerical evaluation of ${\cal V}_{\rm fin}$ we proceed as follows:

• •

After inserting into Eq. (12) the exact one-loop and two-loop one-particle reducible form factors and the high-energy expansion of the remaining parts ${\cal V}_{\rm fin}$ can be written as

[TABLE]

where $\mathcal{V}_{0}$ contains all parts that are exact in $m_{t}$ and $m_{h}$ (i.e., $F^{(0)}_{\textnormal{tri}}$, $F^{(0)}_{\textnormal{box1}}$, $F^{(0)}_{\textnormal{box2}}$, $F^{(1)}_{\textnormal{dt1}}$ and $F^{(1)}_{\textnormal{dt2}}$) and the second term in Eq. (18) contains those parts which involve333Exact results for $F^{(1)}_{\textnormal{tri}}$ are available from Refs. [45, 46, 47]. For simplicity, in the following we nevertheless use our expansions which provide a very good approximation of the exact result [9]. $F^{(1)}_{\textnormal{tri}}$, $F^{(1)}_{\textnormal{box1}}$ and $F^{(1)}_{\textnormal{box2}}$. In Eq. (18) we explicitly show the dependence on $m_{t}$ but suppress dependence on $m_{h}$; note that ${\cal V}_{i}$ contains an expansion up to $m_{h}^{2}$.

• •

At this point we fix all numerical values except the top quark mass, i.e., $\sqrt{s}$, $p_{T}$ and $m_{h}$.

• •

Next we apply the replacements444$\log m_{t}$ terms are not replaced. $m_{t}^{2k}\to m_{t}^{2k}x^{k}$ and $m_{t}^{2k-1}\to m_{t}^{2k-1}x^{k}$ for the odd and even powers of $m_{t}$. We insert the numerical value for $m_{t}$ and consider ${\cal V}_{\rm fin}^{N}$ as an expansion in $x$. In Ref. [9] terms up to order $m_{t}^{16}$ were presented. Since then the expansion has been extended to $m_{t}^{32}$ which implies that ${\cal V}_{\rm fin}^{N}$ is available up to $x^{16}$. The analytic results for the form factors can be obtained from [48].

• •

Next we construct Padé approximants of ${\cal V}_{\rm fin}^{N}$ in the variable $x$ which means that we write Eq. (18) as a rational function of the form

[TABLE]

where $a_{i}$ and $b_{i}$ are determined by comparing the coefficients of $x^{k}$ after expanding the right-hand side of Eq. (19) in $x$.

As an alternative approach one can construct Padé approximations for ${\cal V}_{\rm fin}^{N}-{\cal V}_{0}$, which have $a_{0}=0$ and different values for the remaining coefficients. Both approaches lead to very similar final results, so in our analysis we concentrate on the one outlined in Eq. (19).

• •

For $N=32$, Padé approximations with $n+m=16$ can be constructed. We restrict our analysis to Padé approximants which are close to “diagonal” (where $n=m$). We require $|n-m|\leq 2$. Furthermore, we demand that expansions include at least terms up to order $m_{t}^{30}$. This leads to a list of five Padé approximants $\mathbb{Q}=\left\{[7/8],[8/7],[7/9],[8/8],[9/7]\right\}$.

• •

We aim for an approximation of ${\cal V}_{\rm fin}$ in the two-dimensional $\sqrt{s}$–$p_{T}$ plane where for each point a separate Padé approximant is constructed. Due to the structure of the ansatz (Eq. (19)), the Padé approximants may develop poles in the complex $x$ plane. Poles close to $x=1$ might lead to unphysical results. For this reason we assign a weight to each Padé approximant, which depends on the distance of the closest pole to $x=1$, and use this information to construct for each pair $(\sqrt{s},p_{T})$ a central value and an estimate of the uncertainty. In detail, we proceed as follows

• –

For each phase-space point $(\sqrt{s},p_{T})$ we compute for all Padé approximants in $\mathbb{Q}$ (see above) the value at $x=1$ and the distance of the closest pole which we denote by $\alpha_{i}$ and $\beta_{i}$, respectively.

• –

We introduce a re-weighting function, which reduces the impact of values $\alpha_{i}$ from Padé approximations with poles close to $x=1$. We define

[TABLE]

and assign $\omega_{i}$ to each value $\alpha_{i}$.

• –

We use the values $\alpha_{i}$ and $\omega_{i}$ to compute the central value from the weighted average and the uncertainty from the standard deviation as follows

[TABLE]

This procedure provides for each point $(\sqrt{s},p_{T})$ a result of the form $\alpha\pm\delta_{\alpha}$ which is based on Padé approximation.

4 Padé improved master integrals

In this section we construct $[8/8]$ Padé approximants (see Eq. (19)) at the level of the master integrals, for which numerical results can be obtained using FIESTA [49] and pySecDec [50]. In Fig. 2 we show the real and imaginary parts of the non-planar seven-line master integrals $G_{59}(1,1,1,1,1,1,1,-1,0)$ and $G_{59}(1,1,1,1,1,1,1,-2,0)$ (see Refs. [8, 9] for the notation and graphical representations) as a function of $\sqrt{s}$. In each panel several lines are shown which correspond to different choices of $p_{T}$. For better readability we shift some of the lines such that they are well separated, at least in some parts of the phase space, which leads to arbitrary units on the $y$-axis. Solid lines correspond to the Padé approximant555Similar results are also obtained for other choices. $[8/8]$ and the dots are obtained using pySecDec. One observes an impressive agreement between the Padé-improved and numerical results, even for the lower $p_{T}$ values around 100–200 GeV (the lower, blue-coloured lines). The small spikes visible above $\sqrt{s}=500$ GeV in some of the plots are due to the proximity of poles in the complex plane of the $[8/8]$ Padé approximants. In our final results, such spikes are removed by the re-weighting procedure described at the end of Section 3.

For illustration we show for $p_{T}=350$ GeV the results of the asymptotic expansions up to order $m_{t}^{30}$ and $m_{t}^{32}$ as dashed curves. For $\sqrt{s}\approx 2000$ GeV reasonable agreement is found with the numerical result and the Padé approximation. However, for smaller values of $\sqrt{s}$ one observes that the expansions quickly deviate from the exact result.

We obtain similar results for all non-planar master integrals and are thus confident that the procedure of Section 3 applied to ${\cal V}_{\rm fin}^{N}$ will provide a good approximation, even for relatively small values of $p_{T}$.

5 Numerical results for ${\cal V}_{\rm fin}$

In this section we consider ${\cal V}_{\rm fin}$ as a function of $\sqrt{s}$ and $p_{T}$ and compare to the exact results obtained in [12]. The results of [12] are available from [42] in the form of a grid in the $\sqrt{s}$–$p_{T}$ plane, where an uncertainty from numerical integration is assigned to each data point. For the renormalization scale the value $\mu=\sqrt{s}/2=m_{hh}/2$ has been chosen. Furthermore we use the values $m_{t}=173~{}\mbox{GeV}$ and $m_{h}=125~{}\mbox{GeV}$.

In Fig. 3 we show all data points from [42], normalized to their central values, as a function of $p_{T}$ (the dark blue points with uncertainty bars). Note that in general, different data points belong to different values of $\sqrt{s}$. Fig. 3 also contains Padé results for ${\cal V}_{\rm fin}^{N}$ constructed from $N=30$ and $N=32$ input, again normalized to the central values of the grid points from [42] (the coloured points without uncertainty bars). Additionally, the results of the expansions ${\cal V}_{\rm fin}^{30}$ and ${\cal V}_{\rm fin}^{32}$ are shown as green and light-blue data points, respectively. Note that the data points based on ${\cal V}_{\rm fin}^{N}$ are computed using the same input values as those of the grid points.

As expected, good agreement is found for large values of $p_{T}$ (which implies large values of $\sqrt{s}$). Most of the data points lie within one sigma of the grid points [42]. One also observes that some of the points are outside the one-sigma range, however, still agree within two sigma. The interesting region of Fig. 3 is $p_{T}\lesssim 400$ GeV. Here, the high-energy expansion diverges rather quickly and the agreement with the grid points breaks down. The Padé-improved results, however, follow the dark blue points until $p_{T}\approx 200$ GeV. Some of the Padé approximants reproduce the exact numerical result even down to $p_{T}\approx 150$ GeV with reasonable precision. This behaviour motivates a closer look into the comparison of Padé-improved and numerical results for fixed values of $p_{T}$.

We now fix $p_{T}$ and consider ${\cal V}_{\rm fin}$ as a function of $\sqrt{s}$. For small values of $p_{T}$ and $\sqrt{s}$ the grid points are dense. However, for $p_{T}\gtrsim 300$ GeV and/or $\sqrt{s}\gtrsim 1000$ GeV they become quite sparse. Furthermore, if one wants to perform an analysis for fixed $p_{T}$ one can in principle only use a few data points from the grid which makes a comparison difficult. On the web-page [42] an interpolation routine is provided which allows for an extension of the grid points to the whole phase space. However, we find that in regions where the grid is only sparsely populated this interpolation routine seems to provide unreliable results. In order to separate interpolated points with solid support from nearby grid-points from interpolated points without such support, we enhance an interpolated data point at $P_{0}=\left(\sqrt{s}_{0},p_{T,0}\right)$ by an error estimate as follows:

• •

Define a region around $P_{0}$ as

\Delta=\left\{(\sqrt{s},p_{T})\Big{|}|\sqrt{s}-\sqrt{s_{0}}|\leq 5~{}\mathrm{GeV},|p_{T}-p_{T,0}|\leq 10~{}\mathrm{GeV}\right\}.

• •

$\mathbb{P}$ is the set of data points of the grid [42] which lie in $\Delta$: $\mathbb{P}=\{{\cal V}_{\rm fin}^{\rm grid}(P_{1})\pm\delta_{1},{\cal V}_{\rm fin}^{\rm grid}(P_{2})\pm\delta_{2},\ldots,{\cal V}_{\rm fin}^{\rm grid}(P_{n})\pm\delta_{n}\}$, where $\delta_{i}$ are the corresponding uncertainties.

• •

If $\mathbb{P}$ is empty no uncertainty can be assigned to the interpolated value ${\cal V}_{\rm fin}^{\rm int}(P_{0})$. Note that such a point has no support from the actual grid points.

• •

For non-empty set $\mathbb{P}$ we define $\sigma=\sum_{i=1}^{n}\left|\delta_{i}\right|/n$ as a mean uncertainty assigned to ${\cal V}_{\rm fin}^{\rm int}(P_{0})$.

In Fig. 4 we show ${\cal V}_{\rm fin}$ as a function of $\sqrt{s}$ for four different values of $p_{T}$. The blue dots correspond to the results obtained from the grid [42] using the procedure described above. Blue dots with no uncertainty bar have no nearby grid points from which one can estimate an uncertainty. The other dots correspond to our Padé-improved results which are obtained using the prescription from Section 3. If no uncertainty is visible for these points, it is smaller than the dot size.

Let us start with the discussion of Figs. 4(a) and (b) which correspond to $p_{T}=250$ GeV and $p_{T}=200$ GeV, respectively. For $\sqrt{s}<800$ GeV both the numerical and the Padé results have small uncertainties and agree very well. Between $\sqrt{s}\approx 800$ GeV and $\sqrt{s}\approx 1400$ GeV the Padé results behave smoothly but the (interpolated) numerical results show strong variation which is due to the interpolation procedure used in Ref. [42]. This is also true for $\sqrt{s}>1400$ GeV where the results from [42] show an unphysical constant behaviour. This behaviour suggests that above $\sqrt{s}\approx 800$ GeV one should not trust the results of [42] but rather the approximations obtained from the high-energy expansion [8, 9]. On the other hand, for $\sqrt{s}\lesssim 800$ GeV, the good agreement of the Padé results with the numerical calculation provides confidence regarding the reliability of the Padé procedure.

For $p_{T}=150$ GeV, see Fig. 4(c), the Padé procedure develops uncertainties of about 10% to 20% for $\sqrt{s}\lesssim 800$ GeV. It is nevertheless quite impressive that agreement with the numerical results, which have small uncertainties, is found. For higher values of $\sqrt{s}$ it seems that one can trust the results from [42] up to about $\sqrt{s}=1300$ GeV, above which they again become constant, which is unphysical.

Although it is far from the region of convergence of the high-energy expansion, we show in Fig. 4(d) the results for $p_{T}=100$ GeV. Here, the Padé method develops large uncertainties over the whole range of $\sqrt{s}$. It is, however, interesting to note that the central value shows good agreement with the numerical results for $\sqrt{s}\lesssim 1500$ GeV. In this plot we also show, as a solid red curve, results for the large-$m_{t}$ expansion of ${\cal V}_{\rm fin}$, which is constructed using the large-$m_{t}$ expansion of the form factors, computed to order $1/m_{t}^{12}$ in [22]. We observe agreement with the exact results (blue dots) up to $\sqrt{s}\approx 400$ GeV which constitutes a good consistency check.

The discussion of the plots in Fig. 4 shows that the Padé method provides accurate results even for relatively small values of $p_{T}$. Furthermore, it provides realistic estimates of the uncertainties. In Fig. 5 we show ${\cal V}_{\rm fin}$ as a function of $\sqrt{s}$ for fixed values of $p_{T}$ (shown in different colours, see the plot legend for details). The plot contains the curves for the four $p_{T}$ values of Fig. 4 and a further eight choices of $p_{T}$, with the highest value $p_{T}=650$ GeV. The dots represent the results from [42]. Where available, the uncertainties are explicitly indicated. For $p_{T}\geq 200$ GeV the Padé results are shown as solid lines. Note that in this region of the phase space the uncertainty is below the thickness of the lines. One observes that the solid lines agree with the data points within the indicated uncertainties, which are in general much larger than the Padé uncertainty. For $p_{T}=100$ GeV and $p_{T}=150$ GeV we reproduce in Fig. 5 the curves from Fig. 4 (see black and dark violet data points).

We now define a criterion which provides a prescription for the improvement of the grid [42]. In order to have guidance we show in Fig. 6 the relative uncertainty of the Padé results in the $\sqrt{s}$–$p_{T}$ plane. We also overlay all grid points from [42] and use the same colour scale for their uncertainties. Note that the kinematic boundary is obtained from the requirement that $1-4(m_{h}^{2}+p_{T}^{2})/s$ (see Eq. (11)) is positive.

From Fig. 6 we learn that the uncertainty is below $0.1\%$ for $p_{T}\gtrsim 200$ GeV and then grows towards lower $p_{T}$ relatively quickly. Still, even for $p_{T}\approx 150$ GeV the uncertainty is around a few percent for most values of $\sqrt{s}$. Note that larger relative uncertainties for larger values of $\sqrt{s}$ are observed since in this region $\widetilde{\mathcal{V}}_{\textnormal{fin}}$ is small.

On the basis of this observation we extend the grid provided in [42] as follows:

• •

We increase the number of points computed using the full NLO result from 3398 to 6320. The new points are sampled according to the distribution of unweighted events and, therefore, populate the same kinematic regime as the original points.

• •

For $\sqrt{s}\geq 700$ GeV and $p_{T}\geq 150$ GeV we add points from the Padé approximation.

• •

For $\sqrt{s}<700$ GeV and $p_{T}\geq 200$ GeV we add points from the Padé approximation.

The boundary above which we include points from the Padé approximation is denoted as a yellow line in Fig. 6. We note here that if one reproduces Figs. 5 and 6 using the 6320 points described above the behaviour is qualitatively the same and we therefore refrain from showing them in this paper.

In Fig. 7 we compare the Padé results to the improved version of the grid, which provides precise results in the whole relevant phase space. We note that the wiggly behaviour and the deviation of the grid data points from the Padé approximation for larger values of $\sqrt{s}$ and smaller values of $p_{T}$ could be improved by including further data points from the Padé approximation. This behaviour would then be pushed to higher values of $\sqrt{s}$. We judge the performance of the grid as displayed by Fig. 7 to be sufficient for the phenomenological applications of this paper, and further improvements of the grid not to be necessary. This improved grid can be downloaded from [42].

6 Applications

In the following we discuss differential distributions w.r.t. the Higgs boson pair invariant mass $m_{hh}$ and the “single inclusive” Higgs boson transverse momentum $p_{T,h}$ for hadronic centre-of-mass energies $\sqrt{s_{H}}=14$ TeV and $\sqrt{s_{H}}=100$ TeV. The emphasis of this analysis is the comparison of the current [42] and improved grid introduced in the previous section.

For our analysis we use the parton distribution functions PDF4LHC15_nlo_100_pdfas [51, 52, 53] and adopt the corresponding value for $\alpha_{s}$. For the top quark and Higgs boson masses we use $m_{h}=125$ GeV and $m_{t}=173$ GeV and choose $\mu_{0}=m_{hh}/2$ as the central value for the renormalization ($\mu_{R}$) and factorization ($\mu_{F}$) scales. The uncertainties due to higher-order QCD corrections are estimated using the usual seven-point scale variation around $\mu_{0}$, i.e., for $\mu_{R}$ and $\mu_{F}$ we introduce $\mu_{R,F}=c_{R,F}\mu_{0}$ with $c_{R,F}\in\{0.5,1,2\}$ and omit the extreme choices $(c_{R},c_{F})=(0.5,2)$ and $(c_{R},c_{F})=(2,0.5)$.

In Fig. 8 we show our results for $\sqrt{s_{H}}=14$ TeV. In the upper panels we present the $m_{hh}$ and $p_{T,h}$ differential distributions, and in the lower panels we display the ratio of the NLO corrections to the LO values ($K$ factor). The LO values are shown in black and the coloured curves correspond to different versions of the NLO prediction, all of which contain the full real radiation corrections and only differ in the way that the virtual corrections are implemented. The blue curve, denoted “FTapprox”, incorporates the virtual corrections computed in the infinite top quark mass limit and rescaled by the exact LO prediction. The red curve is based on the grid constructed in Ref. [13] but improved by increasing the number of points from 3398 to 6320 (see discussion above). Finally, the green curve is based on the new grid, the construction of which is described in Section 5. This curve constitutes our best prediction. The grey and green bands around the corresponding curves have been obtained by independent variations of $\mu_{R}$ and $\mu_{F}$ as described above.

It is interesting to note that for small $m_{hh}$ and $p_{T,h}$ there is perfect agreement of the red and green curves, which is expected since in this region the dependence on ${\cal V}_{\rm fin}$ comes primarily from the region in the (partonic) $\sqrt{s}$–$p_{T}$ plane where the support of the old grid was dense. For higher values of $m_{hh}$ and $p_{T,h}$, one observes a difference between the red and the green curves. However, in both cases the red curve lies well within the green uncertainty band.

The $m_{hh}$ and $p_{T,h}$ distributions for $\sqrt{s}_{H}=100$ TeV are shown in Fig. 9, where the same notation is used as in Fig. 8. Note that now a significant difference is observed between the red and green curves; for higher values of $m_{hh}$ and $p_{T,h}$ the red curve lies outside the green uncertainty band. As an example let us consider $p_{T,h}=2000$ GeV. For this value the $K$ factor is reduced from $K\approx 1.7$ to $K\approx 1.5$ after including the high-energy results in the grid.

Let us mention that in Figs. 8 and 9, the same phase-space points have been used for all curves. Thus, the differences between the curves is only due to the different implementations of the virtual corrections.

We should emphasize that one observes no change in the total cross section due to the change from the red to the green curve, since the main contribution to $\sigma_{\rm tot}$ comes from smaller centre-of-mass energies. However, Figs. 8 and 9 show that it is important to use the improved grid for phenomenological analyses, if one wishes to consider large values of $m_{hh}$ or $p_{T,h}$, even for $\sqrt{s_{H}}=14$ TeV. In these regions the predictions based on “FTapprox” deviate significantly from the green curve.

7 Conclusions

We provide optimized predictions for the NLO corrections to Higgs boson pair production by combining the exact numerical results with analytic expressions for the form factors obtained in a high-energy expansion. For the latter the region of convergence is significantly improved by constructing Padé approximants, which are validated at the level of master integrals. Furthermore, we identify regions in the phase space where both the exact numerical evaluations and the Padé results provide precise predictions and find good agreement. We thus combine both approaches and generate a new grid which is available from [42]. The analytic expressions for the high-energy expansion of the form factors are available from [48].

We apply the improved grid to phenomenological studies of the Higgs boson pair invariant mass and Higgs boson transverse momentum distributions at LHC energies and for $\sqrt{s_{H}}=100$ TeV. We show that at high energies the improvements are noticeable and we recommend to use the updated grid for phenomenological studies, even for $\sqrt{s_{H}}=14$ TeV.

Acknowledgements

This research was supported in part by the COST Action CA16201 (‘Particleface’) of the European Union and by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under grant 396021762 — TRR 257 “Particle Physics Phenomenology after the Higgs Discovery”. MK was supported in part by the Swiss National Science Foundation (SNF) under grant number 200020-175595. DW acknowledges the support of the DFG-funded Doctoral School KSETA.