NLO QCD corrections to pseudoscalar Higgs production in the MSSM

We present a calculation of the two-loop quark-squark-gluino contributions to pseudoscalar Higgs boson production via gluon fusion in the MSSM. We regularize the loop integrals using the Pauli-Villars method, and obtain explicit and compact analytic results based on an expansion in the heavy particle masses. Our results - valid when the pseudoscalar Higgs boson is lighter than squarks and gluinos - can be easily implemented in computer codes for an efficient and accurate determination of the pseudoscalar production cross section.


Introduction
With the coming into operation of the Large Hadron Collider (LHC), a new era has begun in the search for the Higgs boson(s). At the LHC the main production mechanism for the Standard Model (SM) Higgs boson, H SM , is the loop-induced gluon fusion mechanism [1], gg → H SM , where the coupling of the gluons to the Higgs is mediated by loops of colored fermions, primarily the top quark. The knowledge of this process in the SM includes the full next-to-leading order (NLO) QCD corrections [2,3,4], the next-to-next-to-leading order (NNLO) QCD corrections [5] including finite top mass effects [6], soft-gluon resummation effects [7], an estimate of the next-to-next-to-next-to-leading order (NNNLO) QCD effects [8] and also the first-order electroweak corrections [9,10,11].
The Higgs sector of the Minimal Supersymmetric extension of the Standard Model (MSSM) consists of two SU (2) doublets, H 1 and H 2 , whose relative contribution to electroweak symmetry breaking is determined by the ratio of vacuum expectation values of their neutral components, tan β ≡ v 1 /v 2 . The spectrum of physical Higgs bosons is richer than in the SM, consisting of two neutral CP-even bosons, h and H, one neutral CP-odd boson, A, and two charged scalars, H ± . The couplings of the MSSM Higgs bosons to matter fermions differ from those of the SM Higgs, and they can be considerably enhanced (or suppressed) depending on tan β. As in the SM case, the gluon-fusion process is one of the most important production mechanisms for the neutral Higgs bosons, whose couplings to the gluons are mediated by top and bottom quarks and their supersymmetric partners, the stop and sbottom squarks.
In the case of the CP-even bosons h and H the gluon-fusion cross section in the MSSM is known at the NLO in QCD. 1 The contributions arising from diagrams with quarks and gluons can be obtained from the corresponding SM results with an appropriate rescaling of the Higgs-quark couplings. The contributions arising from diagrams with squarks and gluons were first computed under the approximation of vanishing Higgs mass in ref. [13]. The complete top/stop contributions, including the effects of stop mixing and of the two-loop diagrams involving gluinos, were computed under the same approximation in ref. [14], and the result was cast in a compact analytic form in ref. [15]. Later calculations aimed at the inclusion of the full Higgs-mass dependence in the squark-gluon contributions, which are now known in a closed analytic form [16,17,18,19].
The approximation of vanishing Higgs mass in the contributions of two-loop diagrams allows for compact analytic results that can be implemented in computer codes for a fast and efficient evaluation of the Higgs production cross section. For what concerns the top-gluon contributions, the effect of such approximation on the result for the cross section has been shown [20,19] to be limited to a few percent, as long as the Higgs mass is below the threshold for creation of the massive particles running in the diagrams (in this case, the top quarks). While this condition may also apply to the two-loop diagrams involving top, stop and gluino, it obviously does not apply to the corresponding diagrams involving the bottom quark, whose contribution can be relevant for large values of tan β. For the latter diagrams the dependence on the Higgs mass should in principle be retained, which has proved a rather daunting task. A calculation of the full quark-squark-gluino contributions via a combination of analytic and numerical methods was presented in ref. [21] (see also ref. [22]), but neither explicit analytic results nor a public computer code have been made available so far. However, ref. [23] presented an evaluation of the bottom-sbottom-gluino diagrams based on an asymptotic expansion in the large supersymmetric masses that is valid up to and including terms of where m φ denotes a Higgs boson mass and M denotes a generic superparticle mass. This expansion should provide a good approximation to the full result, at least comparable to the one obtained for the top-stop-gluino diagrams, as long as the Higgs boson mass is below all the heavy-particle thresholds. An independent calculation of the bottom-sbottom-gluino contributions, restricted to the limit of a degenerate superparticle mass spectrum, was also presented in ref. [24], confirming the results of ref. [23].
In the case of the CP-odd boson A the calculation of the production cross section is somewhat less advanced. Due to the structure of the A-boson coupling to squarks, only loops of top and bottom quarks contribute to the cross section at LO, with the bottom loops being dominant for even moderately large values of tan β. In the limit of vanishing A-boson mass, m A , the contributions from diagrams with quarks and gluons were computed at NLO in ref. [25] and at NNLO in ref. [26] (see also ref. [27]). For arbitrary values of m A the NLO contributions arising from two-loop diagrams with quarks and gluons, as well as from one-loop diagrams with emission of a real parton, were computed in ref. [3]. Supersymmetric particles contribute to the cross section at NLO through two-loop diagrams involving quarks, squarks and gluinos. The top-stop-gluino contributions were computed in ref. [28] in the limit of vanishing m A . The analytic result for generic values of the stop and gluino masses was deemed too voluminous to be explicitly displayed in ref. [28], and was instead made available in the fortran code evalcsusy.f [14]. On the other hand, the two-loop bottom-sbottom-gluino contributions, which can be relevant for large values of tan β, have never been directly computed so far.
In this paper we aim to reduce the gap in accuracy between the available NLO calculations of the production cross sections for CP-odd and CP-even Higgs bosons of the MSSM, exploiting the techniques we developed for computing the top-stop-gluino [15] and bottom-sbottom-gluino [23] contributions in the CP-even case. In particular, we present an evaluation of the two-loop top-stopgluino contributions to the pseudoscalar production cross section valid up to and including terms of . We show how the terms of order zero in m 2 A can be cast in an extremely compact analytic form, fully equivalent to the result of ref. [28], and we investigate the effect of the first-order terms. We also evaluate the same contributions via an asymptotic expansion in the large superparticle masses, valid up to and including terms of O(m 2 A /M 2 ) and O(m 2 t /M 2 ). While the latter result is valid for m t , m A ≪ M but does not assume a hierarchy between m t and m A , the former is expected to provide a better approximation in the region with m A < m t and relatively light superparticles, M ≃ m t . As a byproduct, we also obtain a result for the bottom-sbottom-gluino contributions valid up to and including terms of O(m 2 b /m 2 A ) and O(m b /M ). Finally, we compare our results for the bottom-sbottom-gluino contributions to both CP-even and CP-odd Higgs production cross sections with those obtained in the effective-Lagrangian approximation of refs. [29,30]. A non-trivial technical issue that arises in the calculation of the pseudoscalar production cross section is the treatment of the Dirac matrix γ 5 -an intrinsically four-dimensional object -within regularization methods defined in a number of dimensions n d = 4 − 2ǫ. The original calculation of the two-loop quark-gluon contributions of ref. [3] was performed in Dimensional Regularization (DREG), employing the 't Hooft-Veltman (HV) prescription [31] for the γ 5 matrix and introducing a finite multiplicative renormalization factor [32] to restore the Ward identities. In ref. [28] the calculation of the top-gluon and top-stop-gluino contributions to the Wilson coefficient in the relevant effective Lagrangian was performed both in DREG and in Dimensional Reduction (DRED), which, differently from DREG, preserves supersymmetry (SUSY). The latter method does not require the introduction of finite renormalization factors, but it involves additional subtleties concerning the treatment of the Levi-Civita symbol ε µνρσ .
In our calculation of the quark-squark-gluino contributions we avoided all problems related to the treatment of γ 5 by employing the Pauli-Villars regularization (PVREG) method. Being defined in four dimensions, the PVREG method respects both SUSY and the chiral symmetry, therefore no symmetry-restoring renormalization factors need to be introduced. We tested our implementation of PVREG by computing the top-gluon contributions via an asymptotic expansion in the top quark mass, and recovering the result obtained in DREG in refs. [3,17]. As a further cross check, we also computed the quark-squark-gluino contributions using the DREG procedure outlined in ref. [32], and found agreement with the result that we obtained in PVREG.
The paper is organized as follows: in section 2 we summarize general results on the cross section for pseudoscalar Higgs boson production via gluon fusion. In section 3 we outline our implementation of the PVREG method. Section 4 contains our explicit results for the NLO contributions arising from both top-stop-gluino and bottom-sbottom-gluino diagrams, as well as a discussion of suitable renormalization schemes for the bottom contributions and a comparison with the results obtained in the effective-Lagrangian approximation. In section 5 we assess the validity of the expansion in powers of m 2 A in the top contributions, and discuss the numerical relevance of the different NLO contributions. In the last section we present our conclusions. We also include, for completeness, an Appendix in which we present the NLO contributions from one-loop diagrams with emission of a real parton.

