Pseudoscalar MSSM Higgs Production at NLO SUSY--QCD

One of the most important mechanisms at the Large Hadron Collider (LHC) for the production of the pseudoscalar Higgs boson of the Minimal Supersymmetric Standard Model (MSSM) is the loop-induced gluon fusion process $gg\to A$. The higher-order QCD corrections have been obtained a long time ago and turned out to be large. However, the genuine supersymmetric (SUSY--)QCD corrections have been obtained only in the limit of large SUSY particle masses so far. We describe our calculation of the next-to-leading-order (NLO) SUSY--QCD results with full mass dependence and present numerical results for a few representative benchmark points. We also address the treatment of the effective top and bottom Yukawa couplings, in the case of heavy SUSY particles, in terms of effective low-energy theories where the heavy degrees of freedom have been decoupled. Furthermore, we include a discussion of the relation between the SUSY--QCD corrections that we have computed and the Adler--Bardeen theorem for the axial anomaly. In addition, we apply our results to the gluonic and photonic pseudoscalar Higgs decays $A\to gg,\gamma\gamma$ at NLO.


Introduction
The discovery of a Standard-Model-like Higgs boson at the LHC [1] completed the Standard Model (SM) of electroweak and strong interactions. The existence of the Higgs boson [2] is inherently related to the mechanism of spontaneous symmetry breaking while preserving the full gauge symmetry and the renormalizability of the SM [3]. The measured Higgs boson mass of (125.09 ± 0.24) GeV [4] ranks at the weak scale. The existence of the Higgs boson allows the SM particles to be weakly interacting up to high-energy scales [5]. This, however, is only possible for particular Higgs-boson couplings to all other particles, so that the knowledge of the Higgs-boson mass fixes all its properties uniquely. The massive gauge bosons and fermions acquire mass through their interaction with the Higgs field that develops a vacuum ‡ Former academic affiliation 1 expectation value in its ground state. The minimal model requires the introduction of one isospin Higgs doublet and leads after spontaneous symmetry breaking to the existence of one scalar Higgs boson. The SM itself, however, leaves several fundamental questions open as e.g. the nature of Dark Matter, the baryon asymmetry of the universe or the stability of the electroweak against the Planck or grand unification scale. If the SM is extended to a Grand Unified Theory (GUT) scale, radiative corrections to the Higgs-boson mass tend to push it towards the GUT scale, if the Higgs boson couples to particles of that mass order. In order to obtain a Higgs mass at the electroweak scale the Higgs-mass counterterm has to be fine-tuned to cancel these large corrections thus establishing an unnatural situation that asks for a solution. This is known as the hierarchy problem [6]. These open questions call for extensions of the minimal model. To increase the experimental sensitivity to effects beyond the SM (BSM), the SM and BSM parts of measured relevant observables need to be known as precisely as possible in order to allow for a reliable interpretation of potential deviations and effects beyond the SM.
The open problems of the SM motivate extensions of the minimal model which cover e.g. the Two-Higgs-Doublet model (2HDM) [7] or the minimal supersymmetric extension (MSSM) [8,9] as prominent and highly motivated examples. Supersymmetric extensions of the SM provide a solution to the hierarchy problem if the supersymmetric particle masses rank at scales up to a few TeV [10]. Supersymmetry relates fermionic and bosonic degrees of freedom and thus links internal and external symmetries. The MSSM, if embedded in a Grand Unified Theory, predicts a value of the Weinberg angle in excellent agreement with experimental measurements of electroweak precision observables [11]. Moreover, it contains a Dark Matter candidate if R-parity is conserved [12] and allows for generating electroweak symmetry breaking radiatively, since the top mass ranks in the proper region for that mechanism to work [13]. The MSSM introduces two isospin Higgs doublets due to the analyticity of the superpotential, requiring two different doublets for the generation of the up-and down-type fermion masses and the anomaly-freedom with respect to the gauge symmetries [14], since the higgsino states as the supersymmetric partners of the Higgs bosons contribute to the Adler-Bell-Jackiw anomaly [15]. Due to this, the MSSM Higgs sector is a 2HDM of type II at leading order (LO). There are a light (h) and heavy (H) scalar, a pseudoscalar (A) and two charged (H ± ) states as the corresponding mass eigenstates. Since the self-interactions of the Higgs fields, as defined by the corresponding Higgs potential, are entirely fixed by the electroweak gauge couplings, this induces an upper bound on the light scalar Higgs mass that has to be smaller than the Z-boson mass M Z at LO. However, radiative corrections, which are dominated by top-quark-induced contributions, strongly increase this upper bound to about 130 GeV in general [16]. The Higgs sector is uniquely fixed at LO by the value of the pseudoscalar mass M A and the parameter tgβ, defined as the ratio of the two vacuum expectation values of the scalar Higgs fields.
In this work, we will describe the calculation of the full SUSY-QCD corrections at NLO to pseudoscalar Higgs production via the gluon-fusion mechanism gg → A. This process belongs to the dominant MSSM Higgs-boson production processes at the LHC and thus contributes to the present bounds on the so far negative searches for the heavy MSSM Higgs bosons at the LHC. In order to make the predictions for this process reliable, the full NLO corrections within SUSY-QCD have to be computed. The paper is organized as follows. In Section 2, we will summarize the present status of the gluon-fusion cross section. In Section 3, we briefly discuss pseudoscalar Higgs decays to gluons and photons. In Section 4 we will describe our implementation of the stop and sbottom sector followed by the detailed description of our NLO calculation in Section 5. In the latter we also include a discussion of effective Yukawa couplings and the relation of the considered process to the Adler-Bardeen theorem [17]. In Section 6, we discuss numerical results for a few representative benchmark points. We close the paper with our conclusions in Section 7.

