Two-loop QCD corrections to $b + \bar{b} \rightarrow H+H$ amplitude

We present the first results on the two-loop massless QCD corrections to the four-point amplitude $b+\overline{b} \rightarrow H+H$ in the five flavour scheme, treating bottom quarks as massless. This amplitude is sensitive to the trilinear Higgs boson coupling. Our two-loop result for this amplitude constitutes of purely virtual contributions to the next-to-next-to-leading order QCD predictions for the production of a pair of Higgs bosons at the Large Hadron Collider. We have implemented our two-loop results in a numerical code that can be used for further studies.


Introduction
Ever since the discovery of the Standard Model (SM) Higgs boson [1,2], one of the main objectives of the Large Hadron Collider (LHC) physics program has been to understand its properties. This involves the measurements of the Higgs boson couplings to the SM fermions and gauge bosons, its mass (m h ), its CP properties etc. Among these, the Higgs boson self couplings such as the trilinear (λ SM 3 ) and quartic couplings (λ SM 4 ) take prominence, which in the SM, can be unambiguously obtained from the Higgs boson mass. The SM Higgs potential, after the electro-weak symmetry breaking (EWSB), is given by where φ(x) denotes the Higgs field. v ≈ 246 GeV is the vacuum expectation value (vev) of the Higgs field and is fixed by the Fermi constant G F . The Higgs boson mass m h , is found experimentally to be approximately equal to 125 GeV and hence, the SM values for λ SM 3 and λ SM 4 are ∼ 0.13 and ∼ 0.03, respectively. However, presence of beyond the SM (BSM) physics scenarios can modify these couplings, which, in turn, suggests independent measurements of them. Any deviation from the SM values from the experimental measurements, could provide crucial information on the structure of the scalar potential and thus could constrain BSM physics scenarios [3]. Moreover, the measurement of λ SM 3 also provides a way to check that the EWSB follows from the simple Ginzburg-Landau φ 4 potential. The observable that can probe these couplings at the hadron colliders is the production of multiple Higgs bosons [4]. More precisely, the production of a pair of Higgs bosons can probe λ SM 3 but it is difficult to measure due to the smallness of its production cross section and the presence of a large QCD background. However, the study for the high luminosity LHC indicate that the Higgs boson pair production due to gluon fusion can predict λ SM 3 with O(1) accuracy [5][6][7][8].
A pair of Higgs bosons can be produced through several partonic channels, viz gluon fusion, vector boson fusion, associated production with a vector boson or a pair of heavy quarks. Among these channels, the gluon fusion channel is the most dominant one at the LHC. Being a loop-induced channel, gluon fusion gives a minuscule production cross section. Additionally, the large background of this channel makes its measurement experimentally challenging. Hence unless contributions from BSM physics enhance the production cross section, a measurement of this channel will require a considerable integrated luminosity. On the other hand, in such a scenario, the sub-dominant channels in the SM could possibly become interesting as they would receive substantial contributions from new physics. One such channel is the production of a pair of Higgs bosons in bottom quark annihilation. In certain supersymmetric models, namely the Minimal Supersymmetric SM (MSSM) [9], the bottom quark Yukawa coupling is enhanced w.r.t. the top quark Yukawa coupling, in the large tan β region, where tan β is the ratio of vev's of up and down type Higgs fields in the Higgs sector of the MSSM. Hence precise predictions for this channel is of high importance.
The dominant channel for Higgs boson pair production i.e. the gluon fusion channel, is mediated by a top quark loop. This was evaluated at leading order (LO) in perturbative QCD in [10][11][12] decades before. The next-to-leading order (NLO) contributions were obtained in [13] only in the infinite top mass limit, i.e. the top quark loop is integrated out resulting in an effective Lagrangian [13][14][15][16] of gluons and Higgs fields. There are several NLO results [17][18][19][20][21][22] considering finite top quark mass effects which finally led to the full NLO corrections with exact top quark mass dependence [23,24]. In all these works, it has been found that, with an inclusive K-factor close to 2, the QCD corrections at NLO level are as large as that observed for a single Higgs boson production. Hence, the next-tonext-to-leading order (NNLO) corrections were computed in [25], in the effective theory, followed by a soft plus virtual approximated NNLO cross section in [26]. Consecutively, the effect of leading top quark mass corrections also has been included in [27]. Finally a fully differential distribution has been obtained at the NNLO level in [28,29] and also threshold resummation at next-to-next-to-leading logarithmic (NNLL) level in [30,31]. In [32], a reweighting technique has been used to properly account for finite top mass effects at NNLO level. Recently the virtual contributions relevant for next-to-next-to-next-to-leading order (N 3 LO) QCD have also been computed in [33], within the effective theory.
While a plethora of work has been performed to reach ultimate precision for the gluon channel, the sub-dominant channels have not received much attention. Although, as mentioned earlier, in certain BSM physics scenarios, they become consequential. We are particularly interested in the bottom quark annihilation channel where the Higgs boson couples to bottom quarks through the Yukawa coupling (proportional to the mass of the bottom quark), and the bottom quark is massless otherwise [34][35][36]. For single Higgs boson produc-tion through this channel, various work is known up to NNLO [37][38][39][40][41][42] and N 3 LO [43][44][45][46] level in the variable flavour scheme (VFS) [47][48][49][50]. For the production of a pair of Higgs bosons, the NLO correction was first obtained in [51]. Later on, NLO corrections have been obtained for this channel considering several BSM scenarios [52][53][54]. For the latter, the bottom quark annihilation process dominates over the gluon fusion even at LO level. In addition, their NLO QCD corrections are not only sizeable but also larger than the supersymmetric QCD corrections. In order to stabilize the cross section with respect to higher order radiative corrections, NNLO corrections to this channel are desirable. In this paper, as a first step towards the full NNLO QCD corrections, we present the two-loop virtual contributions for the production of a pair of Higgs bosons in bottom quark annihilation channel.
There are two classes of diagrams (we call them Class-A and Class-B, see Sec. 2.4), that contribute at two loops. The vertex type of diagrams which belong to Class-A are already known up to three loops [44]. For the Class-B, the one-loop QCD corrections exist in the literature [51]. Here we compute the two-loop QCD corrections. We have studied the structure of infrared (IR) singularities and found that they are in agreement with the predictions by Catani [55]. The finite results expressed in terms of classical polylogarithms of weight up to 4 is used to study the numerical stability of the amplitudes over a wide range of allowed kinematical variables. The immediate application of our result is to include the universal soft gluon contributions from the real emission diagrams to obtain soft plus virtual contributions up to the NNLO level which will be presented elsewhere.
The paper is organized as follows. In Sec. 2, we discuss the Lagrangian, kinematics and the classes of diagrams that are relevant for our computation. Sec. 3 contains details of the computation, the ultraviolet (UV) renormalization and the structure of IR divergences. We devote Sec. 4 for the numerical evaluation of the amplitude over a wide kinematic region before we conclude in Sec. 5.

