Two-loop amplitudes for di-Higgs and di-pseudo-Higgs productions through quark annihilation in QCD

Through this article, we present the two-loop massless QCD corrections to the production of di-Higgs and di-pseudo-Higgs boson through quark annihilation in the large top quark mass limit. Within dimensional regularisation, we employ the non-anticommuting γ5 and treat it under the Larin prescription. We discover the absence of any additional renormalisation, so-called contact renormalisation, that could arise from the short distance behaviour of two local operators. This finding is in corroboration with the operator product expansion. By examining the results, we discover the lack of similarity in the highest transcendentality weight terms between these finite remainders and that of a pair of half-BPS primary operators in maximally supersymmetric Yang-Mills theory. We need these newly computed finite remainders to calculate observables involving di-Higgs or di-pseudo- Higgs at the next-to-next-to-leading order. We implement the results to a numerical code for further phenomenological studies.


Introduction
The detailed exploration of the discovered scalar resonance [1,2] with a mass of 125.09 ± 0.24 GeV [3] at the Large Hadron Collider (LHC) is underway with spectacular success. After establishing the coupling of the Higgs boson to all electroweak gauge bosons [3], top [4,5] and bottom quarks [6,7], tau lepton [8,9], now LHC has its focus on the measurement of its coupling with muon [10,11]. Reconstructing the Higgs potential experimentally and consequently pinning down the mechanism of electroweak symmetry breaking is one of the most important objectives in current times. This requires the precise determination of Higgs self-coupling which can be accessed directly through the production of Higgs boson pair. While the determination of quartic self-coupling is difficult at the LHC, the probe of trilinear coupling is certainly possible. Through the current run [12][13][14], the limits on the latter is improved significantly. The theoretical value of λ, which dictates the self-coupling, can be determined from the available precise results of the mass and vacuum expectation value of the scalar Higgs boson. However, it is important to verify it independently through experiment.
In the Standard Model (SM), the Higgs boson pair is dominantly produced through gluon fusion containing a top quark loop. The production cross-section at leading order (LO) in QCD has been computed exactly in refs. [15][16][17][18]. Due to the presence of massive JHEP01(2022)189 top quark loop, computing it exactly beyond the LO is a nontrivial task. To reduce the complexity while capturing the dominant contribution, the Higgs Effective Field Theory (HEFT) has been designed by integrating out the top quark loop by treating its mass as infinitely large. In ref. [19], the next-to-LO (NLO) corrections in HEFT were first computed by keeping the LO exact. Later, almost after two decades, the NLO corrections by keeping the exact top quark mass were achieved in refs. [20][21][22][23]. In refs. [24][25][26][27], the parton shower was also included to make the prediction more robust. The threshold resummation effect was investigated in ref. [28].
In HEFT, the NNLO QCD corrections were computed in refs. [29][30][31][32][33]. The relevant Wilson coefficients were presented in ref. [31]. By keeping the exact top quark dependence up to NLO and in the double real radiation, the NNLO correction was further improved in ref. [34] which made use of the result from ref. [33]. These results were further improved by including the soft-gluon resummed contribution in ref. [35]. The differential results at NNLO can be found in refs. [36][37][38]. Quite recently, the focus has moved towards the N 3 LO corrections within HEFT. The computation of the two-loop virtual correction [39] and the four loop matching coefficient for the effective coupling of two Higgs bosons and gluons [40,41] has been achieved. Within HEFT, the N 3 LO results were computed in refs. [42,43], where the result of ref. [43] is NLO-improved employing the results of refs. [24,26]. Recently, in ref. [44], the NNLO cross-section has been improved by computing the higher order terms in top mass expansion.
Although the gluon initiated process for di-Higgs production is dominant, there is a need of incorporating the subdominant channels to make the theoretical predictions as precise as possible. The first step towards this direction was successfully attempted in ref. [45], where the NNLO correction to the process initiated by the bottom quark annihilation was computed. The Higgs boson was assumed to couple to the bottom quark through non-zero Yukawa coupling, otherwise the latter was treated as a massless entity. In this article, we address the computation of two-loop QCD correction to di-Higgs production through quark annihilation in HEFT. On top of its phenomenological relevance, we intend to discover whether there is a need for any additional renormalisation arising from the short distance behaviour of two composite operators. For instance, for the gluon initiated process, its necessity was discovered in refs. [39,46]. In the former reference, it was found out from the operator product expansion (OPE). In this article, we want to find whether there is a need for any additional contact renormalisation arising from short distance behaviour for the quark initiated channel. Moreover, we intend to determine whether any subtle issue shows up owing to working in an effective theory, an example of this can be found out from ref. [47] where one of the authors demonstrated the need for an additional four-point effective vertex.
In addition to the scalar Higgs, the exploration of the Higgs sector at the LHC includes the search for a CP-odd Higgs boson which is often called a pseudo-scalar. In certain beyond the SM (BSM) scenarios, for instance the minimally supersymmetric SM (MSSM), the Higgs sector contains a pseudo-scalar in addition to two neutral and two charged scalars. It is important to have theoretical predictions associated with a CP-odd Higgs boson as precise as with the CP-even counterpart. While the observables associated with the latter JHEP01(2022)189 have been computed to a very high accuracy, the progress involving di-pseudo-scalar is not up to the requirement. One of the reasons is the handling of chiral quantity in dimensional regularisation where space-time dimension is analytically continued to d = 4−2 with being a complex number. In ref. [48], some of us demonstrated how to simplify the calculation of radiative corrections involving a chiral quantity. For a single CP-odd Higgs boson, the production cross-section is available to N 3 LO under soft-virtual approximation [49][50][51][52][53][54]. For the di-pseudo-scalar, the NLO in heavy top limit was obtained in ref. [19]. The first step to go beyond the NLO was attempted in ref. [55] where the two-loop QCD correction for the production of di-pseudo-scalar through gluon fusion was achieved. In this article, we address the computation in the quark annihilation channel in the heavy top limit up to two loops. In the process, we discover the absence of any contact renormalisation, consistent with the analysis through OPE in ref. [56].
In section 2, we discuss the kinematics of the four-point amplitudes and the dimensionless variables which are introduced for future reference. We discuss the effective Lagrangian in the heavy top limit in section 3. A short discussion on handling γ 5 is also included. In section 4, we introduce the form factor decomposition and consequently the projector that is employed to get the form factors up to two-loop in QCD. We briefly discuss the method of computing the loop amplitude that is adopted in this article. The ultraviolet renormalisation is performed in section 5. We discuss the infrared factorisation in terms of universal subtraction operator and some properties of the ultraviolet renormalised form factors in section 6. Finally we perform the numerical analysis of the results in section 7. Then we summarise and draw the conclusion in section 8.

