Two loop QCD amplitudes for di-pseudo scalar production in gluon fusion

We compute the radiative corrections to the four-point amplitude $g+g \rightarrow A+A$ in massless Quantum Chromodynamics (QCD) up to order $a_s^4$ in perturbation theory. We used the effective field theory that describes the coupling of pseudo-scalars to gluons and quarks directly, in the large top quark mass limit. Due to the CP odd nature of the pseudo-scalar Higgs boson, the computation involves careful treatment of chiral quantities in dimensional regularisation. The ultraviolet finite results are shown to be consistent with the universal infrared structure of QCD amplitudes. The infrared finite part of these amplitudes constitutes the important component of any next to next to leading order corrections to observables involving pair of pseudo-scalars at the Large Hadron Collider.

loops compared to those in the full theory. Unlike the case of CP even Higgs boson, inclusive cross section for the production of pseudo scalar Higgs is known [15][16][17] only up to nextto-next-to leading order (NNLO) in pQCD. For N 3 LO predictions, one requires three loop virtual amplitudes and real emission contributions. The computation of virtual corrections is technically challenging [18] as pseudo scalar Higgs boson couples to SM fields through two composite operators that mix under renormalisation. In addition, these operators involve Levi-Civita tensor and γ 5 which are hard to define in dimensional regularisation. The three loop form factor thus obtained was later combined with appropriate soft distribution function [19][20][21] and mass factorisation kernels to obtain soft plus virtual contribution at N 3 LO in QCD [22]. Later, the process dependent resummation constants from the three loop form factors were used to perform threshold resummation in [23] and also make approximate prediction at N 3 LO level. This was possible due to the similarity of the interaction vertices of scalar and pseudo scalar Higgs bosons with the gluons.
Recently, there have been a surge of interest to study the production of pair of Higgs bosons to determine Higgs self coupling, whose strength is a prediction of the SM, if the mass of the Higgs boson is known. Measurement of this coupling will provide an independent test on nature of the Higgs boson. The gluon gluon fusion subprocess producing pair of Higgs bosons through a heavy quark loop [13,24] is the dominant one at the LHC, however the cross section is only few tens of fb, making it very difficult to observe. QCD corrections not only increase the cross section but also stabilise the predictions against renormalisation µ R , and factorisation µ F scales. NLO QCD corrections [14] and later on the top quark mass effects are systematically taken into account in [25][26][27][28][29][30]. Beyond NLO, an EFT where top quark degrees of freedom are integrated out is used. At present, production of pair of Higgs bosons in EFT is known to N 3 LO level [31], for NLO, NNLO, see [25][26][27][28][29][30][32][33][34]. All the two loop virtual amplitudes for g + g → hh that are required for the N 3 LO cross section for the di-Higgs production were obtained in [35]. The production of di-Higgs bosons through bottom quark annihilation was obtained up to NNLO level in [36]. In [37][38][39], the fully differential results at NNLO level are presented. While, there have been flurry of activities in the context of scalar Higgs boson, very little is known for the production of pair of pseudo scalar Higgs bosons at the LHC so far. In [14], LO contribution keeping finite top mass and NLO contributions using EFT framework where top quark degrees of freedom are integrated out have been obtained. Like the production of single pseudo scalar Higgs boson, pair production is also important to understand the nature of the extended Higgs sector. In order to reduce the theoretical uncertainties, it is important to have QCD radiative corrections under control. Due to EFT, it is now possible to go beyond NLO with available tools to make precise as well as stable predictions with respect to the unphysical scales. At NNLO level, we require two loop virtual, one loop single real emission and double real emission amplitudes. In this article, as a first step towards obtaining going beyond NLO QCD corrections, we compute all the one and two loop amplitudes that can contribute to the pure virtual part of the cross section in dimensional regularisation and perform ultraviolet (UV) renormalisation to obtain UV finite results.
The paper is organised as follows: in section-2, we describe how two loop virtual amplitudes are computed. In particular, we introduce the effective Lagrangian, the relevant kinematics, describe how projector method can be applied to obtain the scalar parts of the amplitudes, the subtleties involved in defining the Levi-Civita tensor and γ 5 in dimensional regularisation, ultraviolet renormalisation of strong coupling, over all renormalisations for the composite operators and finite renormalisation for the γ 5 . In section-3, computation of the amplitudes and their infrared (IR) structure are briefly discussed. In section-4, we summarise our results and conclude.
2 Theoretical framework