Pseudoscalar Higgs boson production via gluon fusion at NLO
In this section we recall for completeness some general results on pseudoscalar Higgs boson production via gluon fusion. The hadronic cross section at center-of-mass energy √ s can be written as where τ A = m 2 A /s, µ F is the factorization scale, f a,h i (x, µ F ) the parton density of the colliding hadron h i for the parton of type a (for a = g, q,q), andσ ab the cross section for the partonic subprocess ab → A + X at the center-of-mass energyŝ = x 1 x 2 s = m 2 A /z. The partonic cross section can be written in terms of the LO contribution σ (0) and a coefficient function G ab (z) aŝ The LO term can be written as where G µ is the muon decay constant and α s (µ R ) is the strong gauge coupling expressed in the MS renormalization scheme at the scale µ R . H A is the form factor for the coupling of the pseudoscalar A with two gluons, which we decompose in one-and two-loop parts as Due to the structure of the pseudoscalar coupling to squarks (see section 4), only diagrams involving top or bottom quarks contribute to the one-loop form factor H 1ℓ A . The latter can be decomposed into top and bottom contributions as where T F = 1/2 is a color factor, τ q = 4 m 2 q /m 2 A , and We recall the behavior of K 1ℓ in the limit in which the pseudoscalar mass is much smaller or much larger than twice the mass of the particle running in the loop. In the first case, i.e. τ ≫ 1, which may apply to the top contribution if m A is relatively small, while in the opposite case, i.e. τ ≪ 1, which is relevant for the bottom contribution, The analytic continuation of K 1ℓ (τ ) corresponds to the replacement m 2 A → m 2 A +iǫ , thus the imaginary part of eq. (8) can be recovered via the replacement ln(−4/τ ) → ln(4/τ ) − iπ.
The coefficient function G ab (z) in eq. (2) can be decomposed, up to NLO terms, as with the LO contribution given only by the gluon-fusion channel: The NLO terms include, besides the gg channel, also the one-loop induced gq and qq channels: where the LO Altarelli-Parisi splitting functions are In the equations above C A = N c and C F = (N 2 c − 1)/( 2 N c ) (N c being the number of colors), β 0 = (11 C A − 2 N f )/6 (N f being the number of active flavors) is the one-loop β-function of the strong coupling in the SM, and The two-loop virtual contributions to gg → A, regularized by the infrared-singular part of the contributions from real gluon emission in the one-loop gluon fusion channel, gg → Ag, are displayed in the first line of eq. (11). The second line of that equation contains the non-singular contributions from real gluon emission. Eq. (12) contains the contributions due to the one-loop quark-antiquark annihilation channel, qq → Ag, and to the one-loop quark-gluon scattering channel, gq → qA. General expressions for the functions R gg , R qq , R qg in the case of pseudoscalar production are collected in the Appendix. The two-loop form factor H 2ℓ A receives contributions from diagrams involving quarks and gluons, as well as from diagrams involving quarks, squarks and gluinos. The contributions from two-loop diagrams with quarks and gluons were first computed in ref. [3], and later confirmed in ref. [17]. The contribution to H 2ℓ A arising from top-stop-gluino diagrams was computed in ref. [28] in the limit of vanishing pseudoscalar mass. For what concerns the contribution arising from bottom-sbottom-gluino diagrams, no genuine two-loop calculation has been available so far. In the following sections we present our calculation of both kinds of quark-squark-gluino contributions.