Preliminaries
We consider the production of di-scalar (φ) and di-pseudo-scalar (φ) Higgs boson through quark annihilation where we denote the four-momentum by p i satisfying the on-shell conditions p 2 1 = p 2 2 = 0 and p 2 3 = p 2 4 = m 2 . We represent quark (anti-quark) by q(q). The parameter m denotes the mass of the scalar or pseudo-scalar depending on the process. We define the Mandelstam variables through which, owing to the momentum conservation, satisfy For future requirements, we introduce three dimensionless variables through

JHEP01(2022)189
These partonic channels contribute to the production of di-scalar and di-pseudo-scalar Higgs boson at a hadronic collision. In this article, we calculate the scattering amplitudes within the framework of the heavy-top effective theory which we now turn into.

Theoretical framework under heavy top limit
In the heavy-top effective field theory, that is obtained by integrating out the top quark appearing in the loop, the interacting Lagrangian density describing the one and two scalar Higgs bosons is given by The quantity G µν a denotes the gluon field strength tensor and v = 246 GeV is the vacuum expectation value of the scalar Higgs field. In moving from exact to the effective theory where we work in n f = 5 light flavor, all the top quark mass dependence is encapsulated within the Wilson coefficients C H and C HH . These are determined by matching the effective n f = 5 and full (n f + 1) flavor theories. Furthermore, this matching can be done order by order in perturbation theory which in powers of renormalised strong coupling constant, a s (µ 2 R ) = α s (µ 2 R )/(4π), are found to be [31,40,41,[57][58][59]] In the aforementioned results, m t is the top quark mass in MS scheme determined at the renormalisation scale µ R . We denote the quadratic Casimirs of SU(n c ) group in the fundamental and adjoint representations by C F = (n 2 c − 1)/(2n c ) and C A = n c , respectively. The constant T f equals 1/2. We skip presenting the expression of C HH as it is not required in our current calculation, as discussed in the upcoming section.
In an analogous way, in this effective theory the pseudo-scalar Higgs boson (φ) can also be described through a Lagrangian containing a gluonic (O G ) and a quarkonic (O J ) operator. The effective Lagrangian [60] reads as where the operators are defined in terms of gluon field strength tensor and light quark wave function (ψ) as The transition from exact to effective theory introduces two Wilson coefficients, C G and C J which can be expanded in powers of strong coupling constant:

