Pseudo-scalar Higgs Boson Production at Threshold N$^3$LO and N$^3$LL QCD

We present the first results on the production of pseudo-scalar through gluon fusion at the LHC to N$^3$LO in QCD taking into account only soft gluon effects. We have used the effective Lagrangian that describes the coupling of pseudo-scalar with the gluons in the large top quark mass limit. We have used recently available quantities namely the three loop pseudo-scalar form factor and the third order universal soft function in QCD to achieve this. Along with the fixed order results, we also present the process dependent resummation coefficient for threshold resummation to N$^3$LL in QCD. Phenomenological impact of these threshold N$^3$LO corrections to pseudo-scalar production at the LHC is presented and their role to reduce the renormalisation scale dependence is demonstrated.


Introduction
The spectacular discovery of the Higgs boson [1,2] at the Large Hadron Collider (LHC) has put the Standard Model (SM) of elementary particles in the firm footing. Most importantly, the mystery of the electroweak symmetry breaking [3][4][5][6][7] mechanism can now be solved. The consistency of the measured decay rates of the Higgs boson to a pair of vector bosons namely W + W − , ZZ and fermions bb, τ τ with the precise predictions of the SM for the measured Higgs boson mass of 125 GeV within the experimental uncertainty [8,9] makes this discovery very robust. In addition, there is a strong evidence that the discovered Higgs boson has spin zero and even parity [10,11]. The ongoing 13 TeV run at LHC will indeed provide further scope to study the properties of the Higgs boson in great detail.
While the SM is complete in the sense that all of its predictions have been tested experimentally, the model suffers from various deficiencies as it can not explain baryon asymmetry in the Universe, dark matter, neutrino mass etc. There are several extensions of the SM, motivated to address these issues. The minimal version of Supersymmetric Standard Model (MSSM) [12] is one of the most elegant extensions of the SM and it addresses the above mentioned issues. The Higgs sector of it contains a pair of Higgs doublets which after symmetry breaking gives two CP even Higgs bosons h, H and one CP odd (pseudo-scalar) Higgs boson (A) and two charged Higgs bosons H ± [13][14][15][16][17][18][19][20] . The predicted upper bound on the mass of the lightest Higgs boson (h) up to three loop level is consistent [21][22][23] with the recently observed Higgs boson at the LHC. The efforts to test the predictions of MSSM or its variants have already been underway and the current run at the LHC will shed more light on them. One of them could be to look for CP odd Higgs boson in the gluon fusion through heavy fermions as its coupling is appreciable in the small and moderate tan β, the ratio of vacuum expectation values v i , i = 1, 2. In addition, large gluon flux can boost the cross section.
Since, the leading order production mechanism of the pseudo-scalar of mass m A is through heavy quarks, the cross section is not only proportional to tan β but also square of the strong coupling constant. Like the scalar Higgs boson in SM, the leading order prediction of the pseudo scalar production at the hadron colliders suffers from large theoretical uncertainties due to renormalisation scale µ R that enters in the strong coupling constant and the factorisation scale µ F in the gluon distribution functions of the protons. Predictions based on one loop perturbative Quantum Chromodynamics (pQCD) corrections [24][25][26][27] reduce these uncertainties (in the conventional range with the central scale µ = m A /2 and m A = 200 GeV) from about 48% to 35% while increasing the LO cross section substantially, by as large as 67%. Effective theory approach in the large top quark mass limit provides an opportunity to go beyond NLO. Such an approach [27,28]in the case of scalar Higgs boson production [29][30][31] turned out to be the most successful one as the finite mass effects at NNLO level were found to be within 1% [32][33][34]. NNLO predictions for the production of pseudo-scalar at the hadron colliders are already available [31,35,36]. The NNLO correction increases the NLO cross section by about 15% and reduces the scale uncertainties to about 15%. Due to large gluon flux at the threshold, namely when the mass of A approaches to the partonic centre of mass energy, the cross section is dominated by the presence of soft gluons. These contributions often can spoil the reliability of the predictions based on fixed order perturbative computations. Resummation of large logarithms resulting from soft gluons to all orders in the perturbation theory provides the solution to this problem. The systematic predictions based on the next-to-next-to-leading log (NNLL) resummed result [37][38][39][40][41][42][43][44][45] demonstrate the reliability of the approach and also reduce the scale uncertainties.
A complete calculation at NNLO [29][30][31] and leading logarithms at N 3 LO in the threshold limit [38][39][40][41][42] and NNLL soft gluon resummation [37] for the scalar Higgs boson production are known for more than a decade. Recently there have been series of works on predicting inclusive scalar Higgs boson production beyond this level. The computation of δ(1 − z) contribution at N 3 LO in the threshold limit [46] was the first among them. This was confirmed independently in [47]. Later on the sub-leading collinear logarithms were computed in [48,49]. Spin off of the result presented in [46] is the computation of N 3 LO prediction for the Drell-Yan production [47,50,51] at the hadron colliders in the threshold limit. In addition, one can obtain N 3 LO threshold corrections to the Higgs boson production through bottom quark annihilation [52] and also in association with vector boson [53] at the hadron colliders. Later, along the same direction, rapidity distribution of the Higgs boson in gluon fusion [54], DY [54] and Higgs boson in bottom quark annihilation [55] were obtained at threshold N 3 LO QCD.
A milestone in this direction was achieved by Anastasiou et. al. who have now accom-plished the complete N 3 LO prediction [56] of the scalar Higgs boson production through gluon fusion at the hadron colliders in the effective theory. These third order corrections increase the cross section by a few percent, about 2% and reduce the scale uncertainty by about 2%. Using these predictions, it is now possible to obtain the soft gluon resummation at N 3 LL, see [51,57].
While the next step in the wish list is to obtain the N 3 LO predictions for the pseudoscalar production through gluon fusion, the first task in this direction is to obtain the threshold enhanced cross section at N 3 LO level. One of the crucial ingredients is the form factor of the effective composite operators that couple to pseudo-scalar, computed between partonic states. One and two loop results for them between gluon states were computed for NNLO production cross section [35,36,58], the analytical results up to two loop level can be found in [58]. These were computed in dimensional regularisation where the space time dimension is d = 4 + ǫ. Threshold corrections to pseudo-scalar production at N 3 LO level requires the knowledge of the form factors up to three loop level. We also need to know one and two loop corrections computed to desired accuracy in ǫ, namely up to ǫ 2 for one loop and up to ǫ at two loops. In [59], we obtained the three loop form factors of the effective composite operators between quark and gluon states at three loop level along with the lower order ones to desired accuracies in ǫ. In the present article we will describe how threshold corrections at N 3 LO level can be obtained from the formalism developed in [40,41] using the available information on recently computed three loop form factor of the pseudo scalar Higgs boson [59], the universal soft-collinear distribution [50] and operator renormalisation constant [59][60][61] and the mass factorisation kernels [62,63] known to three loop level. In addition, we compute third order correction to the N -independent part of the resummed cross section [64,65] using our formalism [40,41]. We also present the numerical impact of our findings with a brief conclusion.
The underlying effective theory is discussed in the Sec. 2. This is followed by a short description of the formalism which has been employed to compute the soft-plus-virtual cross section in Sec. 3. We present the analytical results of these findings in the Sec. 4 up to N 3 LO in QCD. In Sec. 5, the N-independent parts of the threshold resummed cross section in Mellin space have been presented up to third order in QCD. Before making concluding remarks, in Sec. 6 we demonstrate the numerical implications of the fixed order soft-plus-virtual cross sections to N 3 LO at LHC.