In our computation of H 2ℓ
A we regularized the loop integrals using the PVREG method. For the purposes of this computation, the main advantage of PVREG is the fact that all the Lorentz indices remain strictly 4-dimensional, thus the γ 5 matrices anticommute with the other gamma matrices and the trace on a string of gamma matrices can be taken using the standard 4-dimensional relations. We recall that in PVREG, given an ultraviolet (UV) divergent integral I(q, m 2 ) where q and m 2 denote collectively the external momenta and masses, its regularized version is constructed as In the equation above the original integral I(q, m 2 ) is combined with a number n of replicas, weighted by coefficients c i , in which some of the masses of the original integral are replaced by the PV mass regulators (m i ), in such a way that the regularized integral is finite if m i are kept finite, but tends to infinity as m i → ∞. The number of added terms, as well as the relation that the coefficients c i should satisfy in order to make I R convergent, depend on the divergent nature of the original integral.
If the latter is only logarithmically divergent, a single subtraction is sufficient to construct I R , i.e., n = 1, For what concerns the infrared (IR) divergences associated to massless particles, in PVREG they are regularized by giving a fictitious mass λ to the massless particle, and later considering the limit λ → 0. All the diagrams contributing to the virtual NLO contributions to pseudoscalar production are at most logarithmically UV-divergent, therefore a single subtraction is sufficient to make them convergent. In this case, PVREG reduces to subtracting from the original diagrams the same diagrams with some of the masses replaced by M P V , and then taking the limit M P V → ∞. In the case of the top-gluon contributions also the limit λ → 0 must be taken on the fictitious gluon mass. In the present calculation, taking the relevant limits for the mass regulators does not introduce additional complications with respect to the same calculation performed in DRED or DREG. This is due to the fact that we are computing the two-loop diagrams via an asymptotic expansion, so that the final result is expressed in terms of two-loop vacuum integrals with different masses and of one-loop integrals. Both kinds of terms are fully known analytically, including all the relevant limits when one or more masses are sent to infinity or to zero. The asymptotic expansion of the relevant diagrams is generated following the procedure described in ref. [23], which amounts to adding to and subtracting from each diagram its IR-divergent part. As discussed in that paper, a diagram minus its IR-divergent part can be evaluated via a Taylor expansion in the external momenta (being this combination IR finite by construction) while its remaining IR-divergent part, which is expressed as a product of two one-loop integrals, must be evaluated exactly.
In order to test our implementation of PVREG we first considered the two-loop top-gluon contributions. These contributions can be split in two parts, one proportional to C F and the other proportional to C A . The latter, which stems from the non-abelian nature of SU (3), is not IR finite but contains a soft and collinear divergence that factorizes on the lowest-order cross-section. In DREG, this IR divergence appears as a 1/ǫ 2 pole multiplying the top contribution to σ (0) . We computed the top-gluon contributions via an asymptotic expansion in the top mass up to and including terms O(m 8 A /m 8 t ). The IR divergences are regularized by giving a mass λ to the gluon, while the UV divergences are regularized by subtracting to any term a replica in which λ is replaced by M P V . The final result is then obtained taking the limits M P V → ∞ and λ → 0. We were able to reproduce in PVREG the known result for the top-gluon contributions obtained in DREG [3,17] once the PVREG IR-divergent term 1/2 log 2 (−m 2 A /λ 2 ) is identified in DREG with 1/ǫ 2 . This is quite non-trivial, because it is known that, in general, regularizing the IR divergences via a fictitious gluon mass does not respect the nonabelian symmetry of SU (3). Thus, one expects to get the correct result only for the part proportional to C F . However, we quantize the Lagrangian employing the Background Field Method (BFM) [33], so that the external background gluons satisfy QED-like Ward identities. Then it is not surprising that PVREG gives the correct results also for the C A part. We also remark that within the BFM the renormalization of the strong gauge coupling is due only to the wave function renormalization of the external background gluons. Thus, the renormalization of α s decouples completely from the rest of the calculation, and can be treated separately in the standard way. As a consequence, even if PVREG is used to regularize the loop integrals, the LO partonic cross section σ (0) can be directly expressed in terms of the running coupling α s (µ R ) as in eq. (3).
In the evaluation of the top-stop-gluino contributions to H 2ℓ A , the two-loop integrals are regularized by subtracting from each of them the same expression with m 2 . We conclude this section with a couple of observations concerning the use of PVREG in the computation of the virtual NLO contributions. First, we recall that in PVREG one obtains directly the correct result without the need of introducing a finite renormalization factor to restore the Ward identities. Second, we note that in PVREG the evaluation of the leading term in the Taylor expansion (i.e., the term corresponding to m A = 0) does not require the computation of counterterm diagrams. This seems natural, because the leading term in the one-loop expression, eq. (7), does not depend on the top mass. However, the same evaluation in DREG or DRED does require the computation of counterterm diagrams. Indeed, in n d dimensions the one-loop leading term in the Taylor expansion contains an O(ǫ) part that depends on the top mass, so that the counterterm diagrams give rise to a non-vanishing contribution. Figure 1: Examples of two-loop quark-gluon diagrams (a), and of two-loop quark-squark-gluino diagrams involving (b) the pseudoscalar-quark coupling or (c) the pseudoscalar-squark coupling. Here, q = t, b and i = 1, 2.
4 Two-loop contributions to the form factor H A To fix our notation, we write down the Lagrangian for the interactions of the MSSM pseudoscalar A with quarks and squarks: 2 where: h t and h b are the top and bottom Yukawa couplings A t and A b are the soft SUSY-breaking Higgs-squark-squark couplings; µ is the Higgs mass term in the MSSM superpotential. Our convention for the sign of µ is such that, e.g., the stop and sbottom left-right mixing angles θ t and θ b obey the relations The fact that the pseudoscalar only couples to two different squark mass eigenstates, while gluons only couple to two equal eigenstates, implies that the form factor H A receives neither one-loop contributions from diagrams with squarks nor two-loop contributions from diagrams with squarks and gluons. However, contributions to H 2ℓ A do arise from two-loop diagrams with quarks and gluons, as well as from two-loop diagrams with quarks, squarks and gluinos. Examples of such diagrams, involving either the pseudoscalar-quark coupling or the pseudoscalar-squark coupling, are given in figure 1.
The two-loop form factor for pseudoscalar production can be decomposed as where K 2ℓ qg denotes the quark-gluon contributions (q = t, b), and K 2ℓ qqg denotes the quark-squark-gluino contributions. In the following we discuss separately the two-loop contributions arising from quarkgluon, top-stop-gluino and bottom-sbottom-gluino diagrams.