JHEP01(2022)189
The absence of higher-order terms beyond one-loop in C G is guaranteed by the Adler-Bardeen theorem [61]. We denote the Fermi constant by G F and mixing angle in a generic two-Higgs doublet model by β. The superscript of C J indicates the corresponding perturbative order.
In the upcoming section, we discuss how we compute the di-scalar and di-pseudo-scalar amplitudes within this framework of effective theory. The presence of the chiral quantity, Dirac's γ 5 (and Levi-Civita symbol µνρσ ), which are inherently defined in 4-dimensions faces the issue of translating properly into d-dimensions in dimensional regularisation (DR) [62,63]. One of the most popular ways of handling it to define it following the 't Hooft-Veltman [62] and Breitenlohner-Maison [64] scheme as This definition has many profound implications which we need to address. Firstly, it breaks the hermiticity property of the operator which is cured by replacing [65,66] instead of directly using the definition (3.6). The contraction of Levi-Civita symbols is done in d-dimensions using where the indices of the metric tensor on the right-hand side are treated in d dimensions [67].
The anti-symmetrisation of the Lorentz indices is denoted through [ ].

Four-point amplitudes
Performing the Lorentz covariant decomposition, we can write down the general form of the amplitude for these scattering processes as where h can either be φ orφ. By applying on-shell Dirac equation and keeping in mind the conservation of chirality along the massless quark line along with the parity even final states, it is straightforward to guess the form of Γ h as The form factor F h can be calculated order by order in perturbation theory. At the leading order, both of these processes are one-loop induced in HEFT since these start at two-loop in exact theory. Our goal is to compute the two-loop QCD contributions to F h in the HEFT.
To compute the form factor (FF), we can apply an appropriate projector on A ij,h , which in this case comes out to be   We start the computation by generating the set of Feynman diagrams in heavy quark effective theory using QGRAF [68] that provides the expression in a symbolic form. For the production of di-Higgs through quark annihilation, we categorise the diagrams into class-A and class-B depending on whether those are proportional to C HH or C 2 H , respectively. The diagrams at the LO belong to class-A, as shown in figure 1. The class of diagrams belonging to this category vanish to all orders since this part of the amplitude is proportional to the mass of external quarks which is zero in massless QCD and hence only the class-B diagrams contribute. At NLO, there are 5 diagrams, and at NNLO there are 143. We have shown some sample diagrams of this category in figure 2. Before moving to pseudo-scalar, let us mention that we do not need to compute the set of Feynman diagrams where a scalar Higgs boson decays to two other Higgs through the triple Higgs coupling since this set of diagrams give a vanishing contribution, similar to class-A.
For the processφφ at LO, there are 2 diagrams (see figure 1) containing the four-point vertex ggφφ. These diagrams also vanish to all orders for the same reason as for the di-scalar. We categorise these as class-A. For higher orders, we classify the diagrams into three types.  we find the appropriate loop momentum shifts to recast each Feynman diagram beyond the tree level to a form belonging to one of the integral families. We have 3 and 6 integral families at one-and two-loop, respectively. The resulting expressions contain Feynman loop integrals which are reduced to a set of master integrals through integration-by-parts (IBP) [72] identities. To perform the IBP reductions, we use LiteRed [73,74] 1 in conjunction with Mint [75]. The latter assists to determine the number of master integrals. We get 10 and 149 master integrals at one-and two-loop, respectively. Using the results of the master integrals [76,77] in terms of multiple polylogarithms, we get the bare form factors as Laurent series in dimensional regulator . To do the series expansion efficiently, we make use of the finite field reconstruction through FiniteFlow [78].