Effective Lagrangian
We work with the effective Lagrangian [40] that describes the interaction of the pseudoscalar field Φ A (x) with the gauge field G aµν and the fermion ψ: The pseudo-scalar gluonic (O G (x)) and the light quark (O J (x)) operators are defined as where f abc is the SU(3) structure constant and ǫ µνρσ is the Levi-Civita tensor. The pseudoscalar fermionic operator is the derivative of the flavour singlet axial vector current The effective Lagrangian is obtained after integrating out the top quark fields in the large top mass limit. Hence, the corresponding Wilson coefficients C G and C J depend on the mass of the top quark m t . As a result of the Adler-Bardeen theorem [41], there is no QCD correction to C G beyond one-loop level. On the other hand, C J begins only at second-order in the strong coupling constant a s ≡ g 2 s /16π 2 = α s /4π. The Wilsons coefficients are given by where G F is the Fermi constant, cot β -the ratio of the vacuum expectation values of the two Higgs doublets, in a model where the CP is not spontaneously broken. C F is the quadratic Casimir in the fundamental representation of QCD and µ R is the renormalisation scale at which a s is renormalised. We use the effective lagrangian 2.1 to obtain amplitudes for the production of pair of pseudo-scalar Higgs bosons A of mass m A up to two loop level in perturbative QCD. We restrict ourselves to the dominant gluon fusion subprocess: where p 1 and p 2 are the momenta of the incoming gluons, p 2 1,2 = 0 and p 3 and p 4 are the momenta of the outgoing pseudo-scalar Higgs bosons, p 2 3,4 = m 2 A . The Mandelstam variables for the above process are given by which satisfy s + t + u = 2m 2 A . It is convenient to express these amplitudes in terms of the dimensionless variables x, y and z as which lead to the constraint x −1 + x = y + z.
As in the case of di-Higgs production amplitude via gluon fusion [24], the di-pseudo scalar production amplitude, can also be decomposed in terms of two second rank Lorentz tensors T µν i (i = 1, 2), as follows: where ǫ µ (p i ) are the polarisation vectors of the initial state gluons. The Lorentz scalar functions M i , i = 1, 2 are independently gauge invariant. δ ab indicates that there is no colour flow from initial to final state. The second rank tensors are given by with p 2 T = (tu − m 4 A )/s is the transverse momentum square of the pseudo-scalar Higgs boson expressed in terms of the Mandelstam variables. The tensor T µν 1 depends only on the initial state momenta p 1,2 . Using momentum conservation, it can be seen that T µν 2 is symmetric under the interchange of the two pseudo-scalar Higgs momenta. The scalar functions M 1,2 can be obtained from M µν ab , by using appropriate d-dimensional projectors P µν i,ab with i = 1, 2, respectively and the projectors are given by: 12) where N corresponds to the SU (N ) colour group.
In the following, we briefly discuss on the type of Feynman diagrams that contribute up to order O(a 4 s ) in QCD. To evaluate the 4-point amplitude g +g → A+A to any order in a s , one needs to calculate the contributing diagrams to that particular order and evaluate the scalar functions M 1,2 , using the projectors P µν i,ab , i = 1, 2. Using the effective Lagrangian eq.(2.1), the higher order corrections to g + g → A + A amplitude are calculated in massless QCD. There are two types of diagrams that contribute to this process. We classify them as type-I and type-II. The form factor type diagrams where a pair of gluons annihilate to a single A, which branches into a pair of As belong to type-I and type-II contains t and u channel diagrams where each A is coupled to pair of gluons, or to quarks. In type-I, we have two classes of diagrams: type-Ia ( fig. 1 left panel) which contains only four point AAgg effective vertex and type-Ib ( fig. 1 right panel) containing both AAg and AAA vertices. These diagrams contribute at tree level (O(a s )) and we need to calculate them to O(a 4 s ) i.e., up to 3-loop order. Since these diagrams are related to form factors of O G between gluons states and O J between quark and gluon states, we can readily obtain them from [18,42,43].
The type-II diagrams consist of (a s ). Since, type-I diagrams are known to required order in a s , the results presented in this paper will mainly include the type-II amplitudes up to two loops in massless perturbative QCD i.e. order O(a 4 s ). We use dimensional regularisation (d = 4 + ǫ) to regularise both UV and IR singularities which appear as poles in ǫ in the UV, soft and collinear regions. Since we will have to deal with the Levi-Civita tensor in O G operator and γ 5 in O J operator, both of which are constructs inherently in 4-dimensions, a consistent method to deal with them in 4 + ǫ dimensions is essential. We discuss the details of a consistent and practical prescription to go over to 4 + ǫ and its implications in the next section. Hence, the scalar amplitudes M i can be written as a sum of amplitudes resulting from types-I and II diagrams as and in the following we concentrate only on M II i .

