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 Minimal Supersymmetric Standard Model, including the effect of non-vanishing external momenta in the self-energies. We obtain corrections of [Formula: see text] and [Formula: see text], 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 [Formula: see text] renormalization scheme or a mixed on-shell (OS)-[Formula: see text] 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-[Formula: see text] 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 Large Hadron Collider.

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 high-precision 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 doublets, and , whose relative contribution to electroweak (EW) symmetry breaking is determined by the ratio of vacuum expectation values (VEVs) of their neutral components, . The spectrum of physical Higgs bosons is richer than in the SM, consisting of two neutral scalars, and , one neutral pseudoscalar, , and two charged scalars, . At the tree level, the neutral scalar masses and and the scalar mixing angle can be computed in terms of the -boson mass , the pseudoscalar mass and , and the bound applies. In a significant portion of the parameter space the lightest scalar has SM-like couplings to fermions and gauge bosons, in which case the tree-level bound on 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  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-9] of the importance of the one-loop corrections2 involving top and stop, full one-loop computations of the MSSM Higgs masses have been provided [10-13], leading logarithmic effects at two loops have been included via renormalization-group methods [14-17], and genuine two-loop corrections of  [18-25],  [18, 24, 26], [27, 28] and  [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 have also been obtained, both via renormalization-group methods [38-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-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-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 third-family 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–21, 25–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 renormalization scheme adopted in Refs. [36, 37, 53] for the parameters of the MSSM lagrangian does not match the “mixed on-shell (OS)–” 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 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 masses – as 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 , 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 : 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 with respect to the zero-momentum 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 momentum-dependent 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 and the OS– 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 and 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 corrections to the neutral Higgs masses in the MSSM appeared [56], relying on the results of Ref. [43] for the decomposition of two-loop 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– 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 strong-gauge and top-Yukawa couplings agree numerically with those of Ref. [56], we do not reproduce the overall values of the momentum-dependent 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 .

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 and combine into the pseudoscalar and the neutral would-be Goldstone boson . The CP-even components and combine into the scalars and . The squared physical masses of the latter are the two solutions for of the equationwhere denotes the inverse-propagator matrix in the basis, being the external momentum flowing into the scalar self-energies. We can decompose aswhere denotes the tree-level mass matrix written in terms of renormalized parameters, and 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 . We express in terms of the parameter renormalized in the scheme, and of the physical masses of the pseudoscalar and of the bosonusing (here and thereafter) the shortcuts and for a generic angle . Neglecting terms that do not contribute at or , our choices for the parameters entering lead to the following expressions for the two-loop parts of the individual entries of : In the equations above, and (with ) denote the two-loop parts of tadpoles and self-energies, respectively, for the scalars , while and denote the two-loop parts of the pseudoscalar and -boson self-energies. In addition, (with ) in Eqs. (5)–(7) denote the two-loop parts of the WFR counterterms for the Higgs fields , which we renormalize as follows:where in the expansion of the square root we have again neglected terms that do not contribute at or . We adopt a definition for the , 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 getwhere and are color factors, in dimensional reduction, and the coupling entering the one-loop counterterm is in turn renormalized in the scheme. Finally, in Eqs. (5)–(7) denotes the two-loop part of the counterterm for the parameter . The choice of a definition for 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 and corrections to the Higgs masses, we need to introduce counterterms only for the parameters that are subject to corrections, namely the top mass , the stop masses and , the stop mixing angle , the soft supersymmetry-breaking Higgs-stop coupling and the masses 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 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 and contributions to tadpoles and self-energies assuming that the relevant quark/squark parameters are renormalized either in the or in the OS scheme. Formulas for the –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 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 and 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 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. (2) including in the full one-loop corrections plus the dominant two-loop corrections of , , and computed in the approximation of vanishing external momentum. From FeynHiggs we obtain the scalar masses and , and an effective mixing angle which diagonalizes the loop-corrected scalar mass matrix at vanishing external momentum. Our full results for the scalar masses are then obtained by adding to the results of FeynHiggs the momentum-dependent parts of the corrections and the whole corrections:Concerning the former, we have where we define , and retain only the real and finite part of the contributions to the two-loop self-energies. For what concerns the corrections, they contain all terms from Eqs. (5)–(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 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 counterterm). Indeed, if the two-loop contributions to the inverse-propagator matrix are computed with 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 drop out of the mass corrections . On the other hand, if the loop-corrected scalar masses and 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 definition for the WFR; therefore the terms involving 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 . 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 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 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 -dimensional gluon, (where is the -dimensional metric tensor) in order to include the corresponding -scalar contribution. We are then left with Feynman integrals of the formwhere are positive (or zero) integer exponents and the ’s are defined asIntegrals belonging to the class above are in general not linearly independent of each other. When the scalar products in the numerator are expressed in terms of the denominators, powers of a present in the original integral might cancel against a in the numerator, possibly generating a Feynman integral in which does not appear, i.e. in which the corresponding line has been shrunk to a point. For given ’s and high enough , some ’s may acquire negative exponents. The computation of a Feynman integral of the type in Eq. (16) therefore reduces to the evaluation of a number of integrals of the formwhere the exponents . In the present case, one has to evaluate 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 . 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 . 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 first-order 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 is possible only in some simple cases. In more general cases the MIs are expanded in powers of , 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 one-propagator 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 one-loop 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 self-energies entering Eqs. (5)–(7) require the knowledge of MIs with up to four different masses, the most complicated ones being those involving simultaneously , and .

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 self-energy MIs. The method of Refs. [47-51] is based on the fact that the value at (or the expansion for small 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 . Given that we include the contribution of light quarks to (in the approximation ) and that 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 three-particle thresholds show up, together with the pseudothresholds, as poles in the coefficients of the DEs. TSIL overcomes the numerical instabilities related to such poles by displacing the -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 and , the latter being the typical time for close or equal to the heavy stop pair threshold and to the three-particle (pseudo)thresholds . In Ref. [52] the relative accuracy of TSIL is claimed to be better than 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, are needed in order to reach a relative accuracy of at least close to a threshold.

Numerical examples

In this section we assess the numerical impact of the momentum-dependent part of the corrections and of the whole 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 agrees, within the theoretical uncertainty of 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– scheme described in Sect. 2. The SM input parameters are chosen as the pole top mass , the running bottom mass , the Fermi constant , the strong gauge coupling , and the pole gauge-boson masses and . To compute the scalar masses and the effective mixing angle entering the corrections in Eqs. (11)–(15), we call FeynHiggs version 2.10.2. We use default values for all settings with the exception of runningMT = 0, i.e. the top mass in the radiative corrections is identified with the pole mass (to match the renormalization conditions imposed both in our OS– calculation and in the one of Ref. [56]). By default, the renormalization scale associated to the definition of the Higgs WFR and of is fixed as .

In Fig. 1 we present our predictions for the lightest-scalar mass in the six benchmark scenarios. We choose and , so that the lightest scalar 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 , obtained from FeynHiggs without additional corrections; the middle bar includes the effect of the momentum-dependent part of the corrections, i.e. the defined in Eq. (12); finally, the lower bar represents our final result for , and includes the effects of both the momentum-dependent part of the corrections and the corrections, i.e. the defined in Eq. (14).

Fig. 1

Predictions for the mass of the lightest scalar in the six benchmark scenarios of Ref. [77], for and . For each scenario, the three bars show: the “unperturbed” mass computed with FeynHiggs 2.10.2 (upper), the inclusion of the momentum-dependent part of the corrections (middle) and the additional inclusion of the whole corrections (lower). From top to bottom, the considered scenarios are (red), (blue), (green), light stop (turquoise), light stau (purple), tau-phobic (orange)

Figure 1 shows that, in all considered scenarios, the momentum-dependent part of the corrections and the whole corrections can shift the prediction for by several hundred MeV each (the largest shifts, of about  GeV, occur in the scenario). However, in all of our examples the two classes of corrections happen to be similar to each other in magnitude, and to enter the prediction for with opposite signs. As a result, their combined effect is always fairly small, less than  MeV.

In Fig. 2 we illustrate the impact of the momentum-dependent corrections and of the corrections on the prediction for the lightest-scalar mass as a function of , and in Fig. 3 we do the same for the heaviest-scalar mass. In each figure, the MSSM parameters for the left plot are chosen as in the benchmark scenario of Ref. [77], while for the right plot they are chosen as in the modified light-stop scenario. In each plot, the dashed lines represent the contribution of the sole momentum-dependent part of the corrections, while the solid lines include both the momentum-dependent corrections and the corrections. The red lines were obtained with while the blue lines were obtained with .

Fig. 2

Corrections to the lightest-scalar mass as a function of , for (red) and for (blue). The other MSSM parameters are chosen as in the scenario (left) or as in the light-stop scenario (right). The dashed lines represent the contribution of the sole momentum-dependent part of the corrections, the solid lines include both the momentum-dependent corrections and the corrections

Fig. 3

Same as Fig. 2 for the corrections to the heaviest-scalar mass

Figure 2 shows that the corrections to the lightest-scalar mass are negligible at low values of , but they become larger and essentially independent of as the latter increases. The transition to this “decoupling” regime – where the lightest scalar has SM-like couplings and its mass is insensitive to the value of – is sharper for larger values of . In both the and light-stop scenarios, the momentum-dependent effects decrease by 300–400 MeV at large . However, as already seen in Fig. 1, the effects enter the prediction for with comparable magnitude but opposite sign, significantly reducing the total size of the correction.

Figure 3 shows that for low values of , where the heaviest scalar is the one with SM-like couplings, the corrections to its mass are comparable to the ones that affect the lightest-scalar mass in the decoupling region. On the other hand, for larger values of the pseudoscalar mass – where – 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 in the plot on the left for the scenario. More-pronounced spikes (note the different scale on the axis) are visible in correspondence with , , and in the plot on the right for the light-stop scenario. Analogous spikes would appear at larger values of in the 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 heaviest-scalar mass amount at most to a few hundred MeV, and they decrease in size with increasing . Moreover, the effect of the 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 corrections to are suppressed by one or two powers of , whereas the 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 itself, they are not entirely irrelevant when compared to the difference , 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 and 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 scheme, i.e. we introduce minimal counterterms that, by definition, subtract only powers of multiplied by coefficients that should be polynomial in the -renormalized masses and couplings. In a second step, we convert our results to the mixed OS– scheme adopted in FeynHiggs, replacing the 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 corrections and we compared our results with those in Ref. [25], finding full agreement. We also successfully compared the 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 contributions to the -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 momentum-dependent corrections with those of Ref. [53], where the two-loop calculation of the Higgs masses was performed entirely in the 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 boson was not computed in Ref. [53]). Rotating our scalar self-energies from the basis to the 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 shown in figure 2 of that paper. This provides a full cross-check of the momentum-dependent contribution to the self-energy, as well as a partial check of the 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 , , and 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 subtraction of their divergences are correct in both papers.

Our results for the momentum-dependent corrections in the mixed OS– scheme can in turn be compared with those of Ref. [56]. To start with, we compared our two-loop 1PI contributions to the scalar and pseudoscalar self-energies 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 momentum-dependent 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 scenario discussed in the previous section, we find that for large the lightest-scalar mass is subject to a negative correction of about 350–400 MeV (depending on ; 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 .

In the scheme, the WFR counterterm for each field 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 scheme, Eq. (18) leads to the 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– 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 towhere is a scale-independent coupling extracted from the pole top mass, and is the finite one-loop shift for the top mass given in Eq. (B2) of Ref. [25]. By converting the coupling in the one-loop term of Eq. (19) into the corresponding coupling, it is easy to see that differs from the WFR constant in Eq. (9) by a finite two-loop term:where denotes the part proportional to in the top self-energy regularized with dimensional reduction:In the equation above is the renormalization scale associated to the Higgs WFR and to , while denotes the coefficient of in the expansion of the Passarino–Veltman function . An explicit expression for 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 scheme and then converted to the OS scheme via a finite shift, the two-loop self-energies and tadpoles contain terms proportional to . Such terms would drop out of the final result for the renormalized inverse-propagator matrix if Eq. (20) was used to obtain the correct definition for the WFR constant , and consequently for , 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 – as obtained from Eq. (20) – instead of and then computing a non-minimal counterterm for 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 . 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 is not performed in the scheme as claimed in the paper, but rather in some non-minimal scheme where the WFR counterterms and differ from their counterparts by finite, non-polynomial terms, and neither the Higgs fields nor 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 definitions for the WFR and for . First of all, to account for the difference in , the input value for the -renormalized parameter should be converted to the corresponding value in the non-minimal scheme of Ref. [56], according toHowever, Eqs. (5)–(7) show that the contributions of to the entries of the Higgs mass matrix are suppressed by powers of . Consequently, the effect on the Higgs masses arising from a two-loop difference in is very small already for . 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 loop-corrected masses differ significantly from their tree-level values.

Conclusions

We computed the two-loop corrections to the neutral MSSM Higgs masses of and – i.e., all two-loop corrections that involve the strong gauge coupling when the only non-vanishing Yukawa coupling is – including the effect of non-vanishing external momenta in the self-energies. 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 scheme, as well as results valid in a mixed OS– scheme where the top/stop parameters are renormalized on-shell. Our results for the scalar and pseudoscalar self-energies in the scheme confirm the results of an earlier calculation, Ref. [53], where they overlap. In addition, we obtained new results for the two-loop contributions to the -boson self-energy that involve the strong gauge coupling. The latter, which were not computed in Ref. [53], enter the corrections to the Higgs masses when the tree-level mass matrix of the scalars is expressed in terms of the physical -boson mass.

We also compared our results for the momentum-dependent corrections in the mixed OS– 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 are renormalized in a non-minimal scheme instead of the usual 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 self-energy integrals, can be a thousand times faster than a multi-purpose 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 corrections and the whole corrections can shift the prediction for the lightest-scalar mass by several hundred MeV. However, we noticed that – at least in the considered scenarios – the two classes of corrections enter the prediction for with opposite sign, and they compensate each other to a good extent. For what concerns the heaviest-scalar mass , 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  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.