Two loop QCD corrections for the process Pseudo-scalar Higgs $\rightarrow 3$ partons

We present virtual contributions up to two loop level in perturbative Quantum Chromodynamic (QCD) to the decay of pseudo-scalar Higgs boson ($A$) to three gluons ($g$) and also to quark ($q$), anti-quark ($\overline q$) and a gluon. With appropriate crossing, they are well suited for predicting the differential distribution of $A$ in association with a jet in hadron colliders up to next-to-next-to-leading order (NNLO) in strong coupling constant and also for the subsequent decay of $A$ to hadrons. We use effective field theory approach to integrate out the top quarks in the heavy top limit. The resulting theory involves two pseudo-scalar composite operators describing the interaction of $A$ with gluons as well as with quark and anti-quark. We perform our computation in dimensional regularisation and use minimal subtraction ($\overline{MS}$) scheme to renormalise strong coupling constant as well as the composite operators. The ultraviolet (UV) finite amplitudes contain infrared (IR) divergences that are found to be in agreement with the predictions by Catani. For both the amplitudes namely $A \rightarrow g g g$ and $A \rightarrow q \overline q g$, the leading transcendental terms at one and two loops are found to be identical to those in a three point form factor (FF) of the half-BPS operator in ${\cal N}=4$ Supersymmetric Yang Mills (SYM) theory when the QCD color factors are adjusted in a specific way. We present our results in terms of harmonic polylogs well suited for further numerical study.


