Double-Higgs boson production in the high-energy limit: planar master integrals

We consider the virtual corrections to the process gg → HH at NLO in the high energy limit and compute the corresponding planar master integrals in an expansion for small top quark mass. We provide details on the evaluation of the boundary conditions and present analytic results expressed in terms of harmonic polylogarithms.


Introduction
One of the main aims of particle physics in the coming years is the exploration of the scalar sector of the theory which describes fundamental interactions, be it the Standard Model or an extension.One has to clarify whether the Higgs boson is a fundamental particle and how the particles of the theory obtain their mass.A process which helps to find answers to these questions is the production of Higgs boson pairs, since it is the simplest process which is sensitive to the triple-Higgs boson coupling.Although experimentally quite challenging, there is a chance that double Higgs production will be observed after the the high-luminosity upgrade of the CERN LHC.
The leading order (LO) corrections to Higgs boson pair production have been computed in Refs.[1,2] including the exact dependence on the top quark mass and the Mandelstam variables.At next-to-leading order (NLO), QCD corrections were computed for the first time in Ref. [3] in the infinite top quark mass limit using an effective theory and an independent cross check was provided in [4] by performing an asymptotic expansion in the full theory.In this way a quantitative estimate of the quark mass effects could be provided.Virtual NLO corrections in the large-m t limit have also been computed in Ref. [5], confirming the results of Ref. [4].Finite top quark mass effects have also been considered in Ref. [6], in which the exact real radiation contribution is combined with the effective-theory virtual corrections.Within the effective theory also next-to-next-toleading (NNLO) contributions are available [7,8].The NNLO result was completed in Ref. [9] in which the three-loop matching coefficient of the effective operator for two Higgs bosons and two, three or four gluons was computed.Note that it differs from that of single Higgs boson production.The result of [9] has been complemented by power-suppressed terms in the top quark mass in Ref. [10], where the soft-virtual approximation was constructed.The resummation of threshold-enhanced logarithms to next-to-next-to-leading logarithmic (NNLL) accuracy has been performed in Refs.[11,12] and transverse momentum resummation has been considered in Ref. [13].Differential distributions through NNLO for various observables were computed in Ref. [14] in the heavy-top limit.Finally, exact NLO results became available in Refs.[15,16] using a numerical approach for the computation of the two-loop virtual corrections.More recently these results have been matched to parton showers in Ref. [17].
In this paper we study a class of massive two-loop four-point functions with massless external particles.We describe in detail the methods used for the computation of the amplitudes and in particular the evaluation of the master integrals.We aim to study double Higgs boson production via the process gg → HH.Numerical NLO results are available [15,16], however the calculation of cross sections is computationally expensive and we want to provide an independent cross check in the high-energy region.We wish to provide results in terms of compact analytic expressions which can be used to construct simple approximations or can be used directly in the kinematic region in which they are valid.In this paper we provide the first step towards this goal by considering the part of the amplitude which is expressed in terms of planar master integrals.
We perform our calculation in the limit of vanishing Higgs boson mass which provides, as we will demonstrate in Section 3, a good approximation to the general case where m H = 0. Furthermore, finite Higgs-mass effects can be incorporated by a simple Taylor expansion.Recently the amplitudes for single-Higgs boson plus jet production have been considered in the limit of large Higgs transverse momentum [18].In this reference an expansion for small Higgs boson mass has been performed and thus the underlying integrals are the same as those of our calculation, so part of our findings can be cross checked against Ref. [18].
In the recent literature one can find several calculations where two-loop box integrals are also involved.However, the underlying integral families and/or the kinematics of the external and internal masses are different.For example, in Ref. [19] the amplitude of a Higgs boson and three partons has been considered.In the limit m H → 0 their integrals are also the same as ours.However, this limit cannot been taken since the calculation is performed in the Euclidean region with the assumption m 2 H < s < 0 and the results are expressed in terms of multiple polylogarithms, which can not easily be analytically continued into other regions.Similar arguments apply to other recent calculations such as [20] or [21]; analytic results have been obtained in terms of multiple polylogarithms which can in principle be evaluated numerically, but are very unwieldy.This is a another reason why we have decided to perform an expansion in the high energy limit.Our final results have a simple structure in terms of harmonic polylogarithms and can be evaluated numerically in a fast and reliable manner.
An interesting approach to obtain simple and easy-to-evaluate expressions for gg → HH at NLO has been developed in Ref. [22] where the large top mass expansion has been combined with expansion terms obtained for the top threshold.A good approximation of the exact (purely numerical) result [15,16] has been constructed by combining the different kinematic regions using Padé approximants.Further improvement is expected after incorporating information about the gg → HH amplitude at high energies which is the main purpose of this work.
The remainder of the paper is organized as follows: we introduce our notation in Section 2. In Section 3 we briefly consider the one-loop corrections to gg → HH in the high-energy limit to provide motivation for our calculation, and Section 4 describes the reduction of the amplitude to master integrals.The main part of the paper is Section 5 in which we discuss the calculation of the master integrals.We describe in detail the method we used to compute the boundary values necessary for the solution of differential equations for the master integrals.In this paper we refrain from presenting long formulae, which instead can be found in the ancillary file of this paper [23].
The partonic cross section is obtained from |M| 2 after integration over the phase space and multiplication by the flux factor.