Gluon Fusion
The dominant channels for pseudoscalar production at a hadron collider are given by gluon fusion, gg → A, and production in association with bottom quarks, qq, gg → Abb, with their relative importance depending on the value of tgβ. For large tgβ, Abb production dominates, with the gluon fusion contribution amounting to up to about 30% close to the present exclusion bounds, depending on the region in the M A − tgβ plane [18,19].

Leading Order
The gluon-fusion mechanism [20] pp → gg → A dominates the pseudoscalar MSSM Higgs boson production at the LHC in the phenomenologically relevant Higgs mass ranges for small and moderate values of tgβ. Only for large tgβ the associated Abb production channel develops a larger cross section due to the enhanced Higgs couplings to bottom quarks [21]. The gluon coupling to pseudoscalar Higgs bosons in the MSSM is built up by loops involving top and bottom quarks, see Fig. 1. The partonic A t, b g g Figure 1: Typical diagram contributing to gg → A at lowest order.
cross section is given at lowest order by [22,23]: where G F denotes the Fermi constant, α s the strong coupling, and µ R the renormalization scale. The scaling variables are defined as , andŝ denotes the partonic c.m. energy squared. The amplitudes A A Q (τ Q ) are obtained as and the MSSM coupling factors g A Q are determined as g A t = 1/tgβ, g A b = tgβ. In the narrowwidth approximation the hadronic cross section is given by with the scaling variable τ A = M 2 A /s, where s specifies the total hadronic c.m. energy squared, and the gluon luminosity at the factorization scale µ F . For small tgβ the top-loop contribution is dominant, while for large values of tgβ the bottom-quark contribution is strongly enhanced.

QCD Corrections
The full two-loop QCD corrections to the gluon-fusion cross section were calculated in the past [23,24,25]. In complete analogy to the SM case, they consist of virtual two-loop corrections to the basic gg → A process and real one-loop corrections due to the associated production of the pseudoscalar Higgs boson with massless quarks and gluons. The final result for the hadronic cross section at NLO can be decomposed as The analytical expressions for arbitrary Higgs boson and quark masses at NLO are rather involved [23,25]. As in the SM case, the quark-loop masses have been identified with the pole mass m Q (Q = t, b), while the QCD coupling and the parton distribution functions (PDFs) of the proton are treated in the MS scheme with five active flavours. The axial γ 5 coupling can be regularized in the 't Hooft-Veltman scheme [26] or its extension by Larin [27], which preserve the chiral symmetry in the massless quark limit by the addition of supplementary counterterms and fulfill the non-renormalization theorem [17] of the ABJ anomaly [15] at vanishing momentum transfer. The same result can also be obtained with the scheme of Ref. [28] that gives up the cyclicity of the traces involving Clifford matrices. The next-tonext-to-leading order (NNLO) QCD corrections have been obtained in the limit of heavy top quarks (HTL) [29]. The QCD corrections are positive and large in total, increasing the MSSM Higgs production cross sections at the LHC by up to about 100%. For the top-loop contributions alone, the (moderate) NNLO corrections in the heavy-top limit (HTL) can be used consistently. Electroweak corrections are unknown so far. The leading terms of the relative QCD corrections in the HTL provide a reasonable approximation for small tgβ up to pseudoscalar Higgs masses of ∼ 1 TeV with a maximal deviation of ∼ 25% for tgβ < ∼ 5 at NLO in the intermediate mass range [30]. The genuine SUSY-QCD corrections are only known in the limit of heavy SUSY particles [31,32]. For large values of tgβ they can be large and approximated by the ∆ b terms. This work improves this incomplete status by calculating the full SUSY-QCD corrections with full virtual quark-, squark-and gluino-mass dependence, which will contribute to the virtual corrections as where C A QCD is the virtual part of the pure QCD corrections. We will compare the full results for C A SQCD with the approximate calculations in the following sections. For the SUSY-QCD corrections we implement the stop and sbottom sector at the NLO level, although the squarks do not contribute at LO, and therefore the definition of a renormalization scheme for their parameters is not required. However, to be in line with the treatment of scalar Higgs production in a future work, where stops and sbottoms contribute at LO already, we choose the same framework. The NLO implementation of the stop and sbottom sectors will be discussed in Section 4.
In the opposite limit, where the pseudoscalar Higgs mass is much larger than the quark mass, the analytical results of the relative QCD corrections coincide with the SM expressions at the leading and subleading logarithmic level for both the scalar and pseudoscalar Higgs bosons up to NLO where the results for small quark masses are known [23]. This coincidence is due to the restoration of the chiral symmetry in the massless quark limit. The leading double and subleading logarithms have been resummed recently [33].

Pseudoscalar Higgs Decays
Although pseudoscalar Higgs decays into gluons and photons do not play a prominent role as for the SM-like light scalar Higgs particle, they can still reach sizeable branching ratios for smaller values of tgβ so that they might be accessible at future e + e − colliders.