Introduction
The discovery of a scalar particle in 2012 at the Large Hadron Collider (LHC) [1,2] is the milestone in high energy physics. The detailed study on the nature of this particle and its couplings to other Standard Model (SM) particles is underway. There are already several compelling evidences pointing to the fact that the discovered boson is none other than the Higgs boson of the SM. Although the SM is remarkably successful in explaining most of the observed phenomena at the subatomic level, it has many shortfalls. For example, it does not have dark matter candidate to explain relic abundance in the early universe, similarly the observed baryon asymmetry and the phenomena of neutrino oscillations. There exist several beyond the SM (BSM) scenarios that address these issues and provide plausible explanations. Often they have larger symmetry and contain more particles than the SM. Most of the BSMs contain the scalar sector with additional scalars providing the scope for rich phenomenology at the LHC. The precise measurements of the coupling of the Higgs boson with other SM particles as well as the dedicated direct searches of new scalars at the LHC can constrain the parameters of these scalar sectors in the BSM models. For example, in the Minimal Supersymmetric Standard Model (MSSM) which addresses some of the shortfalls of the SM contains in addition to other fields, two isospin doublets of Higgs field [3][4][5][6][7][8][9][10]. These two doublets preserve the analyticity of the scalar potential and also maintain the anomaly cancellation. They are responsible for masses of up and down type fermions. After spontaneous symmetry breaking, scalar sector contains three neutral (h, H, A) and two charged (H ± ) Higgs bosons, where h and H are CP-even scalars while A is CP-odd pseudo-scalar. The upper bound on the mass of the light scalar Higgs h can be predicted within the theory as the self couplings of the Higgs fields are fixed in terms of the gauge couplings.
Since the coupling of A with the fermions is appreciable in the small and moderate tan β, the ratio of vacuum expectation values v i , i = 1, 2 of the Higgs doublets, it can be searched at the LHC. Note that large gluon flux at the LHC can also boost the cross section. While, the searches for A have been underway at the LHC, predictions for observables involving A, in the theoretical side, are already available. Like the production of CP scalar Higgs boson, the leading order predictions for A severely suffer from theoretical uncertainties resulting from renormalisation (µ R ) and factorisation (µ F ) scales as well as large higher order radiative corrections. These scales enter at leading order through the renormalised strong coupling constant g s (µ R ) and the parton distribution functions defined at the factorisation scale µ F . This necessitates to go beyond leading order in perturbation theory. In [11][12][13][14] the next-to-leading order (NLO) perturbative QCD corrections to inclusive production of A were computed and it was found that the scale dependence reduced from 48% to 35% and the corrections were found to be as large as 67%. The effective field theory approach where the top quark degrees of freedom are integrated out has provided opportunity to go beyond NLO level to further improve the predictions. The results on NNLO production cross section for A can be found in [15][16][17] which further improves the reliability of the predictions.
There have been continued efforts to go beyond NNLO as the results on the form factors (FFs) and soft gluon contributions at third order level have become available. For this purpose, in [18], we computed FFs of the effective composite operators between quark and gluon states at three loop level in QCD along with the lower order. Using the formalism developed in [19,20] and the third order results on the FFs of the pseudo-scalar Higgs boson [18], the universal soft-collinear distribution [21], operator renormalisation constant [18,22,23] and the mass factorisation kernels [24,25], the first results on the threshold N 3 LO corrections for the inclusive production of A were reported in [26]. In addition, the third order corrections to both the N -independent part of the resummed cross section [27,28] and the matching coefficients in soft-collinear effective theory (SCET) were also computed in [18]. In [29], we obtained the approximate N 3 LO contribution using the full N 3 LO results available for the scalar Higgs boson. This along with the threshold effects at next-to-next-to-next-to-leading logarithm (N 3 LL) accuracy using both conventional [27,28] and SCET [30][31][32][33][34][35][36] setups provide the accurate prediction [29] for the inclusive production of A at the LHC.
In order to probe the nature of A and its coupling to other SM particles, one also needs to study exclusive observables namely the differential cross sections. The distributions of transverse momentum and rapidity of the produced A are often very useful for characterising its properties. In addition, observables involving A recoiled against one or two hard jets can help to efficiently reject the background from other sources. There are already several studies in the case of scalar Higgs boson in the literature. For example, the differential cross sections for scalar Higgs boson via gluon fusion are available including its decay to two photons or two weak gauge bosons, see [37][38][39]. It was noted in [40,41] that observables with jet vetos enhance the significance of the signal considerably allowing one to study the properties of the Higgs boson and its coupling to other SM particles. Recently, in [42], a complete NNLO prediction for Higgs-plus-jet final states and for the transverse momentum distribution of the Higgs boson taking into account the experimental definition of the fiducial cross sections has been achieved. For Higgs-plus-jet production at NNLO, see [43][44][45][46][47]. The relevant two loop amplitudes can be found in [48][49][50][51][52][53][54]. Pseudo-scalar Higgs boson being produced dominantly in gluon initiated processes shares a lot with its counter part, namely the scalar Higgs boson. In particular one finds that the theoretical uncertainties and the perturbative corrections are significantly large. For the pseudo-scalar Higgs boson, the results for differential distributions are available only up to NLO level, for example, see [55,56]. Hence it is desirable to have predictions upto NNLO level for observables like A + 1jet, transverse momentum and rapidity distributions of A, taking into account various experimental cuts. At NNLO, one encounters three different production channels that contribute to these observables: pure virtual, virtual-real emissions, pure real emissions. In this paper, we will present all the relevant contributions resulting from pure virtual amplitudes. In particular, we will consider the processes A → ggg and A → qqg in effective theory where top quark has been integrated out up to two loop level. These amplitudes can be used to predict the production of A in hadron collider up to two loop level in QCD after performing straightforward crossing of kinematic variables. In the effective theory for pseudo-scalar, one is left with two effective operators with same quantum numbers and mass dimensions, hence they mix under renormalisation. We will present the results in the M S scheme.
Our paper is organised as follows. In Section 2.1, the effective Lagrangian for the coupling of A with partons along with the relevant Wilson coefficients are presented. After introducing the notation in Section 2.2, we elaborate on how to deal with γ 5 and Levi-Civita tensor in dimensional regularisation and subtleties involved in operator mixing and finite renormalisation in Section 2.3. In Section 2.4, we describe the method of our computation of relevant matrix elements up to two loop level in QCD and Section 2.5 contains the study of IR structure of the results obtained in the context of Catani's predictions. In Section 3, we discuss our results on one and two loops QCD amplitudes in the light of maximum transcendentality principle in N = 4 SYM. Finally we conclude in Section 4.

Effective Lagrangian
In the generic BSM scenarios, the pseudo-scalar Higgs boson couples to heavy quarks through Yukawa interaction. If we restrict ourselves to the dominant contribution coming from top Yukawa coupling, then the interaction Lagrangian is given by where the model dependent constant g t in MSSM is cos β where tan β is ratio of vacuum expectation values of the two Higgs fields, m t is the top quark mass, v the SM vacuum expectation value.
In the limit of infinite top quark mass, the interaction between a pseudo-scalar A and the fields of the remaining SM is encapsulated by an effective Lagrangian [57] which reads as where the operators are defined as The Wilson coefficients C G and C J are obtained by integrating out the loops resulting from top quark. The coefficient C G does not receive any QCD corrections beyond one loop in the perturbation series expanded in strong coupling constant a s ≡ g 2 s /(16π 2 ) = α s /(4π) due to the Adler-Bardeen theorem [58] and C J starts at second order. Their expressions are given below In the above expressions, G µν a and ψ represent gluonic field strength tensor and light quark fields, respectively and G F is the Fermi constant. Here, a s ≡ a s µ 2 R is the renormalised strong coupling constant at the scale µ R which is related to the unrenormalised one,â s ≡ g 2 s /(16π 2 ) throughâ β i are the coefficients of the QCD β-functions [59], given by with the SU(N) QCD color factors n f is the number of active light quark flavours.