Theory
At the LHC, the dominant channel for the production of a pair of Higgs bosons is the gluon fusion. In addition, there are several sub-leading channels that contribute to the production. We consider one of these channels, namely the production through the bottom quark annihilation process. Since the LO and the NLO [51] QCD effects have already been studied in the literature, as a first step towards the computation of the full NNLO QCD corrections, we evaluate two-loop virtual contributions to the production of a pair of Higgs bosons in this channel. Note that we further need to compute contributions from real emission subprocesses to obtain IR safe observables at the NNLO level. These pure virtual corrections contribute to both the inclusive as well as the differential observables. These results along with the process independent soft gluon contributions, can give us the first result at the NNLO level in the threshold limit, i.e., when the invariant mass of the pair of Higgs bosons approaches the partonic centre of mass energy. We reserve these applications for future publication.
We use the regularized version of the QCD Lagrangian throughout. The regularization scheme that we use, is the dimensional regularization (DR), in which all the fields and couplings of the Lagrangian and the loop integrals that appear in the Feynman diagrams are analytically continued to d = 4 + spacetime dimensions. In addition, we perform traces of Dirac γ matrices in d-dimensions.

The Yukawa interaction
We begin by reviewing the theoretical framework for the production of a pair of Higgs bosons via bottom quark annihilation at hadron colliders. The interaction part of the Lagrangian that is responsible for the production is given by, is the bottom quark mass and v the vev of the Higgs field. In the SM, the ratio of the top quark Yukawa coupling (λ t ) and the bottom quark Yukawa coupling (λ b ) is found to be approximately 35 i.e. λ t /λ b ≈ . 35. In addition, the bottom quark flux in the proton-proton collision is much smaller than the gluon flux. Hence, the contribution from this channel is sub-dominant as compared to the gluon fusion channel. However, in the MSSM [9], this ratio depends on the value of tan β which can increase the contribution resulting from the bottom quark annihilation channel. At LO, where h is the SM like light Higgs boson, H and A are the heavy and the pseudoscalar Higgs bosons, respectively. The parameter α is the angle between weak and mass eigenstates of the neutral Higgs bosons h and H. Since, the bottom quark mass is much smaller than the other energy scales that appear at the partonic level, we set the former to zero except in the Yukawa coupling in perturbation theory [34][35][36]. In particular, the finite mass effects from the bottom quarks are found to be suppressed by the inverse power of mass of the Higgs boson. The number of active flavours is taken to be n f = 5 and we work in the Feynman gauge.