One-loop considerations
Before providing details on our two-loop calculation we want to investigate the quality of the high-energy expansion at one-loop order.In the following we consider the differential partonic cross section dσ dθ (s, t) where the scattering angle θ of the Higgs boson in the center-of-mass frame enters via the following relation, In Fig. 1 we study the √ s-dependence of dσ/dθ for fixed scattering angle θ of 90 degrees.
The exact result (see the solid curve for m H = 0 and the short-dashed curve for m H = 0) is compared to various approximation, computed for m H = 0, incorporating high-energy expansions up to m 32 t (see the long-dashed and dash-dotted curves).For comparison we also show the result based on an effective-theory calculation in which the limit of infinite top quark mass is assumed (dotted curve).We observe that, as expected, the high-energy approximations lead to good results for large values of √ s.A systematic improvement is obtained after including higher order expansion terms.space has to be considered when performing the integration over θ.In practical applications this does not constitute a big problem since θ → 0, π corresponds to the forward and backward scattering of the Higgs boson, where no measurement can be performed.Furthermore, the bulk of the cross section is provided by the central region.For example, if we restrict 0.25π < θ < 0.75π in the exact one-loop corrections we cover around 70% of the full cross section for √ s = 1000 GeV.
Apart from providing an independent and analytically simple expression in the high energy region, which can be used as a cross-check of the exact (numerical) results, our expressions also serve as input of the method based on Padé approximants [22] as already mentioned in the Introduction.

Reduction to master integrals
We generate our amplitudes with the program qgraf [24] and use q2e and exp [25,26] to rewrite the output to FORM [27] notation.exp is also used to assign to each Feynman diagram an integral family which is defined according to the topology and mass distribution of the internal lines.For our application we have defined 34 families.We use FORM to express the amplitude for each diagram as a linear combination of scalar integrals of a given family.
We use the C++ version of FIRE [28] for the reduction of all scalar integrals from each family to master integrals, with LiteRed [29,30] providing #rules for FIRE.All families were reduced using the publicly available version FIRE 5.2.Some families were also re-duced using the development versions FIRE 5.5, 5.6 (which use more information from the LiteRed files) in order to check whether the number of master integrals produced was smaller.This was the case, but the number of master integrals was still not minimal.We describe our procedure to obtain a minimal set of master integrals below.The Mathematica-readable tables generated by FIRE are transformed to FORM Fill statements, so that the reduction can be applied to the amplitude in FORM using a TableBase.The reduction rules are heavily manipulated in FORM before creating this TableBase.
After the FIRE reduction, each integral family contains between 7 and 77 master integrals (1395 in total (1+2 loops)).One must minimize the number of master integrals between all families; the use of FIRE's FindRules[] command yields a total of 231 = 10 + 221 master integrals.This does not constitute a minimal set.We use the following procedure, implemented in FORM, to find and eliminate "master integrals" which are in fact a linear combination of other integrals of the set.
1.For each n-line master integral with a dot on one of its lines, generate the n − 1 integrals which instead have a dot on one of the other lines.Append these integrals to the set of integrals from the amplitude.Ensure that each integral in this extended set is present in the reduction tables.