Notation
We consider the decay of pseudo-scalar Higgs boson to both ggg and qqg which can be expressed as The associated Mandelstam variables are defined as which satisfy the relation where M A is the mass of the pseudo-scalar Higgs boson. We also define following dimensionless invariants, which we use to describe our results in terms of HPLs [60] and 2d-HPLs [61,62] as where (x, y, z) lie between 0 and 1 and satisfy the condition

Operator mixing and UV renormalisation
In the following, we will describe the computation of one and two loop matrix elements that are needed to perform NNLO QCD corrections to production of a pseudo-scalar with a jet in hadron colliders or its decay to three jets. Since the amplitudes for the production of pseudo-scalar in association with a jet from gluon gluon or quark anti-quark initiated channels are related to the decay of the same to three gluons or quark anti-quark gluon by crossing symmetry, in the rest of the paper we will only describe the latter in detail. The composite operators present in the effective Lagrangian Eq. (2.2), develop UV divergences that require additional renormalisation. Also these operators mix under renormalisation due to same quantum numbers. In higher order computations involving chiral quantities, the inherently four dimensional objects like γ 5 and ε µνρσ , the Levi-Civita tensor, in d = 4 dimensions pose problems. In this article, we have followed the most practical and self-consistent definition of γ 5 introduced for multi loop calculations in dimensional regularization by 't Hooft and Veltman through [63] (2.14) Here, ε µνρσ is contracted according to the rule where all the Lorentz indices are considered to be d-dimensional [22]. The prescription used here fails to preserve the anti-commutativity of γ 5 with γ µ in arbitrary d-dimensions.
In addition, the Ward identities, which are valid in a 4-dimensional regularization scheme like the one of Pauli-Villars where γ 5 does not pose any problem, are violated as well. Due to this, it is not possible to restore the correct renormalisation of axial current, which is defined as [22,64]. A finite renormalisation of the axial vector current is required to preserve chiral Ward identities and the Adler-Bardeen theorem. The axial current is defined as in dimensional regularization. The chiral Ward identities can be fixed by introducing a finite renormalisation constant Z s 5 [58,65] in addition to the standard overall UV renormalisation constant Z s M S within the M S-scheme: While the evaluation of the appropriate Feynman diagrams explicitly fixes the constants Z s M S , the finite renormalisation constant can not be fixed without using the anomaly equation. In other words, Z s 5 can be determined by demanding the conservation of the one loop character [66] of the operator relation of the axial anomaly in dimensional regularization: The bare operator [O J ] B is renormalised multiplicatively as the axial current J µ 5 through through renormalisation constants Z GG and Z GJ . The above two equations can be written as where Z JG = 0 to all orders in perturbation theory , The expressions for the above mentioned renormalisation constants Z s M S , Z GG , Z GJ up to O a 3 s are given in [22,23] using operator product expansion. In [18] one of us has calculated the same quantities in a completely different way and found exact agreement with the original calculation. In the latter article, authors have used universality of the IR poles of the FF to determine the UV renormalisation constants and also computed Z s 5 up to O(a 2 s ) by demanding the operator relation of the axial anomaly Eq. (2.18). These renormalisation constants up to sufficient order in the perturbation theory appropriate for our calculation are given below: Using these operator renormalisation constants along with strong coupling constant renormalisation through Z as , we obtain UV finite amplitudes.

