Two-loop QCD corrections to the MSSM Higgs masses beyond the effective-potential approximation

We compute the two-loop QCD corrections to the neutral Higgs-boson masses in the MSSM, including the effect of non-vanishing external momenta in the self-energies. We obtain corrections of O(alpha_t*alpha_s) and O(alpha*alpha_s), i.e., all two-loop corrections that involve the strong gauge coupling when the only non-vanishing Yukawa coupling is the top one. We adopt either the DRbar renormalization scheme or a mixed OS-DRbar scheme where the top/stop parameters are renormalized on-shell. We compare our results with those of earlier calculations, pointing out an inconsistency in a recent result obtained in the mixed OS-DRbar scheme. The numerical impact of the new corrections on the prediction for the lightest-scalar mass is moderate, but already comparable to the accuracy of the Higgs-mass measurement at the LHC.


Introduction
The accuracy of the measurement of the Higgs-boson mass by the ATLAS and CMS collaborations at the Large Hadron Collider (LHC) has already reached the level of 300-400 MeV [1,2] and, being still dominated by statistics, is bound to improve further when the LHC restarts operations in 2015. This puts new emphasis on the need for highprecision calculations in those extensions of the Standard Model (SM), such as the Minimal Supersymmetric Standard Model (MSSM), in which the Higgs-boson mass can be predicted as a function of other physical observables.
The Higgs sector of the MSSM consists of two SU (2) doublets, H 1 and H 2 , whose relative contribution to eleca e-mail: degrassi@fis.uniroma3.it b e-mail: divita@mpp.mpg.de c e-mail: slavich@lpthe.jussieu.fr troweak (EW) symmetry breaking is determined by the ratio of vacuum expectation values (VEVs) of their neutral components, tan β ≡ v 2 /v 1 . The spectrum of physical Higgs bosons is richer than in the SM, consisting of two neutral scalars, h and H , one neutral pseudoscalar, A, and two charged scalars, H ± . At the tree level, the neutral scalar masses m h and m H and the scalar mixing angle α can be computed in terms of the Z -boson mass m Z , the pseudoscalar mass m A and tan β, and the bound m h < | cos 2β| m Z applies. In a significant portion of the parameter space the lightest scalar h has SM-like couplings to fermions and gauge bosons, in which case the tree-level bound on m h has long been disproved by the LEP [3,4]. However, radiative corrections can raise the prediction for the lightest-scalar mass up to the value m h ≈ 125 GeV observed at the LHC, and they bring along a dependence on all MSSM parameters. Among the latter, particularly relevant are the masses and mixing of the scalar partners of the third-generation quarks, the stop and sbottom squarks.
Due to the crucial role of radiative corrections in pushing the prediction for the lightest-scalar mass above the tree-level bound, an impressive theoretical effort has been devoted over more than 20 years to the precise determination of the Higgs sector of the MSSM. 1 After the early realization [5][6][7][8][9] of the importance of the one-loop O(α t ) corrections 2 involving top and stop, full one-loop computations of the MSSM Higgs masses have been provided [10][11][12][13], leading logarithmic effects at two loops have been included via renormalization-group methods [14][15][16][17], and genuine two-loop corrections of O(α t α s ) [18][19][20][21][22][23][24][25], O(α 2 t ) [18,24,26], O(α b α s ) [27,28] and O(α t α b + α 2 b ) [29] have been evaluated in the limit of vanishing external momentum in the Higgs self-energies. All of these corrections have been implemented in widely used computer codes for the calculation of the MSSM mass spectrum, such as FeynHiggs [30], SoftSUSY [31,32], SuSpect [33] and SPheno [34,35]. Furthermore, a complete two-loop calculation of the MSSM Higgs masses in the effective potential approach (i.e., at zero external momentum), including also two-loop corrections controlled by the EW gauge couplings, has been presented in Refs. [36,37]. Some of the dominant three-loop corrections to m h have also been obtained, both via renormalizationgroup methods [38][39][40] and by explicit calculation of the Higgs self-energy at zero external momentum [41,42].
Already at the two-loop level, going beyond the approximation of zero external momentum brings significant complications to the calculation of the Higgs self-energies. Different algorithms for expressing all two-loop self-energy integrals with arbitrary external momentum in terms of a minimal set of Master Integrals (MIs) were developed in Refs. [43][44][45][46]. However, explicit analytical formulas for the MIs can be derived only for special values of the masses of the particles circulating in the loops, whereas in the general case a numerical calculation becomes unavoidable. A method to compute all the MIs of Ref. [46] by numerically solving a system of differential equations in the external momentum was developed in Ref. [47], extending earlier results of Refs. [48][49][50][51]. A library of routines for the computation of the MIs with the method of Ref. [47] was then made available in the package TSIL [52].
A calculation of the two-loop contributions to the Higgs self-energies involving the strong gauge coupling or the thirdfamily Yukawa couplings, based on the methods of Refs. [47,52], was presented in Ref. [53]. That calculation goes beyond the two-loop results implemented in public codes [19][20][21][25][26][27][28][29] in that it includes external-momentum effects, as well as contributions involving the D-term-induced EW interactions between Higgs bosons and sfermions. When combined with the effective-potential results of Refs. [36,37], the results of Ref. [53] provide an almost-complete two-loop calculation of the Higgs masses in the MSSM -the only missing contributions being external-momentum effects that involve only the EW gauge couplings. However, no code for the calculation of the MSSM mass spectrum implementing the results of Refs. [36,37,53] was ever made available, and the way those results are organized does not lend itself to a straightforward implementation in the existing public codes. On one hand, the DR renormalization scheme adopted in Refs. [36,37,53] for the parameters of the MSSM lagrangian does not match the "mixed on-shell (OS)-DR" scheme adopted in FeynHiggs. On the other hand, implementation of the results of Refs. [36,37,53] in SoftSUSY, SuSpect and SPheno, which also adopt the DR scheme, is complicated by the fact that in Refs. [36,37,53] the running masses of the Higgs bosons entering the loop corrections are defined by the second derivatives of the tree-level potential. While this choice amounts to a legitimate reshuffling of terms between different perturbative orders, it restricts the applicability of the calculation to rather specific ranges of renormalization scale where none of the running Higgs massesas defined in Refs. [36,37,53] -is tachyonic. Perhaps as a consequence of these complications, a full decade after the publication of Ref. [53] its results have yet to be included in phenomenological analyses of the MSSM Higgs sector.
In this paper we present a new calculation of the momentum-dependent part of the two-loop corrections to the neutral Higgs masses of O(α t α s ), i.e. those involving both the top Yukawa coupling and the strong gauge coupling. We also compute "mixed" two-loop corrections that we denote by O(αα s ): they involve both the strong gauge coupling and the EW gauge couplings, under the approximation that the only non-vanishing Yukawa coupling is the top one. It is natural to consider these two classes of corrections together, because in both of them the dominant terms affecting the lightest-scalar mass are expected to be suppressed by a factor of O(m 2 Z /m t 2 ) with respect to the zero-momentum O(α t α s ) corrections (in practice, we find that both classes of corrections are considerably more suppressed than that, but still comparable to each other in size).
In our calculation we rely on the integration-by-parts (IBP) technique of Refs. [54,55] to express the momentumdependent loop integrals in terms of the MIs of Ref. [46], which we evaluate by means of the package TSIL. We obtain results for both the DR and the OS-DR renormalization schemes, organized in such a way that they can be directly implemented in the existing codes for the computation of the MSSM mass spectrum. We verify that our results are in full agreement with the ones of Ref. [53] where they overlap. After describing our calculation in some detail, we briefly discuss the numerical impact of the momentum-dependent O(α t α s ) and O(αα s ) corrections to the Higgs masses in a set of representative points in the MSSM parameter space.
While our paper was in preparation, an independent calculation of the momentum-dependent O(α t α s ) corrections to the neutral Higgs masses in the MSSM appeared [56], relying on the results of Ref. [43] for the decomposition of twoloop integrals into MIs and on the package SecDec [57,58] for the numerical evaluation of the latter. The results of that calculation are expressed in the OS-DR scheme, and they have been implemented in the latest version of FeynHiggs. Although we have verified that our results for the contributions of genuine two-loop diagrams involving the stronggauge and top-Yukawa couplings agree numerically with those of Ref. [56], we do not reproduce the overall values of the momentum-dependent O(α t α s ) corrections to the Higgs masses. We trace the reason for the discrepancy to an inconsistency in Ref. [56] concerning the definitions of the wave-function-renormalization (WFR) constants for the Higgs fields and of the parameter tan β.