Quark-gluon contributions
We recall for completeness the results of refs. [3,17] for the contributions to H 2ℓ A arising from diagrams with quarks and gluons (see figure 1a). If the corresponding contribution in the one-loop form factor H 1ℓ A is expressed in terms of the physical quark mass, the two-loop contribution for a given quark q reads If the one-loop form factor is instead expressed in terms of the running quark mass, renormalized in the DR scheme at the scale Q, the two-loop contribution becomes Expressions for the functions denoted here as F 1 (τ ), F 2 (τ ) and F 3 (τ ), valid for arbitrary values of τ , can be found in ref. [17]. They correspond to the functions E (4/τ ) in eq. (4.7), and K (2ℓ,C A ) t (4/τ ) in eq. (4.12) of that paper, respectively. Their limiting behaviors for heavy and light quark are (τ ≪ 1) :

Top-stop-gluino contributions
While a fully analytic computation of the top-stop-gluino contributions to H 2ℓ A valid for arbitrary values of all the relevant particle masses is currently beyond our reach, it is possible to derive approximate analytic results valid in different phenomenologically relevant limits.
To start with, we computed the term K 2ℓ ttg in eq. gluino masses. Such expansion should give a reasonable approximation to the full result when m A is small compared to the other masses, and is anyway restricted to values of m A below the lowest threshold encountered in the diagrams (this usually means m A < 2 m t ). In the limit of vanishing m A we find that our result for K 2ℓ ttg can be cast in an extremely compact form: where the function Φ(m 2 g , m 2 t , m 2 t i ) is given, e.g., in appendix A of ref. [34], and we introduced the shortcut ) . As appears from eqs. (5) and (7), in the limit of vanishing m A the one-loop top contribution to H A reduces to − cot β, i.e., it does not actually depend on any parameter subject to O(α s ) corrections. Therefore, the results in eqs. (27) and (28) do not depend on the renormalization scheme in which the calculation is performed. The contributions to K 2ℓ ttg of the first order in the Taylor expansion in m 2 A are too lengthy to be printed here, but in section 5 we will discuss their relevance in a representative region of the MSSM parameter space.
The two terms between parentheses in eq. (27) come from the diagrams with pseudoscalar-top and pseudoscalar-stop couplings in figures 1b and 1c, respectively. Inserting the explicit expressions for s 2θt and Y t we find Even when the superparticles are much heavier than the pseudoscalar, the validity of the result for K 2ℓ ttg obtained via a Taylor expansion in m 2 A becomes questionable if m A is close to or even larger than m t . To cover this region of the parameter space we performed an asymptotic expansion of K 2ℓ ttg in the large superparticle masses. More specifically, we consider the case (m A , m t ) ≪ M without assuming any hierarchy between m A and m t , and retain terms up to O(m 2 A /M 2 ) and O(m 2 t /M 2 ) in the expansion. Assuming that the top contribution to H 1ℓ A in eqs. (5) and (6) is expressed in terms of the pole top mass, we find where x i = m 2 t i /m 2 g , the one-loop function K 1ℓ (τ ) was defined in eq. (6), and the terms R i collect contributions suppressed by m t /M or m 2 t /M 2 : In the equations above, B denotes the finite part of the Passarino-Veltman function B 0 (m 2 A , m 2 t , m 2 t ) computed at the renormalization scale Q 2 = m 2 t . The comparison between the result for K 2ℓ ttg obtained via a Taylor expansion in m 2 A and the corresponding result obtained via an asymptotic expansion in M will be discussed in section 5. . As a result, assuming that the bottom contribution to H 1ℓ A in eqs. (5) and (8) is fully expressed in terms of the pole bottom mass, we again find a rather compact expression for the term K 2ℓ bbg in eq. (18):