Kinematics
We compute all the relevant one-and two-loop amplitudes in perturbative QCD that contribute to the annihilation of bottom quarks into a pair of Higgs bosons. The scattering process is given by which satisfy the relation s + t + u = 2m 2 h . For convenience, we use the dimensionless variables x, y and z defined [56] as follows The variables x, y and z satisfy The final result will be expressed in term of logarithms and classical polylogarithms which are functions of these scaling variables.

General structure of the amplitude
The external states for the process given in Eq. (2.4) involve two fermions and two scalars, hence the most general structure of the amplitude can be parametrized as where the coefficients C m ≡ C m (x, y, z) with m = 1, 2 are scalar functions. In color space, the amplitude is diagonal in the indices (i, j) of the incoming quarks. Since, we are interested in higher order QCD corrections, we have used symmetries such as Lorentz covariance, parity and time reversal invariances to parametrize the amplitude. In addition we have dropped those terms that vanish when the bottom quarks are massless. The coefficients C m , m = 1, 2, can be determined from the amplitude A ij by using appropriate projection operators denoted by P(C m ), i.e., where the sum includes spin, flavours and colours of the external fermions; N is the number of colours in SU(N) gauge theory. In d-spacetime dimensions, the projectors that satisfy P(C m )T m = 1 and P(C m )T n = 0 ∀ m = n, are found to be Since the application of projection operators on the amplitude gives only Lorentz scalar functions, the algebraic manipulations with loop integrals become straightforward. The square of the amplitudes, that contributes to the total cross-section, can now be obtained from the coefficients C 1 and C 2 using Note that these coefficients are in general complex due to the Feynman loop integrals. We expand the amplitude A ij as well as the coefficients C m in powers of the strong coupling constant defined by a s = g 2 s (µ 2 R )/16π 2 , where g s is the renormalized strong coupling constant and µ R is the normalization scale: and consequently A (l) Our next task is to compute these coefficients C

