Pseudo-scalar Higgs boson production at threshold N3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^3$$\end{document}LO and N3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^3$$\end{document}LL QCD

We present the first results on the production of pseudo-scalar Higgs boson through gluon fusion at the LHC to N3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^3$$\end{document}LO in QCD taking into account only soft-gluon effects. We have used the effective Lagrangian that describes the coupling of the pseudo-scalar Higgs boson with the gluons in the large top quark mass limit. We have used quantities that have recently become available, namely the three-loop pseudo-scalar Higgs boson 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 a threshold resummation to N3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^3$$\end{document}LL in QCD. Phenomenological impact of these threshold N3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^3$$\end{document}LO corrections to pseudo-scalar Higgs boson production at the LHC is presented and their role in the reduction of the renormalization 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 on a 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 − , Z Z 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 the study of 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 cannot explain the baryon asymmetry in the Universe, dark matter, the neutrino mass etc. There are several extensions of the SM, motivated to address these issues. The minimal version of the 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 threeloop level is consistent [21][22][23][24] with the recently observed Higgs boson at the LHC. The efforts to test the predictions of MSSM or its variants already are under way and the cur-rent run at the LHC will shed more light on them. One of them could be to look for a CP odd Higgs boson in the gluon fusion through heavy fermions as its coupling is appreciable for small and moderate tan β, the ratio of vacuum expectation values v i , i = 1, 2. In addition, a large gluon flux can boost the cross section.
Since the leading order production mechanism of the pseudo-scalar Higgs boson of mass m A is through heavy quarks, the cross section is not only proportional to tan β but also the square of the strong coupling constant. Like the scalar Higgs boson in the SM, the leading order prediction of the pseudo-scalar Higgs boson production at the hadron colliders suffers from theoretical uncertainties that are large, particularly due to the renormalization scale μ R arising from the strong coupling constant, and are mild due to the factorization scale μ F in the parton distribution functions. Predictions based on one-loop perturbative Quantum Chromodynamics (pQCD) corrections [25][26][27][28] 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 much as 67 %. The effective theory approach in the large top quark mass limit provides an opportunity to go beyond NLO. Such an approach [28,29] in the case of scalar Higgs boson production [30][31][32] turned out to be most successful, as the finite quark mass effects at NNLO level were found to be within 1 % [33][34][35][36][37]. NNLO predictions in the effective theory for the production of pseudo-scalar Higgs boson at the hadron colliders are already available [32,38,39]. The calculations were further performed after considering the finite mass of the top quarks in [40,41]. The NNLO correction increases the NLO cross section by about 15 % and reduces the scale uncertainties to about 15 %. Due to the large gluon flux at the threshold, namely when the mass of pseudoscalar Higgs boson approaches the partonic center 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 perturbation theory provides the solution to this problem. The systematic predictions based on the next-to-next-to-leading log (NNLL) resummed result [42][43][44][45][46][47][48][49][50] demonstrate the reliability of the approach and also reduce the scale uncertainties.
A complete calculation at NNLO [30][31][32] and leading logarithms at N 3 LO in the threshold limit [43][44][45][46][47] and NNLL soft-gluon resummation [42] for the scalar Higgs boson production have been 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 the δ(1 − z) contribution at N 3 LO in the threshold limit [51] was the first among them. This was confirmed independently in [52]. Later on, the sub-leading collinear log-arithms were computed in [53,54]. A spin-off of the result presented in [51] is the computation of the N 3 LO prediction for the Drell-Yan production [52,55,56] 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 [57] and also in association with vector boson [58] at the hadron colliders. Later, along the same direction, the rapidity distribution of the Higgs boson in gluon fusion [59], Drell-Yan [59], and the Higgs boson in bottom quark annihilation [60] were obtained at threshold N 3 LO QCD.
A milestone in this direction was achieved by Anastasiou et al. who have now accomplished the complete N 3 LO prediction [61] 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 [56,62]. In [63], the SM Higgs resummation is performed for the first time in the soft-collinear effective field theory (SCET) approach to N 3 LL. We also note that in the exact theory, including the finite top mass effects, the three-loop virtual corrections are already available in [64] and the full result is yet to be computed.
While the next step in the wish list is to obtain the N 3 LO predictions for the pseudo-scalar Higgs boson 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 Higgs boson, computed between partonic states. One-and two-loop results for them between gluon states were computed for NNLO production cross section [38,39,65], the analytical results up to two-loop level can be found in [65]. These were computed in dimensional regularization where the space-time dimension is d = 4 + . Threshold corrections to pseudo-scalar Higgs boson 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 [66], we obtained the three-loop form factors of the effective composite operators between quark and gluon states along with the lower order ones to the desired accuracy in . In the present article we will describe how threshold corrections at N 3 LO level can be obtained from the formalism developed in [45,46] using the available information on recently computed three-loop form factor of the pseudo-scalar Higgs boson [66], the universal soft-collinear distribution [55], and the operator renormalization constant [66][67][68] and the mass factorization kernels [69,70] known to three-loop level. In addition, we compute a third order correction to the N -independent part of the resummed cross section [71,72] using our formalism [45,46]. We also present the numerical impact of our findings with a brief conclusion.
The underlying effective theory is discussed in the Sect. 2. This is followed by a short description of the formalism which has been employed to compute the soft-plus-virtual cross section in Sect. 3. We present the analytical results of these findings in the Sect. 4 up to N 3 LO in QCD. In Sect. 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 our concluding remarks, in Sect. 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 Higgs boson 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 [73] describing the interaction between pseudoscalar Higgs boson χ A and the QCD particles in the infinitely large top quark mass limit is given by where the two operators are defined as (2. 2) The symbols G μν a and ψ represent the gluonic field strength tensor and the quark field, respectively. The Wilson coefficients C G and C J of the two operators are the consequences of integrating out the heavy quark loop in the effective theory. C G does not receive any QCD corrections beyond one loop because of the Adler-Bardeen theorem [74], whereas C J starts only at second order in the strong coupling constant. These Wilson coefficients are given by [73] β i 's are the coefficients of the QCD β function [75].

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, renormalized at the factorization scale μ F . Here, Δ A ab ( τ y , m 2 A , μ 2 R , μ 2 F ) are the partonic level cross sections, for the sub-process initiated by the partons a and b, computed after performing the overall operator UV renormalization at scale μ R and the mass factorization 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 softgluon contributions to the pseudo-scalar Higgs boson 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) renormalized virtual part and performing the mass factorization using appropriate counter terms. This combination is often called the soft-plus-virtual (SV) cross section, whereas the remaining portion is known as the 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 ) contain only distributions of the 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 [45,46], using 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 Higgs boson. 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 meromorphic functions of the variable N (see [71,72]) 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 renormalization constant Z A g (â s , μ 2 R , μ 2 , ), the softcollinear distribution Φ A g (â s , q 2 , μ 2 , z, ) arising from the real radiations in the partonic sub-processes and the mass factorization kernels Γ gg (â s , μ 2 F , μ 2 , z, ). In terms of the above-mentioned quantities it takes the following form, as presented in [46,55,57]: In the subsequent sections, we will demonstrate the methodology to get these ingredients to compute the SV cross section of pseudo-scalar Higgs boson production at N 3 LO. In particular, in Sect. 3.1, we show how to obtain the relevant form factor that goes into the computation of pseudo-scalar production through gluon fusion, using the form factors obtained in our earlier work [66]. In Sect. 3.2, we calculate the operator renormalization constant from the relevant UV anomalous dimensions. Mass factorization is discussed very briefly in Sect. 3.3. Finally, we explain the relevant soft-collinear distribution in Sect. 3.4.

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 Higgs boson production through gluon fusion, we need to consider two operators O G and O J , defined in Eq. (2.2), which yield in total two form factors. The unrenormalized gluon form factors at O(â n s ) are defined [66] througĥ In terms of these quantities, the full matrix element and the full form factors can be written as series expansions inâ s : 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 the one-loop level.
The form factor for the production of a pseudo-scalar Higgs boson through gluon fusion,F A,(n) g , can be written in terms of the two individual form factors, Eq. (3.11), as follows: (3.12) In the above expression, the quantities Z i j (i, j = G, J ) are the overall operator renormalization constants which are required to introduce in the context of UV renormalization. These are discussed in our recent article [66] 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 [66]. Using those results we obtain the three-loop form factor for the pseudo-scalar Higgs boson production through gluon fusion. In this section, we present the unrenormalized form factorsF A,(n) g up to three loops 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:  (3.14) The results up to two-loop level is consistent with the existing ones [65] and the three-loop result is the new one.