Bottom-sbottom-gluino contributions
Here x i = m 2 b i /m 2 g , and R 1 collects the contributions suppressed by m b /M : As in the case of the top-stop-gluino contribution, the terms proportional to Y b originate from the diagrams that involve the pseudoscalar-sbottom coupling, while the other terms originate from the diagrams that involve the pseudoscalar-bottom coupling. Inserting the expressions for s 2θ b and Y b in the first term in the right-hand side of eq. (35) we obtain Similarly to what found in ref. [23] for the production of CP-even Higgs bosons, if the one-loop contribution to H A is expressed in terms of the pole bottom mass the bottom-sbottom-gluino diagrams induce potentially large two-loop contributions. According to whether or not we insert the explicit expression for s 2θ b in our formulae, such contributions manifest themselves either as terms enhanced by the ratio mg/m b , as in eq. (35), or as terms enhanced by tan β, as in eq. (37). However, such terms cancel out if the pseudoscalar-bottom coupling entering the one-loop contribution to H A is identified with the DR-renormalized mass m b , while the mass of the bottom quark running in the loop is identified with the pole mass M b (this amounts to rescaling by m b /M b the one-loop result fully computed in terms of M b ). As a result, the two-loop form factor in eq. (18) is shifted as with respect to the result obtained when the one-loop bottom contribution is fully expressed in terms of M b . Here Q is the scale at which the running mass m b is renormalized, and (δm b ) SU SY denotes the SUSY contribution to the bottom self-energy, in units of C F α s /π and in the limit of vanishing m b : where While the shift in eq. (38) removes the contributions enhanced by mg/m b (or tan β), it does introduce potentially large logarithms of the ratio between the renormalization scale Q and the masses of the particles running in the loop. Such logarithms cannot be eliminated by a specific scale choice for m b , unless Q is set to a value much smaller than the bottom mass itself. Therefore, as already found in ref. [23] for the CP-even Higgs bosons, the bottom contributions to H 2ℓ A may turn out to be sizable even in the "mixed" renormalization scheme in which the tan β-enhanced contributions are absorbed in a redefinition of the pseudoscalar-bottom coupling entering H 1ℓ A . Finally, if the bottom contribution to H 1ℓ A is fully expressed in terms of the running bottom mass m b the bottom-sbottom-gluino contribution to the form factor in eq. (35) is shifted as In this case H 2ℓ A contains both terms enhanced by mg/m b and potentially large logarithms, the latter arising from (δm b ) SU SY in eq. (41) as well as from the two-loop bottom-gluon contribution in eq. (20).