Classification of Feynman diagrams
At LO, only three Feynman diagrams contribute, out of which one contains single Yukawa and trilinear couplings, and the remaining ones are quadratic in the Yukawa coupling. We denote the former by Class-A and the latter diagrams by Class-B. The same classes of diagrams contribute beyond LO. We elaborate on these classes of diagrams below: • Class-A: It contains diagrams where an off-shell Higgs boson produced in the bottom quark annihilation process decays to a final state containing a pair of on-shell Higgs bosons (H * → HH) and is proportional to λ SM 3 λ b . They are shown in Fig. 1. Note that the decay part of the amplitudes does not get any QCD corrections, however the initial states do get. These corrections are identical to those that contribute to the amplitudes for producing a single Higgs boson in bottom quark annihilation. The latter is known up to three-loop level in QCD [44].

Computational details
It is easy to see from the form of T i in Eq. (2.8), that only the Class-A diagrams contribute to C 1 and the Class-B to C 2 . Note that the Class-A diagrams are already computed to three loops in QCD [44]. Hence in this section, we briefly discuss how the scalar function C 2 in Eq. (2.9) is computed order by order in perturbation theory. As we mentioned, we use the dimensional regularization, in which the spacetime dimensions are taken to be d = 4 + and perform traces of Dirac γ matrices, contraction of Lorentz indices in d-dimensions. For convenience, we work with the bare form of the Lagrangian and evaluate the coefficient C 2 in powers of bare coupling constantâ s , whereâ s =ĝ 2 s /16π 2 ,ĝ s being the dimensionless strong coupling constant. Beyond LO, one-and two-loop amplitudes containing massless quarks, anti-quarks and gluons develop IR divergences in addition to UV ones. There are two types of IR divergences, viz soft and collinear divergences. The soft ones are due to soft gluons and the collinear ones arise due to massless quarks and gluons. Dimensional regularisation regulates both these divergences in addition to UV divergences.
We have used QGRAF [57] to generate the Feynman diagrams at every order in the strong coupling constant. Beyond one-loop, large number of Feynman diagrams contributes to the amplitude. We find that there are 2 diagrams at the Born level, 10 diagrams at one-loop and 153 diagrams at two-loop level. We then multiply these amplitudes with the projection operator P(C 2 ) defined in Eq. (2.10) to obtain the scalar function C 2 . Substitution of Feynman rules and computation of various traces involving Dirac and Gell-Mann matrices, are done using in-house routines that use publicly available packages such as FORM [58] and Mathematica. At this stage we end up with a large number of one-and two-loop Feynman integrals. The projection operators guarantee that all the tensor integrals are converted to scalar integrals. We rearrange all the Feynman integrals into a few chosen integral families through shifting of loop momentum. To achieve this, we use the package Reduze2 [59]. At one-loop, the following three integral families can accommodate all the Feynman integrals where, i takes one of the values {1, 2, 3} whose elements are arranged cyclically. A typical two-loop topology contains at most seven propagators. However, there are nine different Lorentz invariants (k i .k j , k i .p j ) which can appear in the numerator of an integral. Hence, we introduce two auxiliary propagators in each of the two-loop integral families. The following two sets describe the six integral families that we use at two-loops, {P 0 , P 1 , P 2 , P 1:i , P 2:i , P 1:i,i+1 , P 2:i,i+1 P 1:i,i+1,i+2 , P 2:i,i+1,i+2 } , Here, This large number of Feynman integrals belonging to different integral families can be written in terms of a smaller set of integrals, so-called master integrals (MIs). This can be achieved by using the integration-by-parts (IBP) [60,61] and the Lorentz Invariance (LI) [62] identities, which are implemented in the Mathematica based package LiteRed [63]. Finally, we obtain 10 and 149 MIs at one-and two-loops, respectively. The resulting set of MIs is systematically mapped on to those evaluated in [56,64] as Laurent series in up to the required order. Finally, substituting the results of MIs from [56,64], we obtain the two-loop result for the coefficient C 2 . Both UV and IR divergences appear as poles in at every order inâ s . In the next section, we demonstrate how the renormalization of the strong and the Yukawa couplings render these coefficients UV finite leaving only IR divergences.