γ 5 within dimensional regularisation
Due to the axial anomaly, the pseudo-scalar gluonic operator O G = ǫ µνρσ G aµν G aρσ is related to the divergence of the axial vector current O J = ∂ µ (ψγ µ γ 5 ψ). Computation of higher order corrections with chiral quantities, involve inherently d = 4 dimensional objects like γ 5 and the Levi-Civita tensor ǫ µνρσ , and this warrants a prescription in going away from 4-dimension i.e. d = 4 + ǫ. There exist several prescriptions to deal with γ 5 in dimensional regularisation [44,45]. In multi-loop computations that use dimensional regularisation, we use the self-consistent prescription for γ 5 that was proposed by 't Hooft and Veltman [44]. In this prescription, one defines γ 5 as where Levi-Civita tensor is purely 4-dimensional, while the Lorentz indices on the γ µ i are in d = 4 + ǫ dimensions. To maintain the anti-commuting nature of γ 5 with d-dimensional γ µ i , the symmetrical form of the axial current has to be used  this is in concurrence with the above definition of γ 5 in eq. 2.14, and will lead to The O G and O J operators now take the form Contraction of two Levi-Civita tensors that result from either O G operator or the mixing of O G and O J operators is given by the Lorentz indices in this determinant, could now be considered as d-dimensional and the consequence would be, addition of only the inessential O(ǫ) terms to the renormalisated quantity [46]. This prescription though is not without consequence-a finite renormalisation of the axial vector current [47] is required in order to fulfill the chiral Ward identities and the Adler-Bardeen theorem. This will be discussed further in the next section.

