Higgs boson pair production at NNLO with top quark mass effects
Received: March
Higgs boson pair production at NNLO with top quark mass e ects
M. Grazzini 0 1 3 6 7
G. Heinrich 0 1 3 4 7
S. Jones 0 1 3 4 7
S. Kallweit 0 1 3 5 7
M. Kerner 0 1 3 4 7
J.M. Lindert 0 1 2 3 7
J. Mazzitelli 0 1 3 6 7
0 CH1211 Geneva 23 , Switzerland
1 Fohringer Ring 6 , 80805 Munchen , Germany
2 Institute for Particle Physics Phenomenology, Durham University
3 Winterthurerstrasse 190 , CH8057 Zurich , Switzerland
4 Max Planck Institute for Physics
5 TH Division, Physics Department , CERN
6 PhysikInstitut, Universitat Zurich
7 Durham DH1 3LE , U.K
We consider QCD radiative corrections to Higgs boson pair production through gluon fusion in proton collisions. We combine the exact nexttoleading order (NLO) contribution, which features twoloop virtual amplitudes with the full dependence on the top quark mass Mt, with the nexttonexttoleading order (NNLO) corrections computed in the largeMt approximation. The latter are improved with di erent reweighting techniques in order to account for niteMt e ects beyond NLO. Our reference NNLO result is obtained by combining oneloop doublereal corrections with full Mt dependence with suitably imation. We present predictions for inclusive cross sections in pp collisions at p reweighted realvirtual and doublevirtual contributions evaluated in the largeMt approx27 and 100 TeV and we discuss their uncertainties due to missing Mt e ects. Our approximated NNLO corrections increase the NLO result by an amount ranging from +12% at
NLO Computations; QCD Phenomenology