A → gg
The decay of pseudoscalar Higgs bosons into gluons is loop-induced, see Fig. 2. The dominant contributions originate from top and bottom loops, while lighter quarks as e.g. the charm quark yield contributions at the per-cent or sub-per-cent level only. The LO expression of the gluonic pseudoscalar Higgs decay reads [22,23] A t, b g g Figure 2: Typical diagrams contributing to A → gg at lowest order.
where we adopted the same notation as in Eq. (1) using the same quark form factors as given in Eq. (2). The NLO QCD and SUSY-QCD corrections can be cast into the form with the NLO coefficient E A splitting into pure QCD corrections and genuine SUSY-QCD corrections, The QCD part can be expressed as [23,34] where ∆ m denotes finite mass effects at NLO [23], and N F is the number of active light flavors included as final-state quarks as well. For e.g. tgβ = 1 the mass effects amount to ∆ m ≈ 1.3, if the quark masses are defined as pole masses, but are larger for increasing values of tgβ due to the rising significance of the bottom contributions. The expression without ∆ m corresponds to the heavy-quark limit of the relative QCD corrections. The coefficient E A SQCD coincides with the one for the gluon-fusion cross section of Eq. (6),

A → γγ
As for the gluonic pseudoscalar Higgs decay, its decay into photon pairs is a loop-induced process with top and bottom quarks providing the dominant contributions, but also charginos, see Fig. 3. At LO, the pseudoscalar decay width into photon pairs reads [22,23] whereχ ± denotes the two chargino mass eigenstates, N cf is the color factor of the fermions of charge e f contributing to the loops. The LO form factors A A i (τ i ) (i = t, b,χ ± ) follow the with the charge factors Q ii , S ii (i = 1, 2) given in Refs. [9,22]. They are related to the mixing angles between the chargino statesχ ± 1,2 . The NLO QCD and SUSY-QCD corrections can be defined as a shift of the corresponding LO quark-form factors, where the pure QCD corrections D Q,QCD to the quark form factor vanish in the heavy-quark limit due to the Adler-Bardeen [17] theorem for these leading contributions. This means they are induced by pure quark-mass effects [23]. The implementation of the QCD corrections D A Q,QCD follows Ref. [23], i.e. the running quark masseŝ where m Q (µ) denotes the MS mass [35] and K b = 12.4, K t = 10.9 [36], are used for the loopquark masses at the scale µ = M A /2 such that the relations M A = 2m Q (m Q ) = 2m Q (Q = t, b) define the virtual quark thresholds in terms of the quark pole masses m Q . The genuine SUSY-QCD corrections, represented by the coefficient D A Q,SQCD , will be discussed in Section 6.3.

4 Squark Masses and Couplings
In the following the parametrization of the stop and sbottom sectors will be described in detail at LO and at NLO starting from the soft SUSY-breaking parameters, where the extension to NLO requires a dedicated scheme choice for our gluon-fusion calculation. We will follow the set-ups described in Refs. [37,38] with corresponding modifications.

Sfermion Masses and Couplings at LO
Since the scalar sfermion current-eigenstatesf L,R , the super-partners of the the left-and right-handed fermions, mix with each other, the corresponding mass eigenstatesf 1,2 are related to the current eigenstates by a rotation involving the mixing angles θ f , These mixing angles grow with the Yukawa couplings of the corresponding SM fermions, i.e. mixing effects are in general only relevant for the third-generation sfermionst,b,τ . The mass matrix in the current-eigenstate basis is given by where r b = r τ = 1/r t = tgβ. A f is the trilinear sfermion coupling of the soft SUSY-breaking part of the Lagrangian, while µ denotes the higgsino mass parameter and m f the fermion mass. The parametersMf L/R absorb the corresponding D-terms, with e f being the electric charge of the sfermion, and I 3L its third isospin component, θ W denotes the Weinberg angle and Mf L/R are the sfermion mass parameters of the soft SUSYbreaking part of the Lagrangian. Hence, the mixing angles are determined from and the squark-eigenstate masses acquire the form 8 In the current-eigenstate basis, the neutral Higgs couplings to sfermions are given by where the couplings g Φ i (i = 1, . . . , 4) are specified in Table Next, we will discuss the extension of the stop and sbottom sectors to the NLO SUSY-QCD level.f

Stops and Sbottoms at NLO
At NLO, we will introduce the soft SUSY-breaking parameters in the MS scheme, i.e. we will start from the soft supersymmetry-breaking parameters MQ L,R (Q 0 ) and A Q (Q 0 ) at the input scale Q 0 which will in general be the SUSY scale, i.e. the average size of the left-and right-handed soft SUSY-breaking mass parameters. The benchmark scenarios of Ref. [19], however, are defined in the on-shell scheme of all involved input parameters. Thus, we will describe how we are implementing the relation between the MS and the on-shell parameters. The bottom and top masses involved in the sbottom and stop mass matrices have to be chosen such that large higher-order corrections to their entries are avoided. We have chosen the top pole mass and a derived bottom mass for the sbottom mass matrix according to Refs. [38]. At LO, the stop/sbottom mass matrix is then given by (q = t, b) where m t is the top pole mass and m b is the derived bottom mass as will be discussed in the following. The D-terms DQ L/R have again been absorbed in the soft SUSY-breaking parameters, The diagonal and off-diagonal entries of the stop/sbottom mass matrix are corrected at higher orders. We absorb the radiative corrections to the diagonal matrix elements in shifted soft mass parameters, MQ L/R , while the corrections to the off-diagonal entries will be absorbed in shifted soft trilinear couplings, The shifted parameters are related to the radiative corrections to the mixing angles and stop/sbottom masses in order to arrive at simple tree-level like expressions at NLO for the stop/sbottom parameters. On the other hand, these shifted parameters correspond to the on-shell scheme introduced in Refs. [38] and thus have to coincide with the input values of the chosen benchmark scenario.