Comparison with the effective-Lagrangian approximation
It is well known that, in the MSSM, loop diagrams involving superparticles induce interactions between the quarks and the "wrong" Higgs doublets, i.e., interactions that are absent from the tree-level Lagrangian due to the requirement that the superpotential be a holomorphic function of the superfields [36]. Such non-holomorphic, loop-induced Higgs-quark interactions result in tan β-enhanced (or tan β-suppressed) corrections to the MSSM predictions for various physical observables. If all superparticles are considerably heavier than the Higgs bosons they can be integrated out of the Lagrangian, in which case the loop-induced corrections are resummed in effective Higgs-quark couplings. In particular, if g φ b denote the tree-level couplings of a neutral Higgs φ = (h, H, A) to bottom quarks (normalized to the SM value), the corresponding effective couplingsg φ b read [29,30] where α is the mixing angle in the CP-even Higgs sector and, to O(α s ), In the calculation of processes involving the Higgs-bottom couplings, it is often found that the tan β-enhanced corrections can be included to all orders in an expansion in powers of α s tan β by inserting the effective couplings of eq. (42) in the lowest-order result. A comparison with our explicit results for the two-loop form factors allows us to test the validity of that procedure in the case of the production of both CP-even [23] and CP-odd Higgs bosons in gluon fusion. 3 We recall that the bottom-quark contributions H 1ℓ ,b φ to the one-loop form factors for the production of the Higgs boson φ = (h, H, A) read where the function G 1ℓ 1/2 (τ ) is given, e.g., in eq. (12) of ref. [23]. Assuming that H 1ℓ ,b φ are expressed in terms of the pole bottom mass, and that the Higgs-sbottom couplings are renormalized in a way that avoids the introduction of additional tan β-enhanced corrections (see ref. [23]), we find that the two-loop form factors read where the ellipses denote contributions suppressed by m b /M or m 2 Z /M 2 , as well as all of the contributions from diagrams involving top and stop, and In practice, the effective-Lagrangian approximation consists in rescaling the one-loop bottom contributions H 1ℓ ,b φ by the same factors that rescale the Higgs-bottom couplings g φ b in eq. (42). Expanding the rescaling factors to the first order in ∆ b it is easy to see that the effective-Lagrangian approximation does indeed reproduce the two-loop terms proportional to ∆ b in eqs. (45)-(47).
It is also interesting to consider the so-called decoupling limit of the MSSM, m A ≫ m Z , in which cot α → − tan β and the light scalar h has SM-like couplings to fermions and gauge bosons. 4 Eq. (42) shows that in this limit the effective coupling of h to bottom quarks is equal to the tree-level coupling, therefore in the effective-Lagrangian approximation there are no tan β-enhanced contributions to H 2ℓ h . Indeed, for cot α → − tan β the terms proportional to ∆ b drop out of the two-loop form factor in eq. (45). However, eq. (45) also shows that in the decoupling limit H 2ℓ h contains additional tan βenhanced contributions, controlled by the left-right sbottom mixing X b = (A b + µ tan β), which are not reproduced by the effective-Lagrangian approximation. However, when the implicit dependence of the sbottom masses and mixing on the bottom mass is taken into account, such contributions turn out to be partially suppressed by powers of m b . Indeed, taking for illustrative purposes the limit in which the diagonal entries of the sbottom mass matrix as well as the squared gluino mass are all equal to M 2 , and expanding the form factor in powers of m b , we find where the ellipses denote terms further suppressed by powers of m b or m Z , as well as all of the contributions from diagrams involving top and stop. The first term in eq. (49) comes from the expansion of the terms proportional to s 2 2θ b in eq. (45), while the second comes from the expansion of terms not shown in eq. (45). The contributions neglected by the effective-Lagrangian approximation can be relevant for values of X b large enough to compensate for the suppression due to m b . It should however be recalled that in the decoupling limit H 1ℓ ,b h is not further enhanced by tan β, thereforedifferently from what happens in the case of the heavy Higgs bosons -the total form factor for h production can still be dominated by the top/stop contributions even for large values of tan β.