The Effective Lagrangian
A pseudo-scalar couples to gluons only indirectly through a virtual heavy quark loop which can be integrated out in the infinite quark mass limit. The effective Lagrangian [66] describing the interaction between pseudo-scalar χ A and the QCD particles in the infinitely large top quark mass limit is given by where the two operators are defined as The Wilson coefficients C G and C J of the two operators are the consequences of integrating out the heavy quark loop in effective theory. C G does not receive any QCD corrections beyond one loop because of Adler-Bardeen theorem, whereas C J starts only at second order in the strong coupling constant. These Wilson coefficients are given by [66] C G = −a s 2 The symbols G µν a and ψ represent gluonic field strength tensor and quark field, respectively. G F stands for the Fermi constant and cotβ is the mixing angle in the Two-Higgs-Doublet model. m A and m t symbolise the masses of the pseudo-scalar and top quark (heavy quark), respectively. The strong coupling constant a s ≡ a s µ 2 R is renormalised at the mass scale µ R and is related to the unrenormalised one,â s ≡ĝ 2 s /16π 2 , througĥ The coefficient of the QCD β function β i are given by [67] with the SU(N) QCD color factors n f is the number of active light quark flavors.

Threshold Corrections
The inclusive cross-section for the production of a colorless pseudo scalar at the hadron colliders can be computed using where, the born cross section at the parton level including the finite top mass dependence is given by Here τ A = 4m 2 t /m 2 A and the function f (τ A ) is given by while the parton flux is given by where, f a and f b are the parton distribution functions (PDFs) of the initial state partons a and b, renormalised at the factorisation scale µ F . Here, ∆ A ab τ y , m 2 A , µ 2 R , µ 2 F are the partonic level cross sections, for the subprocess initiated by the partons a and b, computed after performing the overall operator UV renormalisation at scale µ R and mass factorisation at a scale µ F . The variable τ is defined as q 2 /s with q 2 = m 2 A . The goal of this article is to study the impact of the soft gluon contributions to the pseudo-scalar production cross section at hadron colliders. The infrared safe contribution is obtained by adding the soft part of the cross section to the ultraviolet (UV) renormalised virtual part and performing the mass factorisation using appropriate counter terms. This combination is often called the soft-plus-virtual (SV) cross section whereas the remaining portion is known as hard part. Thus, we write the partonic cross section as with z ≡ q 2 /ŝ = τ /(x 1 x 2 ) . The threshold contributions ∆ A,SV ab (z, q 2 , µ 2 R , µ 2 F ) contains only the distributions of kind δ(1 − z) and D i , where the latter one is defined through On the other hand, the hard part ∆ A,hard ab contains all the terms regular in z. The SV cross-section in z-space is computed in d = 4 + ǫ dimensions, as formulated for the first time in [40,41], using where, Ψ A g z, q 2 , µ 2 R , µ 2 F , ǫ is a finite distribution and C is the Mellin convolution defined as (3.8) Here ⊗ represents Mellin convolution and f (z) is a distribution of the kind δ(1 − z) and D i . The subscript g signifies the gluon initiated production of the pseudo-scalar. The equivalent formalism of the SV approximation is in the Mellin (or N -moment) space, where instead of distributions in z the dominant contributions come from the continuous functions of the variable N (see [64,65]) and the threshold limit of z → 1 is translated to N → ∞.
The Ψ A g z, q 2 , µ 2 R , µ 2 F , ǫ is constructed from the form factors F A g (â s , Q 2 , µ 2 , ǫ) with Q 2 = −q 2 , the overall operator UV renormalisation constant Z A g (â s , µ 2 R , µ 2 , ǫ), the soft-collinear distribution Φ A g (â s , q 2 , µ 2 , z, ǫ) arising from the real radiations in the partonic subprocesses and the mass factorisation kernels Γ gg (â s , µ 2 F , µ 2 , z, ǫ). In terms of the above-mentioned quantities it takes the following form, as presented in [41,50,52] In the subsequent sections, we will demonstrate the methodology to get these ingredients to compute the SV cross section of pseudo-scalar production at N 3 LO.