Apply relations from
FindRules[] to the extended set of integrals and consider equations of the form: Apply the reduction tables to these equations and discard all trivial equations.One obtains a set of non-trivial equations which relate some integrals in the original set of master integrals.
3. Use these equations to construct reduction relations into a final linearly independent set of master integrals.We solve the equations to obtain the "most complicated" (highest line count) integrals in terms of simpler integrals.
In Step 2, we consider such equations for all integrals I of complexity < 9, where we define the complexity as the sum of the absolute values of the propagator powers.The reduction rules for higher complexity integrals contain coefficients which are too large to efficiently manipulate with FORM's PolyRatFun.Despite this, the set of equations contains many redundancies.That these equations are all satisfied increases our confidence that our final set of master integrals is a minimal set.Additionally, FindRules[] maps many integrals into different integral families, so this procedure shows some consistency between our families.We note that no approximation (except m H = 0) is applied during the reduction procedure.In particular, we retain the exact m t dependence.
Following this procedure we reduce the number of two-loop master integrals from 221 to 161.A list of them, and all 10 one-loop1 master integrals, can be found in Appendix A.
At one-loop order we obtain the minimal set of ten master integrals simply by applying FindRules[] to the master integrals of the three one-loop families.The additional twoloop reduction relations are applied to the FIRE reduction relations before we create the FORM TableBase which we apply to the amplitude and to the right-hand-side of differential equations (see Subsection 5.1).
The computation of these master integrals is described in Section 5.

Calculation of master integrals
For the calculation of the master integrals we use the method of differential equations [31,32].We solve the differential equations using an appropriate ansatz which is described in Subsection 5.1.The boundary conditions (see Subsection 5.2) are fixed by evaluating the master integrals in the limit m t → 0. In some cases it is sufficient to evaluate the integrals in this limit for fixed t = s = −1.

Differential equations
We compute the master integrals in an unphysical region where two Mandelstam variables (s and t) are negative and u is positive.In this region, the integrals which we compute are real valued.We can analytically continue results obtained here into the physical region.
For each master integral we have three differential equations which are obtained by taking derivatives w.r.t.m 2 t , s and t.The derivatives are computed using LiteRed.Note that only two of the three differential equations are needed to construct the result.The third provides a consistency check.The generation of the system of differential equations requires the extension of the FIRE reduction tables.Note, however, that the additional integrals which are required are not difficult to reduce.
Differentiating the vector of master integrals, (MI), w.r.t.x = t, m 2 t and applying the reduction tables to the result leads to systems of equations where K x is a square matrix.
To solve the differential equations we follow two approaches.In the first, we make an ansatz for each master integral which is suitable to describe the solution in the limit m t → 0. Guided by asymptotic expansion we use (see also Refs.[18,33]) where l is the number of loops.To determine the coefficients c of the ansatz we use the following procedure: (1) Use the differential equation for t and determine the coefficients of the leading terms in the m t → 0 expansion.This requires the solution of a system of first-order differential equations for the t-dependent coefficients c.Boundary conditions are needed for one specific value of t, e.g., for t = s = −1.
(2) Use the differential equation for m 2 t to obtain relations between the coefficients of the higher order m 2 t terms and the leading terms determined in (1).Since the m t dependence is explicit in the ansatz one only has to solve a system of linear equations.
(3) The results for the master integrals are inserted into the differential equation for s, which must be satisfied.
At one-loop order the matrix K t in Eq. ( 12) has a triangular structure.Thus, starting from the simplest integral, one can run through the vector of master integrals and solve the system integral-by-integral.
At two-loop order K t in Eq. ( 12) rather has a block-triangular structure.It contains blocks of up to four integrals, which form coupled systems of differential equations.The integrals within these blocks must be determined together.
Using the above approach based on t-independent boundary conditions we were not able to obtain results for all two-loop master integrals.Presumably this is due to an inconvenient choice of our set of master integrals.For this reason we developed a second approach where we determine the master integrals in the limit m 2 t → 0, keeping the full t dependence.Afterwards we only have to solve the m 2 t differential equation which, as mentioned above, reduces to solving a system of linear equations.This approach provides results for all master integrals.Where possible, we compared the results of the two approaches and found complete agreement.
Let us end this subsection by making a brief comment on the possibility to introduce a canonical basis [34] for our master integrals.We made several attempts to produce such a basis using the publicly available programs Fuchsia [35] and CANONICA [36].We were not able to obtain a canonical basis for our master integrals for all sectors.Since we are interested in the small-m t limit we did not insist on obtaining a canonical basis.