Numerical examples
We will now illustrate the effect of the two-loop quark-squark-gluino contributions to the form factor for pseudoscalar Higgs production in a representative region of the MSSM parameter space.
The SM parameters entering our calculation include the Z boson mass m Z = 91.1876 GeV, the W boson mass m W = 80.399 GeV and the strong coupling constant α s (m Z ) = 0.118 [37]. Since the squarks do not contribute to the one-loop amplitude for pseudoscalar production, the only parameters entering H 1ℓ A in addition to the quark masses are tan β and m A . Neither of those parameters is subject to one-loop O(α s ) corrections, therefore we need not specify a renormalization scheme for them (although it is natural to consider m A as the pole pseudoscalar mass). The remaining input parameters are mg, µ, A t , A b and the soft SUSY-breaking mass terms for stop and sbottom squarks, m Q , m U and m D . Since these parameters only enter the two-loop part of the form factor we need not specify a renormalization scheme for them either. For simplicity, in our numerical examples we will set all the SUSY-breaking parameters, as well as the supersymmetric mass parameter µ, to a common value M . Note however that the squark mass eigenstates will differ from M , because of the supersymmetric (F-term and D-term) contributions to the squark mass matrices as well as of the left-right mixing terms.
In figure 2 we show the top-stop-gluino contribution to the two-loop form factor for pseudoscalar production, i.e., the term K 2ℓ ttg entering eq. (18), as a function of the common SUSY mass M , for m A = 150 GeV and tan β = 2. Even for the lowest value of M considered in the plot, M = 100 GeV, the stop and sbottom masses are above the threshold for real-particle production. The dashed line represents the result obtained in the limit of vanishing m A , shown explicitly in eqs. (27) and (28) It can be seen in figure 2 that the two-loop top-stop-gluino contribution K 2ℓ ttg is of non-decoupling nature, i.e., it does not tend to zero when all the superparticle masses become large (note that the superpotential parameter µ increases together with the SUSY-breaking parameters). In addition, the comparison between the solid and dashed lines shows that when the common SUSY mass M is close to m A the combined effect of the terms of O(m 2 A /m 2 t ) and O(m 2 A /M 2 ) can be as large as 20%-25% with respect to the result obtained for vanishing m A . However, when M increases the effect of the terms of O(m 2 A /M 2 ) becomes quickly negligible. The remaining discrepancy between the solid and dashed lines for moderate to large values of M is due to the terms of O(m 2 A /m 2 t ), and it amounts to a modest 6% for the value of m A considered in this example.
To assess the importance of the terms of O(m 2 A /m 2 t ) for larger values of m A , we plot in figure 3 the real part of K 2ℓ ttg as a function of the pseudoscalar mass, up to a value m A = 500 GeV well above the threshold for real top-quark production. The common SUSY mass is set to the relatively large value M = 1 TeV, and tan β = 2. As in figure 3, the dashed and solid lines represent the results obtained at the zeroth and first order of the Taylor expansion in m 2 A , respectively. The comparison between those lines shows that when m A approaches 2 m t the effect of the terms of O(m 2 A /m 2 t ) gets as large as 30% with respect to the result obtained for vanishing m A . However, it is natural to wonder whether a Taylor expansion in m 2 A can give an accurate approximation to K 2ℓ ttg for values of m A close to or larger than m t . To address this question, we show in figure 3 as a dot-dashed line the result of the asymptotic expansion in M , given explicitly in eqs. (30)- (34). This result was derived under the assumption that both m A and m t are much smaller than M , which is indeed the case for M = 1 TeV, but it does not require any specific hierarchy between m A and m t . The comparison between the dot-dashed and solid lines shows that the Taylor expansion at the first order in m 2 A provides a good description of the dependence of K 2ℓ ttg on the ratio m A /m t up to values of m A of the order of 250 GeV. On the other hand, when m A reaches the threshold for real top production (i.e., at the cusp of the dot-dashed line) the result of the asymptotic expansion in M is roughly 80% larger in absolute value than the result at the first order of the Taylor expansion in m 2 A , and a full 140% larger than the result obtained for vanishing m A .
In summary, it appears that the compact result for K 2ℓ ttg given in eqs. (27) and (28), which was derived for m A = 0, can be safely applied only to scenarios in which m A is smaller than m t . While the inclusion of the terms proportional to m 2 A pushes the validity of the Taylor expansion up to larger values of m A , the expansion fails when m A gets close to the threshold for real top production. In that case one can use the result of the asymptotic expansion in M , provided that the latter is still considerably larger than m A .
We are now ready to discuss the relative importance of the various two-loop contributions to the form factor for pseudoscalar production. We will see that, at least in the region of the parameter space that we consider in this example, the results are qualitatively similar to what we found in ref. [23] for the case of the heavy scalar H.
A precise NLO determination of the cross section for pseudoscalar production would require us to take into account the contribution of one-loop diagrams with real parton emission, and to perform an integration over the phase space (see section 2). However, for the purpose of illustrating the relative importance of the various two-loop contributions, we can just define a factor K A that contains the ratio of two-loop to one-loop form factors appearing in eq. (11): In the left panel of figure 4 we plot K A as a function of tan β, for m A = 150 GeV and all SUSY mass parameters equal to M = 500 GeV. The one-loop form factor H 1ℓ A in eq. (50) contains both the top and bottom contributions, computed under the approximations of eqs. (7) and (8)   dominate the two-loop form factor up to values of tan β around 5. For larger values of tan β the contribution of the bottom-sbottom-gluino diagrams (included in the solid line) becomes the dominant one, and K A grows linearly with tan β. This behavior can be understood by recalling that, as can be seen in eq. (16), the Yukawa coupling of the pseudoscalar to bottom quarks is enhanced by tan β with respect to the coupling of the SM Higgs, while the coupling to top quarks is suppressed by tan β. Consequently, for moderate to large values of tan β both the one-loop and the two-loop form factors in K A are dominated by the contribution of the diagrams controlled by the pseudoscalar-bottom coupling, with the result that the coupling itself cancels out in the ratio. However, the dominant contribution from the bottom-sbottom-gluino diagrams in the OS scheme, see eq. (37), contains an additional tan β-enhancement, which explains the linear rise of K A . On the other hand, the proximity between the dotted and dashed lines shows that, in the OS scheme, the contribution to H 2ℓ A of the two-loop diagrams with bottom quarks and gluons is very small. This is due to a partial cancellation among the three terms entering K 2ℓ bg in eq. (19), and to the fact that, in this scheme, the term F 2 (τ b ) is not enhanced by the potentially large logarithm of the ratio between the bottom mass and the renormalization scale, as can be seen by comparing eqs. (19) and (20).
As discussed in section 4.3, all tan β-enhanced terms cancel out in a "mixed" renormalization scheme in which the pseudoscalar-bottom Yukawa coupling in the one-loop part of the result is identified with the DR-renormalized MSSM bottom mass m b (Q), where Q is a reference scale that we take equal to m A , while the mass of the bottom quark running in the loop is identified with the pole mass M b . To determine m b (m A ), we first evolve the MS-renormalized SM mass m b (m b ) up to the scale m A via the NLO-QCD renormalization group equations, then we convert it to the DR-renormalized SM mass m SM b (m A ) via the appropriate shift, and finally we convert it to the MSSM running mass according to where ∆ b is given in eq. (43), and δ b is proportional to the part of (δm b ) SU SY in eq. (39) that is not enhanced by tan β: The "mixed" renormalization prescription is realized by computing the one-loop bottom contribution K 1ℓ (τ b ) in eq. (5) in terms of the pole mass M b , then rescaling it by a factor m b (m A )/M b . The two-loop form factor H 2ℓ A must then be shifted as in eq. (38). In the right panel of fig. 4 we present the result of this manipulation. The input parameters and the meaning of the different lines are the same as for the plot in the left panel. The proximity between the dashed and solid lines, and the flatness of the lines for moderate to large values of tan β, show that the contribution of the two-loop bottomsbottom-gluino diagrams is rather small in this renormalization scheme, and it does not induce an additional tan β-enhancement. However, the comparison between the dotted and dashed lines shows that there is a sizable contribution to K A from the two-loop diagrams involving bottom quarks and gluons. This is due to the fact that the shift in eq. (38) brings back a large logarithm, ln(m 2 b /m 2 A ), which compensates the scale dependence of the running mass m b .