UV renormalisation, operator renormalization and mixing
In dimensional regularisation with d = 4 + ǫ, the bare strong coupling constant denoted bŷ a s is related to its renormalized coupling by a ŝ  20) up to O(a 3 s ). β i are the coefficients of the QCD β function and are given by [48] where n f is the number of flavors and T F = 1/2. As we work in an effective theory obtained after integrating out the top quark fields in the large top quark mass limit, n f = 5. The Casimirs of SU(N) are given by C F and C A : For type-I diagrams which begin to contribute at LO, the Z as up to order O(a 3 s ) will be needed while for type-II diagrams, one order lower is sufficient.
Apart from the renormalisation of strong coupling in the massless QCD, the amplitudes require the renormalisation of vertices resulting from the composite operators O G and O J of the effective Lagrangian eq. (2.1). The renormalised operators are denoted by [ ] parenthesis, while the bare quantities without the parenthesis.
The renormalisation of O J is related to the renormalisation of the singlet axial vector current J µ 5 which needs the standard overall UV renormalisation constant Z s M S and a finite renormalisation constant Z s 5 . The later is necessary in dimensional regularisation in order to ensure the nature of operator relation resulting from axial anomaly [49] [ which is true in Pauli-Villars, a 4-dimensional regularisation. To preserve eq. 2.23 in 4 + ǫ dimensions, the multiplicative finite renormalisation constant Z s 5 is required. The bare operator O J is renormalised multiplicatively, exactly in the same way as the singlet axial vector current J µ 5 , through whereas the bare pseudo-scalar gluon operator O G mixes with fermionic operator O J under the renormalisation through with the corresponding renormalisation constants Z GG and Z GJ . Combining the above two equations in a matrix form, we have where Z JG = 0 to all orders in perturbation theory and Z JJ ≡ Z s 5 Z s M S . The renormalisation constants required for above equation are available up to O(a 3 s ) [46], [50] which was computed using OPE. For earlier works on this, see [51,52]. Using a completely different method the same quantities were calculated by some of us [18] and found to be in full agreement. The UV renormalisation constant of the singlet axial vector current J µ 5 in the M S scheme is 28) and the finite renormalisation constant Z s 5 is The renormalisation constants for O G and O J operators up to two loops are given by  GG,g can be related toM GG,g using eq. 2.20 and eq. 2.30. Expanding the renormalisation constants Z KL in eq. 2.30 as We find that the UV singularities that appear at one-loop and two-loop levels can be taken care of by the coupling constant renormalisation Z as and operator renormalisation Z ij . At this point we would like to stress that there could be additional contact terms required as a result of the behaviour of product of operators O G O G or O G O J at short distances. As shown in [50], we find that there are no contact terms as a result of these product of operators at short distances. For earlier works on this, see [51,52].