Boundary conditions
The main tools which we use to compute the boundary conditions are the method of regions [37,38] and Mellin-Barnes techniques (see, e.g., [38]).Additionally, we make use of the PSLQ algorithm [39] and exploit the anticipated dependence on irrational numbers of our final result to obtain exact expressions.
In the following we provide details of each step of the calculation and give concrete examples for the master integral G 6 (1, 1, 1, 1, 1, 1, 1, 0, 0).See Appendix A for definitions of the integral families.
In this section we assume the scaling where the parameter χ ≪ 1 is introduced for convenience.
To begin, we express the Feynman integral in its alpha representation (using the routines provided in FIESTA [40]).This is a convenient starting point to apply the method of regions.For example, the integral G 6 (1, 1, 1, 1, 1, 1, 1, 0, 0) is expressed as where the functions U and F , the so-called first and second Symanzik polynomials, are given by In Eq. ( 15) we have introduced analytic regularization parameters δ i to regularize collinear divergences, which appear later in the calculation.The original integral is obtained by taking the sequence limit δ i → 0 for all δ i .
To implement the asymptotic expansion for χ → 0 we use the program asy.m [41].It provides scaling rules for the alpha parameters for the various regions which have to be considered.For the integral in Eq. ( 15) there are 13 relevant regions.One corresponds to the so-called hard region, in which all seven alpha parameters scale as "1".There are twelve so-called soft-collinear regions where some of the parameters have the weight "1" and others the weight "χ".For example, for region 2, we have that region 2 : In total only four regions need to be considered.The remaining eight regions can be obtained by simple symmetry considerations.
After the expansion, the original integral is expressed as a sum of homogeneously scaling integrals where the summation n spans the relevant regions and the subscript "(n)" indicates that the polynomials U and F from Eq. ( 16) are specific to the corresponding region.They are expanded to leading order in χ.Note that each integral on the r.h.s. is homogeneous in m 2 t (or χ) but not in s and t, since they are O(1) parameters.In the hard region there is only one soft parameter, m 2 t , and thus a naive Taylor expansion of the integrand has to be performed.This leads to purely massless integrals which are known in the literature [42,43].We have cross-checked these results up to the order in ǫ necessary for our application.Note that the contribution of the hard region is regular in δ i which means that one can take the limit δ i → 0 at the very beginning.
The soft-collinear regions are more involved.In the following we outline the calculation of the contribution from region 2 as an example.The calculation for other regions proceeds analogously.
1. We express each integral in terms of two-dimensional Mellin-Barnes integrals.For our example integral, we find the following form where z 12 = z 1 + z 2 , δ 123 = δ 1 + δ 2 + δ 3 and so on.Note that the integration contours of z 1 and z 2 are chosen to be straight lines parallel to the imaginary axis, satisfying −1 < Re (z 1 ) < Re (z 2 ) < 0.
2. Next we use the package MB.m [44] to analytically continue the regularization parameters δ i and ǫ to zero.As a result we obtain two-dimensional Mellin-Barnes representations which depend only on z 1 , z 2 and possibly on t/s.Note that the poles in δ i cancel among the contributions from the different regions, which provides a good check of our calculations.
3. We now transform the two-dimensional Mellin-Barnes integrals into one-dimensional integrals.In general this step is non-trivial; we provide more details in Appendix B.
For simple cases barnesroutines.m[45] can be used.
4. At this point we arrive at two types of one-dimensional integral: those which are independent of t/s and others which depend on t/s, such as We perform a high-precision numerical evaluation (300 digits) of the t-independent integrals and, after summing the contributions from all regions, apply the PSLQ algorithm [39] to re-construct the analytic result as a rational linear combination of 33 products (up to weight 6) of numbers from the set and use a further 200 digits to verify each result.The Nielsen generalized polylogarithm S 3,3 (−1) is implemented in Mathematica as PolyLog[3,3,-1].
For the t-dependent integrals we make an ansatz containing harmonic polylogarithms (HPLs) [46] up to weight 6 with alphabet ν i ∈ {−1, 0}, which we Taylor-expand in t/s.The series is expressed as a multivariate polynomial in (t/s), log(t/s), log(−m 2 t /s).We obtain the Taylor series of the integrals by taking their residues at z 1 = 0, 1, 2, 3, ... .We then use 50 low-order terms of the Taylor series to fix the coefficients c {ν i },n of the ansatz and check the result using a further 200 higher-order terms of the Taylor series.
Using the above procedure, we obtain the following δ i -independent result for our sample master integral G 6 (1, 1, 1, 1, 1, 1, 1, 0, 0), − 48H(−2, 0, 0, 0; t/s) + 12H(−1, −2, 0, 0; t/s) − 12H(−1, 0, 0, 0, 0; t/s) where l m = log(−m 2 t /s) and l t = log(t/s).To solve the m 2 t -differential equation for most of our master integrals, it is sufficient to obtain just the leading term in the small-m t expansion of the boundary condition.However, for 9 integrals it is also necessary to compute the next-to-leading term in the asymptotic expansion.In most of these cases we can simply apply the method discussed above and expand up to the next-to-leading term in χ.However, for two of the seven-line integrals, Step 3 above is hard to apply at the next-to-leading order.For these integrals, we use the corresponding t-differential equations to obtain the next-to-leading boundary conditions.