Ultraviolet divergences and operator renormalisation
Amplitude beyond LO often suffers from ultraviolet (UV) divergence. Since we are working in massless QCD, the first part of the cure consists of renormalising the strong coupling constant in MS scheme through [79]  the full UV divergences. Since this is a property associated to the operator itself, we call it as operator renormalisation. For di-scalar, the renormalised Lagrangian, [L φ eff ] R , is related to the bare part through with the operator renormalisation constant in MS scheme [80][81][82] Incorporating the aforementioned renormalisations, we get the UV renormalised FF where F φ appearing on the right-hand side is the one that we calculate by applying the projector (4.3) on the set of Feynman diagrams and performing coupling constant renormalisation through (5.1). By rewriting F φ in (5.4) after factoring out the Wilson coefficient, we get corresponds to n-th loop bare form factor obtained by evaluating the Feynman diagrams in figure 2. By performing the required renormalisations and expanding all quantities in powers of a s as we arrive at the fully UV renormalised FF up to NNLO expanded in powers of renormalised a s In the aforementioned equation, the renormalised FFs at NLO and NNLO are denoted by [F (1) φ ] R and [F (2) φ ] R , respectively. We provide the corresponding finite remainders in the supplementary material of this work. The coefficients Z Before moving to the next topic, we point out that unlike the case of di-scalar Higgs boson production through gluon fusion [39,46], we do not need any additional renormalisation arising from the contact terms because of the product of two composite operators at a short distance. Although, in principle, one could expect it to be present.