Stops
Starting from the on-shell parameters the treatment of the stop sector is identical to the LO level discussed before. The relation of the on-shell to the MS parameters, however, is affected by the NLO corrections. At tree-level, the mixing angleθ Q is derived from where the tree-level squark masses mq 1/2 according to Eq. (20) have been used 1 .
1 The standard range for the squark mixing angle is chosen between 0 and π. The masses of the stop/sbottom mass eigenstates acquire radiative corrections, The self-energies Σ 11/22 of the stops/sbottoms can be derived from the diagrams in Fig. 4, where Mg denotes the gluino mass and the scalar one-loop integrals are defined as (n = 4−2ǫ) The scaleμ denotes the 't Hooft mass of dimensional regularization. The mass counterterms δm 2Q 1,2 of Eq. (28) are related to the counterterms of the input parameters, using the tree-level mixing angleθ Q of Eq. (27), and C F = 4/3. The counterterms of the parameters M 2Q L/R (Q 0 ) and A Q (Q 0 ) are defined in the MS scheme, The counterterm of the pole quark mass m q is given by where δ SU SY = 1/3 is a finite counterterm required to restore the supersymmetric relation between the Higgs-boson couplings to quarks and squarks within dimensional regularization [40]. The definition of the mixing angleθ Q in Eq. (27) corresponds to the following counterterm at NLO, However, this mixing angle definition induces artificial singularities in physical observables for stop/sbottom masses mq 1,2 close to each other [41]. To avoid such singularities, the mixing angle of the squark fields has been renormalized via the anti-Hermitian (on-shell) counterterm [41], with the off-diagonal part Σ 12 of the stop/sbottom self-energy (see Fig. 4) describing transitions from the first to the second mass eigenstate or vice versa, that will be absorbed in the shifted A Q value of Eq. (26). This shift defines the relation between the on-shell coupling A Q and the MS one A Q (Q 0 ).
Using the NLO corrected squark pole masses of Eq. (28) and the radiatively corrected mixing angle θ q , the shifted (on-shell) squared soft SUSY-breaking squark mass parameters L/R can be obtained from the sum rules, while the shifted (on-shell) trilinear couplings A Q are derived from the relation In terms of these shifted (on-shell) parameters the radiatively corrected squark masses and mixing angles are given by LO-like expressions, The scale of the strong coupling constants α s in Eqs. (29,31,32,33,36) has been identified with the input scale Q 0 .
These relations have been used for the determination of the MS parametersM 2Q L/R (Q 0 ) and A q (Q 0 ) iteratively until the on-shell parameters agreed with the input value of the chosen benchmark scenario.

Sbottoms
The procedure described for the stops is necessary to obtain the MS parameter Mt L (Q 0 ) that by virtue of the SU(2) gauge symmetry is identified with the MS parameter Mb L (Q 0 ), Due to potentially large tgβ-enhanced contributions in the sbottom sector the procedure has to be modified. This modification addresses the treatment and renormalization of the bottom mass m b and of the trilinear coupling A b . Therefore, the bottom mass is not introduced as the pole mass, but as a derived quantity, since it represents the contribution of the bottom Yukawa coupling to the sbottom sector. To achieve a working scheme, we are starting from Eq. (40) for the mixing angle that at NLO is still defined via the anti-Hermitian counterterm of Eq. (35). The trilinear coupling A b , however, is now defined from the proper Ab 1b2 vertex [38]. This definition avoids large tgβ-enhanced contributions in the renormalization 13 of A b . The bottom mass m b entering the sbottom mixing matrix is then treated as a derived quantity. This leads to the explicit counterterms, where the term F is defined as [42] The derived bottom massm b is then determined aŝ where m b (Q 0 ) denotes the MS bottom mass at the input scale Q 0 , δm b the counterterm of Eq. (42) and δm b the MS counterterm of the bottom mass, where δ SU SY = 1/3 is a SUSY-restoring counterterm. This MS counterterm defines the running bottom mass with decoupled SUSY contributions, i.e. the running bottom mass of the SM. The derived bottom massm b is then used for the sbottom mixing matrix throughout.
In the analogous way we determine the MS value A b (Q 0 ) of the trilinear coupling, but this will not be used in our analysis. The shifted (on-shell) sbottom mass parametersMb L/R are finally determined from the corresponding sum rules of Eq. (38). This set-up of the sbottom sector is then used for iteration until the on-shell parameterMb R agrees with the input parameter of the benchmark scenario.
An alternative approach is provided by a purely fixed-order implementation of the difference between Mb L and Mt L , with q = t, b. The counterterms δmq 1/2 are given in Eq. (34), the counterterm δθ q in Eq. (35) and the counterterm δm q in Eq. (33) for the top pole mass m q = m t and in Eq. (42) for the (derived) bottom mass m q =m b . This approach does not require any iteration, since the on-shell parameters of the benchmark scenario can immediately be used to derive the parameters of the sbottom sector. We have compared both approaches and found agreement of the sbottom parameters at the few-per-mille level.

Higgs Couplings to Stops and Sbottoms
The NLO neutral Higgs couplings to squarks in the current-eigenstate basis are given by with the on-shell trilinear couplings A Q and the couplings g Φ i of Table 1. The quark mass m Q denotes either the top pole mass in the stop case or the derived bottom massm b for the sbottom sector. The related couplings to the stop/sbottom mass eigenstatesQ 1,2 are derived by the rotations according to Eq. (22) by the radiatively corrected mixing angle θ Q . For pseudoscalar Higgs bosons, we obtain vanishing diagonal couplings g Ã at the NLO level as at LO.