Solving the differential equations
Using the boundary conditions discussed in the previous subsection we can solve the differential equations of Subsection 5.1 in a straightforward way.All of our results are expressed in terms of HPLs.The expansion depth is limited only by the size of the intermediate expressions which enter the system of linear equations for the coefficients in our ansatz.We have expanded each master integral such that the final result for the amplitude gg → HH is available up to order m 16 t .Our results for the master integrals can be downloaded from [23].
For illustration we show in Fig. 2 the results for two master integrals: G 6 (1, 1, 1, 1, 1, 1, 1, 0, 0), which is used as an example in Subsection 5.2, and G 20 (1, 1, 1, 1, 1, 2, 1, 0, 0).We plot the real and imaginary parts of the ǫ 0 term as a function of √ s and choose t = −s/2, which corresponds to θ = π/2 (cf.Fig. 1), m t = 175 GeV and µ 2 = s.For clarity we multiply each integral by appropriate powers of m t and s such that the leading term starts with m 2 t and is dimensionless.In each case we display the approximations including m 2 t , m 4 t , m 8 t and m 16 t terms.
In the panels showing G 6 (1, 1, 1, 1, 1, 1, 1, 0, 0) we compare the approximations to the exact result, which has been obtained numerically using pySecDec [47].For this integral we observe a rapid convergence.In fact, the m 4 t , m 8 t and m 16 t curves agree with each other and the exact points down to √ s ≈ 600 GeV and the two highest approximation even down to √ s ≈ 500 GeV.It is interesting to note that the m 16 t curve reproduces to high accuracy the turning point at √ s ≈ 400 GeV and the steep rise below that energy.
It can not be expected that all master integrals show such good convergence properties.
In fact there are integrals, in particular some of the non-planar contributions, where the expansion parameter is m 2 t /u instead of m 2 t /s which results in a smaller radius of convergence.
For G 20 (1, 1, 1, 1, 1, 2, 1, 0, 0) we were not able to obtain stable numerical results using pySecDec and thus we only show our approximations.We observe a similar pattern as for the LO cross section shown in Fig. 1: the inclusion of more terms extends the convergence range in √ s down to smaller values.Furthermore, the curves including m 8 t and m 16 t terms agree down to √ s ≈ 900 GeV and it can be expected that above this energy a good approximation to the exact result can be provided.