Ultraviolet renormalization
The scalar function C 2 computed in powers of the bare coupling constantâ s contains both UV and IR divergences. Note that the entire amplitude is proportional to the square of λ b , the bare Yukawa coupling. We use the modified minimal subtraction (M S) scheme to perform the UV renormalization of the amplitudes. In this scheme, the renormalized strong coupling constant a s is related to the bare strong coupling constant,â s , through the renormalization constant Z a s (µ 2 R ), at the renormalization scale µ R aŝ where Z a s (µ 2 R ), up to two-loops is given by Here, S ≡ exp[(γ E − ln 4π) 2 ] is the phase-space factor in d-dimensions, γ E = 0.5772... is the Euler-Mascheroni constant and µ 0 is an arbitrary mass scale introduced to makeĝ s dimensionless in d-dimensions. The constants β 0 and β 1 are the coefficients of beta function which, for n f light quark flavours, are found [65][66][67][68][69] to be are the Casimirs of SU(N) group and T F = 1/2. Similar toâ s , the renormalization of the Yukawa coupling constantλ b leads to renormalized λ b (µ 2 R ) at the renormalization scale µ R througĥ where the coefficients Z (i) λ,j are given by The perturbative expansion of the amplitude for the aforementioned process in terms of the bare strong and Yukawa couplings is given by ij is the l th loop unrenormalized amplitude. Similarly, the coefficient C 2 replicates similar perturbative expansion of the following form, In terms of the renormalised couplings, the coefficient C 2 takes the form We obtain the coefficients C 2 . (3.11) These constants C (l) 2 , l = 0, 1, 2, that result after performing the renormalization of the strong and the Yukawa couplings, are UV finite. However they are sensitive to both soft and collinear divergences which will be the topic of our next section. These soft and collinear divergences show up in terms of poles in .

Infrared divergences and their factorization
The UV finite amplitudes still contain divergences resulting from soft and collinear regions of the loop integrals. They result from soft gluons and massless collinear quarks and gluons in the loops. In the physical observables, the soft and the collinear divergences coming from the final states of the virtual diagrams cancel against those resulting from the phase space integrals of the real emission processes. Due to the Kinoshita-Lee-Nauenberg (KLN) theorem [70,71], the cancellation takes place order by order in perturbation theory. While the soft divergences cancel fully, the collinear divergences resulting from initial massless states, do not cancel at the subprocess level. Thanks to the collinear factorization theorem [72] these initial state collinear divergences can be factored out in a process independent way and absorbed into the bare parton distribution functions. This procedure is called mass factorization which is also a consequence of KLN theorem applied at the hadronic level. While all these IR divergences that appear in the amplitudes do not pose any problem for the physical observables, they provide valuable information about the universal structure of the infrared divergences in the QCD amplitudes. In fact, it can be shown that these divergences systematically factor out from the amplitudes to all orders in perturbation theory [73,74]. These factored IR divergences demonstrate the universal structure in terms of certain soft and collinear anomalous dimensions. An elegant proposal was put forth by Catani who predicted IR pole structure of the amplitudes up to two-loop level in non-abelian gauge theory [55]. He demonstrated that the n-particle QCD amplitudes factorize in terms of the universal IR subtraction operator denoted by I. This I-operator has a dipole structure [55] containing process independent universal cusp and collinear anomalous dimensions. Thanks to the wealth of results from two-loop calculations of the three-parton qqg amplitudes [75] and 2 → 2 scattering amplitudes [76][77][78], that involve nontrivial color structures [78,79], the I-operator is completely known up to two-loop level. In [80], the authors provide further insight on the factorization and resummation properties of QCD amplitudes in the light of Catani's proposal and demonstrate a connection between divergences governed by soft and collinear anomalous dimensions, see also [81,82]. There have been several efforts [83,84] to determine the structure of I-operator beyond two-loop level. Following [55] we express one and two-loop UV renormalized amplitudes in terms of the I-operator as (3.12) The matrix elements of the subtraction operator for the bottom quark, I b are given by [80] given by as follows According to the proposal by Catani, the coefficients C (i),fin 2 ( ) should be free of IR divergences and hence are finite as → 0. Since the resulting expression at two-loops level, C (2),fin 2 ( ) is quiet lengthy, we had to simplify the expression first at the color factor level and then for each color factor, terms of uniform transcendentality were further simplified. We find that our final result is in accordance with Catani's predictions for the IR poles, which serves as an important check on the correctness of our computation.
It is interesting to observe that the singlet contributions which are proportional to the color factor C F n b T F , for n b = 1, develops IR divergences at the intermediate stages of the computation. However at the end, all the IR singularities cancel among themselves contributing only to the IR finite part. This is consistent with the IR pole structure predicted by Catani. The resulting finite constant C (2),fin 2 that results after subtracting the IR divergences using Catani's I-operators is too lengthy to be presented here and hence attached as ancillary file in Mathematica format.