Matrix elements
Our next step is to compute all the relevant matrix elements resulting from virtual amplitudes for the decay of pseudo-scalar to three gluons and also to quark antiquark gluon. They are obtained from the amplitudes |A f , where f = ggg, qqg up to two loop level. They contain two sub-amplitudes namely |M Λ f computed using O G (Λ = G) and O J (Λ = J) operators multiplied by the appropriate Wilson coefficients C Λ : The above UV finite sub-amplitudes |M Λ f can be expressed in terms unrenormalised subamplitudes |M , (2.32) and C (1) (2.33) The computation of S Λ,(i) a for a = g, q; Λ = G, J, GJ, JG, i = 0, 1, 2 closely follows the steps used in [67][68][69][70][71][72][73][74]. We briefly describe the main steps involved. The relevant tree, one and two loop Feynman diagrams are generated using the package QGRAF [75]. While there are only few diagrams at tree and one loop level, at two loops we encounter large number of diagrams. We find that the number of two loop diagrams is 1306 for |M qqg . We set all the external particles on-shell, in other words, the quarks, anti-quarks, gluons are kept massless and the pseudo-scalar has non-zero mass M A . The raw QGRAF output is converted to a format that can be further used to perform the substitution of Feynman rules, contraction of Lorentz and color indices and simplification of Dirac and Gell-Mann matrices using a set of in-house routines written in the symbolic manipulating program FORM [76]. We have included ghost loops in the Feynman gauge. For the external on-shell gluons, we have kept only transversely polarization states in d-dimensions. The next step is to organise all the ddimensional one and two loop integrals such that they can be reduced to a minimum set of scalar integrals, called master integrals (MIs). To do this, we first use Reduze2 [77,78] to shift loop momenta to get suitable integral classes. Each integral belonging to specific class can be expressed in terms of MIs using a set of integration-by-parts (IBP) [79,80] and Lorentz invariance (LI) [81] identities. The LI identities are not linearly independent from the IBP identities [82], however they help to increase the speed to solve the large number of linear equations resulting from IBP. The method of IBP to reduce certain class of integrals to a set of MIs is achieved using Laporta algorithm, [83]. This has been implemented in various symbolic manipulation packages such as AIR [84], FIRE [85], Reduze2 [77,78] and LiteRed [86,87]. We have used LiteRed [86,87] to perform the reductions of all the integrals to MIs. For one and two loop, we find the number of MIs are 7 and 89 respectively. The next task is to express all the MIs that we obtain after using LiteRed to those computed analytically in [61,62]. Using these MIs in the appropriate kinematic regions, we obtain analytical results for S f for f = ggg, qqg. While these amplitude squares are UV finite, they are sensitive to IR divergences. Thanks to universality of these divergences, our results can be verified up to finite terms.

Universal Structure of IR divergences
Beyond leading order in perturbation theory, the UV renormalised amplitudes in gauge theories suffer from IR divergences due to the presence of massless particles. These divergences are classified into two categories, namely soft and collinear. Soft divergences arise when momentum of a massless particle in the loop goes to zero and the collinear ones result when a massless particle in the loop becomes parallel to one of the external particles. Note that the amplitudes are not observables, instead the cross sections or decay rates made out of these amplitudes are the observables. Thanks to Kinoshita-Lee-Nauenberg (KLN) theorem [88,89], the finiteness of the observables can be guaranteed by summing over the degenerate final states and performing mass factorisation for the initial states. While these divergences go away in the physical observables, the amplitudes demonstrate very rich universal structure in the IR region at every order in perturbation theory. The IR structure of amplitudes in QCD is well understood and in particular, the universal structure of IR poles was predicted in a seminal paper [90] by Catani for n-point two loop amplitudes. The iterative structure of singular part of the UV renormalised amplitudes in QCD was exploited in [90] to predict the subtraction operators that capture the universal IR divergences. Later on Sterman and Tejeda-Yeomans related the predictions by Catani to the factorisation and resummation properties of QCD amplitudes [91]. It is easy to convince oneself that iterative structure is the result of factorisation. The generalization of Catani's proposal for arbitrary number of loops and legs for SU (N ) gauge theory using soft-collinear effective theory was given by Becher and Neubert [92]. A similar result was also presented by Gardi and Magnea [93] by analyzing Wilson lines for hard partons.
We follow here Catani's proposal and express the UV renormalised amplitudes |M (1) The IR singular part of S f ( A f |A f ) computed up to a 5 s agrees perfectly with the predictions based on Catani's proposal. Hence, we present only the finite part of S f for f = ggg, qqg, namely S ggg,f in and S qqg,f in in the appendix A.