The Form Factor
The quark and gluon form factors represent the QCD loop corrections to the transition matrix element from an on-shell quark-antiquark pair or two gluons to a color-neutral operator O. For the pseudo-scalar production through gluon fusion, we need to consider two operators O G and O J , defined in Eq. In terms of these quantities, the full matrix element and the full form factors can be written as a series expansion inâ s as where Q 2 = −2 p 1 .p 2 = −q 2 and p i (p 2 i = 0) are the momenta of the external on-shell gluons. Note that |M J,(n) g starts at n = 1 i.e. from one loop level. The form factor for the production of a pseudo-scalar through gluon fusion,F A,(n) g , can be written in terms of the two individual form factors, Eq. (3.11), as follows: In the above expression, the quantities Z ij (i, j = G, J) are the overall operator renormalisation constants which are required to introduce in the context of UV renormalisation. These are are discussed in our recent article [59] in great detail. The ingredients of the form factor F A g , namely, F G g and F J g have been calculated up to three loop level by some of us and presented in the same article [59]. Using those results we obtain the three loop form factor for the pseudo-scalar production through gluon fusion. In this section, we present the unrenormalized form factorsF A,(n) g up to three loop where the components are defined through the expansion We present the unrenormalized results for the choice of the scale µ 2 R = µ 2 F = q 2 as follows: The results up to two loop level is consistent with the existing ones [58] and the three loop result is the new one. These are required in the context of computing SV cross-section which is discussed below. The form factor F A g (â s , Q 2 , µ 2 , ǫ) satisfies the KG-differential equation [68][69][70][71][72] which is a direct consequence of the facts that QCD amplitudes exhibit factorisation property, gauge and renormalisation group (RG) invariances: In the above expression, all the poles in dimensional regularisation parameter ǫ are captured in the Q 2 independent function K A g and the quantities which are finite as ǫ → 0 are encapsulated in G A g . The solutions of the KG equation in the desired form is given in [40] as (see also [50,52]) A A g 's are called the cusp anomalous dimensions. The constants G A g,i 's are the coefficients of a i s in the following expansions: However, the solutions of the logarithm of the form factor involves the unknown functions G A g,i which are observed to fulfil [58,73] the following decomposition in terms of collinear (B A g ), soft (f A g ) and UV (γ A g ) anomalous dimensions: where, the constants C A g,i are given by [41] C (3.20) In the above expressions, X A g,i with X = A, B, f and γ A g,i are defined through the series expansion in powers of a s : f A g 's are introduced for the first time in the article [58] where it is shown to fulfil the maximally non-Abelian property up to two loop level whose validity is reconfirmed in [73] at three loop level. Moreover, due to universality of the quantities denoted by X, these are independent of the operator insertion. These are only dependent on the initial state partons of any process. Hence, being a process of gluon fusion, we can make use of the existing results up to three loop: f g can be found in [58,73], A g,i in [62,63,73,74] and B g,i in [62,73] up to three loop level. Utilising the results of these known quantities and comparing the above expansion of G A g,i (ǫ), Eq. (3.19), with the results of the logarithm of the form factors, we extract the relevant g A,k g,i and γ A g,i 's up to three loop. For soft-virtual cross section at N 3 LO we need g A,1 g.3 in addition to the quantities arising from one and two loop. The form factors for the pseudo-scalar production up to two loop can be found in [58] and the three loop one is calculated very recently in the article [59] by some of us. However, in this computation of SV cross section at N 3 LO, we need the form factor in a particular form which is little bit different than the ones presented in our recent article [59], though the required one can be extracted from the results provided there. For readers' convenience, we have presented the form factors F A g up to three loop in the beginning of this section which have been employed to extract the required g A,k g,i 's using Eq. (3.16), (3.17) and (3.19). Below we present our finding of the relevant g A,k g,i 's up to three loop level: As a matter of emphasising the fact, note that the γ A g,i 's are found to satisfy is the usual QCD β-function. For more elaborate discussion on this, see recent article [59] (also see [60,61]).