SUSY-QCD corrections at NLO
The genuine SUSY-QCD corrections at NLO are determined by the Feynman diagrams shown in Fig. 5 that displays only the non-vanishing graphs. Additional permutations of the external gluons have to be added. The matrix element for the LO expression and the SUSY-QCD corrections can be parametrized as where q 1 , q 2 denote the two incoming momenta of the gluons and ǫ µ (q i ) their polarization vectors, δ ab the Kronecker symbol of the adjoint SU(3) c color space and ǫ µναβ the fourdimensional Levi-Civita tensor. In this paper we will mainly use the γ 5 prescription of Larin [27], where the product of Levi-Civita tensors is replaced by the determinant of ndimensional metric tensors in n = 4−2ǫ dimensions. Within this framework we can construct a projector on the anticipated form factors A A LO/SQCD , so that In order to set up a simple notation in close connection to the QCD corrections of Eq. (5) we will normalize the genuine SUSY-QCD corrections to the individual form factors at LO, where C A Q,SQCD depends on all ratios of the pseudoscalar Higgs, quark, squark and gluino masses. The LO form factor has been defined in terms of the expressions of Eq. (2). In the following we will describe the technical details for the numerical integration to determine the complex coefficient C Q,SQCD by exemplifying our method for the first diagram of Fig. 5. In order to regularize virtual thresholds we have added a small imaginary part to the quark and squark masses 2 , with a positive regulatorǭ > 0, which defines the analytical continuation of our two-loop amplitudes. We work with a small but finite value ofǭ that is small enough to achieve results in the narrow-width approximation. For the parametrization of the two-loop diagrams, we follow the same procedure used and described in Refs. [43] for Higgs-boson pair production and adopted in earlier works [44].

Feynman Parametrization
The parametrization of the first two-loop diagram of Fig. 5 reads where we sum over l, m ∈ {1, 2} in T µν 1 , k, q are the loop momenta that are integrated over and the chiral coupling factors IP j (j = 1, 2) are defined as and IP j emerges from IP j by the replacement γ 5 → −γ 5 . After applying the contraction with the projector P µν onto the contribution to the virtual form factor, we introduce Feynman parameters x 3 , x 4 , x 1 , x 2 for the second to fifth propagator (in this ordering) and 1 − j x j for the first one, (k 2 − m 2 Q ). With the substitutions we obtain a four-dimensional Feynman-parameter integral over x, y, z, v with integration boundaries from 0 to 1. The shift in both the numerator and denominator symmetrizes the k-integration that is performed in a simple and systematic way for the emerging integral. The residual q-dependent denominator after the k-integration is treated as a propagator for the q-integration after extracting all coefficients in front of the term q 2 . We introduce a fifth Feynman parameter r for this propagator and 1 − r for the last purely q 2 -dependent propagator of Eq. (54). Applying the second shift both in the numerator and denominator we perform the symmetric q-integration. In this way, we finally arrive at an integral of the type with x = (x, y, z, v, r) and d 5 x = dx dy dz dv dr. The term H( x) denotes the full numerator and includes singular and higher powers of the dimensional regulator ǫ. N( x) is the final denominator, where the ratios are defined as ρ Q = m 2 This denominator is maximally a second-order polynomial in all Feynman parameters we have introduced. The poles of H( x) in ǫ originate from powers of k 2 and q 2 in the numerators of the k-and q-integrals. We have chosen the convention to normalize all mass parameters to the gluino mass Mg. In order to cope with the LO form factor in an easier way, we have rewritten the coefficients of all integrals as The factors ρ ǫ Q (1 + 2ǫ)(1 + ǫ 2 ζ 2 ) are added to the integrands before expansion in ǫ. The final contribution to the coefficient C A Q,SQCD is then given by The final integral is finite for this diagram. For the other diagrams, we follow the same procedure accordingly. All diagrams are infrared finite, since all virtual particles are massive, but the residual Feynman integrals contain end-point singularities in several cases that are subtracted in the usual way according to the description of Ref. [43]. The integration of the subtracted part yields the corresponding UV singularities.

Integration by Parts
In our numerical analysis, we cross the virtual bb, tt thresholds and for large pseudoscalar masses theb 1b * 2 ,t 1t * 2 thresholds as well. The parametrization of the integrals discussed so far is not sufficiently stable above these thresholds due to the high power of the denominator N( x) that becomes small in the Feynman-parameter regions in the vicinity of the virtual thresholds. We need to adopt imaginary regulatorsǭ < ∼ 10 −3 in order to obtain numbers independent of this regulator. The small size of this regulator makes the integral numerically unstable. A stabilization of the integration can be achieved by an integration by parts (IBP) to reduce the power of the denominator. In general, for this purpose, one can write where N is the dominator of the integral, p 0 and p i are polynomials and ∆ is constant in the variables x i . For simplicity we drop the arguments x everywhere. The polynomials p 0 and p i can be found by constructing the Gröbner basis of the set {N, ∂N ∂x i }. We find that These equations can be applied iteratively to reduce the power of the denominator further. Not every choice of the parameters for the integration by parts will yield a stable result. Potential issues can arise from singularities in the boundary terms as well as singularities that arise when ∆ = 0, which can happen when N = 0 and all ∂N ∂x i = 0. Choosing only a subset of the Feynman parameters yields shorter expressions that can be evaluated faster. For practical purposes, it is thus usually best to find a parametrization where using a single Feynman parameter for the integration by parts is sufficient to stabilize the numerical integration. We exemplify the two examples encountered in our calculation. If N is linear in the Feynman parameter x 1 , i.e.
there are two possible choices for the polynomials A linear combination of these two solutions is also valid. If N is quadratic in the Feynman parameter i.e.
the polynomials are given by For example in the first diagram we have achieved stabilization for x 1 = v. The denominator is linear in this parameter and we have With this explicit parametrization at hand, the following manipulation can be performed, according to Eq. (64). Since the powers of all denominators are reduced and the original denominator N( x) appears in the argument of a logarithm in the last integral the numerical integration appears to be stable for the imaginary regulator down toǭ < ∼ 10 −4 which is sufficient for the narrow-width limit.
In cases of a Feynman parameter entering the denominator in second order, and making make use of the identities of Eq. (69) (we drop the arguments of N) we arrive at the special situation that the derivative appears in second power. This allows us to perform two IBPs of the original integral [43], where for simplicity we dropped the arguments x everywhere.