Neutral Higgs boson masses in the MSSM
We outline here our calculation of the two-loop corrections to the masses of the neutral Higgs bosons in the MSSM with real parameters (we do not consider the possibility of CP violation in the Higgs sector). We decompose the neutral components of the two Higgs doublets into their VEVs plus their CP-even and CP-odd fluctuations as follows: The CP-odd components P 1 and P 2 combine into the pseudoscalar A and the neutral would-be Goldstone boson G 0 . The CP-even components S 1 and S 2 combine into the scalars h and H . The squared physical masses of the latter are the two solutions for p 2 of the equation where S ( p 2 ) denotes the 2 × 2 inverse-propagator matrix in the (S 1 , S 2 ) basis, p being the external momentum flowing into the scalar self-energies. We can decompose S ( p 2 ) as where M 2 0 denotes the tree-level mass matrix written in terms of renormalized parameters, and M 2 ( p 2 ) collectively denotes the radiative corrections. At each order in the perturbative expansion, the latter include both the contributions of one-particle-irreducible (1PI) diagrams and non-1PI counterterm contributions arising from the renormalization of parameters that enter the lower-order parts of S ( p 2 ). We express M 2 0 in terms of the parameter tan β renormalized in the DR scheme, and of the physical masses of the pseudoscalar and of the Z boson using (here and thereafter) the shortcuts c φ ≡ cos φ and s φ ≡ sin φ for a generic angle φ. Neglecting terms that do not contribute at O(α t α s ) or O(αα s ), our choices for the parameters entering M 2 0 lead to the following expressions for the two-loop parts of the individual entries of M 2 ( p 2 ): δ tan β (2) tan β +δZ In the equations above, T (2) i and (2) i j (with i, j = 1, 2) denote the two-loop parts of tadpoles and self-energies, respectively, for the scalars S i , while (2) AA and (2) Z Z denote the two-loop parts of the pseudoscalar and Z -boson selfenergies. In addition, δZ (2) i (with i = 1, 2) in Eqs. (5)-(7) denote the two-loop parts of the WFR counterterms for the Higgs fields H 0 i , which we renormalize as follows: where in the expansion of the square root we have again neglected terms that do not contribute at O(α t α s ) or O(αα s ).
We adopt a DR definition for the Z i , which can then be determined from the anomalous dimensions of the Higgs fields and from the β functions of the couplings entering the anomalous dimensions. Taking from the general formulas of Refs. [59,60] only the terms relevant to our approximation, we get where N c = 3 and C F = 4/3 are color factors, = (4−d)/2 in dimensional reduction, and the coupling α t entering the one-loop counterterm δZ (1) 2 is in turn renormalized in the DR scheme. Finally, δ tan β (2) in Eqs. (5)-(7) denotes the two-loop part of the counterterm for the parameter tan β. The choice of a DR definition for tan β implies that, in our approximation, its counterterm can be expressed via the WFR counterterms: All tadpoles and self-energies in Eqs. (5)-(7) include both 1PI two-loop contributions and non-1PI contributions arising from the renormalization of the parameters entering their one-loop counterparts. Since we are focusing on the O(α t α s ) and O(αα s ) corrections to the Higgs masses, we need to introduce counterterms only for the parameters that are subject to O(α s ) corrections, namely the top mass m t , the stop masses mt 1 and mt 2 , the stop mixing angle θ t , the soft supersymmetry-breaking Higgs-stop coupling A t and the masses mq i of all squarks other than the stops. The latter enter the one-loop tadpoles and self-energies of the Higgs bosons via D-term-induced EW couplings, and the one-loop self-energy of the Z boson via the gauge interaction. In our calculation we neglect all Yukawa couplings (and hence quark masses) other than the top one, 3 therefore none of the other squarks mix. We obtained results for the O(α t α s ) and O(αα s ) contributions to tadpoles and self-energies assuming that the relevant quark/squark parameters are renormalized either in the DR or in the OS scheme. Formulas for the DR-OS shift of the parameters in the top/stop sector can be found, e.g., in appendix B of Ref. [25], while the shifts for the remaining squark masses can be obtained by setting m t = θ t = 0 in the corresponding formulas for the stop masses. We remark that the right-hand sides of Eqs. (5)- (7) are constructed to give finite entries in the inverse-propagator matrix of the scalars. Indeed we have explicitly verified that -after summing all 1PI and counterterm contributions -the 1/ 2 and 1/ poles in the right-hand sides of Eqs. (5)-(7) cancel out.
In principle, the two-loop contributions to the Higgs inverse propagator given in Eqs. (5)-(7) must be combined with a full calculation of the corresponding one-loop contributions, and used to determine the physical Higgs masses by solving directly Eq. (2). However, as will be discussed in 3 The corrections to the Higgs masses involving the bottom Yukawa coupling could become relevant for large values of tan β. When all parameters in the bottom/sbottom sector are renormalized in the DR scheme, the O(α b α s ) corrections can be obtained from the corresponding results for the O(α t α s ) corrections via trivial replacements. On the other hand, an OS renormalization of the bottom/sbottom parameters would entail additional complications, as discussed in Refs. [27][28][29]. Anyway, the regions of the MSSM parameter space where the O(α b α s ) corrections to the Higgs masses are most relevant are being severely constrained by direct searches of Higgs bosons decaying to tau leptons at the LHC [61][62][63][64]. the next section, the computing times required for the evaluation of momentum-dependent two-loop integrals are not negligible. A numerical search for the solutions of Eq. (2) could significantly slow down the codes for the calculation of the Higgs masses, making them unsuitable for extensive phenomenological analyses of the MSSM parameter space. It is therefore convenient to compute the Higgs masses in two steps, with a procedure similar to the one discussed in Refs. [10,26]. We first call FeynHiggs, which solves Eq.
Concerning the former, we have where we define  7): where again we take the real part of all two-loop self-energies.
We remark that Ref. [56] proposes an alternative two-step procedure to include the momentum-dependent parts of the O(α t α s ) corrections in FeynHiggs, differing from the one outlined above only by higher-order effects.
Finally, a comment is in order about the dependence of the corrections to the Higgs masses on the WFR constants. 4 In principle, the predictions for the physical Higgs masses at a given order in the perturbative expansion should not depend directly on the WFR constants (although they could still depend indirectly on them via the tan β counterterm). Indeed, if the two-loop contributions to the inverse-propagator matrix are computed with p 2 equal to the tree-level scalar masses and then rotated to the mass-eigenstate basis via the tree-level mixing angle, so that the computation is performed strictly at the two-loop level, the terms in Eqs. (5)-(7) that depend explicitly on δZ (2) i drop out of the mass corrections m 2 h,H . On the other hand, if the loop-corrected scalar masses m 2 h and m 2 H and the effective mixing angle α are used, as in Eqs. (12)-(15) above, or if the zeroes of the inverse-propagator matrix are determined numerically, the corrections to the scalar masses retain a dependence on the WFR counterterms. In our calculation we adopt a DR definition for the WFR; therefore the terms involving δZ (2) i are purely divergent and cancel out against other divergent terms in the individual entries of the inverse-propagator matrix, hence they do not show up in Eqs. (12) and (13). If, however, one adopts a non-minimal definition of the WFR, Higgs-mass corrections computed as in Eqs. (14) and (15) will contain non-vanishing terms that depend on the finite part of δZ (2) i . Albeit formally of higher order in the perturbative expansion, these terms can be numerically relevant when the loop-corrected scalar masses differ significantly from their tree-level values (as is the case for a SM-like scalar h with mass around 125 GeV).