Operator Renormalisation Constant
The strong coupling constant renormalisation through Z as is not sufficient to make the form factor F A g completely UV finite, one needs to perform additional renormalisation to remove the residual UV divergences which is reflected through the presence of non-zero γ A g in Eq. (3.19). This additional renormalisation is called the overall operator renormalisation which is performed through the constant Z A g . This is determined by solving the underlying RG equation: Using the results of γ A g,i from Eq. (3.24) and solving the above RG equation, we obtain the overall renormalisation constant up to three loop level given by We emphasise that Z A g = Z GG which is introduced in Eq. (3.12) has been discussed in great detail in [59]. The complete UV finite form factor [F A g ] R in terms of this Z A g is This is presented in our recent article [59] up to three loops in the form of hard matching coefficients of soft-collinear effective theory.

Mass Factorisation Kernel
The UV finite form factor contains additional divergences arising from the soft and collinear regions of the loop momenta. In this section, we address the issue of collinear divergences and describe a prescription to remove them. The collinear singularities that arise in the massless limit of partons are removed in the M S scheme using mass factorisation kernel Γ (â s , µ 2 , µ 2 F , z, ǫ). The kernel satisfies the following RG equation : where, P z, µ 2 F are Altarelli-Parisi splitting functions (matrix valued). Expanding P z, µ 2 F and Γ (z, µ 2 F , ǫ) in powers of the strong coupling constant we get and The RG equation of Γ (z, µ 2 F , ǫ), Eq. (3.29), can be solved in dimensional regularisation in powers ofâ s . In the M S scheme, the kernel contains only the poles in ǫ. The solutions up to the required order Γ (3) (z, ǫ) in terms of P (i) (z) can be found in Eq. (33) of [40]. The relevant ones up to three loop, P (0) (z), P (1) (z) and P (2) (z) are computed in the articles [62,63]. For the SV cross section only the diagonal parts of the splitting functions P (i) gg (z) and kernels Γ (i) gg (z, ǫ) contribute.