JHEP01(2022)189
Pseudo-scalar Higgs. In the case of di-pseudo scalar Higgs boson, we also have to perform coupling constant and operator renormalisation in MS scheme, for example. However, for the operator O J , which is constructed out of the singlet axial vector current ψγ µ γ 5 ψ, there is another phenomenon takes place due to the fact that our scheme of chiral quantity depicted through (3.6) violates the defining anti-commuting property of γ 5 in the evanescent 4 − d = 2 dimensional space. Consequently, this gives rise to additional terms of UV-divergence origin in the dimensionally regularised amplitudes which need to be removed [83,84] manually. This is manifested as a loss of correct chiral Ward identity. To cure this phenomenon, we need to introduce additional UV renormalisation constants [65,[84][85][86] in addition to the pure MS renormalisation constant. Therefore, the renormalised operator [O J ] R can be written as where the corresponding bare quantity O J is defined in (3.4). We denote the total operator renormalisation constant of the axial vector current by Z J consisting of a pure MS renormalisation factor Z s MS and an additional finite renormalisation part Z s 5 . The latter is needed to restore the correct form of the axial Ward identity which is known to be anomalous.
The Z GG is available to O(a 4 s ) in refs. [52,85,87,89,91] which coincides with the strong coupling constant renormalisation Z a , a result claimed to hold to all orders in ref. [90] and proved explicitly in ref. [92]. The other constant Z GJ which starts at O(a s ) is computed to O(a 3 s ) in refs. [52,56,85]. For the di-pseudo-scalar production, effectively, we have to renormalise three different combinations of operators corresponding to three categories of Feynman diagrams which are Inserting these into external quark states, we arrive at the following form of the UV renormalised FFs: where these individual components are given by

JHEP01(2022)189
The FF, Fφ ,G 1 G 2 , appearing on the right-hand side of the above equations are obtained by applying the projector ( 14) The notations are exactly similar to the ones in (5.6). The lower limit in each of these expansions is written as per their leading order expressions. By incorporating these, we arrive at the following fully UV renormalised FF, introduced in (5.11) and (5.12): By adding these three contributions order-by-order, we get the final UV renormalised FF at NLO and NNLO, [F (1) φ ] R and [F (2) φ ] R , respectively, with

JHEP01(2022)189
We provide the finite remainder of FF at NLO and NNLO as supplementary material. Interestingly, similar to the case of di-scalar Higgs boson, we find that we do not need any additional renormalisation arising from the contact terms owing to the product of operators at a short distance. For the gluonic channel, its absence is shown both from the operator product expansion [56] and direct two-loop calculation [55]. In the next section, we talk about the infrared structure of the UV renormalised FF and subsequently the behaviour of finite remainders.

Infrared factorisation and finite remainders
The UV renormalised FF still contains divergences originated from the soft and collinear configurations of loop momentum, which can be parameterised in terms of a universal subtraction operator I (1) ( ) as [93,94] [F h,fin ] R at → 0 is the finite remainder. The IR subtraction operator for the processes under consideration in conventional dimensional regularisation (CDR) is given by We check that the renormalised FF do exhibit this predicted infrared structure. The double pole is checked analytically, whereas the single one is done numerically due to the complexity of the algebraic expression. We check it on several set of kinematic points. This serves as a very strong check on the correctness of our computation. By subtracting the IR singular terms from the UV renormalised FF, we obtain the finite remainders for all the FF. These finite remainders will contribute to observables, say the inclusive cross-section, at the NNLO for the process of di-Higgs or di-pseudo-Higgs boson productions through quark annihilation at hadron collider. In the upcoming subsections, we discuss the kinematic symmetry associated with the interchange of two final state particles and the behaviour of leading transcendental weight terms.

Permutation symmetry
The UV renormalised FF or equivalently the finite remainder of the FF for both processes should exhibit the kinematic symmetry associated with the interchange of two final state particles. Due to the presence of an initial state quark and anti-quark pair, the FF are expected to be antisymmetric upon interchanging either the momenta of two incoming (p 1 ↔ p 2 ) or two outgoing (p 3 ↔ p 4 ) particles. This gets translated to an interchange of the Mandelstam variables, t ↔ u, or the scaleless parameters, y ↔ z, as introduced in (2.4). This implies We check that the resulting FF for both the processes do exhibit the expected kinematic symmetry by verifying the above equation at one and two-loop. Although this checking is

JHEP01(2022)189
quite simple for one-loop, this is quite involved at two-loop owing to the presence of special functions such as Li 2,2 . With the help of a C++ code [95], we evaluate this special function numerically and verify the kinematic symmetry. This serves as a strong check on the finite part of the FF.

Behaviour of leading transcendentality weight terms
Understanding the connection among various gauge theories through studying the analytic structures of scattering amplitudes is an active area of investigation. In particular, the connection between maximally supersymmetric Yang-Mills theory (mSYM) and QCD is of fundamental importance. In refs. [96][97][98][99], the anomalous dimensions of leading twisttwo-operator in mSYM are found to be identical to the highest transcendental (HT) weight counterparts in QCD, and consequently, the principle of maximal transcendentality (PMT) is conjectured. This principle, albeit an observational fact, states that algebraically the HT weight terms of certain quantities in mSYM and QCD are identical upon changing the representation of fermions in QCD from fundamental to adjoint through C A = C F = 2n f T F = n c with T F being the normalisation factor in fundamental representation. The HT weight terms of quark and gluon two-point FF with a single operator insertion in QCD are found to be identical to the scalar FF of half-BPS primary in mSYM to three loops [100] (up to some normalisation factor of 2 L with L being the loop). The diagonal terms of the two-point pseudo-scalar [52] and stress energy tensorial FF [101,102] also obey the conjecture. Moreover, the three-point scalar and pseudo-scalar FF are also found to respect the PMT [103][104][105][106][107][108][109][110][111][112]. In ref. [113], the four loop collinear anomalous dimension in the planar mSYM is computed employing this conjecture. Beyond form factors, even some observables, such as the asymptotic limit of energy-energy correlator [114] and soft function [99] are also observed to be consistent with PMT.
As the domain of validity of this conjecture is still not well established, the natural question arises whether the PMT carries over to more general class of FF and correlation functions. Already some instances of the violation of this conjecture are found. The three-point quark [115] and gluon [116] FF with a single stress energy tensorial operator insertion [117], and Regge limit of amplitudes [118] do not obey the principle. In refs. [119,120], it has been shown that the conjecture also fails to hold true for two-point FF with two operators insertion. In view of this, we investigate whether there is any similarity between the HT weight terms of the scalar and pseudo-scalar FF as computed in this article and that of mSYM theory. We find they do not share the same HT weight terms at one-as well as two-loop. For instance, the one-loop finite remainders contain three weight 2 polylogarithms, Li 2 (−x), Li 2 (1/(1 + y)) and Li 2 (1/(1 + z)), with the coefficients , F cotβ for the pseudo-scalar Higgs boson production. The corresponding one-loop two-point scalar FF in mSYM [120] with two operator insertion of half-BPS primary contain different kinematic coefficients.

Numerical analysis of the results
In this section, we study the behaviour of the finite remainders [F (2) φ,fin ] R and [F (2) φ,fin ] R at two loop for the production of di-Higgs and di-pseudo-Higgs boson through quark annihilation respectively in the large top quark mass limit. The results are expressed in terms of the classical polylogarithms with the scaling variables x and y being their arguments. In figure 4 (5), we plot the real and imaginary parts of [F with respect to the square of the corresponding Wilson coefficient C F cotβ. It is evident from the plots that the amplitudes are anti-symmetric under cos(θ) → − cos(θ), as expected for a purely fermionic amplitude. As a result, the amplitude vanishes when cos(θ) = 0. Since this anti-symmetry property has not been used in our calculation, it, in fact, serves as a strong check on our results. We observe that the numerical values of real as well as imaginary parts of both [F (2) φ,fin ] R and [F (2) φ,fin ] R and corresponding to different values of cos(θ) start to converge as x → 1.

Conclusions and outlook
In this article, we present the two-loop four-point amplitudes for the productions of two scalar and two pseudo-scalar Higgs bosons through quark annihilation in Higgs effective field theory. The HEFT is obtained by treating the top quark as infinitely heavy. We perform the calculation under dimensional regularisation. The chiral quantities, Dirac's γ 5 and ε µνρσ , are treated through Larin's prescription. The method of projection is employed JHEP01(2022)189 to perform the calculation. From the results, one can extract the helicity amplitudes in 't Hooft-Veltman scheme. The ultraviolet renormalisation includes the strong coupling constant as well as operator renormalisations. For the pseudo-scalar, the latter involves operator mixing.
In any Greeen's function containing multiple local operators, one could expect the presence of an additional divergence that could arise from the short distance behaviour of the two operators. However, we find no additional divergence of this kind in both the amplitudes. Consequently, there is no contact renormalisation that is needed. This is in corroboration with the analysis through operator product expansion.
The fully ultraviolet renormalised amplitudes exhibit the infrared poles in dimensional regulator as dictated by the prediction in terms of Catani's universal subtraction operator. We check that the finite remainders exhibit expected kinematic symmetry. We also implement the finite remainders in a numerical code that can be used for further phenomenological studies. The finite remainders are provided as ancillary files. Given the availability of these results, we are now one step closer to make the prediction of an observable for the production of di-Higgs or di-pseudo-Higgs arising from a subdominant channel which is important at this time of very high precision physics.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.