Calculation of two-loop diagrams with nonzero momentum
The computation of the two-loop corrections to the neutral MSSM Higgs masses considered in this paper requires the knowledge of the tadpole and self-energy diagrams entering Eqs. (5)- (7). While the strategy for the computation in the zero-momentum approximation is well known, the evaluation of the self-energies with arbitrary external momentum is more involved. We illustrate in this section the details of our calculation, which we performed in a fully automated way. The relevant diagrams are generated with FeynArts [67], using a modified version of the original MSSM model file [68] that implements the QCD interactions in the background field gauge. The diagrams contributing to the vacuum polarization of the Z boson are contracted with a suitable projector in order to single out their transverse part. The color factors are simplified with a private package and the Dirac algebra is handled by FORM [69]. The computation is performed in dimensional reduction, which we can implement in this case by generating the diagrams in dimensional regularization and replacing, in each diagram involving an internal d-dimensional gluon, g μν → g μν + gμν (where gμν is the 2 -dimensional metric tensor) in order to include the corresponding -scalar contribution. We are then left with Feynman integrals of the form where α, . . . , η, a 1 , . . . , a 5 are positive (or zero) integer exponents and the D i 's are defined as Integrals where the exponents n i ∈ Z. In the present case, one has to evaluate O(300) different Feynman integrals. There exists a convenient procedure for dealing with such large numbers of different Feynman integrals in a more efficient way than their direct evaluation. Dimensionally regularized integrals, at arbitrary loop order and with arbitrary number of external legs, satisfy identities of the IBP type [54,55]. These identities are linear relations that connect integrals with different sets of exponents {n 1 , . . . , n 5 }. After a set of independent integrals, the MIs, has been identified, all the remaining integrals can then be reduced to linear combinations of the MIs, the coefficients being rational functions of the masses, the kinematic invariants and the space-time dimension d. One practical advantage of such a procedure is its "divide and conquer" spirit. On the one hand few MIs encode the analyticity properties (singularities, thresholds, branch cuts) of the problem under consideration. On the other hand, the evaluation of the large number of different Feynman integrals entering a computation is reduced to a problem of linear algebra, which can easily be automated, if the MIs are known. We perform the reduction to MIs with the public code REDUZE [70,71], which implements the Laporta algorithm [72] and produces the IBP identities relevant to our case.
The evaluation of the MIs is in general a complicated problem and can proceed via different techniques, like the integration in the Feynman, Schwinger or Mellin-Barnes representations. A remarkable consequence of the aforementioned IBP relations is that the MIs obey linear systems of firstorder differential equations (DEs) in the kinematic invariants, which provide an alternative means for their computation [73,74]. Finding the analytic solution of the DEs for arbitrary d is possible only in some simple cases. In more general cases the MIs are expanded in powers of = (4 − d)/2, giving rise to a (generally coupled) system of DEs for the expansion coefficients. In the limit of vanishing external momentum, two-loop self-energies become two-loop vacuum diagrams, which reduce via IBP to linear combinations of only one genuine two-loop MI and products of the one-loop onepropagator MI [75]. Two-loop self-energies with arbitrary external momentum, in the general case with five different masses in the loops, reduce via IBP to linear combinations of 30 MIs [46]. The finite part of such MIs can be expressed in terms of four functions, in addition to the well-known oneloop MIs. 5 As already mentioned, analytic solutions for such functions have been derived only for special patterns of up to two internal masses (only one function is known in a particular case with three different masses). On the other hand, the diagrams that in our approximation contribute to the selfenergies entering Eqs. , and m 2 g . In our computation we rely on the package TSIL [52], which implements (besides all the analytically known cases) the numerical solution of the DEs for the two-loop selfenergy MIs. The method of Refs. [47][48][49][50][51] is based on the fact that the value at p 2 = 0 (or the expansion for small p 2 in the case of logarithmically divergent integrals) is known for each function and can be used to build the set of initial conditions needed for the solution of the DEs. In the computation of the self-energies entering Eqs. (5)-(7) we need to evaluate the corresponding MIs at p 2 = m 2 Z , m 2 A , m 2 h , m 2 H . Given that we include the contribution of light quarks to Z Z (in the approximation m q = 0) and that m A is a free parameter, it is clear that the way thresholds are handled in the numerical evaluation of the MIs is of particular relevance. In the DEs approach, the physical two-and threeparticle thresholds show up, together with the pseudothresholds, as poles in the coefficients of the DEs. TSIL overcomes 5 In the presence of infrared divergences associated to loops of massless quarks, the IBP reduction of the considered diagrams to MIs requires two additional functions. Their expression in terms of logarithms and polylogarithms can be obtained from the results of Ref. [76]. the numerical instabilities related to such poles by displacing the p 2 -integration contour in the upper half-plane when the momentum is above the smallest (pseudo)threshold. The evaluation at (or very close to) the (pseudo)thresholds is performed through a variant of the algorithm, which is slightly less efficient but ensures reliable results in such critical cases. As an example, the time needed on an Intel Core i7-4650U CPU for the evaluation of the complete set of MIs, for any of the mass patterns entering the self-energies, ranges between 5 × 10 −4 s and 8 × 10 −2 s, the latter being the typical time for p 2 close or equal to the heavy stop pair threshold and to the three-particle (pseudo)thresholds mt i + mg ± m t . In Ref. [52] the relative accuracy of TSIL is claimed to be better than 10 −10 for generic cases, or worse in cases with large mass hierarchies. Being TSIL a package dedicated to the evaluation of the MIs for two-loop self-energy diagrams, it is not surprising that its speed and accuracy prove much better than those quoted in Ref. [56], where the general-purpose package SecDec is used and, in the most complicated case, 100 s are needed in order to reach a relative accuracy of at least 10 −5 close to a threshold.

