Higgs boson production at hadron colliders at N3LO in QCD
JHE
Higgs boson production at hadron colliders at N3LO
Bernhard Mistlberger 0
Theory Division 0
0 CH1211 , Geneva 23 , Switzerland
We present the Higgs boson production cross section at Hadron colliders in the gluon fusion production mode through N3LO in perturbative QCD. Speci cally, we work in an e ective theory where the top quark is assumed to be in nitely heavy and all other quarks are considered to be massless. Our result is the rst exact formula for a partonic hadron collider cross section at N3LO in perturbative QCD. Furthermore, our result is an analytic computation of a hadron collider cross section involving elliptic integrals. We derive numerical predictions for the Higgs boson cross section at the LHC. Previously this result was approximated by an expansion of the cross section around the production threshold of the Higgs boson and we compare our ndings. Finally, we study the impact of our new result on the state of the art prediction for the Higgs boson cross section at the LHC.
Higgs Physics; Perturbative QCD

HJEP05(218)
1 Introduction SetUp 2 3
4
5
6
Calculation of coe cient functions
3.1
3.2
Computation of matrix elements
An elliptic integral in Higgs production
3.3 Iterated integrals 3.4
Analytic solution for partonic coe cient functions
Results
Conclusions Comparison with results based on a threshold expansion
A The elliptic integral
B
Various ingredients for Higgs boson production
With the discovery of the Higgs boson [1, 2] at the Large Hadron Collider (LHC) at CERN
we have entered a new era of particle physics phenomenology. With conclusive evidence for
the existence of the Higgs boson the Standard Model (SM) of particle physics has become
a self consistent theory. It explains the mechanism of electroweak symmetry breaking,
the origin of elementary particle masses and it allows to derive concise predictions to
energies far beyond current experimental reach. The SM is however a phenomenologically
incomplete theory and needs to be extended to obtain a satisfying description of all known
physics. Higgs boson measurements will provide a unique window to deepen our
understanding of fundamental interactions and to stringently test possible extensions of our
current knowledge.
The inclusive cross section for the production of a Higgs boson represents a prototypical
example of experimental and theoretical synergy. Its role in the extraction of fundamental
coupling constants is key and it provides an invaluable tool to discover potential deviations
from the SM. Experimentally it can be determined at the LHC to astounding precision. In
order to exploit the full potential of LHC phenomenology experimental precision must be
matched by equally precise theoretical prediction.
{ 1 {
The dominant production mechanism of a Higgs boson at the LHC is gluon fusion. In
comparison with other processes perturbative QCD corrections to the gluon fusion cross
section are large. In order to match current and future experimental precision this simple
fact demands computation of this process to very high order in perturbation theory.
Nexttoleading oder (NLO) [3{6] corrections to this process are available since more than two
decades. Corrections at nexttonextto leading order (NNLO) were computed in refs [7{9]
in an e ective theory (EFT) of QCD where the top quark is considered to have in nite
mass and all other quarks are massless [10{13]. In ref. [14] nexttonexttonextto leading
order (N3LO) corrections were computed in terms of an expansion around the production
threshold of the Higgs boson. This result marked the rst computation of a Hadron
collider observable to this order in perturbation theory. At the desired level of precision the
inclusion of many subdominant e ects, such as electroweak corrections and quark mass
e ects, in a prediction for the hadron collider observable are essential. Furthermore, a
critical assessment of all sources of uncertainties is required. A comprehensive study achieving
this goal was performed in ref. [15] and resulted in the state of the art prediction for LHC
measurements (see also refs. [16, 17]).
In this article we go beyond the previous approximation of the N3LO corrections to
the Higgs boson gluon fusion cross section in the EFT in terms of a threshold expansion
and compute it exactly.
Our calculation strongly relies on various ingredients already
entering the computation of ref. [14]. Speci cally, we require matrix elements integrated
over phase space for the production of the Higgs boson in association with up to three
partons and involving up to three loops. Purely virtual corrections were computed in
refs. [18, 19]. Contributions with one parton in the nal state and two loops were calculated
in refs. [20{24]. Matrix elements involving two
nal state partons and one loop (RRV)
or tree level matrix elements with three
nal state partons (RRR) were computed for
the purposes of refs. [14, 25{27] in terms of a threshold expansion. Furthermore, our
result relies on infrared subtraction terms formed out of convolutions [28, 29] of splitting
functions [30, 31] and an ultraviolet counter term based on lower loop amplitudes [32].
Both were already computed for the purpose of ref. [14].
In order to obtain our result we compute N3LO corrections to the partonic cross
section due to RRV and RRR matrix elements. The integration over the loop and
nal
state momenta involves complicated, highdimensional integrals. In order to facilitate
our computation we employ the framework of reverse unitarity [7, 33{36] that allows to
relate phase space integrals to cuts of loop integrals. Subsequently, we employ powerful
loop integration techniques to actually compute our phase space integrals. In particular,
we make use of integrationbypart identities [37, 38] in order to express our integrated
matrix elements in terms of a limited set of master integrals. We then proceed to compute
these master integrals using the framework of di erential equations [39{41]. The solution
of di erential equations requires the calculation of one boundary condition per master
integral. To obtain these boundary conditions we perform an expansion of every master
integral in terms of a threshold expansion. We then match the individual terms in the
expansion to socalled soft master integrals that were explicitly computed in refs. [42, 43].
{ 2 {
When solving di erential equations for RRR master integrals we encounter an
obstruction in the form of 2
2 systems of di erential equations that cannot be solved by
conventional means. The solution to these systems is given in terms of elliptic integrals.
The appearance of elliptic integrals in the computation of Feynman integrals is well
established [44{52] but still poses a considerable challenge. The majority of known analytic
results for Feynman integrals can be expressed in terms of iterated integrals referred to
as generalised poly logarithms [53]. A profound understanding of their analytic
properties [53{57] has been key to the success of higher order perturbation theory. The quest
for a similar understanding of iterated integrals involving elliptic functions is subject of
ongoing research and has already produced vast literature [58{78]. In particular, methods
to nd solutions for di erential equations, the understanding of functional relations among
such integrals and their analytic continuation from one kinematic regime to another are
of importance. In this article we present our pragmatic solution to the problem at hand
and produce a hadron collider cross section that involves the analytic treatment of an
Having obtained analytic results for all required matrix elements with di erent parton
multiplicity in the nal states we combine them to form the exact correction to the partonic
Higgs boson production cross section at N3LO. We then convolute our newly obtained result
and all required lower order cross section with parton distribution function in order to derive
physical predictions for hadron collider cross sections. We study in detail the deviations
of our results from the previous approximation of the N3LO cross section [14, 15]. Our
computation allows us to remove one source of uncertainty due to the truncation of the
threshold expansion from the state of the art prediction for the Higgs boson production
cross section [15] and we update the previous result.
This article is structured as follows: in section 2 we setup the notation for our
computation of the inclusive Higgs boson production cross section. Next, we discuss in detail the
analytic computation of the missing RRV and RRR coe cient functions in section 3. We
outline the general computational framework in section 3.1. We discuss the treatment of
elliptic integrals found when solving di erential equations in section 3.2. In section 3.3 we
introduce a class of iterated integrals that serve as the main building blocks for our nal
result. Next, we describe the structure of our analytic results in section 3.4. In section 4
we present numerical results for the EFT Higgs boson cross section through N3LO in QCD
perturbation theory. We compare our new results to previous predictions obtained with a
threshold expansion in section 5. Finally, we draw our conclusions in section 6.
2
SetUp
Higgs boson. In this article we consider scattering processes of two protons that produce at least a
Proton(P1) + Proton(P2) ! H(ph) + X;
(2.1)
P1 and P2 are the momenta of the colliding protons and ph the momentum of the Higgs
{ 3 {
boson. The master formula for the inclusive Higgs boson production cross section is given by
z ^ij (z; m2h):
Here, we employed the parton model and factorization of long and short range interactions
into parton distribution functions fi(x) and partonic cross sections. The momenta of the
colliding partons are related to the proton momenta by p1 = x1P1 and p2 = x2P2 = x1z P2.
We de ne
=
z =
m2h
S
m2h
s
;
;
S = (P1 + P2)2:
s = (p1 + p2)2:
The sum over i and j ranges over all contributing partons. Furthermore, we de ne the
variable z = 1
z. The partonic Higgs cross section is given by ^ij (z; m2h).
In this article we compute the partonic cross section through N3LO in perturbative
QCD in an e ective theory where the top quark is in nitely heavy and has been integrated
out [10{13]. In this theory the Higgs boson is coupled directly to gluons via an e ective
operator of dimension ve [79{82],
Le = LSM;5
4
SM Lagrangian with nf = 5 massless quark
avours. The Wilson coe cient C0 is
obtained by matching the e ective theory to the full SM in the limit where the top quark is
in nitely heavy.
Within the e ective theory, we can write the partonic cross section as
1
z ^ij (z; m2h) = (C0)2 ^0 ij (z) = (C0)2 ^0
1
X
n=0
0 n
S
i(jn)(z):
Dividing out the Born cross section,
we can write the bare partonic coe cient functions as,
The initial state dependent prefactors Nij are given by
Ngg =
4(1
Ngq = Nqg =
Nqq = Nqq = NqQ =
1)2
1
1
)(nc2
4nc2
;
:
1)nc
;
1
)2(nc2
4(1
{ 4 {
(2.2)
(2.3)
(2.4)
(2.5)
(2.6)
(2.7)
(2.8)
Here, g, q, q and Q indicate that the initial state parton is a gluon, quark, antiquark
or quark of di erent avour than q respectively. d H+m is the phase space measure for
the production of a Higgs boson and m partons and is explained in more detail below.
Sn in the coupling constant expansion of the modulus
squared of all amplitudes for partons i and j producing a
nal state Higgs boson and
m partons summed over polarizations and colors. To compute the nth order partonic
coe cient functions we require all combinations lloop matrix elements with m external
particles such that m + l = n.
The occurring loop amplitudes are plagued by ultraviolet divergencies which we
regulate using dimensional regularisation and work in d = 4
2 spacetime dimensions. We
matrix elements with xed parton multiplicity in the nal state are separately infrared
divergent. These infrared divergences are canceled by summing over all contributing squared
matrix elements and performing a suitable rede nition of the parton distribution functions.
The resulting partonic cross section is free of divergencies and we refer to the
corresponding partonic coe cient function as ij (z). Various de nitions regarding renormalisation
and mass factorisation can be found in appendix B. The cross section, eq. (2.2), can be
written in terms of nite partonic coe cient functions and physical parton distribution
functions fiR as
(2.9)
(2.10)
The partonic coe cient functions can be split into two contributions
i(jn)(z) = ij
(n); SV(z) + i(jn); reg.(z):
The term ij
(n); SV(z) is comprised of distributions that act on parton distribution functions.
The superscript SV signi es that this term represents socalled softvirtual contributions
that arise from kinematic con gurations where any parton produced in conjunction with
Higgs boson is soft. The coe cient ij
(
3
); SV(z) was computed in ref. [26] and con rmed by
ij
ref. [27]. The coe cient functions i(j3); reg.(z) represent the socalled regular contributions.
Their functional form was approximated with a power series in 1
z in refs. [14, 25, 36].
The main result of this article is the complete computation of the coe cient functions
(
3
); reg.(z). We supply this result in a machine readable format in an ancillary le together
with the arXiv submission of this article.
3
Calculation of coe
cient functions
In order to obtain the partonic coe cient functions i(j3)(z) we require contributions arising
from matrix elements with up to three loops (l
3) and up to three partons (m
3) in
the
nal state such that 3 = l + m. The purely virtual matrix elements were computed
in refs. [18, 19].
Matrix elements with two loops and one emission were computed in
refs. [20{24]. Matrix elements with two real emissions and one loop (RRV) and three
{ 5 {
d H+m =
where
real emissions (RRR) are sofar publicly only available in terms of the rst two expansion
terms in the expansion around the production threshold of the Higgs bosons [42, 43]. In
this article we complete the computation of the N3LO coe cient functions. We start by
outlining the strategy involved in this computation. Next, we explain the treatment of
an ellitpic integral that is part of the RRR coe cient functions. We introduce a class of
iterated integrals that serve as building blocks of our partonic coe cient functions. Finally,
we obtain the N3LO coe cient functions and describe their structure.
3.1
Computation of matrix elements
In order to obtain RRV and RRR coe cient functions we start by generating all required
Feynman diagrams with QGRAF [83]. Next, we perform spinor and colour algebra in a
private c++ code based on GiNaC [84]. With this we obtain the loop and phase space
integrand for our partonic coe cient functions.
Next, we want to perform the inclusive integral of our integrands over all loop momenta
and nal state parton momenta. The phase space measure for producing a Higgs boson
and m partons is given by
(d2dp)hd (2 ) +(p2h
m2h)(2 )d d
p1 + p2 + ph +
m+2
X pi
i=3
i=3
! m+2 (d2dp)id (2 ) +(pi2);
Y
They satisfy the condition
We can now apply integrationbypart (IBP) identities [37, 38] on our combined loop and
phasespace integrands. A private c++ implementation of the Laporta algorithm [85]
allows us to express our partonic coe cient functions in terms of a limited set of master
integrals. To compute these master integrals we work with the method of di erential
{ 6 {
+(p2
m2) = ( p
0
m) (p2
m2):
All nal state momenta are chosen ingoing such that the energy component in the above
equation appears with a minus sign. In order to perform the loop and phase space
integration we employ the framework of reverse unitarity [7, 33{36] that allows to treat phase space
and loop integrals on equal footing. In particular, we represent the onshell constraints in
terms of cut propagators.
The subscript c serves as a reminder that this propagator is cut. Cut propagators can be
di erentiated just like normal propagators.
+(p2
m2) !
p2
1
1
m2 c
d
1
1
a
dx f (x) c
=
a
f (x) c
a+1 df (x)
dx
:
f (x) c
a
f (x)b = <
:
8 h 1 ia b
f(x) c
0
; if a > b
; if b
a
:
(3.1)
(3.2)
(3.3)
(3.4)
(3.5)
(3.6)
(3.7)
equations [39{41]. This method allows to derive a system of partial di erential equations
for a vector of our master integrals I~(z) of the form
Here, I~(z) is a vector of n master integrals and A(z; ) is a n
n matrix with ratios
of polynomials in z and
as entries. In order to have a complete system of di erential
equations we de ne 550 and 362 master integrals for RRR and RRV respectively.
The commonly used strategy to solve such di erential equations is to nd a n
n
A0(z; )I~0(z):
T:
Here, A0(z; ) is holomorphic in
as
! 0. Having obtained such a form the solution for
our master integrals can be easily expressed in terms of a Laurent series in
by
"
I~0(z) = I +
Z z
dz0A0(z0; ) + 2
Z z
Z z0
dz0
dz00A0(z0; )A0(z00; ) + : : : I~0 :
0
(3.8)
#
Here, I~00 represents a vector of boundary conditions that has to be determined by other
means. For the RRV and RRR master integrals such a boundary condition is easily obtained
by matching the full solution obtained in eq. (3.8) to an expansion of the required integrals
I~(z) around the point z = 1 that can be performed by means presented in refs. [42, 43].
The art in solving di erential equations rests in
nding an adequate transformation
matrix T . For certain di erential equations in a single parameter an algorithmic solution
exists [86{89] and was nicely formulated in ref. [89]. This method applies if a transformation
matrix can be found that is comprised of ratios of polynomials in the parameters z and
. For a large subset of integrals in our vector of master integrals I~ such transformations
can be found and we rely on a private implementation of the algorithm outlined in ref. [89]
to do so.
For another large class of master integrals it is necessary to
nd a transformation
matrix that contains square roots of polynomials of our parameter z. For these cases we
can
nd the desired transformation by
nding suitable algebraic variable transformations
that rationalises the square roots involved. Once the roots are rationalised we can again
employ the aforementioned algorithm. We point out that this procedure is not particularly
algorithmic but leads to a desired solution fairly easily.
We encounter a further obstruction when solving di erential equations for the system
of RRR master integrals. This obstruction involves the presence of elliptic integrals and
we elaborate on our solution in the following section.
{ 7 {
p1
(3.11)
! 0.
(a)
(b)
N3LO. The computation of these integrals involves elliptic integrals. Solid lines represent Feynman
propagators. Solid lines crossed by the dashed line correspond to cutpropagators. The doubled
line represents the onshell constraint of the Higgs boson.
HJEP05(218)
3.2
An elliptic integral in Higgs production
When solving di erential equations for master integrals contributing to the triple real
coe cient functions of Higgs boson production at N3LO we encounter two coupled 4
4
systems of di erential equations that we could not decouple order by order in the
dimensional regulator by conventional means. In this section we discuss the di erential equations
in question and present our solution.
In
gure 1 we display two scalar phase space integrals. Let us choose four master
integrals with the same propagators as the scalar integral in gure 1(b).
Ei =
Z
These four integrals satisfy a system of di erential equations of the form
The vector ~y(z) represents the inhomogeneous part of the di erential equation. The matrix
A1(z; ) in the homogeneous part of the di erential equation is holomorphic in
as
The homogeneous part of the di erential equation that does not decouple as
! 0 is given
{ 8 {
by the matrix
A0(z) = BBB
0
11 2z 3 z
z2 11z 1 0 0 z2 11z 1
0
0
0
1
C
C :
C
A
0
0
1
z
= 0 two of the master integrals decouple and we are left with a 2
2
system for the homogeneous solution of the di erential equation.
E0 !
4
E0
1
= AT :
= 0:
(3.12)
(3.13)
We show in appendix A that the functions tij (z) can be written in terms of complete
elliptic integrals and prefactors. However, this solution is quite unwieldy and we choose
another approach here. For all practical purposes it is su cient to simply de ne the
functions tij (z) to be the solution to the di erential equation eq. (3.14). The homogeneous
di erential equations for the master integrals E10 and E40, de ned by
E1 = t22E10 + t21E40;
E4 = t11E40 + t12E10;
are decoupled as we send
! 0. The inhomogeneity can then be decoupled order by
order by in
by standard techniques. A general solution for the di erential equations can
subsequently be found as illustrated by eq. (3.8).
The second set of master integrals that have the same propagators as the scalar
integral depicted in
gure 1(a) can be chosen in such a way that the homogeneous part of
their di erential equations takes identically the same form as the one already discussed.
Therefor we can apply the same transformation matrix to decouple the system order by
order in . With this we found a transformation matrix T that allows us to express the
di erential equations for all master integrals required for RRV and RRR contributions to
Higgs production at N3LO in the desired form, eq. (3.7).
In order to derive numerical results for the functions tij we can solve the di erential
equations eq. (3.14) in terms of a generalised power series ansatz using the Frobenius
method. Consider for example an ansatz for the solution of the system of di erential
{ 9 {
equations as an expansion around z = 0 and z = 1.
tij (z) =
tij (z) =
1
X znbi(jn):
n=0
1
n=0
X znci(jn) + log(z) X di(jn)zn:
1
n=0
(3.16)
(3.18)
(3.19)
We derived the required structure of our ansatz by regarding the asymptotic solution of
the di erential equations around the considered expansion points.
Here, ti0j and ti1j are some numerical boundary constants.
Inserting the ansaetze into the system of di erential equations we nd the following
b(1n1+2) =
By comparison to the asymptotic solution given in eq. (3.17) we nd all starting conditions
for the solution to the recurrence relations. Speci cally, we nd the conditions bi(jn) = ci(jn) =
di(jn) = 0 if n < 0 and d(201) = 0 and d(101) = c(201). Furthermore, the general solution for t22 is
identical to the solution for t21 and the one for t12 is identical to the solution for t11 up to
the choice of boundary constants.
Any choice of boundary conditions will lead to a transformation matrix that satis es
the di erential equations eq. (3.14). The only restriction we impose is that the
transformation has to be invertible, i.e. that det(TE ) 6= 0. In accordance with this criterium we make
z =
z =
1
2
1
2
p
Consequently, the power series of the functions tij (z) around the point z = 1 has a radius
of convergence r1 = 1. Similarly, the power series around the point z = 0 has radius of
conp
vergence r0 = j 12 11
on the interval z 2 (0; j 12 11
5 5 j. The domains of convergence for the two power series overlap
5 5 j). In order to determine the boundary constants ti0j
in terms of the ti1j we rst compute the truncated power series around both limits under
consideration for each tij (z). Next, we evaluate both series for each tij (z) at a point within
the interval z 2 (0; j 12 11
5 5 j). Equating the results allows us to establish a relation
among the constants ti0j and ti1j up to a small, numerical remainder. The remainder can be
systematically improved upon by increasing the truncation order of the power series.
Let us brie y introduce a simple method of estimating the size of the remainder of the
truncated series. Suppose a function f (x) is given by the convergent series
the following choice for the asymptotic solution of the di erential equation: t111 = t122 = 0
and t112 = t121 = 1. We nd that with this choice the determinant of the transformation
matrix is given to all orders in z by
t11t22
Fixing the asymptotic behaviour of the functions tij (z) in one limit automatically
determines their behaviour at any other point. Computing the explicit values for ti0j explicitly
given the choice we made for ti1j is however a nontrivial task. At this point it is useful to
re ect on the practical aim of our computation. We desire a solution that is numerically
su ciently precise to determine the complete N3LO coe cient functions for values of z in
the interval [0; 1] as required for cross section predictions. In this light our solution for the
tij (z) should allow for the desired precision and should be improvable if necessary. This
can be achieved by computing an approximation based on a truncated power series.
The regular singular points of our 2
2 system of di erential equations (3.14) are located at the following values.
jai+mj = rimjaij:
If we truncate the series before order N its remainder would be given by
Suppose that asymptotically the ratio of to consecutive coe cients remains constant.
(3.21)
(3.22)
(3.23)
(3.24)
Under this assumption we can estimate the modulus of the remainder to be bounded by
jR(f (x); x; N )j
jaN jx
1
i=0
N X(rN x)i =
aN xN
1
rN x
= Rest(f (x); x; N ):
(3.25)
Note, that the series converges for jrN xj < 1.
In order to obtain su ciently high precision for our coe cient functions we perform
an expansion of the functions tij around the expansion points z = 0, z = 1 and z = 12
.
For each expansion we compute several hundred terms and match the boundary conditions
within the overlaps of the domains of convergence. Estimating the remainder of the power
series expansion at our matching points suggests that we can easily determine the boundary
values with a relative accuracy of 10 42 or better if needed. In addition to estimating the
remainder as described above we evaluate the di erent power series for the same tij for
several points in the intervals where all series converge and only observe relative deviations
at levels smaller than 10 42.
In order to further study the convergence of our power series approximation we may
regard the asymptotic behaviour of the recurrence relations given in eq. (3.18) and eq. (3.19)
as n ! 1.
We see that asymptotically b(1n1) approaches a constant and c(1n1) and d(1n1) tend towards zero.
For the other coe cients we nd the asymptotic solutions
b(2n1) =
c(2n1) =
d(2n1) =
9
22
11
2
11
2
are smaller than one. The number
5 5 is larger than one. From this we again
p
22
5 5 and
p
z < 1=
in eq. (3.25).
draw the conclusion that the power series around the expansion point z = 1 is convergent
everywhere within the unit interval.
The power series around z = 0 is convergent if
of the procedure to estimate the remainder of the power series truncated at order N de ned
5 5 . This asymptotic analysis also supports the validity
Z z
0
In this section we brie y introduce a class of iterated integrals [90] that is particularly
convenient to express the solution of di erential equations as in eq. (3.8). We de ne
HJEP05(218)
J (!~; z) = J (!n(z); : : : ; !1(z); z) =
dz0!n(z0)J (!n 1(z0); : : : ; !1(z0); z0);
(3.28)
with J (z) = 1: We refer to one !i(z) as a letter and to an ordered set of letters, f!n(z); : : : ;
!1(z)g that de nes an iterated integral as a word.
Many well known classes of iterated integrals, such as harmonic poly logarithms
(HPLs) [91] or generalised poly logarithms (GPLs) [53], that are widely used in
particle physics, are subclasses of this type of iterated integrals. For example the GPLs are
The presence of the elliptic integrals tij (z) in the solution of our di erential equations does
not allow for a solution purely in terms of GPLs. For this reason it becomes necessary to
de ne an extension of GPLs in this article. Already several generalisations of GPLs to
accommodate elliptic functions exist in the literature (see for example [62, 73, 74, 92{94]). In
the following we review several properties of iterated integrals (see for example [55{57, 95]).
Iterated integrals form a socalled shu e algebra.
J (!n(z); : : : ; !1(z); z) J (!n+m(z); : : : ; !n+1(z); z) =
J (! (n+m)(z); : : : ; ! (1)(z); z);
X
2 (n;m)
where
(n; m) denotes the set of all shu es of n + m elements, i.e., the subset of the
symmetric group Sn+m de ned by
(n; m) = f 2 Sn+mj
1(n) < : : : <
1(1) and
1(n + m) < : : : <
(3.30)
1(n + 1)g :
(3.31)
For example, consider the product of two iterated integrals with two integrations each.
J (a; b; z)J (c; d; z) = J (a; b; c; d; z) + J (a; c; b; d; z) + J (a; c; d; b; z)
+J (c; a; b; d; z) + J (c; a; d; b; z) + J (c; d; a; b; z):
(3.32)
Here the letters a, b, c and d may be generic functions of z. Special care needs to be taken if the integrand of our iterated integrals diverges at the value of the lower integration bound. In this article we only consider simple poles of the { 13 {
integrand at the end points since they simply are the only type of divergence that appears
in the computation we are interested in. Speci cally, we de ne the case where all letters
of a word of lenght n are given by !(z) = z1 then
1
z
1
z
J
; : : : ; ; z
=
1
n!
logn(z):
If the letter z1 appears in the rightmost entry of the word of an iterated integral we de ne
it in a way that is consistent with the shu e algebra. Consider the shu e relation
; z J (!n(z); : : : ; !1(z); z) = J
!n(z); : : : ; !1(z); ; z
+ J
!n(z); : : : ; ; !1(z); z
+ : : : :
1
z
1
z
Z 1 z
0
Z z
1
1
z
f (z)
f (0)
z
; z
1
z
(3.33)
(3.34)
(3.35)
Here, the ellipsis indicates all other terms arising from the shu e product. Assuming that
all !i(z) in the above equation are holomorphic as z ! 0 the only iterated integral with
an endpoint divergence is the rst on the right hand side of the equation. We de ne our
iterated integrals to be regulated in such cases such that the above equation holds true.
Solving for the iterated integral in question we nd
J !n(z); : : : ; !1(z); ; z = log(z)J (!n(z); : : : ; !1(z); z)
J !n(z); : : : ; ; !1(z); z + : : :
If the rightmost letter is divergent as z ! 0 but has the form f(zz) , with f(z) being
holomorphic around z = 0, then we may regularise our function by writing it as
f (z)
z
J
!n(z); : : : ; !1(z);
; z
The last line of the above equation is then regulated as discussed above. If several
rightmost letters have poles at the lower end point of the integration we simply iterate the above
procedure.
We want to be able to rewrite an elliptic integral with argument z in terms of iterated
integrals with argument z = 1
by regarding a transformation from z to z.
z or w = 12
z. Let us illustrate how this can be achieved
J (!n(z); : : : ; !1(z); z) =
dz0!n(z0)J (!n 1(z0); : : : ; !1(z0); z)
dz0!n(1
z0)J (!n 1(1
z0); : : : ; !1(1
dz0!n(1
z0)J (!n 1(1
z0); : : : ; !1(1
z0); 1
z0); 1
z0)
z0)
+
dz0!n(1
z0)J (!n 1(1
z0); : : : ; !1(1
z0); 1
z0): (3.37)
The last line in the above equation is a numerical constant. In order to write the integral
in the penultimate line in terms of an iterated integral with upper integration bound z we
have to rst rewrite the iterated integral in the integrand with an upper integration bound
z0. To do this we simply apply the above equation iteratively to the integrand. Notice,
that the above procedure may be ill de ned if the integrand we are considering is divergent
at any of the end points. This case is easily avoided by shu e regulating both end points
prior to applying eq. (3.37). Let us demonstrate this step with a well known example.
Consider the iterated integral
J
1
;
z 1
1
z
; z
= J
=
1
z
; z J
log(z) log(1
1
z
z)
; z
J
J
1
1
1
1
1
; ; z
z z
1
; ; z
z z
In the above equation we employed a shu e identity such that right most letter of the
function is regular at the new lower integration point z = 1 and that the left most letter is
regular at the new end point z = 0. We now may write
1
1
1
z0
z0
; z0
; z0
2
6
:
Z z
0
Z z
0
=
=
= J
1
dz0
1
1
z0
dz0 J
1
;
z 1
1
z0
J
z
1
; z
=
dz0
Z z
0
0
log(z0)
1
1
z0
dz0 J
z0
1
log(z) log(1
z) +
1
z0
; z0
2
6
:
Combining the the results of eq. (3.38) and eq. (3.39) we nd the famous diLogarithm
identity.
J
1
;
z 1
1
z
; z
=
J
1
;
z 1
1
z
; z
In this example it was possible to determine the integration constant to be 62
analytically. If this is not possible the constant can also be determined numerically with
nite
precision by simply evaluating the function under consideration before and after variable
transformation numerically for any value of z.
The iterated integral representation of eq. (3.28) allows to easily compute truncated
power series expansions for the iterated integrals. For example
J
1
1
z
; z
=
Z z dz0
z
=
By proceeding iteratively we can easily compute the power series in z for any iterated
integral to arbitrary power.
In order to obtain compact expressions for our analytic results it is of importance to be
able to derive functional relations among our iterated integrals. One of the big advantages
of GPLs is that their functional relations are well studied (see for example [53, 55{57, 96]).
The case of generic iterated integrals is not understood at the same level. In ref. [62] it was
(3.38)
(3.39)
(3.40)
(3.41)
outlined how relations among iterated integrals involving elliptic functions can be found
using IBP identities. Here, we proceed di erently.
First, note that our nal analytic result will be a linear combination of iterated integrals
functions.
relation
c1J
7
60
we nd
can be satis ed for some ci 6= 0 for arbitrary values of z. The coe cients ai(z) and
corresponding iterated integrals J (!~i; z) are understood to be identical to those appearing
in our nal result. In order to determine the unknown coe cients ci we expand eq. (3.43)
in z. Every coe cient of every power in z has to vanish separately in order for the equation
to be satis ed. This allows us to build a system of equations that is large enough to solve
for the unknown coe cients ci. If we nd a certain linear combination of iterated integrals
and coe cients that cannot be constrained with this procedure we found a relation of
Let us illustrate the procedure with a trivial example. Consider the simple shu e
1
;
1
z 1 + z
1
; z
+ c2J
1
1 + z 1
;
1
z
; z
+ c3J
1
1
z
; z J
1
1 + z
; z
= 0; (3.44)
and let us pretend we do not know already know the coe cients ci. After expanding in z
and prefactors ai(z),
then the equation
If there are relations among di erent iterated integrals appearing in this linear combination
X ai(z)J (!~i; z):
i
X ciai(z)J (!~i; z) = 0;
i
; ci 2 Q;
(3.42)
(3.43)
We can now create a system of equations by regarding each coe cient in z separately.
Technically, we want to nd the kernel of the system of equations. We nd that the kernel
for our example is spanned by the vector fc1; c1; c1gT . This means we found the shu e
identity
1
1
J
;
1
z z + 1
; z
+ J
1
z + 1 1
;
1
z
; z
J
1
1
z
1
z + 1
; z J
; z
= 0:
(3.47)
Of course this procedure only guarantees that the sofound relations are satis ed up to the
order in z at which we truncate our power series. However, we may convince ourselves that
the relations are correct by computing as many higher order terms as are to our liking. A
more involved example of such an identity is given by
J t11; 1
Analytic solution for partonic coe cient functions
In the previous sections we described how we derive di erential equations for all
master integrals required for RRR and RRV partonic coe cient functions. Furthermore, we
outlined how we nd a suitable transformation matrix that transforms the di erential
equations into the form of eq. (3.7). Once, this form is obtained the solution to the di erential
equations can be conveniently written as in equation eq. (3.8). Iterated integrals as given in
eq. (3.28) are particularly suited to represent this solution. Once we calculated all master
integrals and computed all boundary conditions we simply insert the master integrals into
our IBP reduced matrix elements and obtain the desired result for the partonic coe cient
functions. In this section we describe the structure of our nal result for the partonic
coe cient functions.
The set of all letters, the socalled alphabet, that appear in the iterated integrals that
constitute the Higgs boson cross section at N3LO is given by
1;
1
1
1
1
1
z ; z ; z + 1 ; pz ; p
4
1
p ;
z z 1
1
z ; p p
z z + 4
; p
t11; t12; t21; t22; 1
t11 ; t11 ;
t11 ;
t12 ; t12 ;
t12 ; t21 ; t22 :
z z z + 1 1
z z z + 1 z
z
p
z
p
z
z + 4
1
; p4z + 1
;
p4z + 1
z
;
(3.49)
Note, that the alphabet required to describe all our master integrals individually contains
additional letters that drop out in the nal expression.
The partonic coe cient functions are comprised of iterated integrals with up to ve
letters. Typically we nd that there are several thousand di erent iterated integrals in each
partonic coe cient function. Applying the procedure outlined in the previous section to
nd functional identities among these integrals we nd that we can express them in terms
of only 365 di erent iterated integrals that cannot be rewritten as GPLs in a straight
forward fashion. Out of those 188 have letters containing elliptic integrals tij . For the
remaining ones a representation in terms of GPLs may exist.
Having derived moderately compact expressions for our coe cient functions we want to
nd a method to evaluate them numerically. The conceptually simplest way to evaluate the
iterated integrals is to perform every integral numerically. The fact that all our integrals are
real valued and are nite renders this approach straight forward. Integrating 5 dimensional
integrals numerically is however not particularly fast if a certain level of precision is desired.
As an alternative, we want to represent the entire partonic coe cient functions in terms
of power series expansions.
Let us rst investigate for which values of z we can perform a convergent series
expansion. In order to extract this information we regard all singularities and branch points
that occur in our alphabet and the algebraic factors of our coe cient functions. We nd
that they are located at values of z of
1
2
p
Here, we included the regular singular points of the di erential equations of the elliptic
sector, eq. (3.13). In order to evaluate our functions to high precision within the physical
interval, z 2 [0; 1], we decide to perform a power series expansion around the points z = 1,
z = 12 and z = 0. The associated radii of convergence are then r1 = 1, r 1 = 12 and
2
To obtain a series expansion around our three di erent expansion points we perform
an expansion of all iterated integrals as outlined in the previous section. As the default
upper bound for our iterated integrals is the parameter z the expansion around the point
z = 1 can be carried out simply by expanding the iterated integrals at the integrand
level and integrating subsequently as demonstrated in eq. (3.41). In order to obtain an
expansion around z = 0 and z = 12 we rst reexpress our iterated integrals in terms of
iterated integrals with upper integration bound z and 1=2
z respectively. As outlined in
section 3.3 this procedure requires us to determine certain integration constants which we
obtain numerically by matching series expansions around di erent expansion points. To
ensure that the numerical error introduced by truncating series expansions is su ciently
small we estimate it as explained in eq. (3.25).
We expand the coe cient of every iterated integral in the partonic coe cient function
separately around each of our three expansion points and combine the result with the
expansion for the iterated integrals. In order to obtain numerical values for the coe cient
z = 1 at O((1
z)50), the expansion around z = 12 at O (
1
2
functions within the unit interval of z we evaluate the expansion around z = 1 in the
interval z 2 [0:75; 1], the expansion around z = 12 in the interval z 2 113 ; 0:75 and the
expansion around z = 0 in the interval z 2 0; 113 . We truncate the expansion around
z)200 and the expansion
around z = 0 at O z100 . Using the estimator introduced in eq. (3.25) we nd that this
approximates the coe cient functions at any point in the unit interval to a relative
numerical precision of 10 10 or better. This is supported by evaluating the di erent expansions
for several points within the overlaps of their respective domains of convergence and
calculating their di erence. The numerical precision may of course be improved arbitrarily by
simply including more terms in the respective series expansions.
1500
500
10000.
z
function of the parameter z. The gg coe cient function was rescaled uniformly by a factor of 10 2.
4
In the previous sections we calculated analytic results for the partonic coe cient functions
i(j3)(z). Our analytic results agree with the power series around z = 1 for the same functions
obtained in refs. [14, 15]. The leading behaviour of the coe cient functions as z ! 0 was
correctly predicted in ref. [97]. The coe cient function
q(3Q) was calculated already in
ref. [98] and agrees with our result. We derived a representation of the coe cient functions
in terms of power series expansions that is particularly useful for numerical evaluation.
In this section we present numerical results for the Higgs boson production cross section
through N3LO.
Let us start by regarding the functional dependence of our coe cient functions. In
gure 2 we display the shape of the regular coe cient functions for each distinct partonic
initial state. The quark  gluon and and gluongluon initial state coe cient functions
behave as
log5(1
z) as we approach the value z = 1. The coe cient functions with
two quarks in the initial state are tending towards zero in this limit. The limit z ! 0 is
characterised by a power divergence and all coe cient functions behave as
log5(z) .
z
In order to derive physical predictions for hadron collider phenomenology we need
to convolute our partonic coe cient functions with parton distribution functions (PDF).
Throughout this article we will use the PDF sets PDF4LHC15 [99]. We choose a Higgs
boson mass of 125 GeV and a top quark mass of mt(mt) = 162:7 GeV. We choose a value for
the strong coupling constant of
S(mZ = 91:1876 GeV) = 0:118. If not stated otherwise
we derive numerical predictions for protonproton collider with a center of mass energy of
The gure displays the contribution of N3LO coe cient function to the Higgs boson
production cross section for the gg (red), qg (green), qq (orange), qq(blue) and qQ(purple) initial
state as a function of the perturbative scale . The gg and qg coe cient function were rescaled
uniformly by a factor of 10 2.
13 TeV. We use a private c++ code to perform the numerical convolutions of PDFs and
partonic coe cient functions.
In gure 3 we display the impact of the N3LO corrections on the hadronic cross section
for di erent initial states as a function of the perturbative scale . The gluongluon (red)
and quarkgluon (blue) initial state contributions were rescaled by a factor of 10 2 in
order to t nicely. We observe that the numerical impact of these two channels is clearly
dominant over all other initial state con gurations. The nominally smallest corrections for
each channel can be found in an interval of
2 [40; 90] GeV.
In gure 4 we combine the contribution from all partonic coe cient functions and
evaluate their contribution to the hadronic cross section including lower orders in perturbation
theory as a function of the perturbative scale . We show LO, NLO, NNLO and N3LO
predictions in green, orange, blue and red respectively. We observe that the dependence on
the perturbative scale is greatly reduced at N3LO compared to lower orders. Furthermore,
NNLO and N3LO predictions overlap within the interval of
2
m4h ; mh .
To derive a concrete numerical prediction we choose the value of the cross section at
= m2h . We vary the perturbative scale in the interval
m4h ; mh in order to estimate the
e ect of missing higher order corrections at N4LO and beyond. As can be seen from
gure 4
this procedure is not conservative enough at leading and nexttoleading order. Regarding
the progression of the series from NLO onward we observe convergent behaviour. The
nominal size of the corrections is greatly reduced at each successive order. Uncertainty
estimates based on scale variation overlap at NNLO and N3LO.
LHC 13 TeV
PDF4LHC15.0
P P > H+X
LO
NNLO
NLO
N3LO
HJEP05(218)
30
50
70
90
110
130
150
170
190
210
230
250
μ [GeV]
. The green, orange, blue and red lines correspond to a prediction made by
truncating the perturbative series at LO, NLO, NNLO and N3LO respectively.
Our prediction for the Higgs boson production cross section at the LHC based on a
computation in perturbative QCD in the large top quark mass limit through N3LO of
P P !H+X = 45:18
0:13:148 pb = 45:18 pb
0:369:34%:
(4.1)
5
Comparison with results based on a threshold expansion
In ref. [14] N3LO corrections to the Higgs boson production cross section were computed
using an approximation based on a power series around the point z = 1 truncated at
O((1
z)30). The expansion around z = 1 exploits a kinematic enhancement of the gluon
luminosity in the collision of protons for lower values of partonic center of mass energy to
yield reliable predictions. The point z = 1 represents the production threshold for a Higgs
boson, i.e. the lowest possible amount of energy required to produce a Higgs boson. In
ref. [15] seven additional terms in the power series were added. The quality of a threshold
expansion for N3LO corrections was furthermore studied in refs. [25, 26, 100]. Having now
the complete coe cient functions at our disposal we want to re ect on previous estimates
and compare our exact analytical ndings to the approximate results.
Using the same setup as in the previous section to derive numerical predictions we
nd that the hadronic cross section through N3LO in perturbative QCD in the in nite top
quark mass limit based on thirty terms in the threshold limit is given by
P P !H+X
Threshold30 = 45:07
0:12:643 pb = 45:07 pb
0:358:23%:
(5.1)
We observe a di erence of 0:11 pb with respect to our new prediction, eq. (4.1). The scale
variation interval in eq. (4.1) is slightly larger. In ref. [15] it was estimated the e ect of
missing higher order terms in the threshold expansion are less than 0:18 pb. We now see
that this estimate was su ciently conservative.
In the remainder of this section we want to study the behaviour of N3LO corrections as
a function of the order where the threshold expansion is truncated. In particular we want
to investigate its performance for contributions arising from di erent partonic initial states.
In
gure 5 we show the N3LO correction due to di erent initial sate partons based on a
threshold expansion (red) as a function of the order at which the expansion is truncated. In
blue we also display our new result to all orders in the threshold expansion as a reference.
We observe that the rst four terms show particularly large changes in the derived
prediction. Starting from the fth term we observe slow asymptotic improvement towards the
full result. The nominally largest gluongluon and quarkgluon channels are approximated
better than their purely quark initiated counter parts. The sum of all channels can be seen
in gure 5(f).
In order to see more clearly the quality of the threshold expansion for each channel
we show in
gure 6 the impact of N3LO corrections on the hadronic cross section due to
di erent partonic initial states. The predictions in red are now based on a threshold
expansion normalised to the respective all order result. The xaxis shows the order at which
the threshold expansion is truncated. The line in blue at one serves as a reference. We
observe that contributions originating from the gluongluon channel are approximated within
several permille including only a few terms in the expansion. Similarly the quarkgluon
initiated contributions are approximated reasonably well below a level of ten percent. All
other contributions are considerably di erent from the exact result and receive corrections
of the order of 100% even with thirty terms in the expansion. Their nominal e ect on
the inclusive cross section is however negligible. The fact that the threshold expansion
works best for gluonic initial states can be explained by the fact that the probability to
extract gluons from a proton is peaked towards lower momentum fractions, i.e. closer to
the production threshold. For quarks this enhancement is not as large. The relatively
slow improvement towards the exact result of the predictions as more and more terms in
the threshold expansion are included can be understood from the high energy behaviour
of the partonic coe cient functions. As we displayed in
gure 2 the coe cient functions
have a power like divergence
log5(z) as z ! 0. While the threshold expansion is formerly
z
convergent within the entire physical interval a relatively slow convergence to capture the
high energy behaviour can be expected.
6
Conclusions
In this article we present an exact computation of the Higgs boson production cross section
at hadron colliders through N3LO in perturbative QCD in the in nite top quark mass limit.
The main result of this article are analytic formulae for N3LO corrections to the regular
partonic coe cient functions. We provide these functions in an ancillary le together with
the arXiv submission of this article.
3 ■ ●■ ■ ■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■
10
30
40
■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■
O
L
3 g
N g
σ1
0
0.01
0.008
0. ● ●
0
●
●
●
●
■ Ful
O
L
3 g
N q
σ0.6
0.8
O
L
3
Nσ1
0
●
0
10
30
40
Truncation Order
■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■
● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●
● Expansion
■ Ful
10
30
40
Truncation Order
■ Ful
■ Ful
PDF4LHC15.0
The gure shows in red the contribution of the partonic coe cient function to the N3LO
correction of the Higgs boson cross section approximated by a threshold expansion. The xaxis labels
the order at which the expansion is truncated. The line in blue represents the contribution to all
orders in the threshold expansion and is displayed as a reference. Figures (a), (b), (c), (d), (e) and
(f ) show the contribution due to the gg, qg, qq, qq, qQ initial state and the sum of all channels
respectively.
To obtain our result we compute matrix elements for the production of a Higgs boson in
association with three partons at tree level and with two partons at the oneloop level. In
order to perform required phase space integrals we employ the framework of reverse unitarity
and make use of loop integration techniques such as IBP identities and master integrals.
We compute all required master integrals using the framework of di erential equations.
■ ■ ■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■
20
(b)
10
30
40
●
■ ■ ■ ■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■ ●■
10
30
40
20
(d)
Truncation Order
LHC 13 TeV
PDF4LHC15.0
μ=125 GeV
● Expansion
■ Ful
10
30
40
20
(f)
Truncation Order
0.99
0.9
3LO_ ,lFuq0.8
N q
0.9
l 0.8
u
LO,F
3NqQ
O
L
3NqQ
/σ0.7
σ0.6
0.5 ● ●
●
10
Truncation Order
Truncation Order
Truncation Order
PDF4LHC15.0
■ Ful
● Expansion
■ Ful
LHC 13 TeV
PDF4LHC15.0
3 g
LO
3Nqg
N/σq1.04
σ1.02
σ
0.96 ●
0.95
■ Ful
● Expansion
■ Ful
LHC 13 TeV
PDF4LHC15.0
Truncation Order
(d)
10
20
30
40
Truncation Order
(f)
■ Ful
/σ0.6
O
L
3 q
Nσq0.5
0.4 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●
40
10
20
30
40
1. ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■
1. ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■
● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●
● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●
correction of the Higgs boson cross section approximated by a threshold expansion normalised to
the all order result. The xaxis labels the order at which the expansion is truncated. The line
in blue represents the contribution to all orders in the threshold expansion and is displayed as a
reference. Figures (a), (b), (c), (d), (e) and (f) show the contribution due to the gg, qg, qq, qq, qQ
initial state and the sum of all channels respectively.
When solving di erential equations we encounter elliptic integrals in the solution for triple
real radiation master integrals. We
nd that an analytic solution for our master integrals
can be easily found by embedding the solution of our di erential equations in a fairly
general class of iterated integrals. We discuss in detail how we
nd relations among iterated
integrals involving elliptic functions and how we evaluate them e ciently numerically.
2 TeV
7 TeV
8 TeV
13 TeV
14 TeV
28 TeV
100 TeV
+0:05pb
0:09pb
+0:70pb
1:14pb
+0:90pb
1:43pb
+2:07pb
3:16pb
+2:34pb
3:54pb
+7:02pb
9:93pb
+44:53pb
56:95pb
+4:17%
8:02%
+4:17%
6:76%
+4:18%
6:69%
+4:26%
6:48%
+4:28%
6:46%
+4:54%
6:42%
+5:51%
7:05%
(PDF)
0:03 pb ( 3:17%)
0:31 pb ( 1:89%)
0:40 pb ( 1:87%)
0:89 pb ( 1:85%)
1:00 pb ( 1:86%)
2:98 pb ( 1:96%)
19:98 pb ( 2:51%)
+0:04pb
0:04pb
+0:44pb
0:45pb
+0:56pb
0:56pb
+1:25pb
1:26pb
+1:40pb
1:42pb
+4:10pb
4:03pb
+24:89pb
21:71pb
+3:69%
3:36%
+2:66%
2:68%
+2:63%
2:66%
+2:59%
2:62%
+2:60%
2:62%
+2:70%
2:65%
+3:12%
2:72%
Having obtained analytic expressions for all required partonic cross sections we embed
them in a numerical code and derive predictions for hadron collider cross sections. We nd
that N3LO corrections are small compared to corrections at previous orders and that the
dependence on the perturbative scale is greatly reduced. We perform a detailed
comparison with a previous approximation of N3LO corrections based on an expansion around the
production threshold of the Higgs boson including 37 terms [14, 15]. We observe that our
new results are in excellent agreement with this approximation. Dominant contributions
due to gluon initiated partonic cross sections are approximated rather well by the
threshold expansion. Quark initiated contributions on the other hand are approximated rather
poorly. The estimate of missing higher orders in the threshold expansion in refs. [14, 15]
was su ciently conservative to cover the di erence to the exact result.
To derive precise predictions for hadron collider phenomenology many e ects beyond
the e ective theory cross section considered in this article have to be take into account.
The
niteness of quark masses and neglected electroweak e ects play an important role.
It is particularly important to critically asses all nonnegligible sources of uncertainty. A
detailed study of the inclusive production cross section for the Higgs boson considering all
such e ects was conducted in ref. [15]. Repeating this analysis is beyond the scope of this
article. However, it easily possible to modify the nal predictions for hadron collider cross
sections of ref. [15] such that the results of this article are taken into account. Speci cally,
we include the exact contributions to the cross section at N3LO in the EFT and remove
uncertainties due to the truncation of the threshold expansion. Otherwise, we can simply
use the results of ref. [15] that are neatly combined in a new numerical code iHixs 2 [101].
In table 1 we show updated predictions for the gluon fusion Higgs boson production cross
section at the LHC as in ref. [101].
HJEP05(218)
Acknowledgments
I would like to thank Babis Anastasiou, Claude Duhr, Falko Dulat and Franz Herzog for
many invaluable discussions. I would like to thank Claude Duhr and Babis Anastasiou for
useful comments on the manuscript. I am grateful to Lorenzo Tancredi and Stefan Weinzierl
for useful discussions about elliptic integrals appearing in this computation. Furthermore,
I would like to thank Achilleas Lazopoulos and Falko Dulat for comparisons of
numerical results using the latest version of iHixs 2 [101]. My research was supported by the
European Comission through the ERC Consolidator Grant HICCUP (No. 614577).
A
The elliptic integral
In section 3.2 we discuss a coupled system of two di erential equations that describes the
homogeneous solution to master integrals appearing in triple real radiation matrix elements
when integrated over phase space. The particular system is given by
HJEP05(218)
E0 !
4
E0
1
=
0
3 z
E0 !
4
E0
1
:
Equivalently, we can say that E40 satis es a second order di erential equation.
3z2
z (z2
22z
11z
3)
11z
1) E40 = 0:
E10 = z
(A.1)
(A.2)
First, a solution to this di erential equation was found by Stefan Weinzierl in terms of an
elliptic integral.
The homogeneous part of a di erential equation for a Feynman integral has to be
satis ed by the maximum cut of the corresponding Feynman integral. In refs. [75, 76] it
was demonstrated how one can nd a solution for a coupled homogeneous part of a system
di erential equations for Feynman integrals. In ref. [41] it was even proposed that it is
su cient to normalise the leading singularities of Feynman integrals to constants in order
to decouple their di erential equations order by order in the dimensional regulator. For
this to hold true the physical linear combinations of leading singularities themselves must
satisfy the homogeneous di erential equation for
= 0. Computing the leading singularity
of E4 we nd
Leading Singularity (E4)
Z
z) x
3
z) (x3
x2z + 2x2 + 2xz + x
x2z + 2x2 + 2xz + x
z
z)
:
(A.3)
We can rewrite the quartic polynomial under the square root as
x2z + 2x2 + 2xz + x
z = (x
r1)(x
r2)(x
r3)(x
r4):
(A.4)
I1 =
I2 =
=
Z r3
r2
Z r4
r3
p(r4
K(1
m):
r2)(r3
r1)
r1)(x
r2)(x
r3)(x
r4)
r1)(x
r2)(x
r3)(x
r4)
Following the prescription of ref. [59] we de ne two integrals 4
2
+
1
1
Here, K(m) is the complete elliptic integral of the rst kind. We nd that both integrals
I1 and I2 are solutions to our second order di erential equation eq. (A.2). In principle
we could now follow a procedure outlined in ref. [59] to construct a transformation matrix
TE that allows us to decouple the system of di erential equations order by order in .
Speci cally, we nd that the functions tij (z) de ned in section 3.2 are given by linear
combinations
ci 2 C:
The derivatives of the functions I1 and I2 with respect to z yield a sum of elliptic integrals
of rst and second kind with algebraic prefactors. We can determine the coe cients ci
analytically by equating the power series expansions of the above equation with the results
obtained in section 3.2. However, any of these analytic expressions is quite unwieldy.
(A.5)
(B.1)
Various ingredients for Higgs boson production
In this appendix we summarise various standard ingredients for the perturbative calculation
of the inclusive Higgs boson production cross section.
In order to perform renormalisation in the MS scheme we substitute the bare coupling
and Wilson coe cient as
S0 =
S( 2)
C0 = CZC :
e E Z :
The renormalisation factors for the strong coupling constant and Wilson coe cient required
for a computation through N3LO [19] are given by
Z
= 1 +
ZC = 1
S
S
0
0
30 +
6 2
2 0 1 +
2
3
2
+ O( S4):
2
+ O( S4):
The coe cients at the various orders in the coupling constant i are given by the QCD
beta function [102{105].
S
S 2
1
1
6 3
Pi(k0)
Pi(k0)
Pk(j0)
Pk(l0)
0Pi(j0)
+
1
2
Pk(j1)
3 0Pi(k0)
Pk(j0) + 2 02Pi(j0)
Pi(k1)
Pk(j0) + 2Pi(k0)
2 0Pi(j1)
2 1Pi(j0)
+
Pi(j2) :
In order to obtain infrared nite cross sections we are required to perform a suitable rede nition of our parton distribution functions.
Z 1
0
fi(x) = fiR
g)(z) =
dxdyf (x)g(y) (xy
z):
(B.3)
The infrared counter term
consists of convolutions [29] of splitting functions Pi(jn) [30, 31]
and can be derived from the DGLAP equation. Its perturbative expansion required for
an N3LO accurate calculation of the di erential Higgs boson production cross section is
+
log
mt2
2
+ O( S4
) :
(B.4)
(B.5)
In the e ective theory with nf light avours and the top quark decoupled from the
running of the strong coupling constant, the MSscheme Wilson coe cient reads [79{82]
C( 2) =
S
3 v
S
11
4
6865
31104
mt2
2
2892659
41472
897943
9216
S 2 2777
288
55
54
3 +
77
1728
log
209
64
log
mt2
2
log2
19
16
+
log
1
18
mt2
2
log2
mt2
2
nf
40291
20736
110779
13824 3 nf
1733
288
log
mt2
2
67
96
n
2
f
mt2
2
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.
[1] ATLAS collaboration, Observation of a new particle in the search for the Standard Model
Higgs boson with the ATLAS detector at the LHC, Phys. Lett. B 716 (2012) 1
[arXiv:1207.7214] [INSPIRE].
[2] CMS collaboration, Observation of a new boson at a mass of 125 GeV with the CMS
experiment at the LHC, Phys. Lett. B 716 (2012) 30 [arXiv:1207.7235] [INSPIRE].
HJEP05(218)
Nucl. Phys. B 453 (1995) 17 [hepph/9504378] [INSPIRE].
QCD corrections, Phys. Lett. B 264 (1991) 440 [INSPIRE].
Nucl. Phys. B 646 (2002) 220 [hepph/0207004] [INSPIRE].
[7] C. Anastasiou and K. Melnikov, Higgs boson production at hadron colliders in NNLO QCD,
[8] R.V. Harlander and W.B. Kilgore, Nexttonexttoleading order Higgs production at hadron
colliders, Phys. Rev. Lett. 88 (2002) 201801 [hepph/0201206] [INSPIRE].
[9] V. Ravindran, J. Smith and W.L. van Neerven, NNLO corrections to the total crosssection
for Higgs boson production in hadron hadron collisions, Nucl. Phys. B 665 (2003) 325
[hepph/0302135] [INSPIRE].
Quarks, Z. Physik C 18 (1983) 69.
[10] T. Inami, T. Kubota and Y. Okada, E ective Gauge Theory and the E ect of Heavy
[11] M.A. Shifman, A.I. Vainshtein and V.I. Zakharov, Remarks on Higgs Boson Interactions
with Nucleons, Phys. Lett. B 78 (1978) 443 [INSPIRE].
[12] V.P. Spiridonov and K.G. Chetyrkin, Nonleading mass corrections and renormalization of
the operators m psibar psi and g2( ), Sov. J. Nucl. Phys. 47 (1988) 522 [INSPIRE].
[13] F. Wilczek, Decays of Heavy Vector Mesons Into Higgs Particles, Phys. Rev. Lett. 39
(1977) 1304 [INSPIRE].
[14] C. Anastasiou, C. Duhr, F. Dulat, F. Herzog and B. Mistlberger, Higgs Boson
GluonFusion Production in QCD at Three Loops, Phys. Rev. Lett. 114 (2015) 212001
[arXiv:1503.06056] [INSPIRE].
[15] 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].
[16] 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].
[17] R.V. Harlander, S. Liebler and H. Mantler, SusHi Bento: Beyond NNLO and the heavytop
limit, Comput. Phys. Commun. 212 (2017) 239 [arXiv:1605.03190] [INSPIRE].
[18] P.A. Baikov, K.G. Chetyrkin, A.V. Smirnov, V.A. Smirnov and M. Steinhauser, Quark and
gluon form factors to three loops, Phys. Rev. Lett. 102 (2009) 212002 [arXiv:0902.3519]
[19] T. Gehrmann, E.W.N. Glover, T. Huber, N. Ikizlerli and C. Studerus, Calculation of the
quark and gluon form factors to three loops in QCD, JHEP 06 (2010) 094
[arXiv:1004.3653] [INSPIRE].
contribution to inclusive Higgs production at N3LO, JHEP 02 (2015) 077
[22] F. Dulat and B. Mistlberger, RealVirtualVirtual contributions to the inclusive Higgs cross
section at N3LO, arXiv:1411.3586 [INSPIRE].
[20] C. Anastasiou, C. Duhr, F. Dulat, F. Herzog and B. Mistlberger, Realvirtual contributions
to the inclusive Higgs crosssection at N 3LO, JHEP 12 (2013) 088 [arXiv:1311.1425]
[23] T. Gehrmann, M. Jaquier, E.W.N. Glover and A. Koukoutsakis, TwoLoop QCD
Corrections to the Helicity Amplitudes for H ! 3 partons, JHEP 02 (2012) 056
QCD, JHEP 03 (2015) 091 [arXiv:1411.3584] [INSPIRE].
Lett. B 737 (2014) 325 [arXiv:1403.4616] [INSPIRE].
[26] C. Anastasiou et al., Higgs boson gluon{fusion production at threshold in N3LO QCD, Phys.
[27] Y. Li, A. von Manteu el, R.M. Schabinger and H.X. Zhu, Softvirtual corrections to Higgs
production at N3LO, Phys. Rev. D 91 (2015) 036008 [arXiv:1412.2771] [INSPIRE].
[28] S. Buehler and A. Lazopoulos, Scale dependence and collinear subtraction terms for Higgs
production in gluon fusion at N3LO, JHEP 10 (2013) 096 [arXiv:1306.2223] [INSPIRE].
[29] M. Hoschele, J. Ho , A. Pak, M. Steinhauser and T. Ueda, Higgs boson production at the
LHC: NNLO partonic cross sections through order
and convolutions with splitting
functions to N3LO, Phys. Lett. B 721 (2013) 244 [arXiv:1211.6559] [INSPIRE].
[30] S. Moch, J.A.M. Vermaseren and A. Vogt, The Three loop splitting functions in QCD: The
Nonsinglet case, Nucl. Phys. B 688 (2004) 101 [hepph/0403192] [INSPIRE].
[31] A. Vogt, S. Moch and J.A.M. Vermaseren, The Threeloop splitting functions in QCD: The
Singlet case, Nucl. Phys. B 691 (2004) 129 [hepph/0404111] [INSPIRE].
[32] C. Anastasiou, S. Buehler, C. Duhr and F. Herzog, NNLO phase space master integrals for
twotoone inclusive cross sections in dimensional regularization, JHEP 11 (2012) 062
[arXiv:1208.3130] [INSPIRE].
[33] C. Anastasiou and K. Melnikov, Pseudoscalar Higgs boson production at hadron colliders in
NNLO QCD, Phys. Rev. D 67 (2003) 037501 [hepph/0208115] [INSPIRE].
[34] C. Anastasiou, L.J. Dixon and K. Melnikov, NLO Higgs boson rapidity distributions at
hadron colliders, Nucl. Phys. Proc. Suppl. 116 (2003) 193 [hepph/0211141] [INSPIRE].
[35] C. Anastasiou, L.J. Dixon, K. Melnikov and F. Petriello, Dilepton rapidity distribution in
the DrellYan process at NNLO in QCD, Phys. Rev. Lett. 91 (2003) 182002
[hepph/0306192] [INSPIRE].
[36] C. Anastasiou, L.J. Dixon, K. Melnikov and F. Petriello, High precision QCD at hadron
colliders: Electroweak gauge boson rapidity distributions at NNLO, Phys. Rev. D 69 (2004)
094008 [hepph/0312266] [INSPIRE].
[37] F.V. Tkachov, A theorem on analytical calculability of four loop renormalization group
functions, Phys. Lett. 100B (1981) 65 [INSPIRE].
[38] K.G. Chetyrkin and F.V. Tkachov, Integration by Parts: The Algorithm to Calculate
[39] A.V. Kotikov, Di erential equations method: New technique for massive Feynman diagrams
calculation, Phys. Lett. B 254 (1991) 158 [INSPIRE].
Nucl. Phys. B 580 (2000) 485 [hepph/9912329] [INSPIRE].
[40] T. Gehrmann and E. Remiddi, Di erential equations for two loop four point functions,
[41] J.M. Henn, Multiloop integrals in dimensional regularization made simple, Phys. Rev. Lett.
110 (2013) 251601 [arXiv:1304.1806] [INSPIRE].
[42] C. Anastasiou, C. Duhr, F. Dulat and B. Mistlberger, Soft triplereal radiation for Higgs
production at N3LO, JHEP 07 (2013) 003 [arXiv:1302.4379] [INSPIRE].
[43] C. Anastasiou, C. Duhr, F. Dulat, E. Furlan, F. Herzog and B. Mistlberger, Soft expansion
of doublerealvirtual corrections to Higgs production at N3LO, JHEP 08 (2015) 051
[arXiv:1505.04110] [INSPIRE].
[44] D.J. Broadhurst, The Master Two Loop Diagram With Masses, Z. Phys. C 47 (1990) 115
massive two loop selfenergy diagrams, Nucl. Phys. B 434 (1995) 383 [hepph/9409388]
[46] M. Ca o, H. Czyz, S. Laporta and E. Remiddi, The Master di erential equations for the
two loop sunrise selfmass amplitudes, Nuovo Cim. A 111 (1998) 365 [hepth/9805118]
[47] R. Bonciani, V. Del Duca, H. Frellesvig, J.M. Henn, F. Moriello and V.A. Smirnov,
Twoloop planar master integrals for Higgs ! 3 partons with full heavyquark mass
dependence, JHEP 12 (2016) 096 [arXiv:1609.06685] [INSPIRE].
[48] A. von Manteu el and L. Tancredi, A nonplanar twoloop threepoint function beyond
multiple polylogarithms, JHEP 06 (2017) 127 [arXiv:1701.05905] [INSPIRE].
[49] S. CaronHuot and K.J. Larsen, Uniqueness of twoloop master contours, JHEP 10 (2012)
026 [arXiv:1205.0801] [INSPIRE].
[50] J.L. Bourjaily, A.J. McLeod, M. Spradlin, M. von Hippel and M. Wilhelm, The Elliptic
DoubleBox Integral: Massless Amplitudes Beyond Polylogarithms, Phys. Rev. Lett. 120
(2018) 121603 [arXiv:1712.02785] [INSPIRE].
[51] L.B. Chen, J. Jiang and C.F. Qiao, TwoLoop integrals for CPeven heavy quarkonium
production and decays: Elliptic Sectors, arXiv:1712.03516 [INSPIRE].
[52] M. Czakon and A. Mitov, Inclusive Heavy Flavor Hadroproduction in NLO QCD: The
Exact Analytic Result, Nucl. Phys. B 824 (2010) 111 [arXiv:0811.4119] [INSPIRE].
[53] A.B. Goncharov, Multiple polylogarithms, cyclotomy and modular complexes, Math. Res.
Lett. 5 (1998) 497 [arXiv:1105.2076] [INSPIRE].
[54] A.B. Goncharov, M. Spradlin, C. Vergu and A. Volovich, Classical Polylogarithms for
Amplitudes and Wilson Loops, Phys. Rev. Lett. 105 (2010) 151605 [arXiv:1006.5703]
JHEP 08 (2012) 043 [arXiv:1203.0454] [INSPIRE].
148 (2015) 328 [arXiv:1309.5865] [INSPIRE].
[55] C. Duhr, H. Gangl and J.R. Rhodes, From polygons and symbols to polylogarithmic
functions, JHEP 10 (2012) 075 [arXiv:1110.0458] [INSPIRE].
[56] E. Panzer, Algorithms for the symbolic integration of hyperlogarithms with applications to
[59] S. Laporta and E. Remiddi, Analytic treatment of the two loop equal mass sunrise graph,
HJEP05(218)
Nucl. Phys. B 704 (2005) 349 [hepph/0406160] [INSPIRE].
[60] L. Adams, C. Bogner and S. Weinzierl, The twoloop sunrise graph in two spacetime
dimensions with arbitrary masses in terms of elliptic dilogarithms, J. Math. Phys. 55 (2014)
102301 [arXiv:1405.5640] [INSPIRE].
[61] L. Adams, C. Bogner and S. Weinzierl, The twoloop sunrise integral around four
spacetime dimensions and generalisations of the Clausen and Glaisher functions towards
the elliptic case, J. Math. Phys. 56 (2015) 072303 [arXiv:1504.03255] [INSPIRE].
[62] E. Remiddi and L. Tancredi, An Elliptic Generalization of Multiple Polylogarithms, Nucl.
Phys. B 925 (2017) 212 [arXiv:1709.03622] [INSPIRE].
[63] E. Remiddi and L. Tancredi, Di erential equations and dispersion relations for Feynman
amplitudes. The twoloop massive sunrise and the kite integral, Nucl. Phys. B 907 (2016)
400 [arXiv:1602.01481] [INSPIRE].
[64] L. Adams, C. Bogner, A. Schweitzer and S. Weinzierl, The kite integral to all orders in
terms of elliptic polylogarithms, J. Math. Phys. 57 (2016) 122302 [arXiv:1607.01571]
[INSPIRE].
77 (2017) 77 [arXiv:1610.06207] [INSPIRE].
[65] G. Passarino, Elliptic Polylogarithms and Basic Hypergeometric Functions, Eur. Phys. J. C
[66] S. Bloch, M. Kerr and P. Vanhove, Local mirror symmetry and the sunset Feynman integral,
Adv. Theor. Math. Phys. 21 (2017) 1373 [arXiv:1601.08181] [INSPIRE].
[67] L. Adams and S. Weinzierl, Feynman integrals and iterated integrals of modular forms,
arXiv:1704.08895 [INSPIRE].
arXiv:1706.01299 [INSPIRE].
[68] J. Ablinger et al., Iterated Elliptic and Hypergeometric Integrals for Feynman Diagrams,
[69] F.C.S. Brown and A. Levin, Multiple elliptic polylogarithms, arXiv:1110.6917 [INSPIRE].
[70] J. Broedel, N. Matthes and O. Schlotterer, Relations between elliptic multiple zeta values
and a special derivation algebra, J. Phys. A 49 (2016) 155203 [arXiv:1507.02254]
[INSPIRE].
[71] M. Hidding and F. Moriello, All orders structure and e cient computation of linearly
reducible elliptic Feynman integrals, arXiv:1712.04441 [INSPIRE].
[72] J.L. Bourjaily, A.J. McLeod, M. Spradlin, M. von Hippel and M. Wilhelm, The Elliptic
DoubleBox Integral: Massless Amplitudes Beyond Polylogarithms, Phys. Rev. Lett. 120
(2018) 121603 [arXiv:1712.02785] [INSPIRE].
[73] J. Broedel, C. Duhr, F. Dulat and L. Tancredi, Elliptic polylogarithms and iterated integrals
on elliptic curves I: general formalism, arXiv:1712.07089 [INSPIRE].
[74] J. Broedel, C. Duhr, F. Dulat and L. Tancredi, Elliptic polylogarithms and iterated integrals
on elliptic curves II: an application to the sunrise integral, arXiv:1712.07095 [INSPIRE].
[75] A. Primo and L. Tancredi, On the maximal cut of Feynman integrals and the solution of
their di erential equations, Nucl. Phys. B 916 (2017) 94 [arXiv:1610.08397] [INSPIRE].
[76] A. Primo and L. Tancredi, Maximal cuts and di erential equations for Feynman integrals.
An application to the threeloop massive banana graph, Nucl. Phys. B 921 (2017) 316
[arXiv:1704.05465] [INSPIRE].
Acta Phys. Polon. B 46 (2015) 2125 [INSPIRE].
[77] L. Tancredi, Simplifying Systems of Di erential Equations. The Case of the Sunrise Graph,
HJEP05(218)
[78] J. Broedel, C. Duhr, F. Dulat, B. Penante and L. Tancredi, Elliptic symbol calculus: from
elliptic polylogarithms to iterated integrals of Eisenstein series, arXiv:1803.10256
[79] K.G. Chetyrkin, B.A. Kniehl and M. Steinhauser, Decoupling relations to O( S3) and their
connection to lowenergy theorems, Nucl. Phys. B 510 (1998) 61 [hepph/9708255]
744 (2006) 121 [hepph/0512060] [INSPIRE].
[80] Y. Schroder and M. Steinhauser, Fourloop decoupling relations for the strong coupling,
[81] K.G. Chetyrkin, J.H. Kuhn and C. Sturm, QCD decoupling at four loops, Nucl. Phys. B
[82] M. Kramer, E. Laenen and M. Spira, Soft gluon radiation in Higgs boson production at the
LHC, Nucl. Phys. B 511 (1998) 523 [hepph/9611272] [INSPIRE].
[83] P. Nogueira, Automatic Feynman graph generation, J. Comput. Phys. 105 (1993) 279
[84] C.W. Bauer, R. Kreckel and A. Frink, Introduction to the GiNaC framework for symbolic
computation within the C++ programming language, J. Symb. Comput. 33 (2000) 1.
[85] S. Laporta, High precision calculation of multiloop Feynman integrals by di erence
equations, Int. J. Mod. Phys. A 15 (2000) 5087 [hepph/0102033] [INSPIRE].
[86] M.A. Barkatou and E. P ugel, Computing superirreducible forms of systems of linear
di erential equations via moserreduction: A new approach, in Proceedings of the 2007
International Symposium on Symbolic and Algebraic Computation (ISSAC '07), New York
U.S.A. (2007).
[87] M.A. Barkatou and E. P ugel, On the moser and superreduction algorithms of systems of
linear di erential equations and their complexity, J. Symb. Comput. 44 (2009) 1017
[88] J. Moser, The order of a singularity in fuchs' theory, Math. Zeitschrift 72 (1959) 379.
[89] R.N. Lee, Reducing di erential equations for multiloop master integrals, JHEP 04 (2015)
108 [arXiv:1411.0911] [INSPIRE].
[90] K.T. Chen, Iterated path integrals, Bull. Am. Math. Soc. 83 (1977) 831 [INSPIRE].
[91] E. Remiddi and J.A.M. Vermaseren, Harmonic polylogarithms, Int. J. Mod. Phys. A 15
(2000) 725 [hepph/9905237] [INSPIRE].
new iterated integrals for massive Feynman diagrams, PoS(LL2014)020 [arXiv:1407.4721]
JHEP 03 (2014) 071 [arXiv:1401.4361] [INSPIRE].
HJEP05(218)
Phys. Commun. 174 (2006) 222 [hepph/0507152] [INSPIRE].
m2H =s much less than 1, Phys. Lett. B 535 (2002) 159 [hepph/0203140] [INSPIRE].
[arXiv:1506.02674] [INSPIRE].
[arXiv:1405.5685] [INSPIRE].
arXiv:1802.00827 [INSPIRE].
Convergence of the Threshold Expansion, in Proceedings of 49th Rencontres de Moriond on
function and anomalous dimensions, Nucl. Phys. B 710
YangMills theory with fermions, JHEP 02 (2017) 090 [arXiv:1701.01404] [INSPIRE].
[3] S. Dawson , Radiative corrections to Higgs boson production , Nucl. Phys. B 359 ( 1991 ) 283 [4] D. Graudenz , M. Spira and P.M. Zerwas , QCD corrections to Higgs boson production at proton proton colliders , Phys. Rev. Lett . 70 ( 1993 ) 1372 [INSPIRE].
[5] M. Spira , A. Djouadi , D. Graudenz and P.M. Zerwas , Higgs boson production at the LHC , [6] A. Djouadi , M. Spira and P.M. Zerwas , Production of Higgs bosons in proton colliders: [21] C. Duhr , T. Gehrmann and M. Jaquier , Twoloop splitting amplitudes and the singlereal functions in 4 Loops, Nucl . Phys. B 192 ( 1981 ) 159 [INSPIRE]. Feynman integrals, Comput. Phys. Commun . 188 ( 2015 ) 148 [arXiv: 1403 .3385] [INSPIRE].
[57] C. Duhr , Hopf algebras, coproducts and symbols: an application to Higgs boson amplitudes , [58] S. Bloch and P. Vanhove , The elliptic dilogarithm for the sunset graph , J. Number Theor.
[92] L. Adams , C. Bogner and S. Weinzierl , The sunrise integral and elliptic polylogarithms , PoS(LL2016)033 [arXiv:1606 .09457] [INSPIRE].
[93] J. Broedel , From elliptic iterated integrals to elliptic multiple zeta values , PoS(LL2016 ) 081 .
[94] J. Ablinger , J. Blumlein, C.G. Raab and C. Schneider , Nested (inverse) binomial sums and [96] D. Ma ^tre, HPL, a mathematica implementation of the harmonic polylogarithms , Comput.