Numerical predictions
The finite coefficients, C (i),fin 2 , i = 1, 2, obtained in Eq. (3.12) contain multiple classical polylogarithms, which are functions of the scaling variables x and y. These polylogarithms can be attributed to different transcendental weights, the property that we use to simplify the two loop coefficient, C  perform a numerical evaluation using Mathematica for a wide range of scaling variables. More precisely, in Fig. 4 we plot the real and the imaginary parts of the coefficient C (2),fin 2 as functions of the scaling variable x for different values of cos(θ), where θ denotes the angle between one of the initial state fermions and the Higgs boson in the centre of mass frame of incoming states. We consider m h = 125 GeV and the renormalization scale as µ 2 R = m 2 h /2. In addition, we normalise the coefficient with the factor m 2 h . The amplitude is anti-symmetric under cos(θ) → − cos(θ), as expected for a purely fermionic amplitude. Since this symmetry has not been used in the setup of the calculation, it serves as a strong check on our results. Our expression contains polylogarithms that are multiplied by large rational coefficients, hence we encounter numerical instabilities during the evaluation. To avoid this, we evaluate the polylogarithms at double precision while setting the rational coefficients at higher precision. From the Fig. 4, we observe a stable behaviour for the real and imaginary parts for the range of parameters considered. In addition, the dependence of the coefficient near the phase space boundary x = 0 is displayed in the insets. The simplified analytical results and their numerical implementation are provided as separate ancillary files.

Conclusion
The extraction of the trilinear coupling of the Higgs bosons provides valuable information on the shape of the Higgs potential. One of the most important observables sensitive to this coupling is the production of a pair of Higgs bosons at the LHC. Among various partonic channels that contribute to this process, gluon fusion is the dominant one which is well studied both in effective theory as well as in the full theory. In the effective theory, top quarks are integrated out. As the precision at the hadron collider improves, it is important to incorporate other sub-dominant channels to the production mechanism. In this paper, we have considered one such channel, namely the production of a pair of Higgs bosons in the bottom quark annihilation which is also sensitive to the trilinear coupling. Both LO and NLO QCD contributions exist in the literature and the present paper presents first results on the two-loop virtual contributions at the NNLO level. There are two classes of diagrams that contribute at two-loops, of which, the vertex type of diagrams belong to Class-A are already known up to three-loops and hence we evaluate only the Class-B diagrams. Our results are expressed in terms of classical polylogarithms of weight up to 4. We observe that the infrared poles of the amplitudes are in agreement with the predictions by Catani. We have studied the numerical stability of the coefficient C (2),fin 2 over the range of x and cos(θ) required for further phenomenological studies. The immediate application of our result is to include the universal soft gluon contributions from the real emission diagrams to obtain the soft plus virtual contributions up to the NNLO level which will be presented elsewhere.