Numerical examples
In this section we assess the numerical impact of the momentum-dependent part of the O(α t α s ) corrections and of the whole O(αα s ) corrections on the predictions for the neutral Higgs-boson masses in the MSSM. We focus on six benchmark scenarios introduced in Ref. [77], which identify regions in the MSSM parameter space compatible with the current bounds from SUSY-particle searches and with the requirement that the predicted value of m h agrees, within the theoretical uncertainty of ±3 GeV estimated in Refs. [78,79], with the mass of the SM-like Higgs boson discovered at the LHC. 6 In our numerical examples we adopt the mixed OS-DR scheme described in Sect.  the renormalization conditions imposed both in our OS-DR calculation and in the one of Ref. [56]). By default, the renormalization scale associated to the DR definition of the Higgs WFR and of tan β is fixed as μ R = m t .
In Fig. 1 we present our predictions for the lightestscalar mass m h in the six benchmark scenarios. We choose m A = 500 GeV and tan β = 20, so that the lightest scalar h is SM-like, the bound on its tree-level mass is saturated, and the corrections controlled by the bottom Yukawa coupling, which we do not compute beyond the approximations of FeynHiggs, are not expected to be particularly relevant. For each scenario we show three bars: the upper one represents the "unperturbed" mass m h , obtained from   Figure 2 shows that the corrections to the lightestscalar mass are negligible at low values of m A , but they   Fig. 2 for the corrections to the heaviest-scalar mass become larger and essentially independent of m A as the latter increases. The transition to this "decoupling" regimewhere the lightest scalar has SM-like couplings and its mass is insensitive to the value of m A -is sharper for larger values of tan β. In both the m max h and light-stop scenarios, the momentum-dependent O(α t α s ) effects decrease m h by 300-400 MeV at large m A . However, as already seen in Fig. 1, the O(αα s ) effects enter the prediction for m h with comparable magnitude but opposite sign, significantly reducing the total size of the correction. Figure 3 shows that for low values of m A , where the heaviest scalar is the one with SM-like couplings, the corrections to its mass are comparable to the ones that affect the lightestscalar mass in the decoupling region. On the other hand, for larger values of the pseudoscalar mass -where m H ≈ m A -the corrections to the heaviest-scalar mass show a series of spikes, related to the opening of real-particle thresholds in diagrams that involve a virtual gluon. The first spike is visible in correspondence with m A = 2 m t in the plot on the left for the m max h scenario. More-pronounced spikes (note the different scale on the y axis) are visible in correspondence with m H = 2 mt 1 , m H = mt 1 +mt 2 , and m H = 2 mt 2 in the plot on the right for the light-stop scenario. Analogous spikes would appear at larger values of m A in the m max h scenario, where the stops are heavier. We stress that our results are not reliable in the vicinity of these thresholds: the two-loop correction to the heaviest-scalar mass is actually divergent there, and the height of the spikes in the plots carries no physical meaning. A more sophisticated analysis, taking into account the widths of the virtual particles in the loops as well as non-perturbative QCD effects, would be necessary around the thresholds, but it is beyond the scope of our calculation. Figure 3 also shows that, in the decoupling region and away from the thresholds, the corrections to the heaviestscalar mass amount at most to a few hundred MeV, and they decrease in size with increasing tan β. Moreover, the effect of the O(αα s ) corrections is negligible (the dashed and solid lines are practically overlapping in the plots). Inspection of our analytic formulae shows that, in the decoupling limit, the O(αα s ) corrections to m H are suppressed by one or two powers of tan β, whereas the O(α t α s ) corrections contain unsuppressed contributions proportional to the square of the superpotential Higgs-mass parameter μ. While corrections of this size might be considered negligible in comparison with the value of m H itself, they are not entirely irrelevant when compared to the difference m H − m A , which can be of the order of a few GeV and is the quantity of interest when a large physical mass for the pseudoscalar is taken as input in the calculation.