Operator renormalization constant
The strong coupling constant renormalization through Z a s is not sufficient to make the form factor F A g completely UV finite; one needs to perform additional renormalization to remove the residual UV divergences. This additional renormalization is called the overall operator renormalization; it is performed through the constant Z A g . This is determined by solving the underlying RG equation: In the above expression, the γ A g,i are the UV anomalous dimensions where the components are defined through the expansion in powers of a s : The γ A g,i are determined by explicitly comparing the results of the form factors, Eq. (3.14), against the universal decomposition [65] of the form factors in terms of soft, cusp, collinear and UV anomalous dimensions. The universal decomposition follows from the property of the form factor that it satisfies the K G-differential equation [76][77][78][79][80]. This is a direct consequence of the facts that QCD amplitudes exhibit the factorization property, and the gauge and renormalization group (RG) invariances. The γ A g,i up to three loops (i = 3) are obtained: Let us emphasize and note that the γ A g,i 's are found to satisfy is the usual QCD β-function. For a more elaborate discussion, see the recent article [66] (also see [67,68]). Using the results of γ A g,i from Eq. (3.17) and solving the above RG equation, we obtain the overall renormalization constant up to three-loop level: We emphasize that Z A g = Z GG , which is introduced in Eq. (3.12), has been discussed in great detail in [66]. The complete UV finite form factor This is presented in our recent article [66] up to three loops in the form of hard matching coefficients of soft-collinear effective theory.