Renormalization
In our calculation of the genuine SUSY-QCD corrections we have to renormalize the SUSY-QCD part of the quark mass only, since everything else is already accounted for by the QCD corrections, i.e. the decoupling of all SUSY particles from the evolution of the strong coupling α s and the PDFs that both run with five active flavours in our calculation. The SUSY-QCD part of the on-shell quark-mass counterterm is given by [see Eq. (33)] (75) We renormalize the quark mass on-shell, because the LO form factor A A Q (τ Q ) and the pure QCD corrections are expressed in terms of the quark pole mass 3 . The corresponding counterterm for the gluon-fusion cross section form factor is given by whereÃ A Q (τ Q ) denotes the LO form factor including O(ǫ) terms, with the usual trilogarithms, The derivative is given by where Li 2 denotes the dilogarithm, However, we introduce effective low-energy Yukawa couplings in our calculation, i.e. the Yukawa couplings of a low-energy Two-Higgs-Doublet model (2HDM), where the heavy SUSY particles are integrated out. This implies that the top-and bottom-Yukawa couplings are dressed with ∆ t/b contributions. The SUSY-QCD parts of these contributions are given by .
The expressions for the Yukawa couplings including resummations of the leading cot βenhanced contributions for the top-Yukawa coupling and the tgβ-enhanced terms of the bottom Yukawa coupling can be cast into the form with r Q defined after Eq. (17). These contributions will result in additional terms in the counterterms of our calculation, since the LO form factors A A Q (τ Q ) are proportional to the linear quark-Yukawa coupling. This results in the complete counterterm

Hadronic Cross Section
Our notation can be viewed as a modification of the factor σ A 0 of Eq. (1) as a starting point that can easily be extended to the NLO corrections, whereg A Q (Q = t, b) denote the resummed quark Yukawa couplings of Eq. (82) that absorb ∆ b and ∆ t contributions in the effective Yukawa couplings as the appropriate effective Yukawa couplings in the low-energy effective 2HDM. The factors C A Q,SQCD andĈ A Q,SQCD (Q = t, b) denote the relative SUSY-QCD corrections factors to the individual form factors with and without absorption of the ∆ Q terms, respectively. Within this framework the Yukawa couplings of the QCD corrections will be replaced by these effective Yukawa couplings as well due to the factorizing properties of EFT couplings from the pure QCD corrections. However, the subleading contributions of Eq. (85) involve the LO Yukawa coulings, since ∆ t,b effects only factorize at the leading order of an 1/M 2 SU SY expansion so that the SUSY-QCD remainder does not factorize from the effective Yukawa couplings in general. This will avoid artificial singularities in the scalar MSSM Higgs sector as well [45]. Expressing the LO factor σ A 0 in terms of the effective Yukawa couplings, and referring to Eq. (6), the SUSY-QCD corrections add to the virtual coefficient C A , with the usual QCD-correction coefficient C A QCD and where we are using LO Yukawa couplings g A Q in the numerator, since this contribution constitutes the remainder of the full SUSY-QCD corrections that does not factorize in general terms. In Eq. (85) and for the following discussion of the results, we distinguish between this coefficient for the SUSY-remainder and the corresponding coefficient 4 , that describes the full SUSY-QCD corrections without introducing the effective top and bottom Yukawa couplings, i.e. without absorbing ∆ Q terms in the Yukawa couplings. The contributions ∆A Q,SQCD are given in Eq. (83).

Axial γ 5 Schemes
We have implemented the Larin scheme of Ref. [27] that is a variant of the original 't Hooft-Veltman scheme that has been set-up systematically by Breitenlohner and Maison [26]. We have extracted the Levi-Civita tensor at the pseudoscalar vertex by means of the replacement and just keeping the four γ matrices inside the traces. The diagrams where the pseudoscalar couples to squarks do not have such a vertex. They are however finite, such that a naively anticommuting γ 5 can be used at NLO. The chiral couplings at the quark-squark-gluino vertices are treated fully anticommuting to arrive at traces with one or no γ 5 matrix. Only the contributions with no additional γ 5 matrix contribute after applying the projector of Eq. (49). The projector yields a product of two Levi-Civita tensors that is defined as where the metric tensors inside this determinant are treated as n-dimensional objects. This prescription avoids a splitting of γ matrices and loop momenta into 4-and (n−4)-dimensional components. Since γ 5 as defined in Eq. (90) does not anticommute, an anomalous counterterm has to be added. However the genuine SUSY-QCD contributions to this counterterm vanish. In the 't Hooft-Veltman scheme, the metric tensors in this determinant are defined as strictly 4-dimensional objects so that the numerators of the loop integrals split into 4-and (n−4)-dimensional pieces that have to be treated separately. To avoid additional anomalous counterterms we used anticommuting γ 5 matrices at the QQg-vertices in this scheme as well.
We found full agreement for both schemes. In addition, we have lifted the anti-commuting properties of the γ 5 matrices entering at the QQg-vertices and found mismatches that require anomalous subtractions to restore the chiral properties. Finally, we have implemented the γ 5 scheme of Ref. [28] that gives up the cyclicity of the traces but keeps the full anti-commuting property of the γ 5 matrix. The cyclicity of the trace is equivalent to the arbitrary decision where we start to read the fermion lines. To resolve this, the scheme defines unambiguous reading points in each diagram relative to the external axial couplings. However, since we have no axial vector couplings in our diagrams, the only prescription we have to follow is that the reading point must be outside of subdivergences, e.g. in the fifth diagram of Fig. 5 the reading point must not be at the ggg vertex. We found full agreement with the calculation in the Larin scheme as well. Finally, we have reproduced the limit of large top, stop and gluino masses of Ref. [32] and found full agreement. Ref. [32] worked with Pauli-Villars regularization so that their Clifford algebra is defined in four dimensions strictly resulting in a fully anti-commuting γ 5 . That we have found full agreement with this calculation in the large-mass limit underlines the consistency of our results.