Conclusions
The main focus of this paper is on NLO corrections to double Higgs boson production in the high energy region, where the top quark mass is assumed to be small compared to the kinematic variables s and t.Such considerations complement expansions for large top quark mass and around the threshold which have been considered in the literature, see Refs.[4,22].Furthermore, they provide an indpendent cross check of the exact calculation [15] which relies heavily on numerical methods.
In this paper we perform the reduction of the gg → HH amplitude to master integrals and compute all planar integrals in an expansion for small m 2 t .The expansion depth for each master integral is chosen such that the amplitude includes terms up to order m 16 t .Our analytic results for the master integrals are expressed in terms of HPLs and can be obtained in computer-readable form from [23].

B Reducing the dimensionality of Mellin-Barnes integrals
In this appendix, we consider the following two types of Mellin-Barnes integrals: and choose the contour to be along the imaginary axis.If some of the left poles and the right poles merge, a regularization and a subsequent analytic continuation are required (see the example, discussed below Eq. ( 42)).
Let us first briefly summarize the known properties of Mellin-Barnes integrals.
• If c = a 12 + b 123 , the second Barnes lemma can be applied to Eq. ( 30).However, there is no corresponding expression for Eq. ( 31); differentiation w.r.t. the parameters a i , b j gives relations between expressions which have the form of Eq. ( 31) with different The number of independent relations is smaller than the number of possible choices of X.Thus, no analytic result for Eq. ( 31) can be obtained.
In the following we sketch the derivation of a solution for Eq. ( 30) for the general case, which also yields a solution for Eq. ( 31) after differentiation w.r.t.one of the parameters.
Based on the assumptions about the relation between a i , b j (Eq.( 33)), we can express Eq. ( 30) as where 3 F 2 is the generalized hypergeometric function.The residues at z = a 2 + m are written in a similar manner.At this point the r.h.s. of Eq. ( 36) contains two 3 F 2 .We can transform it into an expression containing only one 3 F 2 using the relation [49], The l.h.s. of Eq. ( 40) is symmetric in a 1 ↔ a 2 and b 1 ↔ b 2 ↔ b 3 , however the symmetry among b j is not obvious on the r.h.s.We can show the symmetry by using the transformation formula of 3 F 2 given in Eq. (39).
In general, the generalized hypergeometric function This condition has to be satisfied when using Eq.(40).If condition ( 41) is violated we perform an analytic continuation to obtain an expression which converges at z = 1.This procedure is well-known [50], so we will not further discuss it here.As an example, let us consider the integral In this case, the right-most left-pole at z 1 = 0 merges with the left-most right-pole.To separate the poles we introduce a regularization parameter δ > 0 as Γ(z 1 ) → Γ(δ + z 1 ) assuming −1 < −δ < Re (z 1 ) < Re (z 2 ) < 0 (43) and analytically continue δ → 0 later.By applying the replacements to Eq. ( 40), we have where the Gauss summation formula has been used.After differentiating w.r.t.c and setting c → 0, we find Finally we can analytically continue δ → 0 as mentioned above; the r.h.s.becomes and therefore Throughout this example, z 2 can assume any value satisfying Eq. (43).It can, in particular, be an integration variable.Thus the two dimensional Mellin-Barnes integral of Eq. ( 42) can be reduced to a one dimensional integral.

Figure 3 :
Figure 3: One-loop master integrals.Solid and dashed lines represent massive and massless scalar propagators, respectively.The external (thin) lines are massless.The four master integrals which are not shown are obtained by crossing.

Figure 4 :
Figure 4: Two-loop planar master integrals.Solid and dashed lines represent massive and massless scalar propagators, respectively.The external (thin) lines are massless.The planar master integrals form (29) which are not shown are obtained by crossing.
The convergence behaviour may change under the replacements b 1 ↔ b 2 ↔ b 3 .By applying the condition (41) to 3 F 2 of Eq. (40), we obtain the condition Re (c − b 3 ) > 0, which is clearly not symmetric under the replacements b 1 ↔ b 2 ↔ b 3 .Thus the convergent domain, in terms of the space spanned by a 1 , a 2 , b 1 , b 2 , b 3 , can differ from expression to expression.