Mass factorization 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 factorization kernel. The kernel contains only the poles in and Altarelli-Parisi splitting functions. For the SV cross section only the diagonal parts of the splitting functions and kernels contribute.

Soft-collinear distribution
The resulting expression from form factor along with operator renormalization constant and mass factorization kernel is not completely finite, it contains some residual divergences which get canceled against the contribution arising from soft-gluon emissions. Hence, the finiteness of Δ A,SV g in the limit → 0 requires that the soft-collinear distribution, Φ A g (â s , q 2 , μ 2 , z, ), has a pole structure in similar to that of residual divergences. In the articles [45,46] it was shown that Φ A g must obey a K G type integro-differential equation, which we call the K G equation [45,46], to remove the residual divergences. 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: In the above expression, Φ g is written in order to emphasize the universality of these quantities i.e. Φ H g can be used for any gluon fusion process, these are independent of the operator insertion. The result up to three-loop level was presented in the article [32,55]. The three-loop one was obtained by using the results of the Higgs boson production cross section at threshold at N 3 LO QCD [51]. This completes all the ingredients required to compute the SV cross section up to N 3 LO that are presented in the next section.

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 Higgs boson 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 [32,38,39]. The result at N 3 LO i.e. Δ A,SV g,3 is the new one computed in this article for the first time.

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 [71,72,81,82]. Alternatively, it is performed in the framework of SCET [83][84][85][86][87][88][89]. 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 [71,72,81,82] The component C A,th g depends on both the initial and the 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 the same for any operator. In addition, it is investigated in the articles [71,72] 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 [46], it was shown how the soft-collinear distribution Φ A g (= Φ g ) 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 as where A A g is the cusp anomalous dimension corresponding to the gluonic operator. All the poles in are contained within K A g and the finite terms are dumped into G A g . The components K A g, j ( ) are defined through the expansion of K A g 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 [46] which reads In the above expression, the superscript A has been omitted to emphasize the universal nature of these quantities. The remaining part of the Eq. 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 ): The above new result of C A,th g,3 along with the universal factor Δ g,N provide the threshold resummed cross section of the pseudo-scalar Higgs boson production at N 3 LL accuracy. A more elaborate discussion of this prescription to perform threshold resummation will be presented elsewhere.