Conclusions
The calculation of the production cross section for the MSSM Higgs bosons is not quite as advanced as in the SM. Indeed, despite valiant efforts [21,22], a full computation of the two-loop quark-squarkgluino contributions, valid for arbitrary values of all the relevant particle masses, has not been made publicly available so far. Approximate analytic results, however, can be derived if the Higgs bosons are somewhat lighter than the squarks and the gluinos. In the MSSM this condition almost certainly applies to the lightest scalar h. Moreover, recent results from SUSY searches at the LHC [40] set preliminary lower bounds on the squark and gluino masses just below the TeV (albeit for specific models of SUSY breaking), suggesting that there might be wide regions of the MSSM parameter space in which the condition also applies to the heavy scalar H and to the pseudoscalar A.
In this paper we presented a calculation of the two-loop quark-squark-gluino contributions to the cross section for pseudoscalar production. We exploited techniques developed in our earlier computations of the production cross section for the CP-even Higgs bosons of the MSSM [15,23] to obtain explicit and compact analytic results based on expansions in the heavy particle masses. We avoided problems related to the definition of the Dirac matrix γ 5 in n d = 4 dimensions, which are specific to the case of pseudoscalar production, by regularizing the loop integrals with the Pauli-Villars method. For what concerns the top-stop-gluino contributions, we provided both the result of a Taylor expansion in the pseudoscalar mass, up to and including terms of O(m 2 A /m 2 t ) and O(m 2 A /M 2 ), and the result of an asymptotic expansion in the superparticle masses, up to and including terms of O(m 2 A /M 2 ) and O(m 2 t /M 2 ). The latter can be easily adapted to the case of the bottom-sbottom-gluino contributions, providing a result valid up to and including terms of O(m 2 b /m 2 A ) and O(m b /M ). We discussed how the tan β-enhanced terms in the bottom-sbottom-gluino contributions can be eliminated via an appropriate choice of renormalization scheme for the parameters entering the one-loop part of the calculation, and compared our results with those obtained in the effective-Lagrangian approximation. All of our results can be easily implemented in computer codes for an efficient and accurate determination of the cross section for pseudoscalar production.
Finally, the results derived in this paper for the production cross section can be straightforwardly adapted to the NLO computation of the gluonic and photonic decay widths of the pseudoscalar Higgs boson in the MSSM, in analogy to what described in section 5 of ref. [15] for the case of the CP-even bosons.