Comparison with earlier calculations
The way we compute the O(α t α s ) and O(αα s ) corrections to the entries of the inverse-propagator matrix for the neutral scalars allows for a relatively easy comparison with earlier calculations. We first renormalize all the relevant parameters in the DR scheme, i.e. we introduce minimal counterterms that, by definition, subtract only powers of 1/ , multiplied by coefficients that should be polynomial in the DRrenormalized masses and couplings. In a second step, we convert our results to the mixed OS-DR scheme adopted in FeynHiggs, replacing the DR top/stop parameters entering the one-loop part of the corrections with the corresponding OS parameters plus the finite one-loop shifts given in Ref. [25].
As a first obvious check, we took the limit of vanishing external momentum in the scalar self-energies entering the O(α t α s ) corrections and we compared our results with those in Ref. [25], finding full agreement. We also successfully compared the O(αα s ) corrections at vanishing external momentum with the results of an independent calculation based on the effective-potential techniques of Ref. [25]. Note, however, that this comparison does not cover the O(αα s ) contributions to the Z -boson self-energy. Concerning the latter, we checked that we can reproduce the result of Ref. [81] for the subset of two-loop diagrams that involve only quarks and a gluon, taking into account the fact that Ref. [81] employed dimensional regularization.
We then compared our results for the momentumdependent corrections with those of Ref. [53], where the two-loop calculation of the Higgs masses was performed entirely in the DR scheme. As mentioned in Sect. 1, the Higgs-mass corrections in Ref. [53] are organized in a different way with respect to our calculation, therefore we could compare only at the level of individual two-loop self-energies for scalars and pseudoscalars (the two-loop self-energy for the Z boson was not computed in Ref. [53]). Rotating our scalar self-energies from the (S 1 , S 2 ) basis to the (h, H ) basis with the tree-level mixing angle defined in Ref. [53], we reproduce perfectly the results for the "top/gluon" and "top/stop/gluino" contributions to hh shown in figure 2 of that paper. This provides a full cross-check of the momentumdependent O(α t α s ) contribution to the self-energy, as well as a partial check of the O(αα s ) contribution, restricted to diagrams involving the stop squarks (the diagrams involving the other squarks are included in the "others" line in the above-mentioned figure). We also checked the analogous contributions to h H , H H , and AA against results provided by the author of Ref. [53], finding again perfect numerical agreement. Although our calculation and the one in Ref. [53] both use TSIL to compute the MIs, and thus cannot be considered entirely independent, the agreement in the results for the self-energies gives us confidence that the computation of two-loop Feynman diagrams in terms of MIs and the DR subtraction of their divergences are correct in both papers.
Our results for the momentum-dependent O(α t α s ) corrections in the mixed OS-DR scheme can in turn be compared with those of Ref. [56]. To start with, we compared our twoloop 1PI contributions to the scalar and pseudoscalar selfenergies with analogous results provided by the authors of Ref. [56], and we found agreement within the accuracy of the sector-decomposition procedure used to compute the loop integrals in that paper. The successful comparison between two sets of self-energies in which the loop integrals were evaluated with TSIL and SecDec, respectively, validates the results for the two-loop MIs, thus reinforcing our cross-check of Ref. [53]. On the other hand, our results for the momentumdependent O(α t α s ) corrections to the Higgs masses, obtained by combining the 1PI diagrams with all the necessary counterterm contributions, differ significantly from the ones in Ref. [56]. Considering for example the m max h scenario discussed in the previous section, we find that for large m A the lightest-scalar mass is subject to a negative correction of about 350-400 MeV (depending on tan β; see the left plot in Fig. 2), whereas the corresponding correction in Ref. [56] is also negative but quite smaller, about 50-60 MeV (see the upper plot in figure 7 of that paper). We traced the reason for this discrepancy to an inconsistency in Ref. [56], related to the renormalization conditions for the Higgs fields and for tan β.
In the DR scheme, the WFR counterterm for each field H 0 i can be related to the divergent part of the derivative of the scalar self-energy with respect to the external momentum: Indeed, when all parameters entering the one-loop part of the scalar self-energies are renormalized in the DR scheme, Eq. (18) leads to the DR WFR counterterms given in Eq. (9), in accordance with the anomalous dimensions of the Higgs fields given in Refs. [59,60]. However, in the mixed OS-DR scheme of Ref. [56] the top/stop parameters in the one-loop self-energies are renormalized OS. In that case, the use of Eq. (18) to determine the WFR counterterms leads to where α OS t is a scale-independent coupling extracted from the pole top mass, and δm t is the finite one-loop shift for the top mass given in Eq. (B2) of Ref. [25]. By converting the coupling α OS t in the one-loop term of Eq. (19) into the corresponding DR coupling, it is easy to see that Z [56] 2 differs from the DR WFR constant in Eq. (9) by a finite two-loop term: where δ m t denotes the part proportional to in the top selfenergy regularized with dimensional reduction: In the equation above μ R is the renormalization scale associated to the Higgs WFR and to tan β, while B (s, x, y) denotes the coefficient of in the expansion of the Passarino-Veltman function B 0 . An explicit expression for B can be found, e.g., in Eq. (2.31) of the TSIL manual [52].
In the calculation of Ref. [56], where the top/stop parameters entering the one-loop part of the corrections are directly renormalized OS instead of being first renormalized in the DR scheme and then converted to the OS scheme via a finite shift, the two-loop self-energies and tadpoles contain terms proportional to δ m t . Such terms would drop out of the final result for the renormalized inverse-propagator matrix if Eq. (20) was used to obtain the correct DR definition for the WFR constant Z DR 2 , and consequently for δ tan β, but they survive if the WFR constant is defined as in Eq. (19). To prove that these terms are indeed at the root of the observed discrepancies, we modified our own calculation, using Z [56] 2 -as obtained from Eq. (20) -instead of Z DR 2 and then computing a non-minimal counterterm for tan β via Eq. (10). We checked that, with this modification, we reproduce exactly the corrections to the renormalized inverse propagator shown in figures 5 and 10 of Ref. [56]. We also reproduce the corrections to the scalar masses shown in figures 7, 8, 12 and 13 of that paper, although small discrepancies persist in the case of the heaviest scalar when its mass is above the threshold m H = 2 m t . These residual discrepancies are formally of higher order in the perturbative expansion, and they result from different approximations in the two-step procedure for the determination of the poles of the propagator (namely, we drop the imaginary parts of the two-loop self-energies, while Ref. [56] keeps them).
In summary, we have found that in Ref. [56] the two-loop renormalization of the Higgs fields and of the parameter tan β is not performed in the DR scheme as claimed in the paper, but rather in some non-minimal scheme where the WFR counterterms and δ tan β differ from their DR counterparts by finite, non-polynomial terms, and neither the Higgs fields nor tan β obey their usual renormalization-group equations (because of the explicit scale dependence of the additional terms). This inconsistency should be taken into account when comparing the results of Ref. [56] with those of calculations that actually employ DR definitions for the WFR and for tan β. First of all, to account for the difference in δ tan β, the input value for the DR-renormalized parameter tan β should be converted to the corresponding value in the non-minimal scheme of Ref. [56], according to tan β [56] However, Eqs. (5)- (7) show that the contributions of δ tan β (2) to the entries of the Higgs mass matrix are suppressed by powers of cos β. Consequently, the effect on the Higgs masses arising from a two-loop difference in δ tan β is very small already for tan β = 5. In fact, the bulk of the numerical discrepancy between our results and those of Ref. [56] is due to higher-order effects that are directly related to the finite shift in the WFR. As discussed at the end of Sect. 2, such effects are included in the Higgs-mass corrections when the latter are computed in terms of loop-corrected Higgs masses and mixing, and can become numerically relevant when the loopcorrected masses differ significantly from their tree-level values.