Numerical impact of SV cross section
In this section, we present our findings on the numerical impact of threshold N 3 LO predictions in QCD for the production of a pseudo-scalar Higgs boson at the LHC and also make comparison with the corresponding results for the SM Higgs boson. While we have all the ingredients up to three-loop level, the one of the Wilson coefficients, namely, C (2) J is known only to two-loop level. Note that C G is exact due to Adler-Bardeen theorem [74]. Due to the unavailability of C (2) J in the literature, we discuss the impact of missing three-loop contribution later on in this section by varying this quantity. As we are interested in quantifying the QCD effects, we assume that pseudo-scalar Higgs boson couples only to top quarks. Hence, the dominant contribution resulting from bottom quark initiated processes can be included in a systematic way in our numerical study but we do not perform it here. Moreover, our predictions are based on the effective theory approach where the top quarks are integrated out and we have only light quarks. Like in the case of predictions for the scalar Higgs boson production in the effective theory, for the pseudo-scalar Higgs boson production we multiply the Born cross section computed using the finite top mass (m t = 172.5 GeV) with higher orders which are obtained in the effective theory. Without loss of generality, we normalize the cross section by cot 2 β. The mass of the pseudo-scalar Higgs boson is taken to be m A = 200 GeV. We use MSTW2008 To estimate the impact of QCD corrections, we define the K-factors as In Fig. 1, for LHC13, we plot the pseudo-scalar Higgs boson 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). We note here that the effective theory (EFT) formalism formally holds only for the pseudo-scalar masses up to top pair threshold. However, from the knowledge of the QCD corrections to the Higgs boson production up to NNLO, we notice that below the top pair threshold, the difference between the results of EFT and the finite top contributions is about 5 % and is even smaller at NNLO, about 1 %. We assumed the same to hold in the QCD corrections for pseudo-scalar production for scalar masses above the top pair threshold. 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. 2, we present similar results but only for pseudo-scalar masses below the top pair threshold, where the effective theory approximation works very well. In Fig. 3, we present the cross sections as a function of the center of mass energy √ S of the incoming protons at the LHC. The increase in the cross sections (left panel) with √ S is simply because of the increase in the corresponding parton fluxes for any given m A . On the contrary, the corresponding K-factors (right panel) increase with decreasing √ S for fixed m A . A similar pattern is shown both in Figs. 1 and 2 where the K-factors increase with m A for a given √ S. The guiding principles for the behavior of the K-factors in these two  cases are the same, namely, as m A approaches √ S, the cross sections are dominated by large soft-gluon effects.
The QCD corrections to pseudo-scalar Higgs boson 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 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 the Higgs and the 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 cases. The difference between the Higgs and the pseudoscalar Higgs boson 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 the 7(14) TeV case.
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-antiquark initiated processes, this contribution vanishes as the quark flavors 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 Higgs boson 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 Higgs boson 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 cases.
We have also studied the impact of missing three-loop contribution to C J i.e. C (2) J . At higher orders starting from N 3 LO SV onwards, the second term, C (2) J , in the Wilson coefficient C J can give a non-zero contribution. To estimate the numerical impact of this term, we assume the following form of C (2) J : and we vary the parameters a, b, and c in the range [−10, 10]. We found that the contribution of such a C (2) J term changes the cross sections only at the third decimal place and hence we ignore its contribution in the rest of our phenomenological study.
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 [91], CT10 [92], MSTW2008 [90] and NNPDF23 [93]. For N 3 LO SV cross sections, however, we use NNLO PDF sets. The corresponding strong coupling constant is directly taken from the LHAPDF [94]. 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 the lowest cross section, while MSTW2008 gives the highest, whereas for NNLO and N 3 LO SV , ABM11 gives the lowest and NNPDF23 gives the highest. The percentage uncertainty arising from PDF sets at any order is defined where σ A max and σ A min are the highest and lowest cross sections at any order obtained from the PDFs considered, respectively. These 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 pseudoscalar Higgs boson production the cross sections are approximately twice the Higgs cross sections, but the percentages 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 that we study. For the SM Higgs or pseudo-scalar Higgs boson with a mass of about 125 GeV, it is far from the threshold region 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 sub-processes present in the FO corrections are expected to be reasonably very high. For a Higgs or pseudo-scalar Higgs boson, 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 [51,59,95], 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 where Δ(z)/G(z) can be decomposed as dominates overΔ hard in such a way that almost the entire NLO and NNLO corrections (Eq. (6.3)) 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 behavior. 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. 4, 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 Higgs boson (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 sub-processes 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 at some subtle hidden aspect and might be useful in the phenomenology.
Motivated by the above observation, one can convolute the perturbative coefficients Δ (3) 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. 4 the SV corrections obtained using G(z) = 1 and G(z) = z 2 for Higgs as well as pseudo-scalar Higgs boson productions.
Next, we present the scale (μ R , μ F ) uncertainties up to N 3 LO SV in Fig. 5 for the choice of m A = 200 GeV. In the rest of our numerical analysis for studying the scale uncertainties, we simply use the modified parton fluxes for the   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 N 3 LO SV cross sections do not show any improvement of the factorization scale uncertainties. However, we observe that with the modified parton fluxes, the factorization scale uncertainties in N 3 LO (SV) get significantly reduced compared to N 3 LO SV . In Fig. 6, we show the combined effect of μ R and μ F scale uncertainties by varying the scale μ between m A /4 and 4m A , where μ = μ R = μ F . Here, the NNLO cross sections show a good improvement over the NLO ones against the scale variations, while the N 3 LO (SV) cross sections are found to be more stable than the NNLO ones. Further, we also study the renormalization and factorization scale variations of both the cross sections for the pro-duction of the SM Higgs boson and the pseudo-scalar Higgs boson for m H = m A = 125 GeV by varying them between m A /4 and 4m A . In Fig. 7, the renormalization scale uncertainties are given for Higgs boson (left panel) and for pseudoscalar Higgs boson (right panel), for μ F = m A = m H . In Fig. 8, we present similar results but for the factorization scale uncertainties keeping μ R = m H = m A . Moreover, in Fig. 9, we present the combined effect by varying μ = μ R = μ F . The pattern of the results for μ R and μ F and the combined variations are similar to the earlier analysis for m A = 200 GeV where the renormalization scale uncertainties get stabilized further after including the third order threshold corrections, while the scale uncertainties due to μ = μ R = μ F variation get significantly improved at N 3 LO (SV) than at NNLO.

Conclusions
In this paper, using the pseudo-scalar Higgs boson form factors that have recently become available 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 a pseudo-scalar Higgs boson at the LHC has been obtained. The computation is performed using the z space representation of the resummed cross section. We have exploited the universal structure of the 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 factorization 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 pseudoscalar Higgs boson 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 renormalization 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 behavior of the QCD effects even though the operators responsible for their interactions with gluons are very different.