Universality of leading transcendental terms
Kotikov and Lipatov [94][95][96] (see also [97,98]) conjectured maximum transcendentality principle which relates transcendental terms in the anomalous dimensions of leading twist two operators in N = 4 SYM with those in the corresponding QCD results [24,25]. In the perturbative calculations, one associates transcendentality weight n to the constants ζ(n), −n and the functions namely HPLs of weight n. In N = 4 SYM, the FFs of the half-BPS type operators show uniform transcendentality. On the other hand QCD results contain terms of all transcendental weights in addition to rational terms (zero transcendentality). From various higher order results that are available in QCD and N = 4 SYM theories, one finds similar connection. In particular, the leading transcendental terms of anomalous dimensions of twist two operators, two and three point FFs in QCD when C A = C F = N and T f n f = N/2 coincide with the corresponding ones of the half-BPS operators in N = 4 SYM. The leading transcendental terms of the QCD amplitude for the Higgs boson decaying to three on-shell gluons [54,99] are also related to the two loop three point MHV FFs of the half-BPS operators [100] in N = 4 SYM. Similar relations are found for two point FFs of quark current, scalar and pseudo-scalar operators constructed out of gluon field strengths, energy momentum tensor of the QCD up to three loops see [18,[101][102][103]. Recently, in [104], we studied several three point FFs of the half-BPS as well as the Konishi operators in N = 4 SYM at two loop level in 't Hooft coupling. We evaluated them between onshell final states gφφ and φλλ where φ, λ, g are scalar, Majorana and gluon states. Unlike the half-BPS, the Konishi operator is not protected by supersymmetry and hence we find non-trivial structure at higher orders. The explicit computation shows that the FF of the Konishi operator does not show uniform transcendentality [101,104,105]. Interestingly, the leading transcendental terms of the Konishi for the states gφφ agree with those of the half-BPS operator [104]. On the other hand, for φλλ states, we did not find any such relation to the corresponding ones for the half-BPS. In the present context, it is tempting to relate our results for FFs to the Konishi operator, rather than the half-BPS owing to the fact that both are not protected. We find that the amplitude A → ggg agrees with the FF of the Konishi operator between gφφ states at the leading transcendental level while A → qqg for O G does not have any relation with that of the Konishi operator computed between φλλ states. Surprisingly, we find that the leading transcendental terms of the amplitude A → qqg multiplied by its born and normalised by its born square contribution agrees with the half-BPS. In summary both the amplitudes namely A → ggg and A → qqg computed using O G operator multiplied by respective born amplitudes and normalised by their born squares have leading transcendental terms identical to those of the half-BPS. Our results with pseudo-scalar operator computed between different on shell states shed more light on the universal structure of higher order terms in gauge theories.

Conclusion
The motivation to compute two loop virtual corrections to decay of a pseudo-scalar to three gluons or quark, anti-quark, gluon is two fold. Firstly, with larger data set that are available at the LHC, constraining BSM scalar sector has become feasible and hence precise theoretical predictions for production of pseudo-scalar and its subsequent decays to jets need to be achieved. Secondly, such a computation of NNLO level for the differential rates is technically challenging due to non-trivial tensorial interaction encountered here. We have successfully computed UV renormalised virtual contributions to A → ggg and A → qqg upto two level in QCD perturbation series using dimensional regularisation. The UV renormalisation is performed in M S scheme. We used an effective field theory approach which in our case involves integrating out top quark degrees of freedom. The resulting effective operators that contribute to the production of pseudo-scalar contain a Levi-Civita tensor and γ 5 matrix. These objects require special treatment in dimensional regularisation as they can not be defined in arbitrary dimensions other than four. Most of prescriptions define them in dimensional regularisation, however all of them fail to preserve the symmetries or Ward identities of the theory. We have used one of them namely the 't Hooft-Veltman prescription and in order to restore the symmetry we performed a nontrivial finite renormalisation on top of other standard UV renormalisations. This way, we could successfully compute UV finite two loop contributions. The resulting expressions do contain soft and collinear divergences due to the presence of massless gluons and quarks. Since these divergences are universal due to factorisation properties of on-shell QCD amplitudes order by order in perturbation theory, we confirm the correctness of our results up to finite terms using Catani's infrared factorisation formula which predict all the poles resulting from soft and collinear configurations up to two loops. Expressing our results using Catani's subtraction operators, the remaining infrared finite terms of the amplitudes are presented in the appendix A. It is well known in the context of scalar Higgs boson that the leading transcendental terms of finite part of H → ggg amplitude up to two loop level are related to the corresponding ones for three particle amplitude with an insertion of the half-BPS operator provided the color factors in QCD are adjusted in a specific way. We find that this is indeed the case for A → ggg, qqg amplitudes presented in this paper. In summary, the results for A → ggg and A → qqg up to two loop level in QCD are presented in this paper. With appropriate analytical continuation, our results can be readily used for the study of production of a pseudo-scalar in association with a jet at two loop level at the hadron colliders.

Acknowledgments
We would like to thank T. Ahmed, G. Das, T. Gehrmann, R. Lee and N. Rana for useful discussions and timely help.