Soft-Collinear Distribution
The resulting expression from form factor along with operator renormalisation constant and mass factorisation kernel is not completely finite, it contains some residual divergences which get cancelled against the contribution arising from soft gluon emissions. Hence, the finiteness of ∆ A,SV g in the limit ǫ → 0 demands that the soft-collinear distribution, Φ A g (â s , q 2 , µ 2 , z, ǫ), has pole structure in ǫ similar to that of residual divergences. In articles [40] and [41] it was shown that Φ A g must obey KG type integro-differential equation, which we call KG equation, to remove that residual divergences: (3.32) K and G play similar roles as those of K and G, respectively. Also, Φ A g (â s , q 2 , µ 2 , z, ǫ) being independent of µ 2 R satisfy the RG equation This RG invariance and the demand of cancellation of all the residual divergences arising from F A g , Z A g and Γ gg against Φ A g implies the solution of the KG equation as [40,41] where, L A g,i (ǫ) are defined in Eq. (3.17). The z-independent constants G A g,i (ǫ) can be obtained by comparing the poles as well as non-pole terms in ǫ ofφ A g,i (ǫ) with those arising from form factor, overall renormalisation constant and splitting functions. We find where, (3.37) However, due to the universality of the soft gluon contribution, Φ A g must be the same as that of the Higgs boson production in gluon fusion: (3.38) In the above expression, Φ g and G k g,i are written in order to emphasise the universality of these quantities i.e. Φ H g and G H,k g,i can be used for any gluon fusion process, these are independent of the operator insertion. The relevant constants G . This was presented in the article [50]. The G H,k g,i 's required to get the SV cross sections up to N 3 LO are listed below:

SV Cross Sections
In this section, we present our findings of the SV cross section at N 3 LO along with the results of previous orders. Expanding the SV cross section ∆ A,SV g , Eq. (3.7), in powers of a s , we obtain where, Here, we present the results of the pseudo-scalar production cross section up to N 3 LO for the choices of the scale µ 2 R = µ 2 F = q 2 for which the Eq. (4.1) reads with the following ∆ A,SV g,i (z, q 2 ): The SV cross section up to NNLO are in agreement with the existing ones, computed in the article [31,35,36].

Threshold Resummation
Despite the spectacular accuracy of the fixed order results which are defined in power series expansions of the strong coupling constant a s , it is necessary, in certain cases, to resum the dominant contributions to all orders in a s to get more reliable predictions and to reduce the scale uncertainties significantly. In case of threshold corrections, due to soft-gluon emission the fixed order pQCD calculation may yield large threshold logarithms of the kind D i , defined in Eq. (3.6), hence we must resum these contributions to all orders in a s . The resummation of these so-called Sudakov logarithms is usually pursued in Mellin space using the formalism developed in [64,65,76,77]. Alternatively, it is performed in the framework of soft-collinear effective field theory (SCET) [78][79][80][81][82][83][84]. Here, we will discuss this in the context of Mellin space formalism.