Calculation of the Amplitude
Our task of computing the amplitude g + g → A + A has reduced to the type-II diagrams up to O(a 4 s ). This involves diagrams with two Agg effective vertices, up to two-loop level in QCD (Type-IIa) and diagrams with one Agg effective vertex and one Aqq effective vertex which involves terms up to one loop in QCD (Type-IIb). Diagrams involving two Aqq effective vertex start at O(a 5 s ) and are not considered here. Applying the projectors P µν i,ab on the amplitudes, we extract the scalar coefficients M i with i = 1, 2 at every order in the perturbation.
All the tree level, one loop and two loop Feynman diagrams in massless QCD are generated using QGRAF [53] where additional vertices resulting from effective lagrangian eq. (2.1) are incorporated. There are two tree level diagrams, 35 one-loop diagrams and 789 two-loop diagrams of type-IIa. For type-IIb which involves effective quark and gluon couplings to pseudo-scalar Higgs, there are no tree level diagrams but 8 diagrams that contribute at one-loop which suffices to generate diagrams up to O(a 4 s ). The raw QGRAF output is converted with the help of in-house codes based on FORM [54] to include appropriate Feynman rules and to perform trace of Dirac matrices, contraction of Lorentz indices and colour indices. At this stage, we encounter huge number of one and scalar 2-loop Feynman integrals, which contain a set of propagator denominators and a combination of scalar products between loop momenta and independent external momenta. These Feynman integrals can be classified in terms of propagator denominators, that they contain. It is hence important to identify the momentum shifts that are required to express each of these diagrams in terms of a standard set of propagators called auxiliary topology. We use REDUZE2 package [55] to achieve this. The auxiliary topologies needed for the present case are same as those found in vector boson pair production [56,57] at two loops.
As expected these large number of scalar integrals are not all independent. To establish the relations, some properties of the Feynman integrals in dimensional regularisation are used. Exploiting the fact that, the total derivative with respect to any loop momenta of these integrals, evaluates to a surface term, which vanishes, leads to integration-by-parts (IBP) identities [58,59]. In addition, the fact that all integrals are Lorentz scalars, gives rise to Lorentz invariance (LI) identities [60]. As a result, these integrals can in turn be expressed in terms of a much smaller set of integrals which are irreducible and appropriately called master integrals (MI). Several automated computer algebra packages are available [55,[61][62][63][64] that use the Laporta algorithm [65] to reduce these Feynman integrals to the MIs. We have used the Mathematica based package LiteRed [64] to perform the reductions of all the integrals to MIs. At one-loop, there are 10 MIs, while at two-loop the number is 149. These two-loop MIs are the same as two-loop four-point functions with two equal mass external legs. The analytical result for the each MI in terms of Laurent series expansion in ǫ is given in [56,57].
At this stage, the renormalisation of the strong coupling constant and of the operators O G and O J , described in section 2.3, removes all the UV singularities. The singularities that still remain are purely of infrared origin and the next section is devoted to it.

Infrared factorization
The UV finite amplitudes that we have computed contain only divergences of infrared origin, which appear as poles in the dimensional regularization parameter ǫ. They are expected to cancel against real emission diagrams for the IR safe observables. While these singularities disappear in the physical observables, the amplitudes beyond leading order show a very rich universal structure in the IR. In [66], Catani predicted the IR poles of twoloop n-point UV finite amplitudes in terms of certain universal IR anomalous dimensions. Later, in [67], factorization and resummation properties of QCD amplitudes were used to understand the IR structure and subsequently the attempts were made to predict the structure of IR poles beyond two loops in [68,69]. Following [66], we obtain g (ǫ) are the IR singularity operators given by

4)
At one loop level, it is straight forward to show analytically that the IR poles are in agreement with the predictions. For the two loop case a fully analytical comparison was possible only for poles ǫ −i with i = 2 − 4. However, due to the large file size for the ǫ −1 pole term, we made a comparison only at the numerical level. We found full agreement with the predictions of Catani up to two loop level for all the IR poles. Having obtained the IR pole cancellations, the finite part defined in eq. 3.1 can be extracted by subtracting the IR poles. Expressions for the finite part are too large to be presented here, however, they are being provided as ancillary files. The finite part of the amplitude corresponding to the projectors 1 and 2 as a function of x for different values of cos θ is plotted normalising appropriately, in fig. 4.

Discussion and Conclusions
In this paper, we have presented the two loop virtual amplitudes that are relevant for studying production of pair of pseudo-scalar Higgs bosons in gluon fusion subprocess at the LHC. This is the dominant sub process that is sensitive to its self coupling. We have done this computation in the EFT where top quark degrees of freedom is integrated out.
In the EFT, the pseudo-scalar Higgs boson directly couples to gluons and light quarks through two local composite operators O G and O J respectively with the strengths proportional to Wilson coefficients that are calculable in perturbative QCD. We used dimensional regularisation to regulate both UV and IR divergences. The composite operators being CP odd, contain Levi-Civita tensor and γ 5 which are inherently four dimensional objects. Hence, a careful treatment was needed to deal with them in d-dimensions. We followed the prescription advocated by Larin. This requires additional renormalisation for the singlet axial vector current up to two loops. In addition, Larin's prescription requires finite renormalistion constant for singlet axial current and is also available. Note that the composite operators mix under UV rerenormalistion. The corresponding renormalisation constants are already known and we use them to obtain UV finite two loop amplitudes. Unlike the amplitudes involving pair of Higgs bosons, we do not need any UV contact counter terms. The UV finite amplitudes thus obtained contain IR divergences due to the presence of massless partons in QCD. We found that these IR poles are in agreement with the predictions by Catani and it provides a test on the correctness of the computation. Our results provide one of the important components relevant for studies related to production of pair of pseudo-scalar Higgs bosons at the LHC up to order O(a 4 s ).