Conclusions
We computed the two-loop corrections to the neutral MSSM Higgs masses of O(α t α s ) and O(αα s ) -i.e., all two-loop corrections that involve the strong gauge coupling when the only non-vanishing Yukawa coupling is h t -including the effect of non-vanishing external momenta in the selfenergies. We relied on an IBPs technique to express the momentum-dependent loop integrals in terms of a minimal set of master integrals, and we used the public code TSIL to evaluate the latter. We obtained results for the Higgs-mass corrections valid when all parameters in the one-loop part of the corrections are renormalized in the DR scheme, as well as results valid in a mixed OS-DR scheme where the top/stop parameters are renormalized on-shell. Our results for the scalar and pseudoscalar self-energies in the DR scheme confirm the results of an earlier calculation, Ref. [53], where they overlap. In addition, we obtained new results for the twoloop contributions to the Z -boson self-energy that involve the strong gauge coupling. The latter, which were not computed in Ref. [53], enter the O(αα s ) corrections to the Higgs masses when the tree-level mass matrix of the scalars is expressed in terms of the physical Z -boson mass.
We also compared our results for the momentumdependent O(α t α s ) corrections in the mixed OS-DR scheme with those of a recent calculation of the same corrections, Ref. [56], and found disagreement. We traced the reason for the discrepancy to the fact that, contrary to what stated in Ref. [56], in that calculation the Higgs fields and the parameter tan β are renormalized in a non-minimal scheme instead of the usual DR scheme. When this difference is taken into account, we reproduce the results of Ref. [56], providing in passing a cross-check of the codes used to evaluate the loop integrals in the two calculations (i.e., TSIL and SecDec, respectively). However, we noticed that TSIL, which implements dedicated algorithms for two-loop selfenergy integrals, can be a thousand times faster than a multipurpose code such as SecDec in the computation of the Higgs-mass corrections. This should be taken into consideration when including the momentum-dependent corrections in phenomenological analyses of the MSSM parameter space.
As to the numerical impact of the corrections computed in this paper, it could at best be described as moderate. We considered six benchmark scenarios compatible with the results of Higgs and SUSY searches at the LHC, and we found that both the momentum-dependent part of the O(α t α s ) correc-tions and the whole O(αα s ) corrections can shift the prediction for the lightest-scalar mass m h by several hundred MeV. However, we noticed that -at least in the considered scenarios -the two classes of corrections enter the prediction for m h with opposite sign, and they compensate each other to a good extent. For what concerns the heaviest-scalar mass m H , the impact of the new corrections is also modest, with the exception of regions around real-particle thresholds where a fixed-order calculation is not reliable anyway.
The predictions for the lightest-scalar mass, as obtained from popular codes for the determination of the MSSM mass spectrum, carry a theoretical uncertainty that has been estimated to be (at least) of the order of ±3 GeV -see, e.g., Refs. [78,79] and the more recent discussion in Ref. [82]. Against this backdrop, the corrections presented in this paper can be considered sub-dominant, and their inclusion in public codes might seem less urgent than, e.g., the inclusion of the dominant three-loop effects [41,42] or the proper resummation of large logarithms in scenarios with multi-TeV stop masses [39,40,82], both of which can shift the prediction for the lightest-scalar mass by several GeV. Nevertheless, one should not forget that the accuracy of the measurement of the Higgs mass at the LHC has already reached the level of a few hundred MeV -i.e., comparable to our sub-dominant corrections -and will improve further when more data become available. If SUSY shows up at last when the LHC operates at 13-14 TeV, the Higgs mass will serve as a precision observable to constrain MSSM parameters that might not be directly accessible by experiment. To this purpose, the accuracy of the theoretical prediction will have to match the experimental one, making a full inclusion of the two-loop corrections to the Higgs masses unavoidable. Our calculation should be regarded as a necessary step in that direction.