Mellin Space Prescription
Under this prescription, the threshold resummation is performed in Mellin-N space where the N -th order Mellin moment is defined with respect to the partonic scaling variable z. In Mellin space, the threshold limit z → 1 corresponds to N → ∞ and the plus distributions D i , Eq. (3.6), take the form ln i−1 N . These logarithmic contributions are evaluated to all orders in perturbation theory by performing the threshold resummation through [64,65,76,77] The component C A,th g depends on both the initial as well as final state particles, though it is independent of the variable N . On the other hand, the remaining part ∆ g,N does not care about the details of the final state particle, it only depends on the initial state partons and the variable N . Being independent of the nature of the final state, ∆ g,N can be considered as a universal quantity which is same for any operator. In addition, it is investigated in the articles [64,65] that it arises solely from the soft parton radiation and it resums all the perturbative contributions a n s ln m N (m ≥ 0) to all orders. Our goal is to calculate the threshold resummation factor C A,th g which encapsulates all the remaining N -independent contributions to the resummed partonic cross section 5.1. Below, we demonstrate the prescription based on our formalism to calculate this quantity C A,th g order by order in perturbation theory.
In the article [41], it was shown how the soft-collinear distribution Φ A g (= Φ g ), Eq. (3.34), captures all the features of the N -space resummation. In this section, we discuss that prescription briefly in the present context. Using the well known identity we can express the soft-collinear distribution 3.34 as where, all the quantities are already introduced in Sec. 3 except K A g,j (ǫ) which is defined through the expansion of K A g , appeared in Eq. (3.32), in powers ofâ s in the following way: The identification of the first plus distribution part of Φ A g , Eq. (5.3), with the factor contributing to the process independent ∆ g,N (q 2 ) has been discussed in the same article [41] which reads In the above expression, the superscript A has been omitted to emphasise the universal nature of these quantities. The remaining part of the Eq. (5.3) along with the other parts, namely, form factor, operator renormalisation constant and mass factorisation kernel in Eq. (3.9) contribute to C A,th g . Expanding this in powers of a s as we determine C A,th g,j up to three loop (j = 3) order which are provided below (with the choice µ 2 R = µ 2 F = q 2 ): by cot 2 β. The mass of the pseudo-scalar is taken to be m A = 200 GeV and the component of the Wilson coefficient C (2) J is considered to be zero due to its unavailability in the literature. We use MSTW2008 [85] parton distribution functions (PDFs) throughout where the LO, NLO and NNLO parton level cross sections are convoluted with the corresponding MSTW2208lo, MSTW2008nlo and MSTW2008nnlo PDFs while for N 3 LO SV cross sections we use MSTW2008nnlo PDFs. The strong coupling constant is provided by the respective PDFs from LHAPDF with α s (m Z ) = 0.1394(LO), 0.12018(NLO) and 0.11707(NNLO).
To estimate the impact of QCD corrections, we define the K-factors as In fig. 1, for LHC13, we plot the pseudo-scalar production cross section as a function of its mass m A . Since we retain the dependence on the m t at the born level, beyond the top pair threshold (τ A > 1), due to change in the functional dependence of τ A one finds a spike at 2m t (left panel). The corresponding K-factors are given in the right panel and are in general found to increase with m A . The NLO correction enhances the LO predictions by as much as 100% for m A = 1 TeV, whereas the NNLO correction adds about an additional 45%. On the other hand the N 3 LO SV correction is found to be about 1.5% of LO for small mass region m A < 300 GeV and for higher m A values the correction at the N 3 LO SV level becomes even smaller, about 0.3% for m A = 1 TeV. In either case, these N 3 LO SV effects show a convergence of the perturbation series. In fig. 3, we present the cross sections as a function of the center of mass energy for fixed m A . A similar pattern is shown both in figs. (1 & 2) where the K-factors increase with m A for a given √ S. The guiding principle for the behaviour of the K-factors in these two cases is the same, namely, as m A approaches √ S, the cross sections are dominated by large soft gluon effects.
Next, we present the scale (µ R , µ F ) uncertainties up to N 3 LO SV in fig. 4 for the choice of m A = 200 GeV. In the left panel, we vary the renormalisation scale µ R between m A /4 and 4m A , keeping µ F = m A fixed. Unlike the Drell-Yan process, for the pseudo-scalar production the renormalisation scale µ R enters even at LO through the strong coupling constant a s . This is identical to the SM Higgs boson production in the gluon fusion channel. This is the main source of large scale uncertainty at LO. It gets significantly reduced when we include NLO and NNLO corrections as expected and it continues to do so at N 3 LO level. In the right panel, we show the factorisation scale uncertainties by varying µ F from m A /4 to 4m A and fixing µ R = m A . Here, the fixed order results show improvement in the reduction of factorisation scale uncertainty from NLO to NNLO. However, due to the lack of parton distribution functions at N 3 LO level and also due to the missing regular contributions from the parton level cross sections, the SV corrections at three loop level do not show any improvement of the factorisation scale uncertainties. In fig. 5 The QCD corrections to pseudo-scalar Higgs production are found to be similar to those of the SM Higgs production due to universal infrared structure of the gluon initiated processes. We give a numerical comparison between their K-factors at various orders. We take m H = m A = 125 GeV and ignore bottom as well as other light quarks and electro weak effects for both the cases. Although the full N 3 LO QCD corrections are already available for the SM Higgs boson, for comparison we take into account only the N 3 LO SV . Table 1 contains the K-factors, defined in Eq. (6.1) up to N 3 LO SV in QCD for both Higgs and pseudo-scalar Higgs boson as a function of √ S. For this mass region, the QCD corrections are positive and hence the K-factors increase with the order in the perturbation theory. Moreover, these K-factors, following the line of argument given before, are found to decrease with √ S but they are identical in both the cases. The difference between the Higgs and the pseudo-scalar cross sections in their respective K-factors is noticed at the second decimal place only. At three loop level, K (3) is found to be around 2.4(2.2) for 7 (14)   The tiny difference between them can be attributed to the presence of an additional operator present in the effective interaction, namely O J which along with the matching coefficient formally enters from NNLO onwards for the gluon initiated processes. For quark anti-quark initiated processes, this contribution vanishes as the quark flavours are massless. The gluon initiated processes involving only O J can contribute at N 4 LO and beyond. However, the interference effects of O G and O J will show up in the gluon initiated processes at NNLO. Thus, the operator O J has non-zero contributions at the lowest order namely at two loop level. However, the presence of such an interference contribution is found to be very small and is the main difference between the SM Higgs and the pseudo-scalar contribution. The QCD corrections through soft and collinear gluon emissions for this interference contribution will be of even higher order and hence will contribute at the three loop level and beyond. In table. 2, we present the Higgs and pseudo-scalar Higgs boson production cross sections up to N 3 LO SV as a function of the scalar mass around 125 GeV. The pseudo-scalar cross section is about twice as big as that of the Higgs boson and the convergence of perturbation series is good and the K-factors are roughly the same for both the cases. Since the predictions are sensitive to the choice of parton density functions, we have estimated the uncertainty resulting from them by choosing the central fit for various well known PDF sets such ABM11 [86] , CT10 [87], MSTW2008 [85] and NNPDF23 [88]. For N 3 LO SV cross sections, however, we use NNLO PDF sets. The corresponding strong coupling constant is directly taken from the LHAPDF. In table. 3, we present the SM Higgs boson and pseudo-scalar Higgs boson production cross sections at NLO, NNLO and N 3 LO SV for LHC13. We find that for NLO, CT10 gives lowest cross section while MSTW2008   gives highest, whereas for NNLO and N 3 LO SV , ABM11 gives lowest and NNPDF23 gives highest. The percentage uncertainty arising from PDF sets at any order is defined as (σ A max − σ A min )/σ A min × 100 where, σ A max and σ A min are the highest and lowest cross sections at any order obtained from the PDFs considered, respectively. This PDF uncertainties in the case of Higgs boson cross sections are about 5.7% at NLO, 8.6% at NNLO and 9.2% at N 3 LO SV . For pseudo-scalar production the cross sections are approximately twice the Higgs cross sections, but the percentage of PDF uncertainties are almost the same.
The SV corrections give a rough estimate of the fixed order (FO) QCD corrections and are often useful in absence of the latter. However, the relative contribution of these SV corrections to the full FO results crucially depends on the kinematic region and in some cases on the process under study. For the SM Higgs or pseudo-scalar Higgs boson with a mass of about 125 GeV, it is far from the threshold region τ = m 2 H /S → 1 for √ S = 13 TeV. Since, the parton fluxes corresponding to this mass region are very high, apart from the threshold logarithms the contributions of the regular terms as well as of other subprocesses present in the FO corrections are expected to be reasonably very high. For Higgs or pseudo- scalar, the prediction at NLO SV level differs from the LO by only a few percent whereas the regular terms at NLO contribute significantly and increase LO prediction by about 70%. Similar is the case even at NNLO. Thus the SV corrections poorly estimate the FO ones, however, if we redefine the hadron level cross sections without affecting the total cross sections in such a way that the parton fluxes peak near the threshold region [46,54,89], then the SV contributions can be shown to dominate over the regular ones. This is due to arbitrariness involved in splitting the parton level cross section in terms of threshold enhanced and regular ones. Using a regular function G(z), we can write the hadronic cross section as σ A (τ ) = σ A,(0) a,b=q,q,g where ∆(z)/G(z) can be decomposed as ∆(z)/G(z) = ∆ SV (z) +∆ hard (z) (6.3) In the above equation the ∆ SV is independent of G(z) (if lim z→1 G(z) → 1) and contains only distributions, whereas the hard part∆ hard is modified due to G(z). Hence the SV part of the cross section at the hadron level depends on the choice of G(z). For the peculiar choice G(z) = z 2 , the ∆ SV dominates over∆ hard in such a way that almost the entire NLO and NNLO corrections (Eq. (6.2)) results from ∆ SV alone. As was noted earlier G(z) = 1 corresponds to the standard SV contribution. Note that the flux Φ ab is modified to Φ mod ab (y) = Φ ab (y)G(τ /y) which is responsible for this behaviour. We may denote the SV cross sections thus obtained with these modified fluxes as NLO (sv) , NNLO (sv) and N 3 LO (sv) while those obtained with the normal fluxes as NLO sv , NNLO sv and N 3 LO sv . In fig. 9, we depict the comparison between the SV cross sections obtained from the modified parton fluxes using G(z) = z 2 and the normal fixed order results that are obtained from the standard parton fluxes, for both the SM Higgs boson (left panel) and the pseudo-scalar (right panel). We notice that the SV results are significantly closer to the corresponding fixed order ones. Incidentally, this agreement is good for NLO as well as for NNLO where different subprocesses appear, and also for several values of √ S where the integration range over the parton fluxes is different. While this could be purely accidental, this good agreement might hint some subtle aspect hidden and might be useful in the phenomenology.
Motivated by the above observation, one can convolute the perturbative coefficients ∆ SV with the modified parton fluxes Φ mod ab (y) for the choice of G(z) = z 2 to get N 3 LO (sv) which could approximate the full N 3 LO result. This way, we present in fig. 9, the SV corrections obtained using G(z) = 1 and G(z) = z 2 for Higgs as well as pseudo-scalar Higgs boson productions.

Conclusions
In this paper, using the recently available pseudo-scalar form factors up to three loops and the third order soft function from the real radiations, a complete N 3 LO threshold correction to the production of pseudo-scalar at the LHC has been obtained. The computation is performed using z space representation of resummed cross section. We have exploited the universal structure of soft function that appears in scalar Higgs boson production at the LHC. We found that the singularities resulting from soft and collinear regions in the virtual diagrams cancel against those from the universal soft functions as well as from mass factorisation kernels. Using our approach, we have also computed the process dependent coefficient that appears in the threshold resummed cross section. This will be useful for resummed predictions at N 3 LL in QCD. Using threshold corrected N 3 LO results, we have presented a detailed phenomenological study of the pseudo-scalar production at the LHC for various center of mass energies as a function of its mass. While the third order corrections are small, they play an important role in reducing the theoretical uncertainty resulting from renormalisation scale. In addition, we have made a detailed comparison against scalar Higgs boson production and found their corrections are very close to each other confirming the universal behaviour of the QCD effects even though the operators responsible for their interactions with gluons are very different.