s = 13 TeV to +7% at p
s = 100 TeV, and the residual uncertainty of the inclusive cross
section from missing Mt e ects is estimated to be at the few percent level. Our calculation
is fully di erential in the Higgs boson pair and the associated jet activity: we also present
predictions for various di erential distributions at p
s = 14 and 100 TeV, and discuss the
size of the missing Mt e ects, which can be larger, especially in the tails of certain
observables. Our results represent the most advanced perturbative prediction available to date
for this process.
1 Introduction 2 Details on the method and approximations
boson production [5{10] or electroweak precision observables [11, 12].
For the gg ! hh production channel, the leading order (LO) calculation was performed
some time ago in refs. [13{15]. The nexttoleadingorder (NLO) corrections with full top
quark mass (Mt) dependence, involving twoloop diagrams with several mass scales, became
available only recently [16, 17], and have been supplemented by softgluon resummation at
small transverse momenta of the Higgs boson pair [18] and parton shower e ects [19, 20].
In the Mt ! 1 limit, also called Higgs E ective Field Theory (HEFT) approxima
tion, pointlike e ective couplings of gluons to Higgs bosons arise. In this limit, the NLO
corrections were rst calculated in ref. [21] and rescaled by a factor BFT=BHEFT, where
BFT denotes the LO oneloop matrix element squared in the full theory. This procedure is
often called \Bornimproved HEFT" approximation.
In refs. [4, 22] an approximation for Higgs boson pair production at NLO, labelled
\FTapprox", was introduced, in which the real radiation matrix elements contain the full
{ 1 {
top quark mass dependence, while the virtual part is calculated at NLO in the HEFT
approximation and rescaled at the event level by the reweighting factor BFT=BHEFT.
At the inlusive cross section level this approximation suggests at the LHC a correction
with respect to the \Bornimproved HEFT" approximation of about
10%, close to the
corresponding correction of
14% later obtained in the full NLO calculation [16, 17].
The nexttonexttoleadingorder (NNLO) QCD corrections in the HEFT
approximation have been computed in refs. [23{26], where ref. [26] provides fully di erential results.
The NNLO HEFT results for the total cross section have been supplemented by an
expansion in 1=Mt2 in ref. [27]. Approximations for the topquark mass dependence of the
twoloop amplitudes in the NLO calculation have been studied in ref. [28] via a Pade
ansatz. Soft gluon resummation has been performed at NLO+NNLL in ref. [29] and at
NNLO+NNLL in ref. [30]. The NNLO+NNLL HEFT results lead to Kfactors of about
1.2 relative to the Bornimproved NLO HEFT result.
In ref. [31], the recommended value for the total gg ! hh cross section was based on the
NNLO+NNLL HEFT results [30], corrected by a factor t accounting for top quark mass
e ects, extracted from ref. [16]. However, this procedure is somewhat ad hoc, and not viable
to study kinematical distributions. In order to account for the NNLO Kfactor in the HEFT
calculation as well as for the correct description of the tt threshold and the highenergy
tails of the distributions, where the top quark loops are resolved, a rst attempt to combine
the two calculations has been made in ref. [17], where the full NLO result for a particular
distribution was reweighted by the NNLO Kfactor obtained from ref. [26] on a binbybin
basis. However, this procedure, called \NLOimproved NNLO" has its drawbacks, as it
needs to be repeated for each observable (and binning) under consideration.
The aim of this paper is to study alternative methods to combine the two results,
i.e. to incorporate top quark mass e ects in the calculation of the production of Higgs
boson pairs at NNLO. One of the studied approximations comprises exact topquark mass
dependence up to NLO and also exact top quark mass dependence in the doublereal
emission contributions to the NNLO cross section at di erential level. The results of this
approximation can be regarded as the most advanced prediction currently available for
Higgs boson pair production in gluon fusion.
This work is organized as follows: in section 2 we describe the technical details of our
calculation, and present the di erent approximations we will consider to incorporate mass
e ects in the NNLO contribution. In section 3 we present our numerical predictions, both
for the total cross section and di erential distributions. Finally, in section 4 we summarise
our results.
2
Details on the method and approximations
We start by presenting the di erent technical ingredients entering our computation, as
well as the de nition of the various approximate ways to include mass e ects in the NNLO
calculation introduced and used in this work. Finally, we also discuss the numerical stability
of our predictions.
{ 2 {
HJEP05(218)9
Our calculation is based on the publicly available computational framework Matrix [32],
which allows the user to perform fully di erential NNLO calculations for a wide class of
processes at hadron colliders. For the purpose of the present work, the public version of the
code has been extended, based on the calculation of ref. [26], to include the production of a
pair of Higgs bosons via gluon fusion. For the calculation of the NNLO corrections the code
implements the qT subtraction formalism [33], in which the genuine NNLO singularities,
located where the transverse momentum of the Higgs boson pair, pT;hh
qT , vanishes, are
explicitly separated from the NLOlike singularities in the hh + jet contribution. The qT
where in particular the contribution d NhhL+Ojet can be evaluated using any available NLO
subtraction procedure to handle and cancel the corresponding infrared (IR) divergencies.1
The remaining qT ! 0 divergence is canceled by the processindependent counterterm
d NCNTLO. The processdependence of the hardcollinear coe cient HNNLO enters only via
hh
the NNLO (HEFT) twoloop virtual corrections [23] through an appropriate subtraction
procedure [37].
The di erence in the square bracket of eq. (2.1) is nite when qT ! 0, but each of
the terms exhibits a logarithmic divergence. Therefore, a technical cut, rcut, needs to be
introduced on qT =Q, where the scale Q is chosen to be the invariant mass of the nalstate
system. More details about the rcut ! 0 extrapolation are provided in section 2.3.
At variance with the calculation of ref. [26], which was strictly done within the HEFT,
this time all the routines needed to compute the full NLO cross section as well as the
di erent NNLO reweightings have been implemented. This includes linking the code to
the NLO twoloop virtual corrections obtained via a grid interpolation [19] and to several
loopinduced amplitudes provided by the OpenLoops amplitude generator [38]. Within
this framework we reproduced the di erential NLO results of refs. [16, 17] at the per
mille level.
The grid for the NLO virtual twoloop amplitudes is based on the calculation
presented in refs. [16, 17], which in turn for the calculation of the twoloop amplitudes relies
on an extension of the program GoSam [39, 40] to two loops [41], using also Reduze2 [42],
SecDec3 [43] and the QuasiMonte Carlo technique as described in ref. [44] for the
numerical integration. These amplitudes (for xed values of the Higgs boson and top quark
masses) are provided in a twodimensional grid together with an interpolation framework,
which allows us to evaluate them at any phase space point without having to perform the
computationally costly twoloop integration. For more details, see refs. [19, 45].
All tree and oneloop amplitudes in the HEFT and also all loopsquared amplitudes
in the full theory as discussed below are obtained via a process independent interface
to OpenLoops [38, 46, 47]. For the latter this comprises loopsquared amplitudes for
1Matrix uses the automated implementation of the CataniSeymour dipole subtraction method [34, 35]
within the Monte Carlo program Munich [36].
{ 3 {
pp ! hh + 1; 2 jets, that need to be evaluated in IR divergent unresolved limits. In
particular the limit qT ! 0 represents a signi cant challenge for the numerical stability
of the hh + 2 jets amplitudes in the full theory. Thanks to the employed algorithms the
numerical stability is under control, as discussed in detail in section 2.3. A major element
of this stability originates from the employed tensor integral reduction library COLLIER [48].
2.2
Approximations for topmass e ects at NNLO
In the following we present three approximations for the NNLO Higgs boson pair production
cross section, which take into account nite top quark mass e ects in di erent ways. In all
cases, we always include the full NLO result when computing the NNLO prediction, and
only apply the di erent approximations to the O( S4) contribution.
NNLONLOi.
The NLOimproved NNLO approximation (NNLONLOi) has already been
presented in ref. [17]. It can be constructed based on an observablelevel multiplicative
approach. In this approximation, for each bin of each histogram we multiply the full NLO
result by the ratio between the HEFT NNLO and NLO predictions for this bin.
NNLOBproj.
A di erent approximation can be obtained by reweighting each NNLO
event by the ratio of the full and HEFT Born squared amplitudes. We denote this procedure
as Bornprojected approximation (NNLOBproj). Of course, in order to do so and due
to the di erent multiplicities involved, an appropriate projection to Bornlike kinematics
is needed; for this purpose we make use of the qT recoil procedure de ned in ref. [49].
Following this prescription, the momenta of the Higgs bosons remain unchanged, and the
new initialstate parton momenta are obtained by absorbing the recoil due to the additional
radiation. Speci cally, denoting the momenta of the incoming partons by p1 and p2, and
the momentum of the Higgs boson pair system by q, the new momentum to be used for
the LO projection k1 (then, k2 = q
k1) is given by
Q2
2
k1T
q p1
k1 = z1 2 q p1 p1 + k1T +
z1 Q2 p1 p2 p2 ;
k1T k1T
=
2
k1T ;
(2.2)
where
z1 =
Q2 + 2 qT k1T + q
(Q2 + 2 qT k1T )2
and k1T is a twodimensional vector in the qT plane which needs to ful ll the condition
k1T ! 0 when qT ! 0, and we set k1T = qT =2 (and therefore k2T = qT =2). This
condition guarantees that the subsequently applied reweighting does not spoil the NNLO
qT cancellation. More details about this procedure can be found in ref. [49].
NNLOFTapprox .
The third approximation we consider is constructed to pro t from the
fact that the doublereal emission contributions to the NNLO cross section require only
oneloop amplitudes in the full theory (FT) and can thus be computed by using OpenLoops.
Of course, the inclusion of these loopinduced amplitudes needs to be done in such a way
{ 4 {
HJEP05(218)9
that the dipole cancellations in the NLO hh + j calculation and the lowqT cancellation for
hh at NNLO are not spoiled.
We will de ne our approximation by using the following procedure: working in the
HEFT, for each nloop squared amplitude that needs to be computed for a given partonic
(n)
subprocess AHEFT(ij ! HH + X), we apply the reweighting
R(ij ! HH + X) =
AFBuolrln(ij ! HH + X)
(0)
AHEFT(ij ! HH + X)
;
where AFBuolrln stands for the lowest order (loopinduced) squared amplitude for the
corresponding partonic subprocess, computed in the full theory.2 We note that, contrary to what
happens in the Bornprojected approach, here the reweighting is de ned using amplitudes
that correspond to the same subprocess under consideration. Therefore, the kinematics is
always preserved and there is no need to de ne a Born projection. Moreover, for amplitudes
that are of treelevel type in the HEFT (as it is the case for the doublereal emission
contributions), this reweighting simply implies using the exact loopinduced amplitudes with
full top mass dependence. The reweighting procedure de ned by eq. (2.4) agrees at NLO
with the socalled FTapprox introduced in ref. [22], therefore we will use the same notation.
Given that the performance of the Bornprojection and FT approximations was already
studied in ref. [17] at NLO, we directly present NNLO predictions in section 3. We point
out that, based on the ingredients entering each of the approximations, the NNLOFTapprox
is expected to be the most advanced prediction for Higgs boson pair production via gluon
fusion. By contrast, the NNLOBproj is expected to be the less accurate, since it is based
on a simple Born level reweighting procedure. Nevertheless, and for comparison purposes,
we always present results for the three approximations described above.
(2.4)
Numerical stability
Before presenting our quantitative predictions, we brie y discuss the numerical stability
of our results. From the computational point of view, the most challenging of the three
approaches to incorporate mass e ects at NNLO is the NNLOFTapprox procedure, as it
involves loopinduced doublereal contributions in the full theory. In particular the dominant
gg ! hhgg amplitude comprises computationally very challenging sixpoint loop integrals
with internal masses. In fact, these contributions have to be evaluated in the numerically
intricate NNLO unresolved limits and to the best of our knowledge, the present calculation
is the rst application of a sixpoint oneloop amplitude integrated over its IR divergent
unresolved limits in an NNLO calculation.
Thanks to the numerical stability of the applied algorithms in OpenLoops together
with Collier, the bulk of the phasespace points remains stable in double precision
when approaching qT ! 0, even close to the dipole singularity, i.e. in the NNLO
doubleunresolved limits. On average the runtime per phase space point for the gg ! hhgg
2Strictly speaking, the reweighting is applied to the
nite part of the loop amplitudes. However, at
oneloop level this procedure reproduces the loop structure of the full theory.
{ 5 {
1.0
0.5
%
[
1
LO 0.0
N
N
σ/
σ
0.5
cut, rcut, normalized with respect to the extrapolated rcut ! 0 result. The dotted lines indicate
the symmetrized uncertainty coming from the extrapolation.
subtraction cut, rcut, for p
amplitude is
1 sec. In principle OpenLoops provides a rescue system, such that
remaining numerically unstable phasespace points can be reevaluated in higher numerical
precision based on reduction with CutTools [50]. However, the runtime of the loopinduced
gg ! hhgg amplitude in OpenLoops is signi cantly increased when CutTools is used
in quadruple precision (to the level of
10 minutes per phasespace point), rendering the
quadruple precision stability system prohibitive for this amplitude for practical purposes.3
Therefore, we restrict the evaluation to double precision and replace potentially unstable
phasespace points close to the dipole singularities, quanti ed by
Li = (pi pj =s^)min,
where the minimum among all potential emitter parton combinations i and j is taken,
with an approximation: below a technical cut
Li, cut we switch from the (loopinduced)
doublereal amplitude in the FT to the (treelevel) doublereal amplitudes in the HEFT,
reweighted at LO. This approach could in principle introduce a bias in the NLO hh+jet
cross section, thereby hampering the lowqT cancellation of the NNLO computation. We
have checked that this is not the case, as detailed in the following.
For the predictions presented in section 3 we use
Li, cut = 10 4 and we varied this
parameter in the range 10 3 to 10 5
, nding results that only di er at the per mille level
or below and which are always compatible within the numerical uncertainties. In
gure 1
we illustrate the resulting dependence of the NNLOFTapprox total cross section on the qT
s = 14 TeV. Due to the previously discussed stability challenges,
we considered values of rcut between 1% and 3:5%, which are larger than the ones typically
used in previous qT subtraction calculations (compared for instance with the default values
in the public Matrix release [32]). Nevertheless our results present a good stability, with
e ects that are below 0:2% in the whole qT =Q range under study, validating this choice.
The rcut ! 0 extrapolation is performed using a linear least
varying the upper bound of the interval (in this case starting from a minimum of 25 points,
which corresponds to an upper bound of rcut = 1:6%, and up to rcut = 3:5%). Then,
the result with the lowest
2=degreesoffreedom value is taken as the best t, and the
2 t. The t is repeated
3Here we want to note that these stability issues will be strongly mitigated in the future based on the
new OpenLoops onthe y reduction method introduced in ref. [47].
{ 6 {
rest is used to estimate the extrapolation uncertainty.4 In the case shown in gure 1 the
extrapolation uncertainty for rcut ! 0, indicated with the dotted lines, is
0:14%.
A further uncertainty arises due to the numerical evaluation of the twoloop integrals
with full topquark mass dependence in the virtual corrections of the NLO contribution.
The error of the numerical integration of the amplitudes is propagated to the total cross
section using Monte Carlo methods, varying the amplitude level results according to the
corresponding error estimates. This leads to changes of the NLO cross section at the per
mille level. Furthermore, we have checked that, within this uncertainty, results based on
the grid for the virtual amplitude are consistent with the ones directly obtained from the
amplitude results calculated in refs. [16, 17]. We want to point out that the uncertainties
can be somewhat larger in di erential results, in particular in the tails of pT and
invariantmass distributions.
This discussion shows that the uncertainties due to the qT subtraction method and
the numerical evaluation of the NLO virtual contribution and grid interpolation are clearly
under su cient control.
3
Results
In this section we present our numerical predictions for inclusive and di erential cross
sections for Higgs boson pair production in pp collisions. We consider centreofmass energies
of 13, 14, 27 and 100 TeV. For the sake of brevity, di erential distributions are presented
only for 14 TeV and 100 TeV. We use the values Mh = 125 GeV for the Higgs boson mass
and Mt = 173 GeV for the pole mass of the top quark.5 We do not consider bottom quark
loops, whose contribution at LO is below 1%. We also neglect top quark width e ects,
which at LO are at the level of 2% for the total cross section [22]. We use the PDF4LHC15
sets [53{58] of parton distribution functions (PDFs), with parton densities and
uated at each corresponding perturbative order (i.e., we use the (k + 1)loop running
at NkLO, with k = 1; 2). As renormalization and factorization scales, we use the central
value 0 = Mhh=2, and we obtain scale uncertainties via the usual 7point scale variation.
S
eval
S
3.1
Inclusive cross sections
In table 1 we present results for the total cross sections at NLO and NNLO in the various
approximations. At NLO we report the exact result, including the full Mt dependence,
and also the FTapprox result. By comparing the two NLO predictions, we see that the
FT approximation overestimates the exact NLO result by 4% (6%) at 14 (100) TeV. At
NNLO the largest prediction is obtained in the NNLOBproj approximation, resulting in an
4We note that, in the current Matrix release, the rcut ! 0 extrapolation and the ensuing uncertainty
estimation is only performed for inclusive ( ducial) cross sections. However, no signi cant e ects have been
observed for kinematic distributions in various dedicated studies (see for instance ref. [51]).
5We note that it is not yet possible to fully assess the e ect of the top mass scheme in the calculation
as the top quark mass is xed in the existing NLO twoloop virtual amplitudes. In the case of single Higgs
boson production, using an MS mass instead of the pole mass for the top quark changes the NLO cross
section by less than one per mille (see e.g. [52]). In the case of Higgs boson pair production, where the
sensitivity to the top quark mass is much stronger, we expect the e ect to be larger but not to exceed the
other uncertainties a ecting our calculation.
{ 7 {
NLO [fb]
NLOFTapprox [fb]
NNLONLO i [fb]
NNLOB proj [fb]
NNLOFTapprox [fb]
Mt unc. NNLOFTapprox
NNLOFTapprox/NLO
13 TeV
1147 +190:9:7%%
1220 +1110::96%%
1337 +45::14%%
1406 +02::58%%
1224 +03::92%%
4:6%
1:067
predictions is also presented. The uncertainties due to the qT subtraction and the numerical
evaluation of the virtual NLO contribution are both at the per mille level.
increase with respect to the exact NLO result of about 20% at 14 TeV. For this collider
energy, the increase within the NNLONLOi approach (which is computed based on the
Mhh distribution) is smaller, being about 18%. Finally, the NNLOFTapprox prediction is
the lowest one, with a 12% increase with respect to the NLO cross section at 14 TeV. For all
the considered approximations and collider energies the scale uncertainties are signi cantly
reduced when including the O( S4
) NNLO corrections. This reduction is largest for the
NNLOB proj and NNLOFTapprox approximations.6 For instance at 14 TeV, the total scale
uncertainty is reduced from about
13% at NLO to +2%
5% at NNLOFTapprox, i.e. by
about a factor of three. This reduction of the scale uncertainties is stronger as we increase
the collider energy, being close to a factor of ve at 100 TeV.
As is well known, scale uncertainties can only provide a lower limit on the true
perturbative uncertainties. In particular, from table 1 we see that the di erence between
the NNLO and NLO central predictions is always larger than the NNLO scale
uncertainties (although within the NLO uncertainty bands). In any case, the strong reduction of
scale uncertainties, together with the moderate impact of NNLO corrections, suggests a
signi cant improvement in the perturbative convergence as we move from NLO to NNLO.
It is also worth mentioning that the three approximations have a di erent behaviour
with ps. For instance at 100 TeV, the increase with respect to the NLO prediction for the
NNLOBproj and NNLONLOi approaches is 23% and 17%, respectively, values that are close
to the ones for 14 TeV (20% and 18%, respectively). By contrast, the NNLOFTapprox result
increases the NLO prediction by 7% at 100 TeV, i.e. the correction is smaller by almost a
factor of two than at 14 TeV (12%), which also means a larger separation with respect to the
other two NNLO approximations. The smaller size of the NNLO corrections in the FTapprox
at higher energies is also consistent with the observed reduction of scale uncertainties.
6The scale uncertainty of the NNLONLOi prediction is de ned as the relative uncertainty of the
HEFT result.
{ 8 {
As was mentioned already in section 2.2, the NNLOFTapprox result is expected to be
the most accurate one among the approximations studied in this work, and therefore it
is considered to be our best prediction. In order to estimate the remaining uncertainty
associated with nite top quark mass e ects at NNLO, we start by considering the accuracy
of the FTapprox approximation at NLO. At 14 TeV the NLO FTapprox result (see table 1)
overestimates the full NLO total cross section by only about 4%, or equivalently by about
11% of the pure O( S3) contribution. If we assume that FTapprox performs analogously at
one order higher, we obtain a
11% uncertainty on the O( S4
) contribution.7 Given that
the relative weight of the O( S4) contributions to the total NNLO cross section is de nitely
smaller than the weight of the O( S3) contributions to the NLO cross section, we obtain a
signi cantly smaller overall uncertainty, in this case of
1:2%. In order to be conservative,
we can increase this estimate by a factor of two. The relative di erence between the
FTapprox and the full NLO result slightly increases with the collider energy. However, at
the same time the relative size of the O( S4) correction decreases. The NNLO uncertainty
obtained with this procedure ranges from
2:3% at 13 TeV to
3:1% at 100 TeV.
We can repeat the above procedure to estimate the uncertainty of the NNLOB proj
approximation, which displays the largest di erences with respect to the NNLOFTapprox
result. Similarly to what we do for FTapprox, we can assign an uncertainty to the NNLOB proj
result by relying on the accuracy of the same approximation at NLO, and conservatively
to
36% at p
multiplying by a factor of two. The ensuing uncertainties range from
s = 100 TeV. We nd that the NNLOFTapprox prediction (always evaluated at
R =
F =
0) is fully contained in the NNLOB proj uncertainty band. Actually, there is
a large overlap between the two approximations, which includes in all the cases the central
value of the NNLOFTapprox, even when the conservative factor of two is not included. This
can be regarded as a nontrivial consistency check for our procedure. We may be tempted
to conclude our discussion by adopting the above procedure for the uncertainty estimate
14% at p
s = 13 TeV
of our NNLOFTapprox result.
However, we have already pointed out that, as ps increases, the di erence between the
5:2% at p
s = 13 TeV, and it becomes 9:2% at p
NNLOFTapprox and the other approximations increases. In particular, the di erence
between the NNLOFTapprox result and our \nexttobest" NNLO prediction, NNLONLO i, is
s = 100 TeV. The signi cant increase of
this di erence with the collider energy suggests us a more conservative approach. Our nal
estimate for the nite top quark mass uncertainty of our NNLOFTapprox result is de ned as
half the di erence between the NNLOFTapprox and the NNLONLO i approximations, and is
reported in table 1 for the di erent values of ps. At p
s = 13 and 14 TeV these
uncertainties are
2:6% and
2:7%, and thus very similar to the ones obtained with the method
discussed above. At p
s = 100 TeV, however, the uncertainty increases to
4:6%, which
appears to be more conservative than the
3:1% obtained with the previous procedure.
7We point out that in order to obtain the pure O( S4) corrections, we have subtracted the lower order
contributions computed with NNLO parton distributions and strong coupling. The corresponding numbers
are a few percent lower than the ones given in table 1 for the NLO results.
{ 9 {
6
5
)4
V
e
G
/
d
/
σd2
1
0
1.4
1.3
LO1.2
d
/
together with the NLO prediction, at 14 TeV (left) and 100 TeV (right). The lower panels show the
ratio with respect to the NLO prediction, and the lled areas indicate the NLO and NNLOFTapprox
scale uncertainties.
3.2
Di erential distributions
In this section we present predictions for di erential Higgs boson pair production at 14 TeV
and 100 TeV. We consider the following kinematical distributions: the invariant mass (Mhh,
gure 2) and rapidity (yhh, gure 3) of the Higgs boson pair, the transverse momenta of
the Higgs boson pair and the leading jet (pT;hh and pT;j1, gures 4 and 5), the transverse
momenta of the harder and the softer Higgs boson (pT;h1 and pT;h2, gures 6 and 7),
and the azimuthal separation between the two Higgs bosons (
hh, gure 8). For the
sake of clarity, we only show the scale uncertainty bands corresponding to the NLO and
NNLOFTapprox predictions.
We start our discussion from the invariantmass distribution of the Higgs boson pair,
reported in
gure 2. We observe that the NNLOBproj and NNLONLOi approximations
predict a similar shape, with very small corrections at threshold, an approximately
constant Kfactor for larger invariant masses, and only a small di erence in the normalization
between them, which increases in the 100 TeV case. The NNLOFTapprox, on the other
hand, presents a di erent shape, in particular with larger corrections for lower invariant
masses, a minimum in the size of the corrections close to the region where the maximum
of the distribution is located, and a slow increase towards the tail. The di erent behavior
of the NNLOFTapprox in the region close to threshold is more evident at 100 TeV, where
the increase is about 30% in the rst bin. Naively we could expect that if this region is
dominated by soft parton(s) recoiling against the Higgs bosons, the Born projection and
FTapprox should provide similar results. We have investigated the origin of this di erence,
and we
nd that in the region Mhh
2Mh the cross section is actually dominated by
events with relatively hard radiation recoiling against the Higgs boson pair (for example,
NNLOBproj
NNLONLOi
NNLOFTapprox
NLO
0
yhh
10
5
b
(f
hh
gether with the NLO prediction, at 14 TeV (left) and 100 TeV (right).
1.4
1.3
LO1.2
s = 100 TeV, the average transverse momentum of the Higgs boson pair in the rst
100 GeV at NLO). In this region the exact loop amplitudes behave
rather di erently as compared to the amplitudes evaluated in the HEFT: as the
production threshold is approached, they go to zero faster than in the massdependent case, thus
explaining the di erences we nd. Within the NNLOFTapprox, the corrections to the Mhh
spectrum range between 10% and 20% at 14 TeV. The scale uncertainty is substantially
reduced in the NNLOFTapprox, and this reduction is particularly strong for large
invariant masses. As observed at the inclusive level, the NNLOFTapprox corrections are smaller
at 100 TeV (except only for the rst bin) and the di erence with respect to the other
approximations is larger.
Next we move to the rapidity distribution of the Higgs boson pair, reported in gure 3.
The NNLO results are similar for all three approximations. This is not unexpected as the
shape of the rapidity distribution is mainly driven by the PDFs. Besides the obvious
difference in the normalization, the largest e ect in the shape of the NNLONLO i distribution
is observed in the central region, which is particularly evident in the 100 TeV case. Again
we observe a clear reduction of scale uncertainties over the whole range under study.
More signi cant di erences between the three approximations are obtained in the pT;hh
distribution, reported in gure 4. The NNLOBproj approximation predicts huge corrections
for large transverse momentum, the result being almost an order of magnitude larger than
the NLO prediction and the other approximations for pT;hh
500 GeV. This behavior is
hardly surprising since already at NLO the Bornprojected result deviates from the exact
NLO prediction in this way [17]. In fact, given that the pT;hh distribution is not de ned
at LO, the NNLOBproj corrections cannot inherit any information about the (full)
lowestorder prediction for this distribution. This is of course not the case for the other two
approximations, which in fact make an almost identical prediction at large pT;hh, with
)
NNLOBproj
NNLONLOi
NNLOFTapprox
NLO
d
/
3.0
2.5
LO2.0
toN1.5
d
/
3.0
2.5
LO2.0
toN1.5
d
/
σd0.10
d
/
σd0.10
0.05
0.01
3.0
2.5
LO2.0
toN1.5
Here jets are clustered with the antikT algorithm [59] with R = 0:4 and pT;j1 > 30 GeV
and j jj
4:4.
large corrections that can be well above 50%, and sizable uncertainties at the level of 30%{
40%, re ecting the NLOnature of this observable. At lower transverse momenta, however,
the NNLONLO i and NNLOFTapprox deviate from each other, and the latter approaches
the NNLOB proj prediction. Once again, the di erent behavior of these approximations is
more pronounced in the 100 TeV distribution, for which the central NNLONLO i curve lies
outside the NNLOFTapprox uncertainty band below pT;hh
200 GeV. Of course, in order to
obtain reliable results in the lowpT;hh region, the corresponding logarithmically enhanced
contributions need to be properly resummed to all orders in the strong coupling constant.
NNLOBproj
NNLONLOi
NNLOFTapprox
NLO
6
5
)
V
e
/G4
b
(f
ph1T
, 3
d
/
NNLOBproj
NNLONLOi
NNLOFTapprox
NLO
HJEP05(218)9
d
/
0.05
0.00
The transverse momentum distribution of the leading jet pT;j1, reported in gure 5,
has similar features as the pT;hh distribution. Again we observe the unphysical excess
predicted by the NNLOB proj approximation, which can be understood using the same
arguments as presented for the pT;hh distribution, and the agreement between NNLOB proj and
NNLOFTapprox at low pT;j1. The di erence between the NNLONLO i and NNLOFTapprox
results is more pronounced here, with the FTapprox predicting a softer spectrum for this
observable, and small corrections that are almost always contained in the NLO scale
uncertainty band.
The transversemomentum distributions of the harder and the softer Higgs boson are
reported in gures 6 and 7, respectively. As can be expected from the pT;hh spectrum, the
NNLOBproj result for pT;h1 features very large corrections as pT;h1 increases. The e ect,
however, is less severe than the one observed in pT;hh because the pT;h1 observable is
already well de ned at LO. The NNLONLOi curve is overall in good agreement with the
NNLOFTapprox prediction: it shows moderate corrections with respect to the NLO result
which increase as pT;h1 increases, while the scale uncertainties are about
15%. At very
small pT;h1 the higherorder corrections become perturbatively unstable as the available
phase space for the real radiation is severely restricted in this regime yielding large
logarithms that should be resummed in order to get a reliable prediction, see also the discussion
in section 3.4 of ref. [19]. For the transverse momentum of the softer Higgs boson, pT;h2,
the NNLO e ect is rather uniform in all three approximations, especially at 14 TeV. The
NNLOFTapprox predicts small corrections of order 10%, while the other two approximations
show larger corrections with a similar shape. In the tail of the distribution the scale
uncertainty at NNLO is larger than at NLO, most likely due to an accidentally small size of
the NLO scale variation (in fact, in this region the NLO corrections almost vanish).
Finally, the distribution in the azimuthal angle between the two Higgs bosons,
hh,
is shown in
gure 8. At LO we have
hh =
, due to the backtoback production of
)
V
V
b
b
NNLOBproj
NNLONLOi
NNLOFTapprox
NLO
d
/
1
0.1
3.0
2.5
LNO2.0
too1.5
i
tra1.0
0.5
0.00.0
pT,h2 (GeV)
d
/
2
0
d
/
100
10
3.0
2.5
LNO2.0
Azimuthal angular separation between the two Higgs bosons at 14 TeV (left) and
the two Higgs bosons at Born level. Real contributions allow
hh to be smaller than ,
and again we observe that the NNLOBproj approximation predicts larger corrections in
the region dominated by hard radiation compared to the other two results, which again
are in good agreement with each other in that region, whereas they start to deviate for
larger angles. For values of
hh close to , this observable receives large corrections from
softgluon emission, and the corresponding large logarithms should be resummed in order
to get a reliable prediction.
We conclude this section by adding a few comments on the niteMt uncertainties
at NNLO for the various di erential distributions. The analysis that was performed for
the total cross section cannot be easily extended to di erential distributions. On one
hand, any accidental agreement between the FTapprox and the full result at NLO in a
given phasespace region would likely lead to an underestimation of the top quark mass
e ects; on the other hand, the regions in which the NLO corrections are very small due
to cancellations between di erent contributions can present very large relative di erences
in the O( S3) contribution of the NLOFTapprox and NLO results, thus leading to arti cially
large uncertainties at NNLO. In addition, there are observables that are by de nition
reproduced in an exact way by the FTapprox at NLO (in our case pT;hh, pT;j1 and
hh),
and the uncertainty estimate procedure that we de ned for the inclusive case is therefore not
applicable. Despite these facts, and based on the performance of the FTapprox at NLO [17]
as well as on the observed di erences between our NNLO approximations, we can try
to assess the order of magnitude of the expected missing Mt e ects for the distributions
presented above.
In the Higgs boson pair invariantmass distribution, for values of Mhh below 500 GeV
the level of accuracy of the FTapprox at NLO is similar to the inclusive case, and therefore
the Mt uncertainty at NNLO is expected to be of a comparable size. In the tail of the
distribution, however, the quality of the FTapprox decreases (see gure 5 of ref. [17]), and
we thus expect the nite top quark mass e ects to be of O(10%) in this region.
The shape of the rapidity distribution of the Higgs boson pair is correctly described
by the FTapprox at NLO (see
gure 8 of ref. [17]), and the di erence to the full result is
only the overall normalization. Based on this, the estimated top quark mass uncertainty
for the NNLOFTapprox result is constant in the whole yhh range and of the same size as for
the inclusive cross section.
The transverse momentum of the harder Higgs boson is very well described at NLO by
the FTapprox (see gure 7 of ref. [17]), being always within the NLO scale uncertainty band.
This fact, together with the close agreement between the NNLOFTapprox and NNLONLO i
predictions, suggests that the missing top quark mass e ects at NNLO are probably of
moderate size.
The same holds true for the transversemomentum distribution of the
softer Higgs boson, except for the tail where at NLO the FTapprox overestimates the full
NLO corrections, which in fact almost vanish in this region.
The remaining distributions, which are either not de ned or trivial at LO, are by
definition reproduced in an exact way by the FTapprox at NLO, and this makes the estimate
of the missing top quark mass e ects at NNLO more di cult. In this case, a possible
approach can be to use the di erence between the NNLOFTapprox and NNLONLO i prediction
as an estimate of the uncertainty (as discussed before, the NNLOB proj prediction is not
expected to be reliable in the regions dominated by hard real radiation, where it largely
deviates from the other two approximations). This procedure would imply relatively low
top quark mass uncertainties for the pT;hh and
hh distributions, except for the low pT;hh
and the
hh
regions, typically below the size of the scale uncertainties, and larger
uncertainties for the leadingjet transverse momentum, for which the di erence between
the two approximations is larger.
In this work we considered Higgs boson pair production through gluon fusion in proton
collisions. We presented new QCD predictions for inclusive and di erential cross sections,
which include the full NLO contribution and also account for nite top quark mass e ects
at NNLO. Our best prediction, denoted NNLOFTapprox, retains the full top quark mass
dependence in the doublereal emission amplitudes, while the remaining realvirtual and
twoloop virtual HEFT amplitudes are treated via a suitable reweighting for the
corresponding subprocesses with a given
nalstate multiplicity. This approximation represents
the most advanced prediction available to date for this process.
The numerical results we obtained for the NNLOFTapprox are quantitatively di erent
from the results obtained in previous combinations. In particular, as far as the total
cross section is concerned, the corrections turn out to be smaller than previous estimates,
increasing the NLO result by about 12% at 13 TeV and 7% at 100 TeV. The reduction of
the scale uncertainties is signi cant, by about a factor of three for LHC energies. Given that
our NNLOFTapprox prediction includes top quark mass e ects in an approximated way, it is
important to assess the corresponding uncertainty. We carefully examined the performance
of our approximations at both the inclusive and di erential levels. The uncertainty on our
reference inclusive NNLOFTapprox prediction is estimated to be about
2:7% at 14 TeV,
increasing with the collider energy to reach
4:6% at 100 TeV.
Regarding di erential distributions, in most of the cases we can observe clear
qualitative di erences with respect to the binbybin reweighting procedure introduced in ref. [17],
in the shape and/or the normalization. For some of the distributions, however, speci cally
the tails of the pT;hh and pT;h1 spectra, both approximations are in very good agreement.
We discussed an estimate of the uncertainty associated with top quark mass e ects at
NNLO at the di erential level, and we found that in most of the cases its magnitude is
comparable to the size of the scale uncertainties, except for the tails of some distributions
where the uncertainty from missing Mt e ects can be dominant.
Acknowledgments
We thank Stefano Catani, Daniel de Florian, Ramona Grober, Andreas Maier, Stefano
Pozzorini and Michael Spira for valuable discussions and comments on the manuscript.
This research was supported in part by the Swiss National Science Foundation (SNF)
under contracts CRSII2141847, 200021156585, and by the Research Executive Agency
(REA) of the European Union under the Grant Agreement number PITN{GA{2012{316704
(HiggsTools).
Open Access.
This article is distributed under the terms of the Creative Commons
Attribution License (CCBY 4.0), which permits any use, distribution and reproduction in
any medium, provided the original author(s) and source are credited.
[2] ATLAS collaboration, Study of the double Higgs production channel H(! bb)H(!
) with
the ATLAS experiment at the HLLHC, ATLPHYSPUB2017001 (2017).
measurement of the Higgs selfcoupling at the LHC: theoretical status, JHEP 04 (2013) 151
[4] R. Frederix et al., Higgs pair production at the LHC with NLO and partonshower e ects,
[5] M. McCullough, An Indirect ModelDependent Probe of the Higgs SelfCoupling, Phys. Rev.
[6] M. Gorbahn and U. Haisch, Indirect probes of the trilinear Higgs coupling: gg ! h and
h !
, JHEP 10 (2016) 094 [arXiv:1607.03773] [INSPIRE].
[7] G. Degrassi, P.P. Giardino, F. Maltoni and D. Pagani, Probing the Higgs self coupling via
single Higgs production at the LHC, JHEP 12 (2016) 080 [arXiv:1607.04251] [INSPIRE].
[8] W. Bizon, M. Gorbahn, U. Haisch and G. Zanderighi, Constraints on the trilinear Higgs
coupling from vector boson fusion and associated Higgs production at the LHC, JHEP 07
(2017) 083 [arXiv:1610.05771] [INSPIRE].
[9] S. Di Vita, C. Grojean, G. Panico, M. Riembau and T. Vantalon, A global view on the Higgs
selfcoupling, JHEP 09 (2017) 069 [arXiv:1704.01953] [INSPIRE].
[10] F. Maltoni, D. Pagani, A. Shivaji and X. Zhao, Trilinear Higgs coupling determination via
singleHiggs di erential measurements at the LHC, Eur. Phys. J. C 77 (2017) 887
[arXiv:1709.08649] [INSPIRE].
[11] G. Degrassi, M. Fedele and P.P. Giardino, Constraints on the trilinear Higgs self coupling
from precision observables, JHEP 04 (2017) 155 [arXiv:1702.01737] [INSPIRE].
[12] G.D. Kribs, A. Maier, H. Rzehak, M. Spannowsky and P. Waite, Electroweak oblique
parameters as a probe of the trilinear Higgs boson selfinteraction, Phys. Rev. D 95 (2017)
093004 [arXiv:1702.07678] [INSPIRE].
[13] O.J.P. Eboli, G.C. Marques, S.F. Novaes and A.A. Natale, Twin Higgs boson production,
Phys. Lett. B 197 (1987) 269 [INSPIRE].
Phys. B 309 (1988) 282 [INSPIRE].
[14] E.W.N. Glover and J.J. van der Bij, Higgs boson pair production via gluon fusion, Nucl.
[15] T. Plehn, M. Spira and P.M. Zerwas, Pair production of neutral Higgs particles in
gluongluon collisions, Nucl. Phys. B 479 (1996) 46 [Erratum ibid. B 531 (1998) 655]
[hepph/9603205] [INSPIRE].
[16] S. Borowka et al., Higgs Boson Pair Production in Gluon Fusion at NexttoLeading Order
with Full TopQuark Mass Dependence, Phys. Rev. Lett. 117 (2016) 012001
[arXiv:1604.06447] [INSPIRE].
[17] S. Borowka et al., Full top quark mass dependence in Higgs boson pair production at NLO,
JHEP 10 (2016) 107 [arXiv:1608.04798] [INSPIRE].
[18] G. Ferrera and J. Pires, Transversemomentum resummation for Higgs boson pair production
at the LHC with topquark mass e ects, JHEP 02 (2017) 139 [arXiv:1609.01691] [INSPIRE].
boson pair production with full top quark mass dependence matched to parton showers, JHEP
08 (2017) 088 [arXiv:1703.09252] [INSPIRE].
[20] S. Jones and S. Kuttimalai, Parton Shower and NLOMatching uncertainties in Higgs Boson
Pair Production, JHEP 02 (2018) 176 [arXiv:1711.03319] [INSPIRE].
[21] S. Dawson, S. Dittmaier and M. Spira, Neutral Higgs boson pair production at hadron
colliders: QCD corrections, Phys. Rev. D 58 (1998) 115012 [hepph/9805244] [INSPIRE].
[22] F. Maltoni, E. Vryonidou and M. Zaro, Topquark mass e ects in double and triple Higgs
production in gluongluon fusion at NLO, JHEP 11 (2014) 079 [arXiv:1408.6542] [INSPIRE].
[23] D. de Florian and J. Mazzitelli, Twoloop virtual corrections to Higgs pair production, Phys.
Lett. B 724 (2013) 306 [arXiv:1305.5206] [INSPIRE].
[24] D. de Florian and J. Mazzitelli, Higgs Boson Pair Production at NexttoNexttoLeading
Order in QCD, Phys. Rev. Lett. 111 (2013) 201801 [arXiv:1309.6594] [INSPIRE].
[25] J. Grigo, K. Melnikov and M. Steinhauser, Virtual corrections to Higgs boson pair production
in the large top quark mass limit, Nucl. Phys. B 888 (2014) 17 [arXiv:1408.2422] [INSPIRE].
[26] D. de Florian et al., Di erential Higgs Boson Pair Production at NexttoNexttoLeading
Order in QCD, JHEP 09 (2016) 151 [arXiv:1606.09519] [INSPIRE].
[27] J. Grigo, J. Ho and M. Steinhauser, Higgs boson pair production: top quark mass e ects at
NLO and NNLO, Nucl. Phys. B 900 (2015) 412 [arXiv:1508.00909] [INSPIRE].
[28] R. Grober, A. Maier and T. Rauh, Reconstruction of topquark mass e ects in Higgs pair
production and other gluonfusion processes, JHEP 03 (2018) 020 [arXiv:1709.07799]
[INSPIRE].
[29] D.Y. Shao, C.S. Li, H.T. Li and J. Wang, Threshold resummation e ects in Higgs boson pair
production at the LHC, JHEP 07 (2013) 169 [arXiv:1301.1245] [INSPIRE].
[30] D. de Florian and J. Mazzitelli, Higgs pair production at nexttonexttoleading logarithmic
accuracy at the LHC, JHEP 09 (2015) 053 [arXiv:1505.07122] [INSPIRE].
[31] LHC Higgs Cross Section Working Group collaboration, D. de Florian et al.,
Handbook of LHC Higgs Cross Sections: 4. Deciphering the Nature of the Higgs Sector,
arXiv:1610.07922 [INSPIRE].
[32] M. Grazzini, S. Kallweit and M. Wiesemann, Fully di erential NNLO computations with
MATRIX, arXiv:1711.06631 [INSPIRE].
[33] S. Catani and M. Grazzini, An NNLO subtraction formalism in hadron collisions and its
application to Higgs boson production at the LHC, Phys. Rev. Lett. 98 (2007) 222002
[hepph/0703012] [INSPIRE].
[34] S. Catani and M.H. Seymour, The dipole formalism for the calculation of QCD jet
crosssections at nexttoleading order, Phys. Lett. B 378 (1996) 287 [hepph/9602277]
[INSPIRE].
[INSPIRE].
[35] S. Catani and M.H. Seymour, A general algorithm for calculating jet crosssections in NLO
QCD, Nucl. Phys. B 485 (1997) 291 [Erratum ibid. B 510 (1998) 503] [hepph/9605323]
http://openloops.hepforge.org.
1889 [arXiv:1111.2034] [INSPIRE].
[arXiv:1608.03846] [INSPIRE].
arXiv:1201.4330 [INSPIRE].
[36] S. Kallweit, Munich is the abbreviation of \MUltichaNnel Integrator at Swiss (CH)
precision"  an automated parton level NLO generator, in preparation.
transversemomentum resummation and hard factors at the NNLO, Nucl. Phys. B 881
(2014) 414 [arXiv:1311.1654] [INSPIRE].
HJEP05(218)9
[42] A. von Manteu el and C. Studerus, Reduze 2  Distributed Feynman Integral Reduction,
numerical evaluation of multiscale integrals beyond one loop, Comput. Phys. Commun. 196
(2015) 470 [arXiv:1502.06595] [INSPIRE].
[44] Z. Li, J. Wang, Q.S. Yan and X. Zhao, E cient numerical evaluation of Feynman integrals,
Chin. Phys. C 40 (2016) 033103 [arXiv:1508.02512] [INSPIRE].
[45] https://github.com/mppmu/hhgrid.
[46] F. Cascioli, P. Maierhofer and S. Pozzorini, Scattering Amplitudes with Open Loops, Phys.
Rev. Lett. 108 (2012) 111601 [arXiv:1111.5206] [INSPIRE].
[47] F. Buccioni, S. Pozzorini and M. Zoller, Onthe y reduction of open loops, Eur. Phys. J. C
78 (2018) 70 [arXiv:1710.11452] [INSPIRE].
[48] A. Denner, S. Dittmaier and L. Hofer, Collier: a fortranbased Complex OneLoop LIbrary in
Extended Regularizations, Comput. Phys. Commun. 212 (2017) 220 [arXiv:1604.06792]
[INSPIRE].
[INSPIRE].
[arXiv:1507.06937] [INSPIRE].
[49] S. Catani, D. de Florian, G. Ferrera and M. Grazzini, Vector boson production at hadron
colliders: transversemomentum resummation and leptonic decay, JHEP 12 (2015) 047
[50] G. Ossola, C.G. Papadopoulos and R. Pittau, CutTools: A program implementing the OPP
reduction method to compute oneloop amplitudes, JHEP 03 (2008) 042 [arXiv:0711.3596]
[51] M. Grazzini, S. Kallweit, D. Rathlev and M. Wiesemann, W Z production at the LHC:
ducial cross sections and distributions in NNLO QCD, JHEP 05 (2017) 139
[arXiv:1703.09065] [INSPIRE].
[52] C. Anastasiou et al., High precision determination of the gluon fusion Higgs boson
crosssection at the LHC, JHEP 05 (2016) 058 [arXiv:1602.00695] [INSPIRE].
[53] J. Butterworth et al., PDF4LHC recommendations for LHC Run II, J. Phys. G 43 (2016)
023001 [arXiv:1510.03865] [INSPIRE].
chromodynamics, Phys. Rev. D 93 (2016) 033006 [arXiv:1506.07443] [INSPIRE].
035 [arXiv:1401.0013] [INSPIRE].
HJEP05(218)9
Representation for Monte Carlo PDFs, Eur. Phys. J. C 75 (2015) 369 [arXiv:1505.06736]
[1] CMS collaboration, Projected performance of Higgs analyses at the HLLHC for ECFA 2016 , CMS PASFTR 16 002 [INSPIRE].
[3] J. Baglio , A. Djouadi , R. Grober, M.M. Muhlleitner, J. Quevillon and M. Spira , The [19] G. Heinrich , S.P. Jones , M. Kerner , G. Luisoni and E. Vryonidou , NLO predictions for Higgs [37] S. Catani , L. Cieri , D. de Florian, G. Ferrera and M. Grazzini , Universality of [38] F. Cascioli , J. Lindert , P. Maierh ofer and S. Pozzorini, The OpenLoops oneloop generator , [39] G. Cullen et al., Automated OneLoop Calculations with GoSam , Eur. Phys. J. C 72 ( 2012 ) [40] G. Cullen et al., GoSam2 . 0: a tool for automated oneloop calculations within the Standard Model and beyond , Eur. Phys. J. C 74 ( 2014 ) 3001 [arXiv: 1404 .7096] [INSPIRE]. [41] S.P. Jones , Automation of 2loop Amplitude Calculations, PoS(LL2016) 069 [43] S. Borowka , G. Heinrich, S.P. Jones , M. Kerner , J. Schlenk and T. Zirke , SecDec3 .0: [54] NNPDF collaboration , R.D. Ball et al., Parton distributions for the LHC Run II , JHEP 04 [56] L.A. HarlandLang , A.D. Martin , P. Motylinski and R.S. Thorne , Parton distributions in the LHC era: MMHT 2014 PDFs, Eur . Phys. J. C 75 ( 2015 ) 204 [arXiv: 1412 .3989] [INSPIRE]. [57] J. Gao and P. Nadolsky , A metaanalysis of parton distribution functions , JHEP 07 ( 2014 ) [58] S. Carrazza , S. Forte , Z. Kassabov , J.I. Latorre and J. Rojo , An Unbiased Hessian