Adler-Bardeen Theorem
According to the analytical results of Ref. [32] the SUSY-QCD coefficient in the large SUSYmass limit (keeping the quark mass small) is given bŷ with ρ i as defined after Eq. (60). In this expression, we have focused just on the leading terms of the large-mass expansion, since this is the relevant contribution of the matching to a low-energy 2HDM. Moreover, in the expression above the ∆ Q terms are not subtracted, i.e. this is the result in terms of the LO Higgs coupling g A Q without ∆ Q -dressing. The coupling Y Q is related to the squark coupling, Inserting the explicit expressions for s 2θ Q and Y Q one arrives at Since the Att operator mixes with the Att * operator the non-decoupling ∆ t contributions to the effective top Yukawa coupling are induced. Working with properly matched low-energy parameters, i.e. effective Yukawa couplings with ∆ Q contributions as in Eq. (82), this term is absorbed in the Yukawa couplings exactly so that the radiative corrections in the low-energy 2HDM with properly defined low-energy parameters are vanishing for the leading O(M 0 This is because in contrast to the MSSM, the chiral symmetry ψ Q → e iαγ 5 ψ Q is only broken by the quark mass term in the effective 2HDM so that only the higher-order corrections to the proper matching of the low-energy 2HDM to the full MSSM contribute. Thus, in the low-energy limit the Adler-Bardeen theorem [17] is fulfilled 5 . Since radiative corrections still arise due to the higher-order corrections to the matching, the Adler-Bardeen theorem [17] builds a deep connection between the explicit structure of the radiative corrections in the full MSSM in the low-energy limit and the radiative corrections in the low-energy EFT. This result is in line with the result of Ref. [46] that the QCD corrections to the effective ggA Lagrangian in the HTL are vanishing if the strong coupling is chosen as the 5-flavour one, i.e. properly decoupling the top-quark contribution from the running of α s or in other words using the properly matched low-energy α s within pure 5-flavour QCD. This also implies that in the large SUSY-mass limit (keeping the top mass small in comparison) no effective ggA operator is generated in the low-energy 2HDM at the dimension-5 level by integrating out the SUSY particles. The same is true as well for the bottom/sbottom contributions so that the SUSY particles do not generate a sbottom-induced effective ggA operator at leading O(M 0 SU SY ) at all. Another situation arises when the top quark is integrated out, i.e. assumed to be much heavier than the pseudoscalar A as well and not assumed to be much lighter than the other SUSY particles. In this case a dimension-5 operator contribution is generated on top of the HTL at LO due to the non-decoupling nature of the top quark as has already been observed for the leading O(G F m 2 t ) corrections to the effective Agg coupling [47]. Since the stops couple to the pseudoscalar in terms of the top Yukawa coupling as well, a new genuine dimension-5 contribution to the Agg coupling, on top of the contribution from the effective Yukawa coupling, emerges starting at NLO. This contribution can be related to the violation of the global Peccei-Quinn symmetry of the MSSM Lagrangian by the µ term [31,32]. This leads to an extension of the related operator identity of the divergence of the axial-vector current by an additional operator involving the stop fields thus destroying the one-to-one correspondence between the pseudoscalar top-Yukawa coupling and the ABJanomaly operator and in this way the translation of the Adler-Bardeen theorem to the Agg operator. It follows that both the ∆ t terms and the genuine radiative corrections to the Agg coupling scale with the µ parameter [48]. Within the EFT view this has to be considered as higher-order corrections to the effective Agg operator in the combined HTL and large-SUSYmass limit, i.e. higher-order corrections to the corresponding matching conditions that scale with µ.

Results
We are now in the position to present and discuss the final results of the NLO SUSY-QCD corrections to pseudoscalar gg → A production, but also to the pseudoscalar decays A → gg and A → γγ. For the numerical analysis we have adopted the M 125 h benchmark scenario [19] that is defined by the following on-shell parameters,

28
As a starting point, the perturbative NLO coefficients C A Q,SQCD are displayed in Figs. 6 and 7 as a function of the pseudoscalar mass M A for tgβ values of 10 and 40, respectively. In these figures, we show the approximate calculations of Ref. [32] as well, i.e. for the stop contribution both approximations of the combined heavy-top/SUSY limit ('approx heavy ') and the pure large SUSY-mass limit ('approx mt '), while for the sbottom contribution only the large SUSYmass limit ('approx') is phenomenologically relevant and shown. The full calculation agrees well with the former approximate calculations for smaller pseudoscalar masses in both the stop and sbottom cases. However, we observe sizeable and increasingly relevant deviations for pseudoscalar masses approaching or exceeding the virtual squark threshold 6 . Moreover, we display the results of the NLO coefficients for the two cases of absorbing the ∆ t/b terms in the corresponding Yukawa couplings and the opposite. It is clearly visible that the ∆ t/b terms approximate the full results quite well for smaller pseudoscalar masses M A so that the results after subtracting them turn out to be quite small. These subtracted results represent the SUSY-remainder, i.e. the contributions beyond the leading parts corresponding to the effective top and bottom Yukawa couplings. It is obvious that the absorption of these contributions leads to a much better perturbative behaviour thus corroborating the effective Yukawa-coupling approach. This is further underlined by the tgβ dependence of the stop and sbottom contributions shown in Fig. 8 for a pseudoscalar mass M A = 1.5 TeV. The description of the SUSY-QCD corrected cross section in terms of the effective low-energy top-and bottom-Yukawa couplings leads to a moderate SUSY-remainder at NLO as long as the pseudoscalar Higgs mass does not approach the virtual stop/sbottom thresholds. At and beyond these virtual thresholds, the SUSY-remainders turn out to be sizeable.
As the next step, we analyze the SUSY-QCD corrections to the hadronic cross section of pseudoscalar Higgs-boson production via gluon fusion. The effect of the corrections on the K-factor at the hadronic level, which is defined as the ratio between the NLO and LO cross sections, is discussed first. We adopt the MSHT20nlo as118 parton density functions and perform the analysis for a c.m. energy of 13 TeV at the LHC. Fig. 9 exhibits the K-factor for tgβ = 10, 40 with effective top-and bottom-Yukawa couplings for the QCD part of the cross section and for the corresponding results of the previous approximate calculations. The QCD part of the K-factors shows the usual sizeable NLO corrections of about 30-50%, while the additional SUSY-QCD remainder turns out to be small or moderate. The comparison implies that effects beyond the approximation become relevant when approaching the virtual stop/sbottom thresholds and above as expected.
These K-factors can be translated to the hadronic production cross sections of pseudoscalar Higgs bosons via gluon fusion as shown in Fig. 10 for two values of tgβ = 10, 40. For the effective bottom-Yukawa couplings, we include the full set of NNLO corrections [51] to lift the accuracy of the factorizing and dominant contributions to the NNLO level, while for the effective top-Yukawa coupling we use the NLO expression in the effective field-theory framework. Here, we present the QCD-corrected cross sections without the effective topand bottom-Yukawa couplings, i.e. without any genuine SUSY-QCD corrections and the approximate and full SUSY-QCD corrected cross sections with the effective Yukawa couplings  as discussed in the previous section. The comparison of the full QCD-corrected cross section (blue line) and the full QCD + SUSY-QCD corrected cross section (red line) supports the high relevance of the SUSY-QCD corrections in total, while the SUSY-remainder plays a role close or above the virtual stop-and sbottom thresholds.

The Gluonic Decay A → gg
The same virtual coefficient as for gg → A contributes to the genuine SUSY-QCD corrections of the gluonic pseudoscalar Higgs decay A → gg according to Eq. (11). The relative QCD and SUSY-QCD corrections to the gluonic decay width are shown in Fig. 11 with the use of effective top and bottom Yukawa couplings. It is clearly visible that the bulk of the genuine SUSY-QCD corrections can be absorbed by the effective top and bottom Yukawa couplings including ∆ t,b contributions. The SUSY-remainder is relevant in regions where finite squark-mass effects become relevant, i.e. close or above the related virtual thresholds. The corresponding partial decay widths Γ(A → gg) are shown in Fig. 12 for tgβ = 10, 40, using effective top and bottom Yukawa couplings for the SUSY-QCD-corrected decay widths, but LO couplings without ∆ t.b terms for the LO and QCD-corrected decay widths. The SUSY-QCD corrections are treated in the same way as for the production cross sections, i.e. ∆ b terms at two-loop order and ∆ t contributions at one-loop level. The main effect of the genuine SUSY-QCD corrections emerges from the factorizing ∆ b,t corrections to the effective top and bottom Yukawa couplings leaving a sizeable SUSY-remainder in regions only where squark-mass effects become relevant. The extended peaking structure around a pseudoscalar mass of 2 TeV originates from the two chargino thresholds that are not affected by corrections due to strong interactions. It should be noted that the same corrections are valid for the reverse process γγ → A as well, which could be probed at a potential future high-energy γγ-collider.

Conclusions
We have calculated the full SUSY-QCD corrections to pseudoscalar Higgs-boson production via gluon fusion gg → A within the MSSM at hadron colliders. We implemented the virtual stop and sbottom sector at the NLO level to be in line with the necessities for the corresponding scalar Higgs-boson production cross sections via gluon fusion gg → h, H. We have analyzed pseudoscalar Higgs-boson production with respect to the introduction of effective low-energy top and bottom Yukawa couplings, i.e. the couplings within the lowenergy 2HDM after integrating out the strongly interacting SUSY particles (stops, sbottoms and gluinos). We found that the bulk of the NLO corrections can be absorbed in these effective Yukawa couplings, while the SUSY-remainder is of moderate size, being significant close or above virtual squark thresholds. We have analyzed the corrections in the context of the Adler-Bardeen theorem and found that this theorem is fulfilled in the large SUSY-mass limit, if the observable is expressed in terms of properly matched low-energy parameters, i.e. top-and bottom-Yukawa couplings. The analogous results have also been obtained for the related rare pseudoscalar Higgs-boson decays A → gg, γγ that, however, only play a minor role in phenomenological analyses at hadron colliders. This work completes the full NLO QCD calculation for pseudoscalar MSSM Higgs production and decay into gluonic and photonic final states and thus serves as a basis for the corresponding theoretical predictions.